dr*1000\dz

更改了间隔dr。

cz.py

@@ -4,7 +4,7 @@ warnings.filterwarnings('ignore')
 import argparse, json
 from common import NoIndent, MyEncoder
 import os, sys
-sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))+'/pyat')
+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)) + '/pyat')
 
 def _argparse():
     parser = argparse.ArgumentParser(description='Average Noise Spectral Level')
@@ -18,13 +18,13 @@ def main():
     if argfile[-5:] != '.json':
         print('输入的参数文件名有误')
         return
-    with open(argfile,'r') as jsonFile:
+    with open(argfile, 'r') as jsonFile:
         paramData = json.load(jsonFile)
         if ("ETOPO1Path" not in paramData) or ("STPath" not in paramData) or \
-            ("freq" not in paramData) or ("sourceDepth" not in paramData) or ("receiveDepth" not in paramData) or\
-            ("lat" not in paramData) or ("lon" not in paramData) or\
-            ("timeType" not in paramData) or\
-            ("outParam" not in paramData):
+                ("freq" not in paramData) or ("sourceDepth" not in paramData) or ("receiveDepth" not in paramData) or \
+                ("lat" not in paramData) or ("lon" not in paramData) or \
+                ("timeType" not in paramData) or \
+                ("outParam" not in paramData):
             print(f'{argfile}中的参数有误，请核实！')
             return
         ETOPO1Path = paramData["ETOPO1Path"]
@@ -42,7 +42,7 @@ def main():
         outParam = paramData["outParam"]
         if "outPath" not in outParam:
             outParam["outPath"] = ""
-    
+
         from DataBase import DataBase
         from czcore import call_cz
         import numpy as np
@@ -50,7 +50,7 @@ def main():
         # 定义时间段
         timeIndex = timeType
         # 创建类
-        database = DataBase(ETOPO1Path, STPath,dataSource="")
+        database = DataBase(ETOPO1Path, STPath, dataSource="")
         # 设置时间段
         database.set_timePeriod(timeIndex)
 
@@ -60,11 +60,11 @@ def main():
         # czdata.TL是损失、czdata.bathym是地形数据
         # czdata.Rays是声线数据
         # czdata.CzDetail是会聚区数据，包括会聚区宽度和距离
-        res = call_cz(database, freq, lat, lon, azi, sourceDepth, receiveDepth)
+        res = call_cz(database, freq, lat, lon, azi, sourceDepth)
         # call_cz(25, 18.0, 114.0, 10, 3500)
 
         # write data to file
-        #输出json文件
+        # 输出json文件
         outData = {}
         outData["param"] = paramData
         if len(res) == 0:
@@ -72,10 +72,10 @@ def main():
             outData["data"] = []
         else:
             outData["result"] = "有会聚区"
-            outDataData, outBathymData,outTlData,outRayData = [],[],[],[]
+            outDataData, outBathymData, outTlData, outRayData = [], [], [], []
             json.encoder.FLOAT_REPR = lambda x: format(x, '.3f')
             for index in range(len(res)):
-                bathymData,tlData,rayData = {},{},{}
+                bathymData, tlData, rayData = {}, {}, {}
                 # 保存地形数据
                 res[index].bathym.r = np.around(res[index].bathym.r, decimals=0)
                 res[index].bathym.d = np.around(res[index].bathym.d, decimals=0)
@@ -84,18 +84,20 @@ def main():
                     "azi": f'{round(res[index].azi)}',
                     "count": len(res[index].bathym.d),
                     "descript": "距离、深度",
-                    "data": NoIndent(list(zip(res[index].bathym.r,res[index].bathym.d)))
+                    "data": NoIndent(list(zip(res[index].bathym.r, res[index].bathym.d)))
                 }
                 # 保存传播损失数据
                 tlpoint = []
                 res[index].depths = np.around(res[index].depths, decimals=0)
                 res[index].ranges = np.around(res[index].ranges, decimals=0)
                 res[index].TL = np.around(res[index].TL, decimals=2)
-                tlpoint = list(zip(np.tile(res[index].depths.reshape(len(res[index].depths),1),len(res[index].ranges)).flatten(), np.tile(res[index].ranges,(len(res[index].depths),1)).flatten(), res[index].TL.flatten()))
+                tlpoint = list(
+                    zip(np.tile(res[index].depths.reshape(len(res[index].depths), 1), len(res[index].ranges)).flatten(),
+                        np.tile(res[index].ranges, (len(res[index].depths), 1)).flatten(), res[index].TL.flatten()))
                 tlData = {
                     "param": paramData,
                     "azi": f'{round(res[index].azi)}',
-                    "count": len(res[index].ranges)*len(res[index].depths),
+                    "count": len(res[index].ranges) * len(res[index].depths),
                     "descript": "距离、深度、损失",
                     "data": NoIndent(tlpoint)
                 }
