dr*1000\dz

更改了间隔dr。

RANDI.py

@@ -3,6 +3,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')
 
 def _argparse():
     parser = argparse.ArgumentParser(description='tl client')
@@ -34,7 +36,7 @@ def main():
         NThreads = int(paramData["NThreads"])
         outParam = paramData["outParam"]
     
-        from pyat.RANDIcore import call_RANDI
+        from RANDIcore import call_RANDI
         NL = call_RANDI(ETOPO1Path, STPath, outParam["outPath"], freq, receiveDepth,
                         lon, lat, timeType, NThreads)
         print(round(NL,2))

cs.py

@@ -3,6 +3,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')
 
 def _argparse():
     parser = argparse.ArgumentParser(description='cs client')
@@ -28,7 +30,7 @@ def main():
         c = float(paramData["c"])
         outParam = paramData["outParam"]
     
-        from pyat.cscore import np, dis_mats, pn_funcs, ag_funcs
+        from cscore import np, dis_mats, pn_funcs, ag_funcs
         
         
         # 输出：

cz.py

@@ -3,6 +3,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')
 
 def _argparse():
     parser = argparse.ArgumentParser(description='Average Noise Spectral Level')
@@ -16,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"]
@@ -40,14 +42,15 @@ def main():
         outParam = paramData["outParam"]
         if "outPath" not in outParam:
             outParam["outPath"] = ""
-    
-        from pyat.czcore import call_cz
+
+        from DataBase import DataBase
+        from czcore import call_cz
         import numpy as np
 
         # 定义时间段
         timeIndex = timeType
         # 创建类
-        database = DataBase(ETOPO1Path, STPath,datasource="")
+        database = DataBase(ETOPO1Path, STPath, dataSource="")
         # 设置时间段
         database.set_timePeriod(timeIndex)
 
@@ -57,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:
@@ -69,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)
@@ -81,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)
                 }
@@ -103,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 = {
@@ -121,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
@@ -133,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
@@ -143,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))

oanp.py

@@ -3,6 +3,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')
 
 def _argparse():
     parser = argparse.ArgumentParser(description='Average Noise Spectral Level')
@@ -29,7 +31,7 @@ def main():
         png = paramData["png"]
         outParam = paramData["outParam"]
     
-        from pyat.OANPd import call_oanp
+        from OANPd import call_oanp
         from matplotlib import pyplot as plt
         import numpy as np
         

profile.py

@@ -1,8 +1,9 @@
 # -*- coding: utf-8 -*-
 import warnings
 warnings.filterwarnings('ignore')
-
 import argparse, json
+import os, sys
+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))+'/pyat')
 
 def _argparse():
     parser = argparse.ArgumentParser(description='剖面分析工具, 根据标准化剖面，得到声道, 跃层分析结果')
@@ -39,9 +40,9 @@ def main():
             if "outPath" not in outParam:
                 outParam["outPath"] = ""
         
-            from pyat.DataBase import DataBase
+            from DataBase import DataBase
             import matplotlib.pyplot as plt
-            from pyat.ProfileTool import ProfileTool
+            import ProfileTool
             import numpy as np
             # 定义时间段
             timeIndex = timeType

pyat/AcousticPara.py

@@ -3,7 +3,7 @@ import math
 from scipy.fftpack import fft
 from matplotlib import pyplot as plt
 from scipy.interpolate import interp1d
-from pyat.DataBase import DataBase
+from DataBase import DataBase
 
 class BathyM:
     def __init__(self, ran=None, depth=None):

pyat/DataBase.py

@@ -1,7 +1,7 @@
 import numpy as np
 from scipy.interpolate import Akima1DInterpolator, interp1d
 from netCDF4 import Dataset