@@ -106,7 +108,7 @@ def main():
                     ray.xy[1] = np.around(ray.xy[1], decimals=2)
                     tempdata = {
                         "pointCount": len(ray.xy[0]),
-                        "data": NoIndent(list(zip(ray.xy[0],ray.xy[1])))
+                        "data": NoIndent(list(zip(ray.xy[0], ray.xy[1])))
                     }
                     raypoint.append(tempdata)
                 rayData = {
@@ -124,9 +126,9 @@ def main():
                     directczdetaillist = []
                     for directcz in res[index].CzDetail.directCz:
                         directczdetail = {
-                            "minRange": f'{round(directcz.minRange,1)}',
-                            "maxRange": f'{round(directcz.maxRange,1)}',
-                            "width": f'{round(directcz.width,1)}',
+                            "minRange": f'{round(directcz.minRange, 1)}',
+                            "maxRange": f'{round(directcz.maxRange, 1)}',
+                            "width": f'{round(directcz.width, 1)}',
                         }
                         directczdetaillist.append(directczdetail)
                     czDeatilData["directCz"] = directczdetaillist
@@ -136,9 +138,9 @@ def main():
                     reflectczdetaillist = []
                     for reflectcz in res[index].CzDetail.reflectCz:
                         reflectczdetail = {
-                            "minRange": f'{round(reflectcz.minRange,1)}',
-                            "maxRange": f'{round(reflectcz.maxRange,1)}',
-                            "width": f'{round(reflectcz.width,1)}',
+                            "minRange": f'{round(reflectcz.minRange, 1)}',
+                            "maxRange": f'{round(reflectcz.maxRange, 1)}',
+                            "width": f'{round(reflectcz.width, 1)}',
                         }
                         reflectczdetaillist.append(reflectczdetail)
                     czDeatilData["reflectCz"] = reflectczdetaillist
@@ -146,19 +148,26 @@ def main():
                 outBathymData.append(bathymData)
                 outTlData.append(tlData)
                 outRayData.append(rayData)
-            with open(("" if outParam["outPath"]=="" else outParam["outPath"] + "/")+"czOut_bathym.json", 'w') as jsonFile:
+            with open(("" if outParam["outPath"] == "" else outParam["outPath"] + "/") + "czOut_bathym.json",
+                      'w') as jsonFile:
                 json_data = json.dumps(outBathymData, cls=MyEncoder, indent=2, ensure_ascii=False)
                 jsonFile.write(json_data)
-            with open(("" if outParam["outPath"]=="" else outParam["outPath"] + "/")+"czOut_tl.json", 'w') as jsonFile:
+            with open(("" if outParam["outPath"] == "" else outParam["outPath"] + "/") + "czOut_tl.json",
+                      'w') as jsonFile:
                 json_data = json.dumps(outTlData, cls=MyEncoder, indent=2, ensure_ascii=False)
                 jsonFile.write(json_data)
-            with open(("" if outParam["outPath"]=="" else outParam["outPath"] + "/")+"czOut_ray.json", 'w') as jsonFile:
+            with open(("" if outParam["outPath"] == "" else outParam["outPath"] + "/") + "czOut_ray.json",
+                      'w') as jsonFile:
                 json_data = json.dumps(outRayData, cls=MyEncoder, indent=2, ensure_ascii=False)
                 jsonFile.write(json_data)
             outData["data"] = outDataData
-            
-        with open(("" if outParam["outPath"]=="" else outParam["outPath"] + "/")+"czOut.json", 'w') as jsonFile:
+
+        with open(("" if outParam["outPath"] == "" else outParam["outPath"] + "/") + "czOut.json", 'w') as jsonFile:
             json.dump(outData, jsonFile, indent=4, ensure_ascii=False)
-    
+
+# import time
 if __name__ == '__main__':
-    main()
+    # start = time.perf_counter()
+    main()
+    # end = time.perf_counter()
+    # print('time = ',(end-start))

pyat/bellhop.py

@@ -3,6 +3,9 @@ from env import Pos, Source, Dom, cInt, SSPraw, SSP, HS, BotBndry, TopBndry, Bnd
 from os import system
 import numpy as np
 from config import AtBinPath
+from multiprocessing import Pool
+from scipy.stats import norm
+import math
 
 
 def write_ssp(sspfile, ssp, bathm, NZmax):
@@ -43,15 +46,18 @@ def calZmax(ssp):
     Zmax = ssp[ssp_idx].z[-1]
     return NZmax, Zmax, ssp_idx
 
-def call_Bellhop(freq, sd, rd, bathm, ssp, sed, base):
+def call_Bellhop_p(filename):
+    system(AtBinPath + "/bellhop " + filename)
+
+def call_Bellhop(freq, sd, bathm, ssp, sed, base, dr=0.1):
     # Source and receiving position
     filename = 'envB'
     Rmax = bathm.r[-1]  # unit: km
-    dr = (1500 / freq) / 1000  # wavelength / 10, unit: km
+    # dr = 0.1  # wavelength / 10, unit: km
     ran = np.linspace(0, Rmax, round(Rmax / dr))  # Receiving range (unit: km)
     NZmax, Zmax, ssp_idx = calZmax(ssp)
     RD = np.linspace(0, Zmax, round(Zmax / 10) + 1)
-    pos = Pos(Source(sd), Dom(ran, rd))
+    pos = Pos(Source(sd), Dom(ran, RD))
     # Range of phase velocity
     cInt.Low = 1400
     cInt.High = 15000
@@ -96,20 +102,35 @@ def call_Bellhop(freq, sd, rd, bathm, ssp, sed, base):
     bdy = Bndry(top, bottom)
     # Beam params
     run_type = 'C'  # incoherently sum beams, see AT docs for more info
-    nbeams = 3000
-    alpha = np.linspace(-45, 45, 2)  # min and max launch angle -20 degrees to 20 degrees
+
+    NAlphaRange = 12
+    Alpha = alphadiv(NAlphaRange,Rmax) # min and max launch angle -20 degrees to 20 degrees
+    totalBeams = beamsnumber(freq,Rmax,max(bathm.r))
     box = Box(Zmax, max(bathm.r))  # bound the region you let the beams go, depth in meters and range in km
     deltas = 0  # length step of ray trace, 0 means automatically choose
-    beam = Beam(RunType=run_type, Nbeams=nbeams, alpha=alpha, box=box, deltas=deltas)  # package
 
     # Write *.env file
-    write_env(filename + '.env', 'BELLHOP', 'Pekeris profile', freq, sspB, bdy, pos, beam, cInt, Rmax)
-    write_ssp(filename, ssp, bathm, NZmax)
-    write_bathy(filename, bathm)
+    ialpha = 0
+    Filenames = []
+    for iAlphaRange in range(len(Alpha)-1):
+        alpha = np.array([Alpha[iAlphaRange], Alpha[iAlphaRange + 1]])
+        nbeams = int(totalBeams * (alpha[1]-alpha[0])/180)
+
+        beam = Beam(RunType=run_type, Nbeams=nbeams, alpha=alpha, box=box, deltas=deltas)  # package
+        filenameI = filename+ str(iAlphaRange)
+        write_env(filenameI + '.env', 'BELLHOP', 'Pekeris profile', freq, sspB, bdy, pos, beam, cInt, Rmax)
+        write_ssp(filenameI, ssp, bathm, NZmax)
+        write_bathy(filenameI, bathm)
+        Filenames.append(filenameI)
+
+    pool = Pool(NAlphaRange)
+    pool.map(call_Bellhop_p, Filenames)
 
-    system(AtBinPath + "/bellhop " + filename)
     # Read sound field
-    [x, x, x, x, Pos1, pressure] = read_shd(filename + '.shd')
+    [x, x, x, x, Pos1, pressure] = read_shd(filename + '0.shd')
+    for iAlphaRange in range(NAlphaRange-1):
+        [x, x, x, x, Pos1, pressure1] = read_shd(filename + str(iAlphaRange) + '.shd')
+        pressure = pressure + pressure1
     TL = -10 * np.log10(np.square(np.abs(pressure)))  # Transmission Loss
 
     return Pos1, TL
@@ -197,6 +218,37 @@ def find_cvgcRays(rays_total):
             Rays.append(ray)
     return Rays
 
+def alphadiv(NalphaRange,Rmax):
+    angle = [-90,90]
+    if Rmax >100 and Rmax<=1000:
+        sigma = 20
+    elif Rmax >=1000:
+        sigma = 15
+    else :
+        sigma = 25
+    for i in range(1,NalphaRange):
+        angle.append(norm.ppf(i/NalphaRange,loc=0, scale=sigma))
+
+    angle.sort()
+    return angle
+
+
+
+def beamsnumber(freq,Rmax,depth):
+    if Rmax >700 :
+        nbeams = max(int(0.01 * Rmax*1000 * freq /1500),300)
+    elif Rmax >300 and Rmax<=700:
+        nbeams = max(int(0.05 * Rmax*1000 * freq / 1500), 300)
+    elif Rmax >100 and Rmax<=300:
+        nbeams = max(int(0.1 * Rmax*1000 * freq / 1500), 300)
+    else:
+        nbeams = max(int(0.3 * Rmax*1000 * freq / 1500), 300)
+
+    d_theta_recommended =math.atan(depth/(10.0 * Rmax*1000))
+    nbeams = min(int(math.pi/d_theta_recommended),nbeams,3000)
+
+    return nbeams
+
 # if __name__=="__main__":
     ''' ********************************************** '''
     ''' ******************** Test ******************** '''

pyat/czcore.py

@@ -11,8 +11,6 @@ from env import Pos, Source, Dom, cInt, SSPraw, SSP, HS, BotBndry, TopBndry, Bnd
 import matplotlib.pyplot as plt
 from config import AtBinPath
 
-file_cz_data = open('cz_data.txt', 'w')
-
 
 class CzRangeWidth:
     def __init__(self, minRange, maxRange, width):
@@ -48,11 +46,11 @@ def calZmax(ssp, sed):
     return Zmax, ssp_idx
 
 
-def call_Bellhop(freq, sd, bathm, ssp, sed, base,azi):
+def call_Bellhop(freq, sd, dr, bathm, ssp, sed, base,azi):
     # Source and receiving position
     filename = 'cz'+str(azi)
     Rmax = bathm.r[-1]  # unit: km
-    dr = (1500 / freq) / 1000  # wavelength / 10, unit: km
+    # dr = 0.1 # wavelength / 10, unit: km
     ran = np.linspace(0, Rmax, round(Rmax / dr))  # Receiving range (unit: km)
     Zmax, ssp_idx = calZmax(ssp, sed)
     """RD 已修改 8.26"""
@@ -100,9 +98,9 @@ def call_Bellhop(freq, sd, bathm, ssp, sed, base,azi):
     bdy = Bndry(top, bottom)
     # Beam params
     run_type = 'I'  # incoherently sum beams, see AT docs for more info
-    nbeams = 3000
+    nbeams = 1000
     alpha = np.linspace(-90, 90, 2)  # min and max launch angle -20 degrees to 20 degrees
-    box = Box(Zmax, 100)  # bound the region you let the beams go, depth in meters and range in km
+    box = Box(Zmax, Rmax+10)  # bound the region you let the beams go, depth in meters and range in km
     deltas = 0  # length step of ray trace, 0 means automatically choose
     beam = Beam(RunType=run_type, Nbeams=nbeams, alpha=alpha, box=box, deltas=deltas)  # package
     """删除ssp 已修改 8.26"""
@@ -118,11 +116,11 @@ def call_Bellhop(freq, sd, bathm, ssp, sed, base,azi):
     return Pos1, TL
 
 
-def call_Bellhop_Rays(freq, sd, bathm, ssp, sed, base,azi):
+def call_Bellhop_Rays(freq, sd, dr, bathm, ssp, sed, base,azi):
     # Source and receiving position
     filename = 'cz_ray'+str(azi)
     Rmax = bathm.r[-1]  # unit: km
-    dr = (1500 / freq) / 1000  # wavelength / 10, unit: km
+    # dr = 0.1  # wavelength / 10, unit: km
     ran = np.linspace(0, Rmax, round(Rmax / dr))  # Receiving range (unit: km)
     Zmax, ssp_idx = calZmax(ssp, sed)
     """RD 已修改 8.26"""
@@ -172,7 +170,7 @@ def call_Bellhop_Rays(freq, sd, bathm, ssp, sed, base,azi):
     run_type = 'R'  # incoherently sum beams, see AT docs for more info
     nbeams = 301
     alpha = np.linspace(-15, 15, 2)  # min and max launch angle -20 degrees to 20 degrees
-    box = Box(Zmax, 100)  # bound the region you let the beams go, depth in meters and range in km
+    box = Box(Zmax, Rmax+10)  # bound the region you let the beams go, depth in meters and range in km
     deltas = 0  # length step of ray trace, 0 means automatically choose
     beam = Beam(RunType=run_type, Nbeams=nbeams, alpha=alpha, box=box, deltas=deltas)  # package
     """删除ssp 已修改 8.26"""
@@ -227,8 +225,6 @@ def call_hjq(r, TL):
                     m = r[out[0]]
                     n = r[out[len(out) - 1]]
                     l = n - m
-                    # file_cz_data.write('\n'"会聚区距离：%.1fm"%m +"-%.1fm" %n)
-                    # file_cz_data.write('\n'"会聚区宽度：%.1fm\n"%l )
                     czrangewidth_list.append(CzRangeWidth(m, n, l))
                     out.clear()
                     signal = 0
@@ -236,13 +232,10 @@ def call_hjq(r, TL):
                 m = r[out[0]]
                 n = r[out[len(out) - 1]]
                 l = n - m
-                # file_cz_data.write('\n'"会聚区距离：%.1fm" % m + "-%.1fm" % n)
-                # file_cz_data.write('\n'"会聚区宽度：%.1fm\n" % l )
                 czrangewidth_list.append(CzRangeWidth(m, n, l))
                 out.clear()
                 signal = 0
     # if signal==1:
-    # file_cz_data.write('\n'" 没有出现会聚区\n")
     # 滤波器效果图：
     # plt.plot(r, TL)
     # plt.plot(r, TL_R)
@@ -266,18 +259,13 @@ def call_hjq(r, TL):
     return czrangewidth_list
 
 
-def call_fangwei(freq, sd, bathm, ssp, sed, base, azi, isPlot=False):
+def call_fangwei(freq, sd, dr, bathm, ssp, sed, base, azi, isPlot=False):
     czdetail = CzDetail("", [], [])
     last = len(ssp[0].c) - 1
-    # if ssp[0].c[0] > ssp[0].c[last]:
-    # file_cz_data.write('\n'"直达会聚区：无")
-    # file_cz_data.write('\n'"反射会聚区：")
-    # else:
-    # file_cz_data.write('\n'"直达会聚区：")
     sed = None
-    rays_total = call_Bellhop_Rays(freq, sd, bathm, ssp, sed, base, azi)
+    rays_total = call_Bellhop_Rays(freq, sd, dr, bathm, ssp, sed, base, azi)
     Rays = find_cvgcRays(rays_total)
-    [Pos, TL] = call_Bellhop(freq, sd, bathm, ssp, sed, base, azi)
+    [Pos, TL] = call_Bellhop(freq, sd, dr, bathm, ssp, sed, base, azi)
     [x, x, x, x, Pos, pressure] = read_shd('cz'+str(azi)+'.shd')
     TL = -10 * np.log10(np.square(np.abs(pressure)))  # Transmission Loss
     # Plot
@@ -297,20 +285,12 @@ def call_fangwei(freq, sd, bathm, ssp, sed, base, azi, isPlot=False):
         czdetail.reflectCz = []
     return Pos.r.range, Pos.r.depth, TL[0, 0, :, :], bathm, Rays, czdetail
 
-def child(database, freq, lat, lon, sd, rd, i):
-    descripts = ["正北方向", "东北方向", "正东方向", "东南方向", "正南方向", "西南方向", "正西方向", "西北方向"]
-    theta = np.arange(0, 359, 45)
-    [bathm, ssp, sed, base] = GetEnv(database, 25, Coord(lon, lat), R=100, azi=theta[i])
-    if np.mean(bathm.d) < 2000:
-        # 不会出现会聚区
-        ddd=[]
-        return None
-    else:
-        ranges, depths, TL, bathym, Rays, czdetail = call_fangwei(freq, sd, rd, bathm, ssp, sed, base, azi=theta[i],isPlot=True)
-        czdetail.descript = descripts[i]
-        return CzData(theta[i], ranges, depths, TL, bathym, Rays, czdetail)
+def child(freq, sd, dr, isPlot, bathm, ssp, sed, base, azi, descript):
+    ranges, depths, TL, bathym, Rays, czdetail = call_fangwei(freq, sd, dr, bathm, ssp, sed, base, azi, isPlot=isPlot)
+    czdetail.descript = descript
+    return CzData(azi, ranges, depths, TL, bathym, Rays, czdetail)
 
-def call_cz(database, freq, lat, lon, azi, sourceDepth=10,  isPlot=False):
+def call_cz(database, freq, lat, lon, azi, sourceDepth=10, R=100, dr=0.1, isPlot=False):
     sd = np.array([sourceDepth])  # Source depth
 
     results = []
@@ -321,30 +301,35 @@ def call_cz(database, freq, lat, lon, azi, sourceDepth=10,  isPlot=False):
     #     return czdata
     ''' 计算单个方向的会聚区 '''
     if azi != None:
-        [bathm, ssp, sed, base] = GetEnv(database, 25, Coord(lon, lat), R=100, azi=azi)
-        if np.mean(bathm.d) < 2000:
-            return None
-        else:
-            ranges, depths, TL, bathym, Rays, czdetail = call_fangwei(freq, sd, rd, bathm, ssp, sed, base, azi,
+        [bathm, ssp, sed, base] = GetEnv(database, 25, Coord(lon, lat), R=R, azi=azi)
+        if np.mean(bathm.d) >= 2000:
+            ranges, depths, TL, bathym, Rays, czdetail = call_fangwei(freq, sd, dr, bathm, ssp, sed, base, azi,
                                                                       isPlot=isPlot)
             czdetail.descript = "{:.1f}°方向".format(azi)
             results.append(CzData(0, ranges, depths, TL, bathym, Rays, czdetail))
-            return results
+        return results
 
     ''' 计算8个方向的会聚区 '''
-    [bathm, ssp, sed, base] = GetEnv(database, 25, Coord(lon, lat), R=100, azi=0)
-    if bathm.d[0] < 2000:
-        # file_cz_data.write("各个方位都不会出现会聚区")
-        return results
-    else:
-        pool = Pool(processes=8)
-        applyResults = []
-        for pi in range(8):
-            applyResults.append(pool.apply_async(partial(child, database, freq, lat, lon, sd, rd), (pi,)))
-            # pi=Process(target=child_0,args=(pi,))
-            # pi.start()
-        pool.close()
-        pool.join()
-        for applyResult in applyResults:
-            results.append(applyResult.get())
-        return results
+    pool = Pool()
+    applyResults = []
+    descripts = ["正北方向", "东北方向", "正东方向", "东南方向", "正南方向", "西南方向", "正西方向", "西北方向"]
+    theta = np.arange(0, 359, 45)
+    for i in range(8):
+        [bathm, ssp, sed, base] = GetEnv(database, 25, Coord(lon, lat), R=R, azi=theta[i])
+        if np.mean(bathm.d) >= 2000:
+            applyResults.append(pool.apply_async(partial(child, freq, sd, dr, isPlot), (bathm, ssp, sed, base, theta[i], descripts[i])))
+        # applyResults.append(pool.apply_async(partial(child, database, freq, lat, lon, sd, isPlot), (pi,)))
+    pool.close()
+    pool.join()
+    for applyResult in applyResults:
+        results.append(applyResult.get())
+    return results
+# timeIndex = '20190101'
+# # 设置时间段
+#
+# STPath = '/home/gan/下载/code'
+# ETOPO1Path = '/home/gan/下载/code/20220920/ETOPO1'
+# database = DataBase(ETOPO1Path, STPath)
+# database.set_timePeriod(timeIndex)
+# res = call_cz(database, 300, 13, 115, None, 60, 150, 0.1, True)
+# a = 1

pyat/krakenc.py

@@ -15,7 +15,7 @@ def calZmax(ssp, sed):
     return Zmax
 
 
-def call_krakenc(freq, sd, rd, bathm, ssp, sed, base):
+def call_krakenc(freq, sd, bathm, ssp, sed, base, dr=0.1):
     filename = "envK"
     envFilename = "envK.env"
     model = "KRAKEN"
@@ -23,12 +23,13 @@ def call_krakenc(freq, sd, rd, bathm, ssp, sed, base):
 
     # Source and receiving position
     Rmax = bathm.r[-1]       # unit: km
-    dr = (1500 / freq / 3) / 1000  # wavelength / 10, unit: km
+    # dr = (1500 / freq / 3) / 1000  # wavelength / 10, unit: km
     ran = np.linspace(0, Rmax, round(Rmax / dr))  # Receiving range (unit: km)
     Zmax = calZmax(ssp, sed)
     RD = np.linspace(0, Zmax, round(Zmax) + 1)
     pos = Pos(Source(sd), Dom(ran, RD))
 
+
     # Range of phase velocity
     cInt.Low = 1400
     cInt.High = 15000
@@ -98,7 +99,7 @@ def call_krakenc(freq, sd, rd, bathm, ssp, sed, base):
     system(AtBinPath + "/krakenc " + filename)
     # Write *.flp file
     Opt_flp = 'RC'
-    pos.r.depth = rd
+    pos.r.depth = RD
     write_fieldflp(filename, Opt_flp, pos, NProf=len(ssp), rProf=bathm.r)
     # Call field.exe
     system(AtBinPath + "/field " + filename)
@@ -108,6 +109,7 @@ def call_krakenc(freq, sd, rd, bathm, ssp, sed, base):
 
     return Pos1, TL
 
+
 # if __name__=="__main__":
     ''' ********************************************************************************* '''
     ''' ***************************** Input Parameters ********************************** '''

pyat/ram.py

@@ -3,7 +3,17 @@ from os import system
 import numpy as np
 from config import AtBinPath
 
-def call_ram(freq, sd, rd, bathm, ssp, sed):
+def calZmax(ssp, sed):
+    d = np.array([])
+    for it in range(len(ssp)):
+        if sed == None:
+            d = np.append(d, ssp[it].z[-1])
+        else:
+            d = np.append(d, ssp[it].z[-1] + sed[it].z[-1])
+    Zmax = max(d)
+    return Zmax
+
+def call_ram(freq, sd, bathm, ssp, sed, dr=0.1):
     ''' ********************************************************************************* '''
     ''' ************************** Sound Field Calculation ****************************** '''
     ''' ********************************************************************************* '''
@@ -13,8 +23,10 @@ def call_ram(freq, sd, rd, bathm, ssp, sed):
     rmax = max(bathm.r)             # Maximum receiving range (unit: m)
     zmax = round(max(bathm.d) * 2)  # Maximum receiving depth (uint: m)
     zmplt = round(max(bathm.d)*1.1)     # Maximum depth of plot (uint: m)
-    dz = 1500 / freq / 15           # Step in depth direction: wavelength / 15
-    dr = 10 * dz                     # Step in range direction
+    dz = 10          # Step in depth direction: wavelength / 15
+    dr = 1000* dr                     # Step in range direction
+    Zmax = calZmax(ssp, sed)
+    rd = Zmax
     pos = PosRAM(sd, rd, rmax, zmax, dr, dz, zmplt)
 
     # Sound speed profile for RAM
@@ -40,6 +52,7 @@ def call_ram(freq, sd, rd, bathm, ssp, sed):
     # Read tl.grid
     return readTL_RAM()
 
+
 # if __name__=="__main__":
     # ''' ********************************************************************************* '''
     # ''' ***************************** Input Parameters ********************************** '''

pyat/tlcore.py

@@ -6,8 +6,8 @@ from krakenc import call_krakenc
 from ram import call_ram
 import matplotlib.pyplot as plt
 
-def call_tl(database, freq, sourceDepth, sourceLat, sourceLon, receiveDepth,
- receiveLat=None, receiveLon=None, range=None, azi=None, modelType='kraken', isPlot=False):
+def call_tl(database, freq, sourceDepth, sourceLat, sourceLon,
+ receiveLat=None, receiveLon=None, range=None, azi=None, dr=0.1, modelType='kraken', isPlot=False):
     """ Input:
     freq - float, Frequency
     sourceDepth - float, Source depth
@@ -40,8 +40,7 @@ def call_tl(database, freq, sourceDepth, sourceLat, sourceLon, receiveDepth,
         return [], [], [], bathym
     if modelType == 'bellhop':
         sd = np.array([sourceDepth])           # Source depth
-        rd = np.linspace(0, receiveDepth, int(receiveDepth/10))   # Receiving depth
-        [Pos, TL] = call_Bellhop(freq, sd, rd, bathym, ssp, sed, base)
+        [Pos, TL] = call_Bellhop(freq, sd, bathym, ssp, sed, base, dr)
         # Plot
         if isPlot:
             plt.pcolor(Pos.r.range, Pos.r.depth, TL[0,0,:,:], vmin=40, vmax=130, cmap='jet_r', shading='auto')
@@ -53,8 +52,7 @@ def call_tl(database, freq, sourceDepth, sourceLat, sourceLon, receiveDepth,
         return Pos.r.range, Pos.r.depth, TL[0,0,:,:], bathym
     elif modelType == 'kraken':
         sd = np.array([sourceDepth])           # Source depth
-        rd = np.linspace(0, receiveDepth, int(receiveDepth/10))   # Receiving depth
-        [Pos, TL] = call_krakenc(freq, sd, rd, bathym, ssp, sed, base)
+        [Pos, TL] = call_krakenc(freq, sd, bathym, ssp, sed, base, dr)
         # Plot
         if isPlot:
             plt.pcolor(Pos.r.range, Pos.r.depth, TL[0,0,:,:], vmin=40, vmax=130, cmap='jet_r', shading='auto')
@@ -66,7 +64,7 @@ def call_tl(database, freq, sourceDepth, sourceLat, sourceLon, receiveDepth,
         return Pos.r.range, Pos.r.depth, TL[0,0,:,:], bathym
     elif modelType == 'ram':
         # 调用ram模块
-        [rd, rr, TL] = call_ram(freq, sourceDepth, receiveDepth, bathym, ssp, sed)
+        [rd, rr, TL] = call_ram(freq, sourceDepth, bathym, ssp, sed, dr)
         # Plot
         if isPlot:
             plt.pcolormesh(rr, rd, TL, vmin=60, vmax=130, cmap='jet_r', shading='auto')
@@ -78,4 +76,16 @@ def call_tl(database, freq, sourceDepth, sourceLat, sourceLon, receiveDepth,
         return rr, rd, TL, bathym
     else:
         return [], [], [], bathym
- 
+
+
+
+
+# timeIndex = '20190101'
+# # 设置时间段
+#
+# STPath = '/home/gan/下载/code'
+# ETOPO1Path = '/home/gan/下载/code/20220920/ETOPO1'
+# database = DataBase(ETOPO1Path, STPath)
+# database.set_timePeriod(timeIndex)
+# range, depth, tl, bathym = call_tl(database=database,freq=200, sourceDepth=60, sourceLat=19, sourceLon=127,
+#                                    receiveLat=17.5, receiveLon=129,  modelType='bellhop', isPlot=True)

tl.py

@@ -1,31 +1,36 @@
 # -*- coding: utf-8 -*-
 import warnings
+
 warnings.filterwarnings('ignore')
 import argparse, json
 from common import NoIndent, MyEncoder
 import os, sys
-sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))+'/pyat')
 
+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)) + '/pyat')
+
+
+# import time
 def _argparse():
     parser = argparse.ArgumentParser(description='tl client')
     parser.add_argument('-a', '--argfile', action='store', dest='argfile', required=True, help='name of arguments file')
     return parser.parse_args()
 
+
 def main():
     parser = _argparse()
     argfile = parser.argfile
     if argfile[-5:] != '.json':
         print('输入的参数文件名有误')
         return
-    with open(argfile,'r') as jsonFile:
+    with open(argfile, 'r') as jsonFile:
         paramData = json.load(jsonFile)
-        if (("ETOPO1Path" not in paramData) or ("STPath" not in paramData) or\
-            "freq" not in paramData) or ("sourceDepth" not in paramData) or\
-            ("receiveDepth" not in paramData) or ("sLat" not in paramData) or\
-            ("sLon" not in paramData) or ("rLat" not in paramData) or\
-            ("rLon" not in paramData) or ("timeType" not in paramData) or\
-            ("modelType" not in paramData) or ("png" not in paramData) or\
-            ("outParam" not in paramData):
+        if (("ETOPO1Path" not in paramData) or ("STPath" not in paramData) or \
+            "freq" not in paramData) or ("sourceDepth" not in paramData) or \
+                ("receiveDepth" not in paramData) or ("sLat" not in paramData) or \
+                ("sLon" not in paramData) or ("rLat" not in paramData) or \
+                ("rLon" not in paramData) or ("timeType" not in paramData) or \
+                ("modelType" not in paramData) or ("png" not in paramData) or \
+                ("outParam" not in paramData):
             print(f'{argfile}中的参数有误，请核实！')
             return
         ETOPO1Path = paramData["ETOPO1Path"]
@@ -41,11 +46,11 @@ def main():
         modelType = paramData["modelType"]
         png = paramData["png"]
         outParam = paramData["outParam"]
-    
+
         from DataBase import DataBase
         from tlcore import call_tl
         import numpy as np
-        
+
         # 定义时间段
         timeIndex = timeType
         # 创建类
@@ -53,10 +58,11 @@ def main():
         # 设置时间段
         database.set_timePeriod(timeIndex)
 
-        [ranges, depths, TL, bathym] = call_tl(database=database, freq=freq, sourceDepth=sourceDepth, sourceLat=sLat, sourceLon=sLon, receiveDepth=receiveDepth,
-                receiveLat=rLat, receiveLon=rLon, modelType=modelType, isPlot=png)
+        [ranges, depths, TL, bathym] = call_tl(database=database, freq=freq, sourceDepth=sourceDepth, sourceLat=sLat,
+                                               sourceLon=sLon,
+                                               receiveLat=rLat, receiveLon=rLon, modelType=modelType, isPlot=png)
 
-        # write data to file                
+        # write data to file
         # 保存地形数据
         bathym.r = np.around(bathym.r, decimals=0)
         bathym.d = np.around(bathym.d, decimals=0)
@@ -66,7 +72,8 @@ def main():
             "descript": "距离、深度",
             "data": NoIndent(list(zip(bathym.r, bathym.d)))
         }
-        with open(("" if outParam["outPath"]=="" else outParam["outPath"] + "/")+"tlOut_bathym.json", 'w') as jsonFile:
+        with open(("" if outParam["outPath"] == "" else outParam["outPath"] + "/") + "tlOut_bathym.json",
+                  'w') as jsonFile:
             json_data = json.dumps(bathymData, cls=MyEncoder, indent=2, ensure_ascii=False)
             jsonFile.write(json_data)
 
@@ -75,16 +82,21 @@ def main():
         depths = np.around(depths, decimals=0)
         ranges = np.around(ranges, decimals=0)
         TL = np.around(TL, decimals=2)
-        tlpoint = list(zip(np.tile(depths.reshape(len(depths),1),len(ranges)).flatten(), np.tile(ranges,(len(depths),1)).flatten(), TL.flatten()))
+        tlpoint = list(zip(np.tile(depths.reshape(len(depths), 1), len(ranges)).flatten(),
+                           np.tile(ranges, (len(depths), 1)).flatten(), TL.flatten()))
         tlData = {
             "param": paramData,
-            "count": len(ranges)*len(depths),
+            "count": len(ranges) * len(depths),
             "descript": "距离、深度、损失",
             "data": NoIndent(tlpoint)
         }
-        with open(("" if outParam["outPath"]=="" else outParam["outPath"] + "/")+"tlOut_tl.json", 'w') as jsonFile:
+        with open(("" if outParam["outPath"] == "" else outParam["outPath"] + "/") + "tlOut_tl.json", 'w') as jsonFile:
             json_data = json.dumps(tlData, cls=MyEncoder, indent=2, ensure_ascii=False)
             jsonFile.write(json_data)
-    
+
+
 if __name__ == '__main__':
-    main()
+    # start = time.perf_counter()
+    main()
+    # end = time.perf_counter()
+    # print('time = ',(end-start))