-from pyat.ProfileDataType import Profile
+from ProfileDataType import Profile
 
 '''
 ************************************************************
@@ -33,7 +33,7 @@ class DataBase():
     """
     获取数据库中的剖面信息
     """
-    def __init__(self,ETOPO1path,STpath,dataSource):
+    def __init__(self,ETOPO1path,STpath,dataSource=None):
         # 打开ETOPO1数据
         fileETOPO1 = ETOPO1path+'/ETOPO1_Bed_g_gmt4.nc'
         dataETOPO1 = Dataset(fileETOPO1)

pyat/DataBaseWOA.py

@@ -0,0 +1,370 @@
+import numpy as np
+from scipy.interpolate import Akima1DInterpolator
+from scipy.optimize import curve_fit
+from netCDF4 import Dataset
+from ProfileDataType import Profile
+
+'''
+************************************************************
+该模块用于从数据库中, 获取剖面信息
+************************************************************
+'''
+
+# 时间段标识
+timeDict = {
+    # 全年平均
+    'Annual': '00',
+    # 1-12月
+    'January': '01',
+    'February': '02',
+    'March': '03',
+    'April': '04',
+    'May': '05',
+    'June': '06',
+    'July': '07',
+    'August': '08',
+    'September': '09',
+    'October': '10',
+    'November': '11',
+    'December': '12',
+    # 北半球四季
+    'Winter': '13',
+    'Spring': '14',
+    'Summer': '15',
+    'Autumn': '16'
+}
+
+
+# def MunkFormula(z, v1, epsilon, z1):
+#     eta = 2*(z-z1)/z1
+#     return v1*(1+epsilon*(np.exp(-eta)+eta-1))
+
+def soundSpeed(tem, sal, dep):
+    '''
+    声速计算公式
+    '''
+    # C = 1449.2 + 4.6*T - 0.055*T^2 + 0.00029*T^3 + (1.34 - 0.01*T)(S-35) + 0.017*D
+    # 式中：T为海⽔测点的温度，单位为℃；S为海⽔测点的盐度，单位为‰；D为海⽔测点的深度，单位为m。
+    vel = 1449.2 + 4.6*tem - 0.055*tem*tem + 0.00029 * \
+        tem*tem*tem + (1.34 - 0.01*tem) * \
+        (sal-35.0) + 0.017*dep
+    return vel
+
+
+class DataBase():
+    """
+    获取数据库中的剖面信息
+    """
+
+    def __init__(self,ETOPO1Path,STpath):
+        print(ETOPO1Path)
+        # 打开ETOPO1数据
+        fileETOPO1 = ETOPO1Path + '/ETOPO1_Bed_g_gmt4.nc'
+        dataETOPO1 = Dataset(fileETOPO1)
+        # 获取ETOPO1网格
+        self.latGrid1 = dataETOPO1.variables['y'][:]
+        self.lonGrid1 = dataETOPO1.variables['x'][:]
+        # 获取海深数据
+        self.botGrid = -dataETOPO1.variables['z'][:]
+        self.STPath=STpath
+        return
+
+    def clearVar(self):
+        """
+        清空变量
+        """
+        self.botDep = None
+        self.salRawPf = None
+        self.temRawPf = None
+        self.velRawPf = None
+        self.velInterpPf = None
+        self.velExtPf = None
+        self.velFullPf = None
+        self.isValid = True
+        return
+
+    def set_timePeriod(self, timeIndex):
+        """
+        设置时间段
+        """
+        timeIndex = '01'
+        self.timeIndex = timeIndex
+        # 打开WOA18数据
+        fileWOA18sal = self.STPath + '/woa18_A5B7_s' + timeIndex + '_04.nc'
+        fileWOA18tem = self.STPath + '/woa18_A5B7_t' + timeIndex + '_04.nc'
+        dataWOA18Sal = Dataset(fileWOA18sal)
+        dataWOA18Tem = Dataset(fileWOA18tem)
+        # 获取WOA18网格
+        self.latGrid2 = dataWOA18Sal.variables['lat'][:]
+        self.lonGrid2 = dataWOA18Sal.variables['lon'][:]
+        # 获取深度数据
+        self.depthGrid = dataWOA18Sal.variables['depth'][:]
+        # 获取盐度数据
+        self.salGrid = np.squeeze(dataWOA18Sal.variables['s_an'][:])
+        # 获取温度数据
+        self.temGrid = np.squeeze(
+            dataWOA18Tem.variables['t_an'][:])
+        self.clearVar()
+        return
+
+    def set_coord(self, lat, lon):
+        """
+        设置坐标
+        """
+        self.lat = lat
+        self.lon = lon
+        self.clearVar()
+        return
+
+    def interpBilinear(self, latId, lonId, dataGrid, dataTpye):
+        """
+        双线性插值
+        """
+        latId0 = int(np.ceil(latId))
+        lonId0 = int(np.ceil(lonId))
+        a = latId - latId0
+        b = lonId - lonId0
+        if (dataTpye == 1):
+            data0 = dataGrid[latId0, lonId0]
+            data1 = dataGrid[latId0+1, lonId0]
+            data2 = dataGrid[latId0+1, lonId0+1]
+            data3 = dataGrid[latId0, lonId0+1]
+        else:
+            data0 = dataGrid[:, latId0, lonId0]
+            data1 = dataGrid[:, latId0+1, lonId0]
+            data2 = dataGrid[:, latId0+1, lonId0+1]
+            data3 = dataGrid[:, latId0, lonId0+1]
+        data01 = (1-a)*data0+a*data1
+        data32 = (1-a)*data3+a*data2
+        data = (1-b)*data01+b*data32
+        return data
+
+    def normCoord(self, latGrid, lonGrid):
+        """
+        坐标单位化
+        """
+        lat0 = latGrid[0]
+        lon0 = lonGrid[0]
+        dlat = latGrid[1]-latGrid[0]
+        dlon = lonGrid[1]-lonGrid[0]
+        latId = (self.lat - lat0)/dlat
+        lonId = (self.lon - lon0)/dlon
+        return (latId, lonId)
+
+    def get_bottomDepth(self):
+        """
+        获得海底深度
+        """
+        if (self.botDep == None):
+            # 单位化经纬度
+            (latId, lonId) = self.normCoord(self.latGrid1, self.lonGrid1)
+            # 获取网格点海深
+            self.botDep = self.interpBilinear(latId, lonId, self.botGrid, 1)
+        return self.botDep
+
+    def get_salRawProfile(self):
+        """
+        获得原始盐度剖面
+        """
+        if (self.salRawPf == None):
+            botDep = self.get_bottomDepth()
+            # 单位化经纬度
+            (latId, lonId) = self.normCoord(self.latGrid2, self.lonGrid2)
+            # 获取盐度剖面信息
+            sal = self.interpBilinear(latId, lonId, self.salGrid, 2)
+            # 获取有效数据段
+            dep = self.depthGrid[sal.mask == False]
+            sal = sal[sal.mask == False]
+            dep = dep.flatten()
+            sal = sal.flatten()
+            # 若为空，则为无效数据
+            if sal.size <= 1:
+                self.isValid = False
+            # 剖面赋值
+            salRawPf = Profile()
+            salRawPf.set_data(dep, sal)
+            salRawPf.set_coord(self.lat, self.lon)
+            salRawPf.set_bottomDepth(botDep)
+            self.salRawPf = salRawPf
+        return self.salRawPf
+
+    def get_temRawProfile(self):
+        """
+        获得原始温度剖面
+        """
+        if (self.temRawPf == None):
+            botDep = self.get_bottomDepth()
+            # 单位化经纬度
+            (latId, lonId) = self.normCoord(self.latGrid2, self.lonGrid2)
+            # 获取温度剖面信息
+            tem = self.interpBilinear(latId, lonId, self.temGrid, 2)
+            # 获取有效数据段
+            dep = self.depthGrid[tem.mask == False]
+            tem = tem[tem.mask == False]
+            dep = dep.flatten()
+            tem = tem.flatten()
+            # 若为空，则为无效数据
+            if tem.size <= 1:
+                self.isValid = False
+            # 剖面赋值
+            temRawPf = Profile()
+            temRawPf.set_data(dep, tem)
+            temRawPf.set_coord(self.lat, self.lon)
+            temRawPf.set_bottomDepth(botDep)
+            self.temRawPf = temRawPf
+        return self.temRawPf
+
+    def get_velRawProfile(self):
+        """
+        获得原始声速剖面
+        """
+        if (self.velRawPf == None):
+            botDep = self.get_bottomDepth()
+            # 获取温盐深数据
+            dep = self.get_salRawProfile().x
+            sal = self.get_salRawProfile().y
+            tem = self.get_temRawProfile().y
+            # 计算声速
+            if self.isValid == True:
+                vel = soundSpeed(tem, sal, dep)
+            else:
+                dep = np.array([0])
+                vel = np.array([0])
+            # 剖面赋值
+            velRawPf = Profile()
+            velRawPf.set_data(dep, vel)
+            velRawPf.set_coord(self.lat, self.lon)
+            velRawPf.set_bottomDepth(botDep)
+            self.velRawPf = velRawPf
+        return self.velRawPf
+
+    def get_velInterpProfile(self):
+        """
+        获得插值后的声速剖面
+        将非等距离采样填充为深度间隔为1m的序列, 并应用Akima插值算法, 对数据进行插值
+        """
+        if (self.velInterpPf == None):
+            botDep = self.get_bottomDepth()
+            # 获得原始声速剖面数据
+            dep = self.get_velRawProfile().x
+            vel = self.get_velRawProfile().y
+            if self.isValid == True:
+                # 对剖面进行垂向插值处理，将非等距离采样填充为深度间隔为1m的序列
+                depInterp = np.arange(0, dep[-1]+1, 1)
+                # 应用Akima插值算法，对数据进行插值
+                func = Akima1DInterpolator(dep, vel)
+                velInterp = func(depInterp)
+            else:
+                depInterp = np.array([0])
+                velInterp = np.array([0])
+            # 剖面赋值
+            velInterpPf = Profile()
+            velInterpPf.set_data(depInterp, velInterp)
+            velInterpPf.set_coord(self.lat, self.lon)
+            velInterpPf.set_bottomDepth(botDep)
+            self.velInterpPf = velInterpPf
+        return self.velInterpPf
+
+    def get_velExtProfile(self):
+        """
+        获得延拓后的声速剖面
+        """
+        if (self.velExtPf == None):
+            botDep = self.get_bottomDepth()
+            if self.isValid == True:
+                # 获得原始声速剖面数据
+                dep = self.get_velRawProfile().x
+                vel = self.get_velRawProfile().y
+                # 获得剖面最低处温盐值
+                tem = self.get_temRawProfile().y[-1]
+                sal = self.get_salRawProfile().y[-1]
+                # 延拓部分深度序列
+                depExt = np.arange(dep[-1]+1, botDep+1)
+                # 按照声速公式延拓，假设温盐不变，仅深度变化
+                velExt = soundSpeed(tem, sal, depExt)
+            else:
+                depExt = np.array([0])
+                velExt = np.array([0])
+            # 剖面赋值
+            velExtPf = Profile()
+            velExtPf.set_data(depExt, velExt)
+            velExtPf.set_coord(self.lat, self.lon)
+            velExtPf.set_bottomDepth(botDep)
+            self.velExtPf = velExtPf
+        return self.velExtPf
+
+    def get_velFullProfile(self):
+        '''
+        获得全海深声速剖面
+        '''
+        if (self.velFullPf == None):
+            botDep = self.get_bottomDepth()
+            if self.isValid == True:
+                # 插值后的速度剖面
+                depInterp = self.get_velInterpProfile().x
+                velInterp = self.get_velInterpProfile().y
+                # 延拓后的速度剖面
+                depExt = self.get_velExtProfile().x
+                velExt = self.get_velExtProfile().y
+                # 全海深声速剖面
+                depFull = np.concatenate((depInterp, depExt))
+                velFull = np.concatenate((velInterp, velExt))
+            else:
+                depFull = np.array([0])
+                velFull = np.array([0])
+            # 剖面赋值
+            velFullPf = Profile()
+            velFullPf.set_data(depFull, velFull)
+            velFullPf.set_coord(self.lat, self.lon)
+            velFullPf.set_bottomDepth(botDep)
+            self.velFullPf = velFullPf
+        return self.velFullPf
+
+    def export_velFullProfile_txt(self, filePath):
+        '''
+        速度剖面导出txt文件
+        '''
+        velPf = self.get_velFullProfile()
+        lat = velPf.lat
+        lon = velPf.lon
+        botDep = velPf.botDep
+        dep = velPf.x
+        vel = velPf.y
+        # 导出文本
+        line1 = 'lat[degree north], {0:>7.4f}, lon[degree east], {1:>7.4f}, BottomDepth[m], {2:>5.0f}\n'.format(
+            lat, lon, botDep)
+        line2 = 'depth[m], velocity[m/s]'
+        np.savetxt(filePath, np.stack((dep, vel), -1),
+                   fmt='%.5e', delimiter=',', comments='', header=line1+line2)
+
+    def export_velFullProfile_nc(self, filePath):
+        '''
+        速度剖面导出nc文件
+        '''
+        velPf = self.get_velFullProfile()
+        lat = velPf.lat
+        lon = velPf.lon
+        botDep = velPf.botDep
+        dep = velPf.x
+        vel = velPf.y
+        
+        data = Dataset(filePath, 'w', format='NETCDF4')
+        data.createDimension('lon', 1)
+        data.createDimension('lat', 1)
+        data.createDimension('botDep', 1)
+        data.createDimension('dep', dep.size)
+        data.createDimension('vel', vel.size)
+
+        data.createVariable('lon', 'f', ('lon'))
+        data.createVariable('lat', 'f', ('lat'))
+        data.createVariable('botDep', 'f', ('botDep'))
+        data.createVariable('dep', 'f', ('dep'))
+        data.createVariable('vel', 'f', ('vel'))
+
+        data.variables['lat'][:] = lat
+        data.variables['lon'][:] = lon
+        data.variables['botDep'][:] = botDep
+        data.variables['dep'][:] = dep
+        data.variables['vel'][:] = vel
+
+        data.close()

pyat/ProfileRead.py

@@ -1,9 +1,9 @@
 import unicodedata
 import numpy as np
-from pyat.ProfileDataType import Profile
+from ProfileDataType import Profile
 from scipy.interpolate import interp1d
 from scipy.signal import butter, filtfilt
-from pyat.DataBase import DataBase
+from DataBase import DataBase
 
 '''
 ************************************************************

pyat/ProfileTool.py

@@ -1,6 +1,6 @@
 import numpy as np
 from scipy.signal import find_peaks
-from pyat.ProfileDataType import Channel, SpringLayer
+from ProfileDataType import Channel, SpringLayer
 
 '''
 ************************************************************

pyat/RANDIcore.py

@@ -1,13 +1,13 @@
 import os
 from multiprocessing import Pool
 from functools import partial
-from time import time
+import time
 from shutil import rmtree
 import numpy as np
 from scipy.interpolate import interp1d
 from copy import deepcopy
-from pyat import ramRW_new as ramRW
-from pyat.DataBase import DataBase,soundSpeed
+import ramRW_new as ramRW
+from DataBase import DataBase
 
 
 class BathyM:
@@ -124,7 +124,7 @@ def RANDI_GetEnv(data, coordS, coordR=None, R=None, azi=None):
         z_ssp = np.append(z_ssp, max(bathym.d))
         sal_m = np.append(sal_m, sal_m[-1])
         tem_m = np.append(tem_m, tem_m[-1])
-    c_ssp = soundSpeed(tem_m, sal_m, z_ssp)
+    c_ssp = DataBase.soundSpeed(tem_m, sal_m, z_ssp)
     ssp = SSP_(0, z_ssp, c_ssp)
     # plt.plot(c_ssp, z_ssp)
     # plt.gca().invert_yaxis()

pyat/WGNPd.py

@@ -1,7 +1,7 @@
 """ Wind Generated Noise Prediction """
 import numpy as np
 from scipy.interpolate import interp1d
-from pyat import DataBase as db
+import DataBase as db
 import matplotlib.pyplot as plt
 
 class MdlPara():

pyat/bellhop.py

@@ -1,8 +1,11 @@
-from pyat.readwrite import write_env, read_shd, get_rays
-from pyat.env import Pos, Source, Dom, cInt, SSPraw, SSP, HS, BotBndry, TopBndry, Bndry, Box, Beam
+from readwrite import write_env, read_shd, get_rays
+from env import Pos, Source, Dom, cInt, SSPraw, SSP, HS, BotBndry, TopBndry, Bndry, Box, Beam
 from os import system
 import numpy as np
-from pyat.config import AtBinPath
+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/config.json

@@ -1,3 +1,3 @@
 {
-    "AtBinPath": "/home/shunli/pro/pyat/at/bin"
+    "AtBinPath": "/home/shunli/pro/at/bin"
 }

pyat/czcore.py

@@ -3,15 +3,13 @@ from multiprocessing import Process, Pool
 from functools import partial
 from os import system
 import numpy as np
-from pyat.AcousticPara import GetEnv, Coord
+from AcousticPara import GetEnv, Coord
 from scipy.signal import savgol_filter, medfilt
-from pyat.DataBase import DataBase
-from pyat.readwrite import write_env, write_bathy, read_shd, get_rays
-from pyat.env import Pos, Source, Dom, cInt, SSPraw, SSP, HS, BotBndry, TopBndry, Bndry, Box, Beam
+from DataBase import DataBase
+from readwrite import write_env, write_bathy, read_shd, get_rays
+from env import Pos, Source, Dom, cInt, SSPraw, SSP, HS, BotBndry, TopBndry, Bndry, Box, Beam
 import matplotlib.pyplot as plt
-from pyat.config import AtBinPath
-
-file_cz_data = open('cz_data.txt', 'w')
+from config import AtBinPath
 
 
 class CzRangeWidth:
@@ -48,11 +46,11 @@ def calZmax(ssp, sed):
     return Zmax, ssp_idx
 
 
-def call_Bellhop(freq, sd, rd, 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, rd, 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, rd, bathm, ssp, sed, base,azi):
     return Pos1, TL
 
 
-def call_Bellhop_Rays(freq, sd, rd, 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, rd, 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, rd, 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, rd, 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, rd, 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,22 +285,13 @@ def call_fangwei(freq, sd, rd, 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, receiveDepth=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
-    rd = np.linspace(0, receiveDepth, 351)  # Receiving depth
 
     results = []
     ranges, depths, TL, bathym, Rays = [], [], [], [], []
@@ -322,30 +301,35 @@ def call_cz(database, freq, lat, lon, azi, sourceDepth=10, receiveDepth=10, isPl
     #     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

@@ -1,8 +1,8 @@
 from os import system
 import numpy as np
-from pyat.readwrite import write_env, write_fieldflp, read_shd
-from pyat.env import Pos, Source, Dom, cInt, SSPraw, SSP, HS, BotBndry, TopBndry, Bndry
-from pyat.config import AtBinPath
+from readwrite import write_env, write_fieldflp, read_shd
+from env import Pos, Source, Dom, cInt, SSPraw, SSP, HS, BotBndry, TopBndry, Bndry
+from config import AtBinPath
 
 def calZmax(ssp, sed):
     d = np.array([])
@@ -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

@@ -1,9 +1,19 @@
-from pyat.ramRW import PosRAM, PROF, SSPRAM, write_env_RAM, readTL_RAM
+from ramRW import PosRAM, PROF, SSPRAM, write_env_RAM, readTL_RAM
 from os import system
 import numpy as np
-from pyat.config import AtBinPath
+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/ramRW.py

@@ -1,4 +1,4 @@
-from pyat.env import *
+from env import *
 import numpy as np
 
 def write_env_RAM(TitleEnv, freq, pos, bathym, ssp, envfil = "ram.in"):

pyat/ramRW_new.py

@@ -1,4 +1,4 @@
-from pyat.env import *
+from env import *
 import numpy as np
 
 

pyat/readwrite.py

@@ -2,7 +2,7 @@ import numpy as np
 from scipy.io import loadmat
 import os
 from struct import unpack
-from pyat.env import Source, Dom, Pos, cInt, Ice, SSPraw, SSP, HS, BotBndry, TopBndry, Bndry, Box, Beam, Modes, Eigenray
+from env import Source, Dom, Pos, cInt, Ice, SSPraw, SSP, HS, BotBndry, TopBndry, Bndry, Box, Beam, Modes, Eigenray
 import re
 from matplotlib import pyplot as plt
 from math import floor

pyat/tlcore.py

@@ -1,13 +1,13 @@
 import numpy as np
-from pyat.AcousticPara import GetEnv, Coord
-from pyat.DataBase import DataBase
-from pyat.bellhop import call_Bellhop
-from pyat.krakenc import call_krakenc
-from pyat.ram import call_ram
+from AcousticPara import GetEnv, Coord
+from DataBase import DataBase
+from bellhop import call_Bellhop
+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)

rl.py

@@ -2,6 +2,8 @@
 import warnings
 warnings.filterwarnings('ignore')
 import argparse, json
+import os, sys
+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))+'/pyat')
 
 def _argparse():
     parser = argparse.ArgumentParser(description='混响级计算')
@@ -33,7 +35,7 @@ def main():
         r = float(paramData["range"])
         outParam = paramData["outParam"]
     
-        from pyat.rlcore import call_rl
+        from rlcore import call_rl
         [rlss, rlsb, rlsv] = call_rl(v,SL,Sv,freq,tau,theta,l,r)
         # [rlss,rlsb] = call_rl(5,170,55,10,9.5,0.3,182.5)
 

tl.py

@@ -1,29 +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')
+
+
+# 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"]
@@ -39,10 +46,11 @@ def main():
         modelType = paramData["modelType"]
         png = paramData["png"]
         outParam = paramData["outParam"]
-    
-        from pyat.tlcore import call_tl
+
+        from DataBase import DataBase
+        from tlcore import call_tl
         import numpy as np
-        
+
         # 定义时间段
         timeIndex = timeType
         # 创建类
@@ -50,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)
@@ -63,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)
 
@@ -72,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))

tl_realtime.py

@@ -3,6 +3,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')
 
 def _argparse():
     parser = argparse.ArgumentParser(description='tl realtime client')
@@ -61,8 +63,8 @@ def main():
 
         from scipy.interpolate import Akima1DInterpolator
         import numpy as np
-        from pyat.AcousticPara import BathyM, Base, Sed, SSP_
-        from pyat.tlcore import call_Bellhop, call_krakenc, call_ram
+        from AcousticPara import BathyM, Base, Sed, SSP_
+        from tlcore import call_Bellhop, call_krakenc, call_ram
         import matplotlib.pyplot as plt
 
         bathym, ssp, sed, base = [], [], [], []

wgnp.py

@@ -2,6 +2,8 @@
 import warnings
 warnings.filterwarnings('ignore')
 import argparse, json
+import os, sys
+sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))+'/pyat')
 
 def _argparse():
     parser = argparse.ArgumentParser(description='Wind Generated Noise Prediction')
@@ -32,8 +34,7 @@ def main():
         png = paramData["png"]
         outParam = paramData["outParam"]
     
-        from pyat.WGNPd import WGNPd
-        from pyat.DataBase import DataBase
+        from WGNPd import WGNPd
         import matplotlib.pyplot as plt
         import numpy as np