diff --git a/.gitignore b/.gitignore index 5d381cc..ab8aff2 100644 --- a/.gitignore +++ b/.gitignore @@ -160,3 +160,8 @@ cython_debug/ # option (not recommended) you can uncomment the following to ignore the entire idea folder. #.idea/ +.vscode/ +packages/ +pyat/at/ +pyqt/ +test/ \ No newline at end of file diff --git a/RANDI.py b/RANDI.py new file mode 100644 index 0000000..fa44df4 --- /dev/null +++ b/RANDI.py @@ -0,0 +1,51 @@ +# -*- coding: utf-8 -*- +import warnings +warnings.filterwarnings('ignore') +import argparse, json +from common import NoIndent, MyEncoder + +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: + paramData = json.load(jsonFile) + if ("ETOPO1Path" not in paramData) or ("STPath" not in paramData) or\ + ("freq" not in paramData) or ("receiveDepth" not in paramData) or\ + ("lat" not in paramData) or ("lon" not in paramData) or\ + ("timeType" not in paramData) or\ + ("NThreads" not in paramData) or ("outParam" not in paramData): + print(f'{argfile}中的参数有误，请核实！') + return + ETOPO1Path = paramData["ETOPO1Path"] + STPath = paramData["STPath"] + freq = float(paramData["freq"]) + receiveDepth = float(paramData["receiveDepth"]) + lat = float(paramData["lat"]) + lon = float(paramData["lon"]) + timeType = paramData["timeType"] + NThreads = int(paramData["NThreads"]) + outParam = paramData["outParam"] + + from pyat.RANDIcore import call_RANDI + NL = call_RANDI(ETOPO1Path, STPath, outParam["outPath"], freq, receiveDepth, + lon, lat, timeType, NThreads) + print(round(NL,2)) + outData = { + "param": paramData, + "noiseLevel": round(NL,2) + } + with open(("" if outParam["outPath"]=="" else outParam["outPath"] + "/")+"RANDIOut.json", 'w') as jsonFile: + json_data = json.dumps(outData, cls=MyEncoder, indent=2, ensure_ascii=False) + jsonFile.write(json_data) + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/common.py b/common.py new file mode 100644 index 0000000..612aac9 --- /dev/null +++ b/common.py @@ -0,0 +1,32 @@ + +# -*- coding: utf-8 -*- +from _ctypes import PyObj_FromPtr +import json, re + +class NoIndent(object): + """ Value wrapper. """ + def __init__(self, value): + self.value = value + +class MyEncoder(json.JSONEncoder): + FORMAT_SPEC = '@@{}@@' + regex = re.compile(FORMAT_SPEC.format(r'(\d+)')) + + def __init__(self, **kwargs): + super(MyEncoder, self).__init__(**kwargs) + + def default(self, obj): + return (self.FORMAT_SPEC.format(id(obj)) if isinstance(obj, NoIndent) + else super(MyEncoder, self).default(obj)) + + def encode(self, obj): + format_spec = self.FORMAT_SPEC # Local var to expedite access. + json_repr = super(MyEncoder, self).encode(obj) # Default JSON. + for match in self.regex.finditer(json_repr): + id = int(match.group(1)) + no_indent = PyObj_FromPtr(id) + json_obj_repr = json.dumps(no_indent.value) + json_repr = json_repr.replace( + '"{}"'.format(format_spec.format(id)), json_obj_repr) + + return json_repr \ No newline at end of file diff --git a/cs.py b/cs.py new file mode 100644 index 0000000..2711cf2 --- /dev/null +++ b/cs.py @@ -0,0 +1,88 @@ +# -*- coding: utf-8 -*- +import warnings +warnings.filterwarnings('ignore') +import argparse, json +from common import NoIndent, MyEncoder + +def _argparse(): + parser = argparse.ArgumentParser(description='cs 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: + paramData = json.load(jsonFile) + if ("m" not in paramData) or ("d" not in paramData) or\ + ("f" not in paramData) or ("c" not in paramData) or\ + ("outParam" not in paramData): + print(f'{argfile}中的参数有误，请核实！') + return + m = int(paramData["m"]) + d = float(paramData["d"]) + f = float(paramData["f"]) + c = float(paramData["c"]) + outParam = paramData["outParam"] + + from pyat.cscore import np, dis_mats, pn_funcs, ag_funcs + import sys + print(sys.path) + + + # 输出： + # 所有阵元的间距的集合 + dis_m = dis_mats(m, d) + print('所有阵元的间距的集合：', dis_m) + + # 相邻阵元的噪声相关性 + pn_res = pn_funcs(f, c, dis_m) + print('相邻阵元的噪声相关性：pn=', pn_res[0, 1]) + + # 阵增益 + ag_res = ag_funcs(m, pn_res) + print('阵增益：ag=', ag_res) + outData = { + "param": paramData, + "dis": np.round(dis_m,2).tolist(), + "pn": np.round(pn_res,2).tolist(), + "ag": np.round(ag_res,2).tolist() + } + with open(("" if outParam["outPath"]=="" else outParam["outPath"] + "/")+"csOut.json", 'w') as jsonFile: + json_data = json.dumps(outData, cls=MyEncoder, indent=2, ensure_ascii=False) + jsonFile.write(json_data) + # # write data to file + # # 保存地形数据 + # bathym.r = np.around(bathym.r, decimals=0) + # bathym.d = np.around(bathym.d, decimals=0) + # bathymData = { + # "param": paramData, + # "count": len(bathym.d), + # "descript": "距离、深度", + # "data": NoIndent(list(zip(bathym.r, bathym.d))) + # } + # 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) + + # # 保存传播损失数据 + # tlpoint = [] + # 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())) + # tlData = { + # "param": paramData, + # "count": len(ranges)*len(depths), + # "descript": "距离、深度、损失", + # "data": NoIndent(tlpoint) + # } + # 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() \ No newline at end of file diff --git a/cz.py b/cz.py new file mode 100644 index 0000000..753f88d --- /dev/null +++ b/cz.py @@ -0,0 +1,154 @@ +# -*- coding: utf-8 -*- +import warnings +warnings.filterwarnings('ignore') +import argparse, json +from common import NoIndent, MyEncoder + +def _argparse(): + parser = argparse.ArgumentParser(description='Average Noise Spectral Level') + 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: + 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): + print(f'{argfile}中的参数有误，请核实！') + return + ETOPO1Path = paramData["ETOPO1Path"] + STPath = paramData["STPath"] + freq = float(paramData["freq"]) + sourceDepth = float(paramData["sourceDepth"]) + receiveDepth = float(paramData["receiveDepth"]) + lat = float(paramData["lat"]) + lon = float(paramData["lon"]) + try: + azi = float(paramData["azi"]) + except: + azi = None + timeType = paramData["timeType"] + outParam = paramData["outParam"] + if "outPath" not in outParam: + outParam["outPath"] = "" + + from pyat.czcore import call_cz + import numpy as np + + # call_cz函数返回汇聚区计算的中间值CzData的列表 + # CzData保存了方位角、距离、深度、损失、地形、会聚区等数据 + # 其中czdata.azi是方位角、czdata.ranges是距离列表、czdata.depths是深度列表、 + # czdata.TL是损失、czdata.bathym是地形数据 + # czdata.Rays是声线数据 + # czdata.CzDetail是会聚区数据，包括会聚区宽度和距离 + res = call_cz(ETOPO1Path, STPath, freq, lat, lon, azi, sourceDepth, receiveDepth, timeType=timeType) + # call_cz(25, 18.0, 114.0, 10, 3500) + + # write data to file + #输出json文件 + outData = {} + outData["param"] = paramData + if len(res) == 0: + outData["result"] = "无会聚区" + outData["data"] = [] + else: + outData["result"] = "有会聚区" + outDataData, outBathymData,outTlData,outRayData = [],[],[],[] + json.encoder.FLOAT_REPR = lambda x: format(x, '.3f') + for index in range(len(res)): + 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) + bathymData = { + "param": paramData, + "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))) + } + # 保存传播损失数据 + 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())) + tlData = { + "param": paramData, + "azi": f'{round(res[index].azi)}', + "count": len(res[index].ranges)*len(res[index].depths), + "descript": "距离、深度、损失", + "data": NoIndent(tlpoint) + } + # 保存声线数据 + raypoint = [] + for ray in res[index].Rays: + ray.xy[0] = np.around(ray.xy[0], decimals=2) + 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]))) + } + raypoint.append(tempdata) + rayData = { + "param": paramData, + "azi": f'{round(res[index].azi)}', + "count": len(res[index].Rays), + "rayData": raypoint + } + # print(len(res[index].Rays)) + czDeatilData = {} + czDeatilData["descript"] = res[index].CzDetail.descript + if len(res[index].CzDetail.directCz) == 0: + czDeatilData["directCz"] = [] + else: + 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)}', + } + directczdetaillist.append(directczdetail) + czDeatilData["directCz"] = directczdetaillist + if len(res[index].CzDetail.reflectCz) == 0: + czDeatilData["reflectCz"] = [] + else: + 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)}', + } + reflectczdetaillist.append(reflectczdetail) + czDeatilData["reflectCz"] = reflectczdetaillist + outDataData.append(czDeatilData) + outBathymData.append(bathymData) + outTlData.append(tlData) + outRayData.append(rayData) + 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: + 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: + 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: + json.dump(outData, jsonFile, indent=4, ensure_ascii=False) + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/oanp.py b/oanp.py new file mode 100644 index 0000000..4b1e1f4 --- /dev/null +++ b/oanp.py @@ -0,0 +1,59 @@ +# -*- coding: utf-8 -*- +import warnings +warnings.filterwarnings('ignore') +import argparse, json +from common import NoIndent, MyEncoder + +def _argparse(): + parser = argparse.ArgumentParser(description='Average Noise Spectral Level') + 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: + paramData = json.load(jsonFile) + if ("soundFile" not in paramData) or ("sense" not in paramData) or\ + ("filter" not in paramData) or ("png" not in paramData) or\ + ("outParam" not in paramData): + print(f'{argfile}中的参数有误，请核实！') + return + soundFile = paramData["soundFile"] + M = float(paramData["sense"]) + A = float(paramData["filter"]) + png = paramData["png"] + outParam = paramData["outParam"] + + from pyat.OANPd import call_oanp + from matplotlib import pyplot as plt + import numpy as np + + [freq, spl, fC, spl_oct] = call_oanp(soundFile=soundFile, M=M , A=A) + # call_oanp(soundFile="./SAW2577_20200923_150000.wav", W=-164.0 , A=0, isPlot=isPlot) + + #输出json文件 + outData = {} + outData["param"] = paramData + outData["descript"] = "frequency and noise level" + outData["count"] = len(fC) + fC = np.around(fC, decimals=0) + spl_oct = np.around(spl_oct, decimals=2) + outData["data"] = NoIndent(list(zip(fC,spl_oct))) + with open(("" if outParam["outPath"]=="" else outParam["outPath"] + "/")+"oanpOut.json", 'w') as jsonFile: + json_data = json.dumps(outData, cls=MyEncoder, indent=2, ensure_ascii=False) + jsonFile.write(json_data) + + if png: + plt.semilogx(freq[20:10000], spl[20:10000]) + plt.semilogx(fC, spl_oct) + plt.xlabel("Frequency / Hz") + plt.ylabel("Noise Level / dB") + plt.show() + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/profile.py b/profile.py new file mode 100755 index 0000000..8bdc357 --- /dev/null +++ b/profile.py @@ -0,0 +1,220 @@ +# -*- coding: utf-8 -*- +import warnings +warnings.filterwarnings('ignore') + +import argparse, json + +def _argparse(): + parser = argparse.ArgumentParser(description='剖面分析工具, 根据标准化剖面，得到声道, 跃层分析结果') + 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: + paramData = json.load(jsonFile) + if ("realtime" not in paramData): + print(f'{argfile}中的参数有误，请核实！') + return + realtime = paramData["realtime"] + if not realtime: + if ("ETOPO1Path" not in paramData) or ("STPath" not in paramData) or\ + ("lat" not in paramData) or ("lon" not in paramData) or\ + ("timeType" not in paramData) or ("png" not in paramData) or\ + ("outParam" not in paramData): + print(f'{argfile}中的参数有误，请核实！') + return + ETOPO1Path = paramData["ETOPO1Path"] + STPath = paramData["STPath"] + lat = float(paramData["lat"]) + lon = float(paramData["lon"]) + timeType = paramData["timeType"] + png = paramData["png"] + outParam = paramData["outParam"] + if "outPath" not in outParam: + outParam["outPath"] = "" + + from pyat.DataBase import DataBase + import matplotlib.pyplot as plt + from pyat.ProfileTool import ProfileTool + import numpy as np + # 定义时间段 + timeIndex = timeType + # 创建类 + data = DataBase(ETOPO1Path, STPath) + # 设置时间段 + data.set_timePeriod(timeIndex) + # 设置坐标 + data.set_coord(lat, lon) + # 全海深速度剖面 + velPf = data.get_velFullProfile() + # 设置剖面分析工具 + velPfTool = ProfileTool() + velPfTool.set_profile(velPf) + # 计算表面声道 + surfCh = velPfTool.get_surfaceChannel() + # 计算深海声道 + deepCh = velPfTool.get_deepChannel() + # 计算跃层信息 + surfSL, mainSL = velPfTool.get_springLayer() + # 导出结果文件 + velPfTool.export_profileParm('velPfParam.txt') + # 导出剖面文件 + # data.export_velFullProfile_txt('velFullProfile.txt') + data.export_velFullProfile_nc('velFullProfile.nc') + # 打印结果 + f = open('velPfParam.txt') + lines=f.readlines() + for line in lines: + print(line) + f.close() + + #输出json文件 + outData = {} + outData["param"] = paramData + outData["lat"] = f'{round(velPf.lat,4)}' + outData["lon"] = f'{round(velPf.lon,4)}' + outData["botDep"] = f'{round(velPf.botDep)}' + properties = { + "表面声道轴深度": f'{round(surfCh.depAxis,2)}' if surfCh.exist else "NONE", + "表面声道声速": f'{round(surfCh.velAxis,4)}' if surfCh.exist else "NONE", + "深海声道轴深度": f'{round(deepCh.depAxis,2)}' if deepCh.exist else "NONE", + "深海声道声速": f'{round(deepCh.velAxis,4)}' if deepCh.exist else "NONE", + "表面跃层区间": f'{round(surfSL.limit[0],2)}, {round(surfSL.limit[1],2)}' if surfSL.exist else "NONE", + "表面跃层强度": f'{round(surfSL.intensity,4)}' if surfSL.exist else "NONE", + "主跃层区间": f'{round(mainSL.limit[0],2)}, {round(mainSL.limit[1],2)}' if mainSL.exist else "NONE", + "主跃层强度": f'{round(mainSL.intensity,4)}' if mainSL.exist else "NONE", + } + velPf.x = np.around(velPf.x, decimals=0) + velPf.y = np.around(velPf.y, decimals=4) + outData["out"] = { + "properties": properties, + "声剖面列表": list(zip(velPf.x,velPf.y)) + } + with open(("" if outParam["outPath"]=="" else outParam["outPath"] + "/")+"profileOut.json", 'w') as jsonFile: + json.dump(outData, jsonFile, indent=4, ensure_ascii=False) + + if png: + # 绘图 + fig, ax = plt.subplots() + # 声速剖面 + ax.plot(velPf.y, velPf.x) + if surfCh.exist: + ax.axhline(surfCh.depAxis, linestyle='--', + color='r', label='surface channel axis') + if deepCh.exist: + ax.axhline(deepCh.depAxis, linestyle='--', + color='g', label='deep channel axis') + if surfSL.exist: + ax.fill_between([1480, 1550], surfSL.limit[0], + surfSL.limit[1], alpha=0.2, color='r') + if mainSL.exist: + ax.fill_between([1480, 1550], mainSL.limit[0], + mainSL.limit[1], alpha=0.2, color='g') + ax.set_xlim(1480, 1550) + ax.invert_yaxis() + ax.set_xlabel('velocity') + ax.set_ylabel('depth') + ax.legend(loc=4) + # 显示绘图 + plt.show() + else: + if ("ETOPO1Path" not in paramData) or ("STPath" not in paramData) or\ + ("profileFile" not in paramData) or ("png" not in paramData) or\ + ("outParam" not in paramData): + print(f'{argfile}中的参数有误，请核实！') + return + ETOPO1Path = paramData["ETOPO1Path"] + STPath = paramData["STPath"] + profileFile = paramData["profileFile"] + png = paramData["png"] + outParam = paramData["outParam"] + if "outPath" not in outParam: + outParam["outPath"] = "" + + import matplotlib.pyplot as plt + from pyat.ProfileTool import ProfileTool + from pyat.ProfileRead import readRawProfile, processProfile + import numpy as np + # 读取实测声速剖面 + velRawPf = readRawProfile(profileFile) + # 对剖面进行后处理 + velPf = processProfile(velRawPf, ETOPO1Path, STPath) + + # 设置剖面分析工具 + velPfTool = ProfileTool() + velPfTool.set_profile(velPf) + # 计算表面声道 + surfCh = velPfTool.get_surfaceChannel() + # 计算深海声道 + deepCh = velPfTool.get_deepChannel() + # 计算跃层信息 + surfSL, mainSL = velPfTool.get_springLayer() + # 导出结果文件 + velPfTool.export_profileParm('velPfParam.txt') + # 打印结果 + f = open('velPfParam.txt') + lines=f.readlines() + for line in lines: + print(line) + f.close() + + #输出json文件 + outData = {} + outData["param"] = paramData + outData["lat"] = f'{round(velPf.lat,4)}' + outData["lon"] = f'{round(velPf.lon,4)}' + outData["botDep"] = f'{round(velPf.botDep)}' + properties = { + "表面声道轴深度": f'{round(surfCh.depAxis,2)}' if surfCh.exist else "NONE", + "表面声道声速": f'{round(surfCh.velAxis,4)}' if surfCh.exist else "NONE", + "深海声道轴深度": f'{round(deepCh.depAxis,2)}' if deepCh.exist else "NONE", + "深海声道声速": f'{round(deepCh.velAxis,4)}' if deepCh.exist else "NONE", + "表面跃层区间": f'{round(surfSL.limit[0],2)}, {round(surfSL.limit[1],2)}' if surfSL.exist else "NONE", + "表面跃层强度": f'{round(surfSL.intensity,4)}' if surfSL.exist else "NONE", + "主跃层区间": f'{round(mainSL.limit[0],2)}, {round(mainSL.limit[1],2)}' if mainSL.exist else "NONE", + "主跃层强度": f'{round(mainSL.intensity,4)}' if mainSL.exist else "NONE", + } + velPf.x = np.around(velPf.x, decimals=0) + velPf.y = np.around(velPf.y, decimals=4) + profileList = list(zip(velPf.x,velPf.y)) + outData["out"] = { + "properties": properties, + "声剖面列表": profileList + } + with open(("" if outParam["outPath"]=="" else outParam["outPath"] + "/")+"profileOut.json", 'w') as jsonFile: + json.dump(outData, jsonFile, indent=4, ensure_ascii=False) + + if png: + # 绘图 + fig, ax = plt.subplots() + # 声速剖面 + ax.plot(velPf.y, velPf.x) + if surfCh.exist: + ax.axhline(surfCh.depAxis, linestyle='--', + color='r', label='surface channel axis') + if deepCh.exist: + ax.axhline(deepCh.depAxis, linestyle='--', + color='g', label='deep channel axis') + if surfSL.exist: + ax.fill_between([1480, 1550], surfSL.limit[0], + surfSL.limit[1], alpha=0.2, color='r') + if mainSL.exist: + ax.fill_between([1480, 1550], mainSL.limit[0], + mainSL.limit[1], alpha=0.2, color='g') + ax.set_xlim(1480, 1550) + ax.invert_yaxis() + ax.set_xlabel('velocity') + ax.set_ylabel('depth') + ax.legend(loc=4) + # 显示绘图 + plt.show() + +if __name__ == '__main__': + main() + diff --git a/pyat/AcousticPara.py b/pyat/AcousticPara.py new file mode 100644 index 0000000..67a3540 --- /dev/null +++ b/pyat/AcousticPara.py @@ -0,0 +1,163 @@ +import numpy as np +import math +from scipy.fftpack import fft +from matplotlib import pyplot as plt +from scipy.interpolate import interp1d +from pyat.DataBase import DataBase + +class BathyM: + def __init__(self, ran=None, depth=None): + self.r = ran + self.d = depth + +class SSP_: + def __init__(self, r, z, c): + self.r = r + self.z = z + self.c = c + +class Sed: + def __init__(self, z, cp, cs, rho, a_p, a_s): + self.z = z + self.cp = cp + self.cs = cs + self.rho = rho + self.a_p = a_p + self.a_s = a_s + +class Base: + def __init__(self, cp, cs, rho, a_p, a_s): + self.cp = cp + self.cs = cs + self.rho = rho + self.a_p = a_p + self.a_s = a_s + +class Coord: + def __init__(self, lon, lat): + self.lon = lon + self.lat = lat + +def GetSeaDepth(data, N, coordS, R, azi): + r = np.linspace(0, R, N) + lon = coordS.lon + r * math.sin(azi) / (111 * math.cos(coordS.lat / 180 * math.pi)) + lat = coordS.lat + r * math.cos(azi) / 111 + z = np.zeros(N) + for it in range(N): + data.set_coord(lat[it], lon[it]) + z[it] = data.get_bottomDepth() + return r, z + + +def GetEnv(data, freq, coordS, coordR=None, R=None, azi=None): + if coordR == None: + azi = azi/180*math.pi + lon = coordS.lon + R * math.sin(azi) / (111 * math.cos(coordS.lat/180*math.pi)) + lat = coordS.lat + R * math.cos(azi) / 111 + coordR = Coord(lon, lat) + + if R == None: + x = (coordR.lon - coordS.lon) * (111 * math.cos(coordS.lat/180*math.pi)) + y = (coordR.lat - coordS.lat) * 111 + R = math.sqrt(x*x + y*y) + if y > 0: + azi = math.atan(x / y) + elif y == 0: + azi = np.sign(x) * math.pi / 2 + else: + azi = math.atan(x / y) + math.pi + + # Get Bathymetry + N = math.ceil(1000 * R / (1500 / freq)) + #print("N = {}".format(N)) + if N == 1: + N1 = 2 + else: + r, z = GetSeaDepth(data, N, coordS, R, azi) + Z = np.square(abs(fft(z-np.mean(z)))) + s1 = np.sum(Z[0:math.ceil(N / 2)-1]) + s2 = 0 + for it in range(math.ceil(N / 2)): + s2 = s2 + Z[it] + if s2 / s1 > 0.98: + break + N1 = max(2 * it, 1) + # plt.plot(r, z) + + #print("N1 = {}".format(N1)) + r, z = GetSeaDepth(data, N1, coordS, R, azi) + bathym = BathyM(r, z) + #plt.plot(r, z) + #plt.gca().invert_yaxis() + #plt.show() + + # Get Sound Speed Profile + + # N2 = math.ceil(1000 * R / (30 * 1500 / freq)) + N2 = N1 +# N2=math.ceil(N2/2) + r = np.linspace(0, R, N2) + lon = coordS.lon + r * math.sin(azi) / (111 * math.cos(coordS.lat / 180 * math.pi)) + lat = coordS.lat + r * math.cos(azi) / 111 + ssp = [] + + z_ssp = np.arange(0, 30, 1) + z_ssp = np.append(z_ssp, np.arange(30, 100, 10)) + z_ssp = np.append(z_ssp, np.arange(100, 2000, 50)) + z_ssp = np.append(z_ssp, np.arange(2000, 10000, 100)) + z_ssp = np.append(z_ssp, np.ceil(bathym.d)) + z_ssp = np.sort(z_ssp) + z_ssp = np.unique(z_ssp) + + + for it in range(N2): + data.set_coord(lat[it], lon[it]) + velFullPf = data.get_velFullProfile() + itp1 = interp1d(velFullPf.x, velFullPf.y) + z_idx = (z_ssp <= velFullPf.x[-1]) & (z_ssp >= 0) + ssp_z = z_ssp[z_idx] + ssp_c = itp1(ssp_z) + ssp.append(SSP_(r[it], ssp_z, ssp_c)) + #ssp.append(SSP_(r[it], velFullPf.x, velFullPf.y)) + #plt.plot(velFullPf.y, velFullPf.x) + + # plt.gca().invert_yaxis() + # plt.show() + + # Sediment + sed = [] + z = np.array([0, 50]) + cp = np.array([1550, 1590]) + cs = np.array([0, 0]) + rho = np.array([1.3, 1.5]) + a_p = np.array([0.1, 0.1]) + a_s = np.array([0, 0]) + for it in range(N2): + sed.append(Sed(z, cp, cs, rho, a_p, a_s)) + + # Base + base = [] + Cb_p = 1550 # Sound speed of P-wave + Cb_s = 0 # Sound speed of S-wave + Rhob = 1.4 # Density of sea bottom + Ab_p = 0.2 # Attenuation of P-wave + Ab_s = 0 # Attenuation of S-wave + for it in range(N2): + base.append(Base(Cb_p, Cb_s, Rhob, Ab_p, Ab_s)) + + return bathym, ssp, sed, base + +# if __name__ == "__main__": +# freq = 25 +# # 定义坐标 +# lat = 18.00 +# lon = 114.00 +# # 定义时间段 +# timeIndex = db.timeDict['Annual'] +# # 创建类 +# data = db.DataBase() +# # 设置时间段 +# data.set_timePeriod(timeIndex) + +# #[bathym, ssp, sed, base] = GetEnv(data, freq, Coord(lon, lat), R=111, azi=-90) +# [bathym, ssp, sed, base] = GetEnv(data, freq, Coord(lon, lat), Coord(113,19)) diff --git a/pyat/DataBase.py b/pyat/DataBase.py new file mode 100644 index 0000000..5c21c5e --- /dev/null +++ b/pyat/DataBase.py @@ -0,0 +1,348 @@ +import numpy as np +from scipy.interpolate import Akima1DInterpolator, interp1d +from netCDF4 import Dataset +from pyat.ProfileDataType import Profile + +''' +************************************************************ +该模块用于从数据库中, 获取剖面信息 +************************************************************ +''' + +# 时间段标识 + + + +# 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): + # 打开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): + """ + 设置时间段 + """ + self.timeIndex = timeIndex + # 打开WOA18数据 + #fileWOA18sal = '.\WOA18\woa18_A5B7_s' + timeIndex + '_04.nc' + #fileWOA18tem = '.\WOA18\woa18_A5B7_t' + timeIndex + '_04.nc' + fileST=self.STpath+'/CA'+timeIndex+'.nc' + dataST = Dataset(fileST) + # 获取WOA18网格 + self.latGrid2 = dataST.variables['lat'][:] + self.lonGrid2 = dataST.variables['lon'][:] + # 获取深度数据 + self.depthGrid = dataST.variables['zlev'][:]*(-1) + # 获取盐度数据 + self.salGrid = np.squeeze(dataST.variables['SALT'][:]) + # 获取温度数据 + self.temGrid = np.squeeze( + dataST.variables['THETA'][:]) + 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 == 0: + 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 == 0: + 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() + if self.isValid == True: + # 获得原始声速剖面数据 + dep = self.get_velRawProfile().x + vel = self.get_velRawProfile().y + # 对剖面进行垂向插值处理，将非等距离采样填充为深度间隔为1m的序列 + depInterp = np.arange(0, dep[-1]+1, 1) + # 应用Akima插值算法，对数据进行插值 + if dep.size == 1: + velInterp = vel + else: + 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(np.ceil(dep[-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() diff --git a/pyat/OANPd.py b/pyat/OANPd.py new file mode 100644 index 0000000..9fef0a5 --- /dev/null +++ b/pyat/OANPd.py @@ -0,0 +1,89 @@ +import numpy as np +from scipy.fft import fft +import soundfile as sf +from numpy import log10 +import math + +# Average Noise Spectral Level +def ANSL(u, fs, M, A=0, N=None, D=None, w=None): + """ + :param u: Voltage data (Unit: V) + :param fs: Sampling rate (Unit: Hz) + :param M: Sensitivity (Unit: dBV re 1uPa) + :param A: Preamplifier (Unit: dB) + :param N: Segment length + :param D: Overlap length + :param w: Window + :return: Spectral level + """ + L = len(u) + if N == None: + if L >= fs: + N = round(fs) + else: + N = L + elif N > L: + raise Exception("Segment length is larger than data length.") + + if D == None: + D = 0 + + if w is None: + w = np.hamming(N) + G = np.mean(np.square(abs(w))) + + K = 1 + math.floor((L-N) / (N-D)) + + T = N / fs + df = 1 / T + X2 = np.zeros(N) + for k in range(K): + X = fft(w * u[round(k*(N-D)):round(k*(N-D)+N)]) + X2 += np.square(abs(X)) + + spl = 10 * log10(X2) - 10 * log10(N*fs*G*K/2) - M - A + freq = np.arange(0, fs/2, df) + spl = spl[0:len(freq)] + + return freq, spl + + +# 1/3 Octave Spectral Level +def OneThirdOct_SPL(freq, spl, f_start, Npoints): + """ + :param freq: Frequency + :param spl: 1 Hz bandwidth spectral level + :param f_start: Starting frequency of the 1/3 octave + :param Npoints: The number of frequency points + :return: 1/3 Octave Spectral Level + """ + fC = f_start * np.power(2, np.arange(Npoints)/3) + fL = f_start * np.power(2, np.arange(Npoints)/3 - 1.0/6) + fH = f_start * np.power(2, np.arange(Npoints)/3 + 1.0/6) + + spl_oct = np.zeros(Npoints) + for it in range(Npoints): + # idx = (freq >= fL[it]) & (freq <= fH[it]) + idx = np.where((freq >= fL[it]) & (freq <= fH[it])) + # print(idx) + # print(spl[idx]) + spl_oct[it] = 10 * log10(np.mean(np.power(10, spl[idx]/10))) + # print(spl_oct[it]) + + return fC, spl_oct + +def call_oanp(soundFile, M, A=0): + [x, fs] = sf.read(soundFile) + u = 3 * x # Voltage data + # M = -164.0 # Sensitivity + # A = 0 # Preamplifier + N = fs # Segment length + D = N/2 # Overlap length + w = np.hanning(N) # Window + + # Average Noise Spectral Level + [freq, spl] = ANSL(u, fs, M, A, N, D, w) + # 1/3 Octave Spectral Level + [fC, spl_oct] = OneThirdOct_SPL(freq, spl, f_start=20, Npoints=28) + + return [freq, spl, fC, spl_oct] \ No newline at end of file diff --git a/pyat/ProfileDataType.py b/pyat/ProfileDataType.py new file mode 100644 index 0000000..7832c66 --- /dev/null +++ b/pyat/ProfileDataType.py @@ -0,0 +1,53 @@ +class Profile(): + """ + 记录剖面数据 + """ + + def __init__(self): + self.x = None + self.y = None + self.lat = None + self.lon = None + self.botDep = None + + def set_data(self, x, y): + ''' + 设置曲线 + ''' + self.x = x + self.y = y + + def set_coord(self, lat, lon): + ''' + 设置坐标 + ''' + self.lat = lat + self.lon = lon + + def set_bottomDepth(self, botDep): + ''' + 设置深度 + ''' + self.botDep = botDep + + +class Channel(): + ''' + 储存声道道信息 + ''' + + def __init__(self): + self.exist = False + self.depAxis = None + self.velAxis = None + + +class SpringLayer(): + ''' + 储存跃层信息 + ''' + + def __init__(self): + self.exist = False + self.limit = None + self.intensity = None diff --git a/pyat/ProfileRead.py b/pyat/ProfileRead.py new file mode 100644 index 0000000..82037db --- /dev/null +++ b/pyat/ProfileRead.py @@ -0,0 +1,168 @@ +import unicodedata +import numpy as np +from pyat.ProfileDataType import Profile +from scipy.interpolate import interp1d +from scipy.signal import butter, filtfilt +from pyat.DataBase import DataBase + +''' +************************************************************ +该模块用读取实测的声速剖面 +************************************************************ +''' + +def is_number(s): + ''' + 判断是否为数字 + ''' + try: + float(s) + return True + except ValueError: + pass + + try: + unicodedata.numeric(s) + return True + except (TypeError, ValueError): + pass + + return False + +def readRawProfile(fname): + # 读取经纬度, 海深等信息 + f = open(fname) + line = f.readline() + line = line.split(',') + if is_number(line[1]): + lat = float(line[1]) + else: + lat = None + if is_number(line[3]): + lon = float(line[3]) + else: + lon = None + if is_number(line[5]): + botDep = float(line[5]) + else: + botDep = None + f.close() + # 读取数据 + data = np.loadtxt(fname, skiprows=2,delimiter=',') + # 剖面赋值 + velPf = Profile() + velPf.set_data(data[:, 0], data[:, 1]) + velPf.set_coord(lat, lon) + velPf.set_bottomDepth(botDep) + return velPf + + +def removeHead(velPf): + ''' + 去掉剖面头部异常数据, 刚入水的 + ''' + lat = velPf.lat + lon = velPf.lon + botDep = velPf.botDep + dep = velPf.x + vel = velPf.y + # 梯度阈值 + gradc = 1.0 + # 去掉前面一段（刚入水）变化计较剧烈的区域 + for i in range(dep.size-1): + grad = np.abs((vel[i+1]-vel[i])/(dep[i+1]-dep[i])) + if grad < gradc: + break + dep = dep[i+1:] + vel = vel[i+1:] + # 剖面赋值 + velPf = Profile() + velPf.set_data(dep, vel) + velPf.set_coord(lat, lon) + velPf.set_bottomDepth(botDep) + return velPf + + +def removeTail(velPf, ETOPO1Path, STPath): + ''' + 去掉剖面尾部异常数据, 触底后的 + ''' + lat = velPf.lat + lon = velPf.lon + botDep = velPf.botDep + dep = velPf.x + vel = velPf.y + # 如果没有给定深度，则直接从数据库中提取 + if botDep == None: + data = DataBase(ETOPO1Path, STPath) + data.set_coord(lat, lon) + botDep = data.get_bottomDepth() + # 过滤掉触底后的曲线 + temp = dep < botDep + dep = dep[temp] + vel = vel[temp] + # 剖面赋值 + velPf = Profile() + velPf.set_data(dep, vel) + velPf.set_coord(lat, lon) + velPf.set_bottomDepth(botDep) + return velPf + + +def resampleProfile(velPf): + ''' + 重采样剖面, 间隔1m采样 + ''' + lat = velPf.lat + lon = velPf.lon + botDep = velPf.botDep + dep = velPf.x + vel = velPf.y + # 得到等间距的采样曲线 + depInterp = np.arange(0, dep[-1]) + f = interp1d(dep, vel, fill_value='extrapolate') + velInterp = f(depInterp) + # 剖面赋值 + velPf = Profile() + velPf.set_data(depInterp, velInterp) + velPf.set_coord(lat, lon) + velPf.set_bottomDepth(botDep) + return velPf + + +def lowPassFilter(velPf): + ''' + 低通滤波器, 使曲线平滑 + ''' + lat = velPf.lat + lon = velPf.lon + botDep = velPf.botDep + dep = velPf.x + vel = velPf.y + # 设计数字滤波器 + freqc = 0.05 + order = 2 + b, a = butter(order, freqc) + # 对曲线进行滤波 + vel = filtfilt(b, a, vel) + # 剖面赋值 + velPf = Profile() + velPf.set_data(dep, vel) + velPf.set_coord(lat, lon) + velPf.set_bottomDepth(botDep) + return velPf + + +def processProfile(velRawPf, ETOPO1Path, STPath): + ''' + 对原始剖面数据进行处理, 得到处理后的剖面 + ''' + # 去掉刚入水，剧烈变化部分 + velPf = removeHead(velRawPf) + # 去掉触底后，无效部分 + velPf = removeTail(velPf, ETOPO1Path, STPath) + # 等间距重采样 + velPf = resampleProfile(velPf) + # 低通滤波 + velPf = lowPassFilter(velPf) + return velPf diff --git a/pyat/ProfileTool.py b/pyat/ProfileTool.py new file mode 100755 index 0000000..01009cb --- /dev/null +++ b/pyat/ProfileTool.py @@ -0,0 +1,237 @@ +import numpy as np +from scipy.signal import find_peaks +from pyat.ProfileDataType import Channel, SpringLayer + +''' +************************************************************ +该模块用于分析剖面, 根据标准化剖面，得到声道, 跃层分析结果 +************************************************************ +''' + + +class ProfileTool(): + ''' + 剖面分析工具, 根据标准化剖面，得到声道, 跃层分析结果 + ''' + + def __init__(self): + self.clearVar() + + def clearVar(self): + """ + 清空函数调用记录 + """ + self.maxTurn = None + self.minTurn = None + self.surfCh = None + self.deepCh = None + self.surfSL = None + self.mainSL = None + return + + def set_profile(self, velPf): + ''' + 设置剖面 + ''' + self.velPf = velPf + self.clearVar() + + def findMaxTurn(self): + ''' + 寻找最大值的拐点 + ''' + if (self.maxTurn == None): + dep = self.velPf.x + vel = self.velPf.y + # 寻找极大值，设定窗口宽度为20m + peaks, properties = find_peaks(vel, width=20) + if (peaks.size != 0): + peaks = peaks[vel[peaks] > vel[0]] + if (peaks.size != 0): + # 定位极大值中最大的那个 + peak = peaks[np.argmax(vel[peaks])] + self.maxTurn = peak + else: + self.maxTurn = -1 + return self.maxTurn + + def findMinTrun(self): + ''' + 寻找最小值的拐点 + ''' + if (self.minTurn == None): + dep = self.velPf.x + vel = self.velPf.y + # 寻找最小值拐点，设定窗口宽度为20m + peaks, properties = find_peaks(-vel, width=20) + if (peaks.size != 0): + peaks = peaks[vel[peaks] < vel[0]] + if (peaks.size != 0): + # 定位极小值中最小的那个 + peak = peaks[np.argmin(vel[peaks])] + self.minTurn = peak + else: + self.minTurn = -1 + return self.minTurn + + def get_surfaceChannel(self): + ''' + 计算表面声道 + ''' + if (self.surfCh == None): + dep = self.velPf.x + vel = self.velPf.y + maxTurn = self.findMaxTurn() + minTurn = self.findMinTrun() + if maxTurn != -1: + self.surfCh = Channel() + self.surfCh.exist = True + self.surfCh.depAxis = dep[0] + self.surfCh.velAxis = vel[0] + else: + self.surfCh = Channel() + self.surfCh.exist = False + return self.surfCh + + def get_deepChannel(self): + ''' + 计算深海声道 + ''' + if (self.deepCh == None): + dep = self.velPf.x + vel = self.velPf.y + maxTurn = self.findMaxTurn() + minTurn = self.findMinTrun() + if minTurn != -1: + self.deepCh = Channel() + self.deepCh.exist = True + self.deepCh.depAxis = dep[minTurn] + self.deepCh.velAxis = vel[minTurn] + else: + self.deepCh = Channel() + self.deepCh.exist = False + return self.deepCh + + def get_springLayer(self): + ''' + 计算跃层 + ''' + if (self.surfSL == None): + dep = self.velPf.x + vel = self.velPf.y + # 寻找拐点 + maxTurn = self.findMaxTurn() + minTurn = self.findMinTrun() + # 初始化跃层 + surfSL = SpringLayer() + surfSL.exist = False + mainSL = SpringLayer() + mainSL.exist = False + if (maxTurn != -1 and minTurn != -1): + # 对应这种形状 + # \ + # \ + # / + # / + # \ + # \ + surfSL.exist = True + surfSL.limit = np.array([dep[0], dep[maxTurn]]) + surfSL.intensity = (vel[maxTurn]-vel[0])/(dep[maxTurn]-dep[0]) + mainSL.exist = True + mainSL.limit = np.array([dep[maxTurn], dep[minTurn]]) + mainSL.intensity = ( + vel[minTurn]-vel[maxTurn])/(dep[minTurn]-dep[maxTurn]) + elif (maxTurn != -1 and minTurn == -1): + # 对应这种形状 + # \ + # \ + # / + # / + # / + surfSL.exist = True + surfSL.limit = np.array([dep[0], dep[maxTurn]]) + surfSL.intensity = (vel[maxTurn]-vel[0])/(dep[maxTurn]-dep[0]) + mainSL.exist = True + mainSL.limit = np.array([dep[maxTurn], dep[-1]]) + mainSL.intensity = (vel[-1]-vel[maxTurn]) / \ + (dep[-1]-dep[maxTurn]) + elif (maxTurn == -1 and minTurn != -1): + # 对应这种形状 + # / + # / + # / + # \ + # \ + # \ + mainSL.exist = True + mainSL.limit = np.array([dep[0], dep[minTurn]]) + mainSL.intensity = (vel[minTurn]-vel[0])/(dep[minTurn]-dep[0]) + elif (maxTurn == -1 and minTurn == -1): + if vel[-1] < vel[0]: + # 对应这种形状 + # / + # / + # / + mainSL.exist = True + mainSL.limit = np.array([dep[0], dep[-1]]) + mainSL.intensity = (vel[-1]-vel[0])/(dep[-1]-dep[0]) + else: + # 对应这种形状 + # \ + # \ + # \ + surfSL.exist = True + surfSL.limit = np.array([dep[0], dep[-1]]) + surfSL.intensity = (vel[-1]-vel[0])/(dep[-1]-dep[0]) + self.surfSL = surfSL + self.mainSL = mainSL + return (self.surfSL, self.mainSL) + + def export_profileParm(self, fileName): + ''' + 计算结果输出文件 + ''' + lat = self.velPf.lat + lon = self.velPf.lon + botDep = self.velPf.botDep + # 计算表面声道 + surfCh = self.get_surfaceChannel() + # 计算深海声道 + deepCh = self.get_deepChannel() + # 计算跃层信息 + surfSL, mainSL = self.get_springLayer() + f = open(fileName, mode='w') + # 输出坐标深度 + line = 'lat[degree north], {0:>7.4f}, lon[degree east], {1:>7.4f}, BottomDepth[m], {2:>5.0f}\n'.format( + lat, lon, botDep) + f.write(line) + f.write('----------------------------\n') + if surfCh.exist: + f.write('表面声道轴深度: {0:>.2f}\n'.format(surfCh.depAxis)) + f.write('表面声道声速: {0:>.4f}\n'.format(surfCh.velAxis)) + else: + f.write('表面声道轴深度: NONE\n') + f.write('表面声道声速: NONE\n') + f.write('----------------------------\n') + if deepCh.exist: + f.write('深海声道轴深度: {0:>.2f}\n'.format(deepCh.depAxis)) + f.write('深海声道声速: {0:>.4f}\n'.format(deepCh.velAxis)) + else: + f.write('深海声道轴深度: NONE\n') + f.write('深海声道声速: NONE\n') + f.write('----------------------------\n') + if surfSL.exist: + f.write('表面跃层区间: {0:>.2f}, {1:>.2f}\n'.format(surfSL.limit[0], surfSL.limit[1])) + f.write('表面跃层强度: {0:>.4f}\n'.format(surfSL.intensity)) + else: + f.write('表面跃层区间: NONE\n') + f.write('表面跃层强度: NONE\n') + f.write('----------------------------\n') + if mainSL.exist: + f.write('主跃层区间: {0:>.2f}, {1:>.2f}\n'.format(mainSL.limit[0], mainSL.limit[1])) + f.write('主跃层强度: {0:>.4f}\n'.format(mainSL.intensity)) + else: + f.write('主跃层区间: NONE\n') + f.write('主跃层强度: NONE\n') + f.close() diff --git a/pyat/RANDIcore.py b/pyat/RANDIcore.py new file mode 100644 index 0000000..13455a0 --- /dev/null +++ b/pyat/RANDIcore.py @@ -0,0 +1,319 @@ +import os +from multiprocessing import Pool +from functools import partial +from 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 + + +class BathyM: + def __init__(self, ran=None, depth=None): + self.r = ran + self.d = depth + + +class SSP_: + def __init__(self, r, z, c): + self.r = r + self.z = z + self.c = c + + +class Sed: + def __init__(self, z, cp, cs, rho, a_p, a_s): + self.z = z + self.cp = cp + self.cs = cs + self.rho = rho + self.a_p = a_p + self.a_s = a_s + + +class Coord: + def __init__(self, lon, lat): + self.lon = lon + self.lat = lat + + +class RANDI_Para: + def __init__(self, baseDir, freq, lon, lat, Range, dRange, dTheta, v_ship, l_ship, theta=0): + self.baseDir = baseDir + self.freq = freq + self.lon = lon + self.lat = lat + self.Range = Range + self.dRange = dRange + self.dTheta = dTheta + self.v_ship = v_ship + self.l_ship = l_ship + self.theta = theta + + def refreshTheta(self, theta): + self.theta = theta + return self + + +def RANDI_GetSeaDepth(data, N, coordS, R, azi): + r = np.linspace(0, R, N) + lon = coordS.lon + r * np.sin(azi) / (111 * np.cos(coordS.lat / 180 * np.pi)) + lat = coordS.lat + r * np.cos(azi) / 111 + + z = np.zeros(N) + for it in range(N): + data.set_coord(lat[it], lon[it]) + z[it] = data.get_bottomDepth() + return r, z + + +def RANDI_GetMeanProfile(data, N, coordS, R, azi): + r = np.linspace(0, R, N) + lon = coordS.lon + r * np.sin(azi) / (111 * np.cos(coordS.lat / 180 * np.pi)) + lat = coordS.lat + r * np.cos(azi) / 111 + sal = np.full([len(data.depthGrid), N], np.nan) + tem = np.full([len(data.depthGrid), N], np.nan) + for i in range(N): + data.set_coord(lat[i], lon[i]) + salPf = data.get_salRawProfile() + temPf = data.get_temRawProfile() + if len(salPf.x) > 1: + sal[0:len(salPf.x), i] = salPf.y + tem[0:len(temPf.x), i] = temPf.y + sal = np.ma.masked_array(sal, np.isnan(sal)) + tem = np.ma.masked_array(tem, np.isnan(tem)) + sal_m = np.mean(sal, 1) + tem_m = np.mean(tem, 1) + z = data.depthGrid[~sal_m.mask] + sal_m = sal_m[~sal_m.mask] + tem_m = tem_m[~tem_m.mask] + return z, sal_m, tem_m + + +def RANDI_GetEnv(data, coordS, coordR=None, R=None, azi=None): + if coordR == None: + azi = azi / 180 * np.pi + lon = coordS.lon + R * np.sin(azi) / (111 * np.cos(coordS.lat / 180 * np.pi)) + lat = coordS.lat + R * np.cos(azi) / 111 + coordR = Coord(lon, lat) + + if R == None: + x = (coordR.lon - coordS.lon) * (111 * np.cos(coordS.lat / 180 * np.pi)) + y = (coordR.lat - coordS.lat) * 111 + R = np.sqrt(x * x + y * y) + if y > 0: + azi = np.arctan(x / y) + elif y == 0: + azi = np.sign(x) * np.pi / 2 + else: + azi = np.arctan(x / y) + np.pi + + # Get Bathymetry + N = int(np.ceil(R / 10) + 1) + r, z = RANDI_GetSeaDepth(data, N, coordS, R, azi) + bathym = BathyM(r, z) + # plt.plot(r, z) + # plt.gca().invert_yaxis() + # plt.show() + + # Get Sound Speed Profile + z_ssp, sal_m, tem_m = RANDI_GetMeanProfile(data, N, coordS, R, azi) + if max(bathym.d) > z_ssp[-1]: + 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) + ssp = SSP_(0, z_ssp, c_ssp) + # plt.plot(c_ssp, z_ssp) + # plt.gca().invert_yaxis() + # plt.show() + + # Sediment + z = np.array([0, 50]) + cp = np.array([1550, 1590]) + cs = np.array([0, 0]) + rho = np.array([1.3, 1.5]) + a_p = np.array([0.1, 0.1]) + a_s = np.array([0, 0]) + sed = Sed(z, cp, cs, rho, a_p, a_s) + + return bathym, ssp, sed + + +def RANDI_writeRAMin(baseDir, data, freq, dTheta, lat, lon, Range, rd, ship_sd, SeaDepthLim): + t1 = time() + btyFile = open("bty", 'w') + for theta in np.arange(0, 360, dTheta): + [bathm, ssp, sed] = RANDI_GetEnv(data, Coord(lon, lat), R=Range, azi=theta) + for dd in bathm.d: + btyFile.write(str(dd) + ' ') + btyFile.write('\n') + idx = np.where(bathm.d < SeaDepthLim) + if len(idx[0]) > 0: + bathm.r = bathm.r[0:idx[0][0]] + bathm.d = bathm.d[0:idx[0][0]] + + title = "Waveguide" + # Source and receiving position + bathm.r = bathm.r * 1000 # km -> m + rmax = max(bathm.r) # Maximum receiving range (unit: m) + zmax = round(max(bathm.d) + 5 * 1500 / freq) # Maximum receiving depth (uint: m) + zmplt = 0 # round(max(bathm.d)*1.05) # Maximum depth of plot (uint: m) + dz = 1500 / freq / 4 # Step in depth direction: wavelength / 15 + dr = 5 * dz # Step in range direction + pos = ramRW.PosRAM(rd, ship_sd, rmax, zmax, dr, dz, zmplt) + + # Sound speed profile for RAM + Z = sed.z + cw_p = ramRW.PROF(ssp.z, ssp.c) # Water sound speed profile + cb_p = ramRW.PROF(Z, sed.cp) # Bottom sound speed profile + rhob_p = ramRW.PROF(Z, sed.rho) # Sea bottom density profile + + Ap = sed.a_p + if Z[-1] < 5 * 1500 / freq: + Z = np.append(Z, 3 * 1500 / freq) + Ap = np.append(Ap, Ap[-1]) + Z = np.append(Z, Z[-1] + 2 * 1500 / freq) + Ap = np.append(Ap, 10) + # Ap = sed[it].a_p + attn_p = ramRW.PROF(Z, Ap) # Sea bottom attenuation profile + + ssp_ram = [ramRW.SSPRAM(0, cw_p, cb_p, rhob_p, attn_p)] + + # Write ram.in file + if not os.path.exists(baseDir + "/TL/temp/theta_{:1d}".format(theta)): + os.makedirs(baseDir + "/TL/temp/theta_{:1d}".format(theta)) + envfil = baseDir + "/TL/temp/theta_{:1d}/ram.in".format(theta) + ramRW.write_env_RAM(title, envfil, freq, pos, bathm, ssp_ram) + t2 = time() + #print(t2 - t1) + btyFile.close() + + +def RANDI_callRAM(baseDir, NThreads, freq, dTheta): + ''' + t1 = time() + for theta in np.arange(0,360,5): + os.chdir(baseDir + "/TL/theta_{:1d}".format(theta)) + os.system(baseDir + "/ram.exe") + + os.chdir(baseDir) + t2 = time() + print(t2-t1) + ''' + t1 = time() + func = partial(callRAM, baseDir) + pool = Pool(NThreads) + pool.map(func, np.arange(0, 360, dTheta)) + t2 = time() + #print(t2 - t1) + + if os.path.exists(baseDir + f"/TL/freq_{round(freq)}"): + rmtree(baseDir + f"/TL/freq_{round(freq)}") + os.rename(baseDir + "/TL/temp", baseDir + f"/TL/freq_{round(freq)}") + + +def callRAM(baseDir, theta): + AtBinPath = os.path.dirname(os.path.abspath(__file__))+'/at/bin' + print(AtBinPath) + os.chdir(baseDir + "/TL/temp/theta_{:1d}".format(theta)) + os.system(AtBinPath + "/ram") + os.chdir(baseDir) + +def RANDI_getShipNum(lat, lon, Range, dRange, theta): + NRange = int(Range / dRange) + #shipNum = np.random.rand(NRange) + shipNum = 0.5*np.ones(NRange) + return shipNum + +def calNoise(randiP): + rr = np.arange(randiP.dRange / 2, randiP.Range, randiP.dRange) + SL = RossShipSL(randiP.freq, randiP.v_ship, randiP.l_ship) + S = np.power(10, SL / 10) + shipNum = RANDI_getShipNum(randiP.lat, randiP.lon, randiP.Range, randiP.dRange, randiP.theta) + filePath = randiP.baseDir + f"/TL/freq_{round(randiP.freq)}/theta_{round(randiP.theta)}" + [r, tl] = ramRW.readTLline_RAM(filePath) + r = r / 1000 + I2 = np.power(10, tl / 10) + I_ = interp1d(r, I2) + I = I_(rr[rr < r[-1]]) + noise = sum(shipNum[rr < r[-1]] * S / I) + return noise + +def RANDI_calNL(randiP, NThreads): + ''' + t1 = time() + theta = np.arange(0, 360, dTheta) + noise1 = np.zeros(len(theta)) + for i in range(len(theta)): + randiP.theta = theta[i] + noise1[i] = calNoise(randiP) + NL1 = 10 * np.log10(sum(noise1)) + t2 = time() + print(t2-t1) + ''' + + t1 = time() + RANDI_P = [] + for theta in np.arange(0, 360, randiP.dTheta): + randiP.theta = theta + RANDI_P.append(deepcopy(randiP)) + pool = Pool(NThreads) + noise = pool.map(calNoise, RANDI_P) + NL = 10 * np.log10(sum(noise)) + t2 = time() + #print(t2-t1) + return NL + + +def RossShipSL(freq, v, l): + if freq < 500: + Lso = -10 * np.log10(np.power(10, -1.06 * np.log10(freq) - 14.34) + + np.power(10, 3.32 * np.log10(freq) - 21.425)) + else: + Lso = 173.2 - 18 * np.log10(freq) + + if freq <= 28.4: + df = 8.1 + elif freq < 191.6: + df = 22.3-9.77*np.log10(freq) + else: + df = 0 + + dl = np.power(l, 1.15) / 3643 + Ls = Lso + 60*np.log10(v/12) + 20*np.log10(l/300) + df*dl + 3 + return Ls + + +def call_RANDI(ETOPO1Path, STPath, baseDir, freq, rd, lon, lat, timeType, NThreads): + # 定义时间段 + timeIndex = timeType + # 创建类 + data = DataBase(ETOPO1Path, STPath) + # 设置时间段 + data.set_timePeriod(timeIndex) + + Range = 500 # Range of sound field calculation (Unit: km) + ship_sd = 6.0 # Ship source depth (Unit: m) + dTheta = 5 # The angle of sector (Unit: °) + dRange = 10 # The + SeaDepthLim = 20 # (Unut: m) + v_ship = 12 # Ship speed (Unit: kt) + l_ship = 300 # Ship length (Unit: ft) + + randiP = RANDI_Para(baseDir, freq, lon, lat, Range, dRange, dTheta, v_ship, l_ship) + + RANDI_writeRAMin(baseDir, data, freq, dTheta, lat, lon, Range, rd, ship_sd, SeaDepthLim) + + RANDI_callRAM(baseDir, NThreads, freq, dTheta) + + NL = RANDI_calNL(randiP, NThreads) + + return NL + + + + + diff --git a/pyat/WGNPd.py b/pyat/WGNPd.py new file mode 100755 index 0000000..6ff987a --- /dev/null +++ b/pyat/WGNPd.py @@ -0,0 +1,143 @@ +""" Wind Generated Noise Prediction """ +import numpy as np +from scipy.interpolate import interp1d +from pyat import DataBase as db +import matplotlib.pyplot as plt + +class MdlPara(): + def __init__(self, title, lat, lon, seaDepth, month, freq, vc, a, b, vs): + self.title = title + self.lat = lat + self.lon = lon + self.seaDepth = seaDepth + self.month = month + self.freq = freq + self.vc = vc + self.a = a + self.b = b + self.vs = vs + +def loadMdlPara(filename): + f = open(filename, mode='r') + + mdlPara = [] + eof = False + while (1): + # Title + while (1): + line = f.readline() + if line == '': + eof = True + break + title = line.strip() + if title != '': + break + if eof: + break + + # Latitude and Longitude + line = f.readline() + line = line.strip().split() + lat = np.double(line[0]) + lon = np.double(line[1]) + + # Sea Depth + line = f.readline() + line = line.strip().split() + seaDepth = np.double(line[0]) + + # Month + line = f.readline() + line = line.strip().split() + month = np.int8(line[0]) + + # Model Parameters + freq = [] + vc = [] + a = [] + b = [] + vs = [] + while (1): + line = f.readline() + line = line.strip().split() + if line[0] == '0': + break + freq.append(np.double(line[0])) + vc.append(np.double(line[1])) + a.append(np.double(line[2])) + b.append(np.double(line[3])) + vs.append(np.double(line[4])) + + mdlPara.append(MdlPara(title, lat, lon, seaDepth, month, freq, vc, a, b, vs)) + + return mdlPara + + +def WGNMdl(mdlPara, windspeed): + f = np.array(mdlPara.freq) + spl = np.zeros(len(f)) + for it in range(len(f)): + if windspeed < mdlPara.vc[it]: + spl[it] = mdlPara.a[it] + mdlPara.b[it] * 10 * np.log10(mdlPara.vc[it]) + elif (mdlPara.vs[it] is not np.NaN) & (windspeed > mdlPara.vs[it]): + spl[it] = mdlPara.a[it] + mdlPara.b[it] * 10 * np.log10(mdlPara.vs[it]) + else: + spl[it] = mdlPara.a[it] + mdlPara.b[it] * 10 * np.log10(windspeed) + + return f, spl + + +def WGNPd(dataFile, seaDepth, windspeed, month, freq): + if min(freq) < 50: + raise UserWarning("The minimum frequency must be larger than 50 Hz") + + mdlPara = loadMdlPara(dataFile) + + z = seaDepth + + Q = [] + for it in range(len(mdlPara)): + q1 = abs(mdlPara[it].seaDepth - z) / (mdlPara[it].seaDepth + z) + q2 = min(abs(month-mdlPara[it].month), 12-abs(month-mdlPara[it].month)) / 6 + alpha = 0.7 + Q.append(alpha * q1 + (1-alpha) * q2) + idx = np.argmin(Q) + # print([z, month, idx]) + + f0, spl0 = WGNMdl(mdlPara[idx], windspeed) + if max(freq) > max(f0): + spl0 = np.append(spl0, spl0[-1] - 5 * np.log2(max(freq)/max(f0))) + f0 = np.append(f0, max(freq)) + if min(f0) > 200: + spl0 = np.insert(spl0, 0, spl0[0] + 5 * np.log2(min(f0)/200)) + f0 = np.insert(f0, 0, 200) + if min(f0) > 50: + spl0 = np.insert(spl0, 0, spl0[0] + 9 * np.log2(min(f0)/50)) + f0 = np.insert(f0, 0, 50) + + itp1 = interp1d(np.log10(f0), spl0) + spl = itp1(np.log10(freq)) + + return freq, spl + + +# if __name__=="__main__": +# lat = 19 +# lon = 116 +# windspeed = 5 # unit: m/s +# month = 3 +# freq = 50 * np.power(2, np.arange(25)/3) + +# # 定义时间段 +# timeIndex = db.timeDict['Annual'] +# # 创建类 +# data = db.DataBase() +# # 设置时间段 +# data.set_timePeriod(timeIndex) + +# f, spl = WGNPd(data, lat, lon, windspeed, month, freq) + +# plt.semilogx(f, spl) +# plt.xlabel('Frequency / Hz') +# plt.ylabel('Noise Level / dB') +# plt.show() \ No newline at end of file diff --git a/pyat/__init__.py b/pyat/__init__.py new file mode 100755 index 0000000..530e036 --- /dev/null +++ b/pyat/__init__.py @@ -0,0 +1,8 @@ +''' +Description: + +Author: Hunter Akins + +''' +import os, sys +sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) \ No newline at end of file diff --git a/pyat/bellhop.py b/pyat/bellhop.py new file mode 100644 index 0000000..924e556 --- /dev/null +++ b/pyat/bellhop.py @@ -0,0 +1,228 @@ +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 os import system +import numpy as np +from pyat.config import AtBinPath + + +def write_ssp(sspfile, ssp, bathm, NZmax): + if sspfile[-3:] != 'ssp': + sspfile += '.ssp' + with open(sspfile, 'w') as f: + f.write(str(len(ssp)) + '\r\n') + for val in bathm.r: + f.write(str(val) + '\t') + f.write('\r\n') + for row in range(NZmax): + for col in range(len(ssp)): + if len(ssp[col].z) <= row: + f.write('{:4.6f}'.format(1800) + '\t') + else: + f.write('{:6f}'.format(ssp[col].c[row]) + '\t') + f.write('\r\n') + return + + +def write_bathy(btyfile, bathm): + if btyfile[-3:] != 'bty': + btyfile += '.bty' + with open(btyfile, 'w') as f: + f.write('CS\r\n') + f.write(str(len(bathm.d)) + '\r\n') + for row in range(len(bathm.r)): + f.write(' ' + str(bathm.r[row]) + ' ' + str(bathm.d[row]) + '\r\n') + return + + +def calZmax(ssp): + d = np.array([], dtype=int) + for it in range(len(ssp)): + d = np.append(d, len(ssp[it].z)) + NZmax =max(d) + ssp_idx = np.argmax(d) + Zmax = ssp[ssp_idx].z[-1] + return NZmax, Zmax, ssp_idx + +def call_Bellhop(freq, sd, rd, bathm, ssp, sed, base): + # Source and receiving position + filename = 'envB' + Rmax = bathm.r[-1] # unit: km + dr = (1500 / freq) / 1000 # 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)) + # Range of phase velocity + cInt.Low = 1400 + cInt.High = 15000 + + # The number of media + NMedia = 1 + ssp_raw = [] + depth = [0] + # Sound speed profile + Z = ssp[ssp_idx].z # Depth + Cp = ssp[ssp_idx].c # Speed of P-wave + Cs = np.zeros(np.shape(Z)) # Speed of S-wave + Rho = np.ones(np.shape(Z)) # Density of the media + Ap = np.zeros(np.shape(Z)) # Attenuation of P-wave + As = np.zeros(np.shape(Z)) # Attenuation of S-wave + ssp_raw.append(SSPraw(Z, Cp, Cs, Rho, Ap, As)) + depth.append(Z[-1]) + + # Sediment is ignored + ''' + if sed != None: + Z = sed[0].z + ssp[0].z[-1] # Depth + Cp = sed[0].cp # Speed of P-wave + Cs = sed[0].cs # Speed of S-wave + Rho = sed[0].rho # Density of the media + Ap = sed[0].a_p # Attenuation of P-wave + As = sed[0].a_s # Attenuation of S-wave + ssp_raw.append(SSPraw(Z, Cp, Cs, Rho, Ap, As)) + depth.append(Z[-1]) + ''' + + Opt_top = 'SVW' + N = np.zeros(NMedia, np.int8) + Sigma = np.zeros(NMedia + 1) + sspB = SSP(ssp_raw, depth, NMedia, Opt_top, N, Sigma) + + # Bottom option + hs = HS(base[0].cp, base[0].cs, base[0].rho, base[0].a_p, base[0].a_s) + Opt_bot = 'A~' + bottom = BotBndry(Opt_bot, hs) + top = TopBndry(Opt_top) + bdy = Bndry(top, bottom) + # Beam params + run_type = 'C' # incoherently sum beams, see AT docs for more info + nbeams = 0 + alpha = np.linspace(-45, 45, 2) # min and max launch angle -20 degrees to 20 degrees + 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) + + system(AtBinPath + "/bellhop " + filename) + # Read sound field + [x, x, x, x, Pos1, pressure] = read_shd(filename + '.shd') + TL = -10 * np.log10(np.square(np.abs(pressure))) # Transmission Loss + + return Pos1, TL + +def call_Bellhop_Rays(freq, sd, rd, bathm, ssp, sed, base): + # Source and receiving position + filename = 'cz' + Rmax = bathm.r[-1] # unit: km + dr = (1500 / freq) / 1000 # wavelength / 10, unit: km + ran = np.linspace(dr, 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)) + # Range of phase velocity + cInt.Low = 1400 + cInt.High = 15000 + + # The number of media + NMedia = 1 + ssp_raw = [] + depth = [0] + # Sound speed profile + Z = ssp[ssp_idx].z # Depth + Cp = ssp[ssp_idx].c # Speed of P-wave + Cs = np.zeros(np.shape(Z)) # Speed of S-wave + Rho = np.ones(np.shape(Z)) # Density of the media + Ap = np.zeros(np.shape(Z)) # Attenuation of P-wave + As = np.zeros(np.shape(Z)) # Attenuation of S-wave + ssp_raw.append(SSPraw(Z, Cp, Cs, Rho, Ap, As)) + depth.append(Z[-1]) + + # Sediment is ignored + ''' + if sed != None: + Z = sed[0].z + ssp[0].z[-1] # Depth + Cp = sed[0].cp # Speed of P-wave + Cs = sed[0].cs # Speed of S-wave + Rho = sed[0].rho # Density of the media + Ap = sed[0].a_p # Attenuation of P-wave + As = sed[0].a_s # Attenuation of S-wave + ssp_raw.append(SSPraw(Z, Cp, Cs, Rho, Ap, As)) + depth.append(Z[-1]) + ''' + + Opt_top = 'SVW' + N = np.zeros(NMedia, np.int8) + Sigma = np.zeros(NMedia + 1) + sspB = SSP(ssp_raw, depth, NMedia, Opt_top, N, Sigma) + + # Bottom option + hs = HS(base[0].cp, base[0].cs, base[0].rho, base[0].a_p, base[0].a_s) + Opt_bot = 'A~' + bottom = BotBndry(Opt_bot, hs) + top = TopBndry(Opt_top) + bdy = Bndry(top, bottom) + # Beam params + nbeams = 301 + alpha = np.linspace(-10, 10, 2) # min and max launch angle -20 degrees to 20 degrees + 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 + run_type = 'R' + 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_bathy(filename, bathm) + + system(AtBinPath + "/bellhop " + filename) + # Read sound field + return get_rays(filename + ".ray") + +def find_cvgcRays(rays_total): + Rays = [] + num_bnc_min = np.Inf + for ray in rays_total[0]: + if ray.num_top_bnc + ray.num_bot_bnc < num_bnc_min and max(ray.xy[1,:]) > 100: + num_bnc_min = ray.num_top_bnc + ray.num_bot_bnc + + for ray in rays_total[0]: + if ray.num_top_bnc + ray.num_bot_bnc > num_bnc_min + 1: + continue + elif max(ray.xy[1,:]) < 100: + continue + else: + Rays.append(ray) + return Rays + +# if __name__=="__main__": + ''' ********************************************** ''' + ''' ******************** Test ******************** ''' + # freq = 25 # Frequency + # sd = np.array([10.0]) # Source depth + # rd = np.linspace(0, 4000, 401) # Receiving depth + # # 定义坐标 + # lat = 19.00 + # lon = 116.00 + # # 定义时间段 + # timeIndex = db.timeDict['Annual'] + # # 创建类 + # data = db.DataBase() + # # 设置时间段 + # data.set_timePeriod(timeIndex) + + # [bathm, ssp, sed, base] = GetEnv(data, freq, Coord(lon, lat), R=100, azi=90) + # [Pos, TL] = call_Bellhop(freq, sd, rd, bathm, ssp, sed, base) + + # #[x, x, x, x, Pos, pressure] = read_shd('envB.shd') + # #TL = -10 * np.log10(np.square(np.abs(pressure))) # Transmission Loss + # # Plot + # plt.pcolor(Pos.r.range, Pos.r.depth, TL[0,0,:,:], vmin=40, vmax=130,cmap='jet_r', shading='auto') + # plt.xlabel("Range/m") + # plt.ylabel("Depth/m") + # plt.gca().invert_yaxis() + # plt.colorbar() + # plt.plot(bathm.r*1000, bathm.d, color='k') + # plt.show() \ No newline at end of file diff --git a/pyat/config.json b/pyat/config.json new file mode 100644 index 0000000..801dbb5 --- /dev/null +++ b/pyat/config.json @@ -0,0 +1,3 @@ +{ + "AtBinPath": "/home/shunli/pro/pyat/at/bin" +} diff --git a/pyat/config.py b/pyat/config.py new file mode 100644 index 0000000..00cbea1 --- /dev/null +++ b/pyat/config.py @@ -0,0 +1,10 @@ +import json +from os import path + +jsonPath = path.dirname(path.abspath(__file__)) + '/config.json' +with open(jsonPath,'r') as jsonFile: + paramData = json.load(jsonFile) + if "AtBinPath" in paramData: + AtBinPath = paramData["AtBinPath"] + else: + AtBinPath = path.dirname(path.abspath(__file__))+'/at/bin' \ No newline at end of file diff --git a/pyat/cscore.py b/pyat/cscore.py new file mode 100644 index 0000000..f1c1476 --- /dev/null +++ b/pyat/cscore.py @@ -0,0 +1,32 @@ +import numpy as np + + +# 函数 +def dis_mats(m, d): # 以线阵为例 + dis = np.zeros((m, m)) + for i in range(m): + D = np.arange(0, m * d, d) - d * int(i) * np.ones(m) + dis[i] = D + dis = np.triu(dis) + dis += dis.T - np.diag(dis.diagonal()) + + return dis # 所有阵元的间距的集合 + + +def pn_funcs(f, c, dis): # 所有阵元的噪声相关性 + k = 2 * np.pi * f / c + # pn = np.sinc(k * dis) + a = np.sin(k * dis) + b = k * dis + + pn = np.divide(a, b, np.ones_like(a), where=b != 0) # sin0/0=1 + + return pn + + +def ag_funcs(m, pn): + ps = np.power(m, 2) + ag_s = ps / np.sum(pn) # + ag = 10 * np.log10(ag_s) + + return ag \ No newline at end of file diff --git a/pyat/czcore.py b/pyat/czcore.py new file mode 100644 index 0000000..c3af077 --- /dev/null +++ b/pyat/czcore.py @@ -0,0 +1,321 @@ +# -*- coding: utf-8 -*- + +from os import system +import numpy as np +from pyat.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 +import matplotlib.pyplot as plt +from pyat.config import AtBinPath + +file_cz_data = open('cz_data.txt', 'w') + +class CzRangeWidth: + def __init__(self, minRange, maxRange, width): + self.minRange = minRange + self.maxRange = maxRange + self.width = width + +class CzDetail: + def __init__(self, descript, directCz, reflectCz): + self.descript = descript + self.directCz = directCz + self.reflectCz = reflectCz + +class CzData: + def __init__(self, azi, ranges, depths, TL, bathym, Rays, CzDetail): + self.azi = azi + self.ranges = ranges + self.depths = depths + self.TL = TL + self.bathym = bathym + self.Rays = Rays + self.CzDetail = CzDetail + +def calZmax(ssp, sed): + d = np.array([],dtype=int) + for it in range(len(ssp)): + d=np.append(d,ssp[it].z[-1]) + Zmax =max(d) + ssp_idx = np.argmax(d) + return Zmax, ssp_idx + +def call_Bellhop(freq, sd, rd, bathm, ssp, sed, base): + # Source and receiving position + filename='cz' + Rmax = bathm.r[-1] # unit: km + dr = (1500 / freq) / 1000 # 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""" + RD = np.linspace(0, Zmax, round(Zmax/10) + 1) + pos = Pos(Source(sd), Dom(ran, RD)) + # Range of phase velocity + cInt.Low = 1400 + cInt.High = 15000 + + # The number of media + NMedia = 1 + ssp_raw = [] + depth = [0] + # Sound speed profile + Z = ssp[ssp_idx].z # Depth + Cp = ssp[ssp_idx].c # Speed of P-wave + Cs = np.zeros(np.shape(Z)) # Speed of S-wave + Rho = np.ones(np.shape(Z)) # Density of the media + Ap = np.zeros(np.shape(Z)) # Attenuation of P-wave + As = np.zeros(np.shape(Z)) # Attenuation of S-wave + ssp_raw.append(SSPraw(Z, Cp, Cs, Rho, Ap, As)) + depth.append(Z[-1]) + + # Sediment + if sed != None: + + Z = sed[0].z + ssp[0].z[-1] # Depth + Cp = sed[0].cp # Speed of P-wave + Cs = sed[0].cs # Speed of S-wave + Rho = sed[0].rho # Density of the media + Ap = sed[0].a_p # Attenuation of P-wave + As = sed[0].a_s # Attenuation of S-wave + ssp_raw.append(SSPraw(Z, Cp, Cs, Rho, Ap, As)) + depth.append(Z[-1]) + """Opt_top 已修改 8.26""" + Opt_top = 'SVW' + N = np.zeros(NMedia, np.int8) + Sigma = np.zeros(NMedia+1) + sspK = SSP(ssp_raw, depth, NMedia, Opt_top, N, Sigma) + + # Bottom option + hs = HS(base[0].cp, base[0].cs, base[0].rho, base[0].a_p, base[0].a_s) + Opt_bot = 'A~' + bottom = BotBndry(Opt_bot, hs) + top = TopBndry(Opt_top) + bdy = Bndry(top, bottom) + # Beam params + run_type = 'C' # incoherently sum beams, see AT docs for more info + nbeams = 0 + 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 + 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""" + # Write *.env file + write_env(filename + '.env', 'BELLHOP', 'Pekeris profile', freq, sspK, bdy, pos, beam, cInt, Rmax) + write_bathy(filename,bathm) + + + system(AtBinPath + "/bellhop " + filename) + # Read sound field + [x, x, x, x, Pos1, pressure] = read_shd(filename + '.shd') + TL = -10 * np.log10(np.square(np.abs(pressure))) # Transmission Loss + + return Pos1, TL + +def call_Bellhop_Rays(freq, sd, rd, bathm, ssp, sed, base): + # Source and receiving position + filename = 'cz_ray' + Rmax = bathm.r[-1] # unit: km + dr = (1500 / freq) / 1000 # 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""" + RD = np.linspace(0, Zmax, round(Zmax / 10) + 1) + pos = Pos(Source(sd), Dom(ran, RD)) + # Range of phase velocity + cInt.Low = 1400 + cInt.High = 15000 + + # The number of media + NMedia = 1 + ssp_raw = [] + depth = [0] + # Sound speed profile + Z = ssp[ssp_idx].z # Depth + Cp = ssp[ssp_idx].c # Speed of P-wave + Cs = np.zeros(np.shape(Z)) # Speed of S-wave + Rho = np.ones(np.shape(Z)) # Density of the media + Ap = np.zeros(np.shape(Z)) # Attenuation of P-wave + As = np.zeros(np.shape(Z)) # Attenuation of S-wave + ssp_raw.append(SSPraw(Z, Cp, Cs, Rho, Ap, As)) + depth.append(Z[-1]) + + # Sediment + if sed != None: + Z = sed[0].z + ssp[0].z[-1] # Depth + Cp = sed[0].cp # Speed of P-wave + Cs = sed[0].cs # Speed of S-wave + Rho = sed[0].rho # Density of the media + Ap = sed[0].a_p # Attenuation of P-wave + As = sed[0].a_s # Attenuation of S-wave + ssp_raw.append(SSPraw(Z, Cp, Cs, Rho, Ap, As)) + depth.append(Z[-1]) + """Opt_top 已修改 8.26""" + Opt_top = 'SVW' + N = np.zeros(NMedia, np.int8) + Sigma = np.zeros(NMedia + 1) + sspK = SSP(ssp_raw, depth, NMedia, Opt_top, N, Sigma) + + # Bottom option + hs = HS(base[0].cp, base[0].cs, base[0].rho, base[0].a_p, base[0].a_s) + Opt_bot = 'A~' + bottom = BotBndry(Opt_bot, hs) + top = TopBndry(Opt_top) + bdy = Bndry(top, bottom) + # Beam params + 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 + 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""" + # Write *.env file + write_env(filename + '.env', 'BELLHOP', 'Pekeris profile', freq, sspK, bdy, pos, beam, cInt, Rmax) + write_bathy(filename, bathm) + + system(AtBinPath + "/bellhop " + filename) + # Read sound field + return get_rays(filename + ".ray") + +def find_cvgcRays(rays_total): + Rays = [] + num_bnc_min = np.Inf + for ray in rays_total[0]: + if ray.num_top_bnc + ray.num_bot_bnc < num_bnc_min and max(ray.xy[1,:]) > 100: + num_bnc_min = ray.num_top_bnc + ray.num_bot_bnc + + for ray in rays_total[0]: + if ray.num_top_bnc + ray.num_bot_bnc > num_bnc_min + 1: + continue + elif max(ray.xy[1,:]) < 100: + continue + else: + Rays.append(ray) + return Rays + +def call_hjq(r, TL): + TL_R = 20 * np.log10(r) + out=[] + + dTL = (TL_R - TL) * (TL <= TL_R) + dTL[np.isnan(dTL)] = 0 + test_fliter_s = savgol_filter(dTL, 103, 3) + test_fliter = medfilt(test_fliter_s,13) + + signal=1 + for r_3 in range(len(r)): + if r[r_3]>30000: + break + + czrangewidth_list = [] + for i in range(r_3,len(test_fliter)): + czrangewidth = CzRangeWidth(0,0,0) + if(test_fliter[i]>=5): + out.append(i) + if i+1 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) + Rays = find_cvgcRays(rays_total) + [Pos, TL] = call_Bellhop(freq, sd, rd, bathm, ssp, sed, base) + [x, x, x, x, Pos, pressure] = read_shd('cz.shd') + TL = -10 * np.log10(np.square(np.abs(pressure))) # Transmission Loss + # Plot + if isPlot: + plt.pcolor(Pos.r.range, Pos.r.depth, TL[0, 0, :, :], vmin=40, vmax=130, cmap='jet_r', shading='auto') + plt.xlabel("Range/m") + plt.ylabel("Depth/m") + plt.gca().invert_yaxis() + plt.colorbar() + plt.show() + czrangewidth_list = call_hjq(Pos.r.range,TL[0, 0, 1,:]) + if ssp[0].c[0] > ssp[0].c[last]: + czdetail.directCz = [] + czdetail.reflectCz = czrangewidth_list + else: + czdetail.directCz = czrangewidth_list + czdetail.reflectCz = [] + return Pos.r.range, Pos.r.depth, TL[0, 0, :, :], bathm, Rays, czdetail + +def call_cz(ETOPO1Path, STPath, freq, lat, lon, azi, sourceDepth=10, receiveDepth=10, timeType='Annual', isPlot=False): + sd = np.array([sourceDepth]) # Source depth + rd = np.linspace(0, receiveDepth, 351) # Receiving depth + # 定义时间段 + timeIndex = timeType + # 创建类 + data = DataBase(ETOPO1Path, STPath) + # 设置时间段 + data.set_timePeriod(timeIndex) + + results = [] + ranges, depths, TL, bathym, Rays = [], [], [], [], [] + + ''' 计算单个方向的会聚区 ''' + if azi != None: + [bathm, ssp, sed, base] = GetEnv(data, 25, Coord(lon, lat), R=100, azi=azi) + if np.mean(bathm.d) < 2000: + return results + else: + ranges, depths, TL, bathym, Rays, czdetail = call_fangwei(freq, sd, rd, bathm, ssp, sed, base, isPlot=isPlot) + czdetail.descript = "{:.1f}°方向".format(azi) + results.append(CzData(0, ranges, depths, TL, bathym, Rays, czdetail)) + return results + + ''' 计算8个方向的会聚区 ''' + [bathm, ssp, sed, base] = GetEnv(data, 25, Coord(lon, lat), R=100, azi=0) + if bathm.d[0] < 2000: + # file_cz_data.write("各个方位都不会出现会聚区") + return results + else: + descripts = ["正北方向", "东北方向", "正东方向", "东南方向", "正南方向", "西南方向", "正西方向", "西北方向"] + theta = np.arange(0, 359, 45) + for i in range(8): + [bathm, ssp, sed, base] = GetEnv(data, 25, Coord(lon, lat), R=100, azi=theta[i]) + if np.mean(bathm.d) < 2000: + # 不会出现会聚区 + continue + ranges, depths, TL, bathym, Rays, czdetail = call_fangwei(freq, sd, rd, bathm, ssp, sed, base, isPlot=isPlot) + czdetail.descript = descripts[i] + results.append(CzData(theta[i], ranges, depths, TL, bathym, Rays, czdetail)) + + return results \ No newline at end of file diff --git a/pyat/env.py b/pyat/env.py new file mode 100755 index 0000000..825ec08 --- /dev/null +++ b/pyat/env.py @@ -0,0 +1,411 @@ +import numpy as np +# from copy import deepcopy +from matplotlib import cm +# from matplotlib.colors import ListedColormap as LCM +from scipy.interpolate import interp1d +import matplotlib.pyplot as plt + + +''' +Objects for interfacing with Michael Porter's models. +Designed in analog with the types of objects called by read_shd, write_fieldflp, write_env, etc. +''' + + +class Source: + def __init__(self, depth): + self.depth = depth + self.x = None + self.y = None + +class Dom: + def __init__(self, ran, depth, offsets=None): + self.range = ran # in km + self.depth = depth + if offsets is not None: + self.offsets = offsets + +class Pos: + def __init__(self, Source, Dom): + self.s = Source + self.r = Dom + +class Ice: + def __init__(self, BumDen, eta, xi): + self.BumDen= BumDen + self.eta = eta + self.xi = xi + +class SSPraw: + def __init__(self,z, alphaR, betaR, rho, alphaI, betaI): + """ + Input: + z - numpy 1d array + Depths the ssp is taken at + alphaR - numpy 1d array + ssp vals + betaR - numpy 1d array + shear speeds + rho - numpy 1d array + density vals at the points + alphaI - numpy 1d array + attenuation (p-wave) + betaI - numpy 1d array + shear attenuation + """ + self.z = z + self.alphaR = alphaR # sound speed in array + self.betaR = betaR # shear wave speed + self.rho = rho + self.alphaI = alphaI # atten. in array (alpha (z)) + self.betaI = betaI # shear attenuation + + def make_sspf(self): + self.sspf = interp1d(self.z, self.alphaR) + + def interp_all(self): + self.betaI_f = interp1d(self.z, self.betaI) + self.betaR_f = interp1d(self.z, self.betaR) + self.rho_f = interp1d(self.z, self.rho) + self.alphaR_f = interp1d(self.z, self.alphaR) + self.alphaI_f = interp1d(self.z, self.alphaI) + return self.alphaR_f, self.betaR_f, self.rho_f, self.alphaI_f, self.betaI_f + +class SSP: + def __init__(self, raw, depth, NMedia, Opt=None, N=None, sigma=0, ranges=0): + self.NMedia = NMedia # number of media layers + self.Opt = Opt # option: not sure what this is for + self.N = N # array with num points in each layer, don't include one for halfpace + self.sigma = sigma # ''interfacial roughness'', one for each layer? see Kuperman and ingenito 1977 for definition, include a point for the halfspace + self.depth = depth # depth array for layers + self.ranges = ranges + self.raw = raw # list of raw ssp profile + self.sspf = None + + def make_sspf(self): + for raw in self.raw: + raw.make_sspf() + if self.NMedia > 1: + self.ssp_vals = np.zeros((len(self.depth), 1)) + if self.NMedia == 2: + layer_depth = self.raw[0].z[-1] + self.sspf = lambda z: np.piecewise(z, [z=layer_depth], [lambda z: self.raw[0].sspf(z), lambda z: self.raw[1].sspf(z)]) + elif self.NMedia == 3: + layer_depth_one = self.raw[0].z[-1] + layer_depth_two = self.raw[1].z[-1] + func1 = lambda z: np.piecewise(z, [z= layer_depth_one], [lambda z: self.raw[0].sspf(z), lambda z: self.raw[1].sspf(z)]) + self.sspf = lambda z: np.piecewise(z, [z= layer_depth_two], [lambda z: func1(z), lambda z: self.raw[2].sspf(z)]) + else: + raise ValueError("too many layers") + else: + self.ssp_vals = self.raw[0].alphaR + self.sspf = self.raw[0].sspf + +class HS: + def __init__(self, alphaR=np.array([]), betaR=np.array([]), rho=np.array([]), alphaI=np.array([]), betaI=np.array([])): + self.alphaR = np.array(alphaR) + self.betaR = np.array(betaR) + self.rho = np.array(rho) + self.alphaI = np.array(alphaI) + self.betaI = np.array(betaI) + +class BotBndry: + def __init__(self, Opt, Hs, depth=[], betaI=[0], ): + self.Opt = Opt # 'A' for analytic or 'CVW' for interpolated ssp + self.hs = Hs + +class TopBndry: + def __init__(self, Opt, depth=[]): + self.Opt = Opt + self.cp = None + self.cs = None + self.rho = None + +class Bndry: + def __init__(self, top, bot): + self.Top = top + self.Bot = bot + +class Box: + def __init__(self, z, r): + self.r = r + self.z = z + +class Beam: + def __init__(self, RunType=None, Type=None,Nbeams=None, Ibeam=None, Nrays=None, alpha=None, deltas=None, box=None, epmult=None, rloop=None, Ibwin=None, Nimage = None): + self.RunType = RunType + self.Type = Type + self.Nbeams = Nbeams + self.Ibeam = Ibeam + self.Nrays = Nrays + self.alpha = alpha + self.deltas = deltas + self.box = box + self.epmult = epmult + self.rloop = rloop + self.Ibwin = Ibwin + self.Nimage = None + +class cInt: + def __init__(self, low, high): + self.Low = low + self.High = high + +class Field: + def __init__(self, **kwargs): + self.pos = kwargs['pos'] + self.greens = kwargs['pressure'] + self.freq = kwargs['freq'] # a list of frequencies + + def contour(self): + shape = np.shape(self.greens) + if len(shape) > 2: + raise ValueError("Field object is too high dimensional, try multi_contour") + pos = self.pos + ax = plt.contourf(pos.r.range, -pos.r.depth, np.log10(abs(self.greens.real))) + return ax + + def multi_contour(self): + print('not implemented, maybe one day...') + +class Modes: + def __init__(self, **kwargs): + self.M = kwargs['M'] + self.k = kwargs['modes_k'] + self.z = np.array(kwargs['z']) + self.phi = kwargs['modes_phi'] + self.top = kwargs['top'] + self.bot = kwargs['bot'] + self.N = kwargs['N'] + self.Nfreq = kwargs['Nfreq'] + self.Nmedia = kwargs['Nmedia'] + self.depth = kwargs['depth'] + self.rho = kwargs['rho'] + self.freqvec = kwargs['freqVec'] + self.init_dict = kwargs + self.num_modes = self.M # easier to remember + + def get_excited_modes(self, sd, threshold): + ''' + return an array of modes that are excited by a source at sd meters depth + threshold + + also populate some structures in moded caled excited_phi and excited_k + + ''' + if sd not in self.z: + raise ValueError("sd not in the depth array, are you sure it's the right depth you're passing?") + depth_ind = [i for i in range(len(self.z)) if self.z[i] == sd][0] + vals = self.phi[depth_ind,:] + const = np.max(abs(vals)) + filtered_inds = [i for i in range(len(self.k)) if abs(self.phi[depth_ind, i]) / const > threshold] + self.excited_phi = self.phi[:, filtered_inds] + self.excited_k = self.k[filtered_inds] + return self.excited_phi, self.excited_k + + def get_source_depth_ind(self, sd): + """ + sd is an int + """ + tol = 1e-2 + sind = [i for i in range(len(self.z)) if abs(self.z[i]-sd) < tol] + if len(sind) == 0 : + raise ValueError("sd not in the depth array, are you sure it's the right depth you're passing?") + else: + self.sind = sind[0] + return self.sind + + def remove_source_pos(self, sd): + """ + Take the source at sd from the mode matrix + Initiate a new attribute to hold the source + modal value + """ + sind = self.get_source_depth_ind(sd) + new_pos_len = len(self.z) - 1 + new_phi = np.zeros((new_pos_len, self.num_modes), dtype=self.phi.dtype) + new_phi[:sind, :] = self.phi[:sind, :] + new_phi[sind:,:] = self.phi[sind+1:,:] + self.phi = new_phi + new_z = np.zeros((new_pos_len), dtype=self.z.dtype) + new_z[:sind] = self.z[:sind] + new_z[sind:] = self.z[sind+1:] + self.z = new_z + self.source_strength = self.phi[sind,:] + return + + def get_source_strength(self, sd): + """ + Get the value of each mode at the source depth sd (meters) + Initialize new attribute for the source strength + """ + sind = self.get_source_depth_ind(sd) + vals = self.phi[sind,:] + self.source_strength = vals + return vals + + def get_receiver_modes(self, zr): + """ + zr is array like + """ + tol = 1e-3 + r_inds = [i for i in range(len(self.z)) if np.min(abs(self.z[i]-zr)) < tol] + receiver_modes = self.phi[r_inds, :] + self.receiver_modes = receiver_modes + return receiver_modes + + def get_source_excitation(self, zs): + """ + For case where there is a track, there may be multiple repeats in zs + """ + tol = 1e-3 + r_inds = [np.argmin(abs(zs[i] - self.z)) for i in range(len(zs))] + strength_modes = self.phi[r_inds, :] + self.strength_modes = strength_modes + return strength_modes + + + def plot(self): + figs = [] + if self.M > 5: + for i in range(self.M): + fig = plt.figure(i) + plt.plot(self.phi[:,i], -self.z) + figs.append(fig) + return figs + + def __repr__(self): + return 'Modes object with ' + str(self.M) + ' distinct modes' + +class Arrival: + def __init__(self, info_list): + self.amp = info_list[0] + self.delay = info_list[1] + self.src_ang = info_list[2] + self.rec_ang = info_list[3] + self.num_top_bnc = info_list[4] + self.num_bot_bnc = info_list[5] + +class Arrivals: + def __init__(self, arrival_list): + self.arrivals = arrival_list + + def plot_cir(self, vals=None,ax=None): + """ + Plot the channel impulse response + Vals provides the option to pass in a set of arrivals to use fo calibrating the y axis + If vals=None, I use the arrivals in the object to calibrate the y-axis + I don't create a figure + I return ymin, ymax, vals so that the values from this arrival set can be used in a more global + scope to set the y-axis limits (useful for creating videos comprised of snapshots of a bunch + of CIR's). For example, I want to plot a bunch of Arrivals on the same plot, or combine a bunch of + CIRs into a video. Then I need to have the same axis limits for all the plots. + I can use default vals (None) for the first plot, then use the ymin, ymax, vals return values + in a larger scope to standardize the axes. Then for the second, third, ... plots, I pass in vals + to make sure they all have the same time axis. + Output - + ymin - float + ymax -float + plus/minus 1.5 times max amplitude value + vals - list + vals[0] is the times + vals[1] is the amplitude vals + """ + arrival_list = self.arrivals + amps = np.array([abs(x.amp) for x in arrival_list]) + phase = np.array([np.angle(x.amp) for x in arrival_list]) + times = np.array([x.delay for x in arrival_list]) + """ + Create the time axis + """ + if type(vals) != type(None): + t = np.linspacea(.9*np.min(vals[0]), 1.1*np.max(vals[0]), 100) + else: + t = np.linspace(.9*np.min(times), 1.1*np.max(times), 100) + zeros = np.zeros(t.size) + if type(ax) == type(None): + fig, ax = plt.subplots(2,1) + ax.scatter(t.real, zeros, s=6) + ax.stem(times.real, amps.real, markerfmt=' ', basefmt=' ',use_line_collection=True) + ax.scatter(times.real, amps,s=16,c='k') + scale = 1.5*np.max(abs(amps)) # set yaxis scale + ymin = - scale + ymax = scale + vals = [times.real, amps] + return ymin, ymax, vals + + def plot_ellipse(self, vals=None, ref_amp =None, src_ang=False, ax=None): + """ + Produce the travel time ellipse with colorbar for amplitude + Input - + self - contains all the arrival info + vals - optional + allows me to use a reference arrival set to calibrate the plot axes + vals[0] = times + vals[0] = amps + """ + arrival_list = self.arrivals + amps = np.array([abs(x.amp) for x in arrival_list]) + times = np.array([x.delay.real for x in arrival_list]) + rec_angles = np.array([x.rec_ang for x in arrival_list]) + src_angles = np.array([x.rec_ang for x in arrival_list]) + """ + Create the time axis + """ + if type(vals) != type(None): + tmin, tmax = np.min(vals[0]), np.max(vals[0]) + else: + tmin, tmax = np.min(times), np.max(times) + if type(ref_amp) == type(None): + amps = amps / np.max(abs(amps)) + else: + amps /= abs(ref_amp) + cvals = 10*np.log10(abs(amps)) + #cvals = np.linspace(-5, 0, len(amps)) + """ Make a nice cmap """ + cmap = cm.get_cmap('Spectral') + """ Configure the x axis """ + + if type(ax) == type(None): + fig, ax = plt.subplots(2,1) + ax.scatter([tmin, tmax], [0,0], s=0) + if src_ang == False: + angles = rec_angles + else: + angles = src_angles + max_db_down = 10 + cvals_normed = abs(cvals) / max_db_down + strong_inds = cvals > -max_db_down + ax.scatter(times[strong_inds], angles[strong_inds], s=25, c=cmap(cvals_normed[strong_inds])) + vals = [times, amps] + ymin, ymax = -1.5*np.max(abs(angles)), 1.5*np.max(abs(angles)) + return ymin, ymax, vals + +class KernInput: + def __init__(self, Field_r, Field_s, env): + self.gr = Field_r.greens_mat + self.num_rcvrs = len(self.gr) + self.gs = np.flip(Field_s.greens_mat, axis=1) # move source off axis + self.env = env + self.Pos = self.Field_r.Pos + + +class Eigenray: + def __init__(self, launch_ang, num_top_bnc, num_bot_bnc, xy_arr): + """ + launch_ang + launch angl in degrees + num_top_bnc + integer + num_bot_bnc + integer + xy_arr - np array + 2 by num steps, first is x position, second is y + (so plt.plot(xy_arr[0,:], xy_arr[1,:]) will plot the ray) + """ + self.src_ang = launch_ang + self.num_top_bnc = num_top_bnc + self.num_bot_bnc = num_bot_bnc + self.xy = xy_arr + + diff --git a/pyat/krakenc.py b/pyat/krakenc.py new file mode 100755 index 0000000..6d3226b --- /dev/null +++ b/pyat/krakenc.py @@ -0,0 +1,139 @@ +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 + +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_krakenc(freq, sd, rd, bathm, ssp, sed, base): + filename = "envK" + envFilename = "envK.env" + model = "KRAKEN" + title = "Krakenc Waveguide" + + # Source and receiving position + Rmax = bathm.r[-1] # 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 + + # Write environmental file + for it in range(len(ssp)): + # The number of media + if sed != None: + NMedia = 2 + else: + NMedia = 1 + + ssp_raw = [] + depth = [0] + # Sound speed profile + Z = ssp[it].z # Depth + Cp = ssp[it].c # Speed of P-wave + Cs = np.zeros(np.shape(Z)) # Speed of S-wave + Rho = np.ones(np.shape(Z)) # Density of the media + Ap = np.zeros(np.shape(Z)) # Attenuation of P-wave + As = np.zeros(np.shape(Z)) # Attenuation of S-wave + ssp_raw.append(SSPraw(Z, Cp, Cs, Rho, Ap, As)) + depth.append(Z[-1]) + + # Sediment + if sed != None: + Z = sed[it].z + ssp[it].z[-1] # Depth + Cp = sed[it].cp # Speed of P-wave + Cs = sed[it].cs # Speed of S-wave + Rho = sed[it].rho # Density of the media + Ap = sed[it].a_p # Attenuation of P-wave + As = sed[it].a_s # Attenuation of S-wave + ssp_raw.append(SSPraw(Z, Cp, Cs, Rho, Ap, As)) + depth.append(Z[-1]) + + # Make sure that the base depth is consistent in all profiles + if depth[-1] != Zmax: + Z = [depth[-1], Zmax] # Depth + Cp = [Cp[-1], base[it].cp] # Speed of P-wave + Cs = [Cs[-1], base[it].cs] # Speed of S-wave + Rho = [Rho[-1], base[it].rho] # Density of the media + Ap = [Ap[-1], base[it].a_p] # Attenuation of P-wave + As = [As[-1], base[it].a_s] # Attenuation of S-wave + ssp_raw.append(SSPraw(Z, Cp, Cs, Rho, Ap, As)) + depth.append(Z[-1]) + NMedia += 1 + + Opt_top = 'SVW .' + N = np.zeros(NMedia, np.int8) + Sigma = np.zeros(NMedia+1) + sspK = SSP(ssp_raw, depth, NMedia, Opt_top, N, Sigma) + + # Bottom option + hs = HS(base[it].cp, base[it].cs, base[it].rho, base[it].a_p, base[it].a_s) + Opt_bot = 'A' + bottom = BotBndry(Opt_bot, hs) + top = TopBndry(Opt_top) + bdy = Bndry(top, bottom) + + # Write *.env file + if it == 0: + write_env(envFilename, model, title, freq, sspK, bdy, pos, [], cInt, Rmax) + else: + write_env(envFilename, model, title, freq, sspK, bdy, pos, [], cInt, Rmax, None, 'a') + + # Call krakenc.exe + system(AtBinPath + "/krakenc " + filename) + # Write *.flp file + Opt_flp = 'RC' + pos.r.depth = rd + write_fieldflp(filename, Opt_flp, pos, NProf=len(ssp), rProf=bathm.r) + # Call field.exe + system(AtBinPath + "/field " + filename) + # Read sound field + [x, x, x, x, Pos1, pressure] = read_shd(filename+'.shd') + TL = -10 * np.log10(np.square(np.abs(pressure))) # Transmission Loss + + return Pos1, TL + +# if __name__=="__main__": + ''' ********************************************************************************* ''' + ''' ***************************** Input Parameters ********************************** ''' + ''' ********************************************************************************* ''' + # freq = 25 # Frequency + # sd = np.array([40.0]) # Source depth + # rd = np.linspace(0, 500, 501) # Receiving depth + + # # 定义坐标 + # lat = 18.00 + # lon = 114.00 + # # 定义时间段 + # timeIndex = db.timeDict['Annual'] + # # 创建类 + # data = db.DataBase() + # # 设置时间段 + # data.set_timePeriod(timeIndex) + + # [bathm, ssp, sed, base] = GetEnv(data, 25, Coord(lon, lat), R=50, azi=0) + + # [Pos, TL] = call_krakenc(freq, sd, rd, bathm, ssp, sed, base) + + # # Plot + # plt.pcolor(Pos.r.range, Pos.r.depth, TL[0,0,:,:], vmin=40, vmax=130, cmap='jet_r', shading='auto') + # plt.xlabel("Range/m") + # plt.ylabel("Depth/m") + # plt.gca().invert_yaxis() + # plt.colorbar() + # plt.show() \ No newline at end of file diff --git a/pyat/ram.py b/pyat/ram.py new file mode 100755 index 0000000..d6f18ce --- /dev/null +++ b/pyat/ram.py @@ -0,0 +1,72 @@ +from pyat.ramRW import PosRAM, PROF, SSPRAM, write_env_RAM, readTL_RAM +from os import system +import numpy as np +from pyat.config import AtBinPath + +def call_ram(freq, sd, rd, bathm, ssp, sed): + ''' ********************************************************************************* ''' + ''' ************************** Sound Field Calculation ****************************** ''' + ''' ********************************************************************************* ''' + title = "Waveguide" + # Source and receiving position + bathm.r = bathm.r * 1000 # km -> m + 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 + pos = PosRAM(sd, rd, rmax, zmax, dr, dz, zmplt) + + # Sound speed profile for RAM + ssp_ram = [] + for it in range(len(ssp)): + rp = bathm.r[it] # Range of this profile + Z = sed[it].z + cw_p = PROF(ssp[it].z, ssp[it].c) # Water sound speed profile + cb_p = PROF(Z, sed[it].cp) # Bottom sound speed profile + rhob_p = PROF(Z, sed[it].rho) # Sea bottom density profile + Z = np.append(Z, Z[-1] + 100) + Ap = np.append(sed[it].a_p, 10) + # Ap = sed[it].a_p + attn_p = PROF(Z, Ap) # Sea bottom attenuation profile + + ssp_ram.append(SSPRAM(rp, cw_p, cb_p, rhob_p, attn_p)) + + + # Write ram.in file + write_env_RAM(title, freq, pos, bathm, ssp_ram) + # Call ram.exe + system(AtBinPath + "/ram") + # Read tl.grid + return readTL_RAM() + +# if __name__=="__main__": + # ''' ********************************************************************************* ''' + # ''' ***************************** Input Parameters ********************************** ''' + # ''' ********************************************************************************* ''' + # freq = 25 # Frequency + # sd = 40.0 # Source depth + # rd = 30.0 # Receiving depth + + # # 定义坐标 + # lat = 18.00 + # lon = 114.00 + # # 定义时间段 + # timeIndex = db.timeDict['Annual'] + # # 创建类 + # data = db.DataBase() + # # 设置时间段 + # data.set_timePeriod(timeIndex) + + # [bathm, ssp, sed, base] = GetEnv(data, 25, Coord(lon, lat), R=100, azi=5) + + # # 调用ram模块 + # [rd, rr, TL] = call_ram(freq, sd, rd, bathm, ssp, sed) + + # # Plot + # plt.pcolormesh(rr, rd, TL, vmin=40, vmax=90, cmap='jet_r', shading='auto') + # plt.xlabel("Range/m") + # plt.ylabel("Depth/m") + # plt.gca().invert_yaxis() + # plt.colorbar() + # plt.show() \ No newline at end of file diff --git a/pyat/ramRW.py b/pyat/ramRW.py new file mode 100755 index 0000000..4e5c97d --- /dev/null +++ b/pyat/ramRW.py @@ -0,0 +1,151 @@ +from pyat.env import * +import numpy as np + +def write_env_RAM(TitleEnv, freq, pos, bathym, ssp, envfil = "ram.in"): + # envfil = "ram.in" + Rmax = pos.rmax + + f = open(envfil, "wt") + + ndr = 1 + ndz = 1 + + zmplt = pos.zmplt + + c0 = np.mean(ssp[0].cw.v) + + nPade = 8 # The number of Pade terms + nStarter = 4 # The number of starter terms + rs = 0.0 # The stability constraints range + + f.write("{:s}".format(TitleEnv) + " ! Title\r\n") + f.write("{:8.2f} ".format(freq) + \ + "{:8.2f} ".format(pos.zs) + \ + "{:8.2f} ".format(pos.zr) + \ + "\t \t \t ! Frequency (Hz) zs zr\r\n") + f.write("{:8.2f} ".format(Rmax) + \ + "{:8.2f} ".format(pos.dr) + \ + "{:3d} ".format(ndr) + \ + "\t \t \t ! rmax dr ndr \r\n") + f.write("{:8.2f} ".format(pos.zmax) + \ + "{:8.2f} ".format(pos.dz) + \ + "{:3d} ".format(ndz) + \ + "{:8.2f} ".format(zmplt) + \ + "\t \t \t ! zmax dz ndz zmplt\r\n") + f.write("{:8.2f} ".format(c0) + \ + "{:1d} ".format(nPade) + \ + "{:1d} ".format(nStarter) + \ + "{:8.2f} ".format(rs) + \ + "\t \t \t ! c0 np ns rs\r\n") + + # Bathymetry + for it in range(len(bathym.r)): + f.write("{:8.2f} ".format(bathym.r[it]) + "{:8.2f} ".format(bathym.d[it]) + \ + " \t \t \t ! rb zb\r\n") + + # Sound speed profile + f.write("-1 -1\r\n") + for it in range(len(ssp)): + # rp + if it != 0: + f.write("{:8.2f} \t \t \t ! rp\r\n".format(ssp[it].rp)) + # cw + for ii in range(len(ssp[it].cw.z)): + f.write("{:8.2f} ".format(ssp[it].cw.z[ii]) + "{:8.2f} ".format(ssp[it].cw.v[ii]) + \ + " \t \t \t ! z cw\r\n") + f.write("-1 -1\r\n") + + # cb + for ii in range(len(ssp[it].cb.z)): + f.write("{:8.2f} ".format(ssp[it].cb.z[ii]) + "{:8.2f} ".format(ssp[it].cb.v[ii]) + \ + " \t \t \t ! z cb\r\n") + f.write("-1 -1\r\n") + + # rhob + for ii in range(len(ssp[it].rhob.z)): + f.write("{:8.2f} ".format(ssp[it].rhob.z[ii]) + "{:8.2f} ".format(ssp[it].rhob.v[ii]) + \ + " \t \t \t ! z rhob\r\n") + f.write("-1 -1\r\n") + + # attn + for ii in range(len(ssp[it].attn.z)): + f.write("{:8.2f} ".format(ssp[it].attn.z[ii]) + "{:8.2f} ".format(ssp[it].attn.v[ii]) + \ + " \t \t \t ! z attn\r\n") + f.write("-1 -1\r\n") + + +def readTL_RAM(filePath=None): + if filePath != None: + f = open(filePath + "\\tl.grid") + else: + f = open("tl.grid") + + rec_marker1 = np.fromfile(f, dtype=np.int32, count=1) + if rec_marker1 == 60: + datatype = np.int32 + else: + rec_maker1 = np.fromfile(f, dtype=np.int32, count=1) + datatype = np.int64 + + + freq = np.fromfile(f, dtype=np.float32, count=1) + zs = np.fromfile(f, dtype=np.float32, count=1) + zr = np.fromfile(f, dtype=np.float32, count=1) + rmax = np.fromfile(f, dtype=np.float32, count=1) + dr = np.fromfile(f, dtype=np.float32, count=1) + ndr = np.fromfile(f, dtype=np.int32, count=1) + zmax = np.fromfile(f, dtype=np.float32, count=1) + dz = np.fromfile(f, dtype=np.float32, count=1) + ndz = np.fromfile(f, dtype=np.int32, count=1) + zmplt = np.fromfile(f, dtype=np.float32, count=1) + c0 = np.fromfile(f, dtype=np.float32, count=1) + nPade = np.fromfile(f, dtype=np.int32, count=1) + nStarter = np.fromfile(f, dtype=np.int32, count=1) + rs = np.fromfile(f, dtype=np.float32, count=1) + lz = np.fromfile(f, dtype=np.int32, count=1) + + rec_marker2 = np.fromfile(f, dtype=datatype, count=1) + + rr = np.arange(ndr*dr, rmax+ndr*dr, ndr*dr) + rd = ndz * dz * (np.arange(lz) + 1) + Nrr = len(rr) + Nrd = len(rd) + + tl = np.zeros([Nrr, Nrd],dtype=np.float32) + for it in range(Nrr): + marker1 = np.fromfile(f, dtype=datatype, count=1) + tl[it][:] = np.fromfile(f, dtype=np.float32, count=Nrd) + marker2 = np.fromfile(f, dtype=datatype, count=1) + + tl = np.transpose(tl) + + return [rd, rr, tl] + + +class PosRAM: + def __init__(self, zs, zr, rmax, zmax, dr, dz, zmplt): + self.zs = zs + self.zr = zr + self.rmax = rmax + self.zmax = zmax + self.dr = dr + self.dz = dz + self.zmplt = zmplt + + +class SSPRAM: + def __init__(self, rp, cw_p, cb_p, rhob_p, attn_p): + self.rp = rp + self.cw = cw_p + self.cb = cb_p + self.rhob = rhob_p + self.attn = attn_p + + +class PROF: + def __init__(self, depth, value): + self.z = depth + self.v = value + + + diff --git a/pyat/ramRW_new.py b/pyat/ramRW_new.py new file mode 100644 index 0000000..ea32dc6 --- /dev/null +++ b/pyat/ramRW_new.py @@ -0,0 +1,168 @@ +from pyat.env import * +import numpy as np + + +def write_env_RAM(TitleEnv, envfil, freq, pos, bathym, ssp): + Rmax = pos.rmax + + f = open(envfil, "wt") + + ndr = 1 + ndz = 1 + + zmplt = pos.zmplt + + c0 = np.mean(ssp[0].cw.v) + + nPade = 3 # The number of Pade terms + nStarter = 4 # The number of starter terms + rs = 0.0 # The stability constraints range + + f.write("{:s}".format(TitleEnv) + " ! Title\r\n") + f.write("{:8.2f} ".format(freq) + \ + "{:8.2f} ".format(pos.zs) + \ + "{:8.2f} ".format(pos.zr) + \ + "\t \t \t ! Frequency (Hz) zs zr\r\n") + f.write("{:8.2f} ".format(Rmax) + \ + "{:8.2f} ".format(pos.dr) + \ + "{:3d} ".format(ndr) + \ + "\t \t \t ! rmax dr ndr \r\n") + f.write("{:8.2f} ".format(pos.zmax) + \ + "{:8.2f} ".format(pos.dz) + \ + "{:3d} ".format(ndz) + \ + "{:8.2f} ".format(zmplt) + \ + "\t \t \t ! zmax dz ndz zmplt\r\n") + f.write("{:8.2f} ".format(c0) + \ + "{:1d} ".format(nPade) + \ + "{:1d} ".format(nStarter) + \ + "{:8.2f} ".format(rs) + \ + "\t \t \t ! c0 np ns rs\r\n") + + # Bathymetry + for it in range(len(bathym.r)): + f.write("{:8.2f} ".format(bathym.r[it]) + "{:8.2f} ".format(bathym.d[it]) + \ + " \t \t \t ! rb zb\r\n") + + # Sound speed profile + f.write("-1 -1\r\n") + for it in range(len(ssp)): + # rp + if it != 0: + f.write("{:8.2f} \t \t \t ! rp\r\n".format(ssp[it].rp)) + # cw + for ii in range(len(ssp[it].cw.z)): + f.write("{:8.2f} ".format(ssp[it].cw.z[ii]) + "{:8.2f} ".format(ssp[it].cw.v[ii]) + \ + " \t \t \t ! z cw\r\n") + f.write("-1 -1\r\n") + + # cb + for ii in range(len(ssp[it].cb.z)): + f.write("{:8.2f} ".format(ssp[it].cb.z[ii]) + "{:8.2f} ".format(ssp[it].cb.v[ii]) + \ + " \t \t \t ! z cb\r\n") + f.write("-1 -1\r\n") + + # rhob + for ii in range(len(ssp[it].rhob.z)): + f.write("{:8.2f} ".format(ssp[it].rhob.z[ii]) + "{:8.2f} ".format(ssp[it].rhob.v[ii]) + \ + " \t \t \t ! z rhob\r\n") + f.write("-1 -1\r\n") + + # attn + for ii in range(len(ssp[it].attn.z)): + f.write("{:8.2f} ".format(ssp[it].attn.z[ii]) + "{:8.2f} ".format(ssp[it].attn.v[ii]) + \ + " \t \t \t ! z attn\r\n") + f.write("-1 -1\r\n") + + +def readTL_RAM(filePath=None): + if filePath != None: + f = open(filePath + "\\tl.grid") + else: + f = open("tl.grid") + + rec_marker1 = np.fromfile(f, dtype=np.int32, count=1) + if rec_marker1 == 60: + datatype = np.int32 + else: + rec_maker1 = np.fromfile(f, dtype=np.int32, count=1) + datatype = np.int64 + + + freq = np.fromfile(f, dtype=np.float32, count=1) + zs = np.fromfile(f, dtype=np.float32, count=1) + zr = np.fromfile(f, dtype=np.float32, count=1) + rmax = np.fromfile(f, dtype=np.float32, count=1) + dr = np.fromfile(f, dtype=np.float32, count=1) + ndr = np.fromfile(f, dtype=np.int32, count=1) + zmax = np.fromfile(f, dtype=np.float32, count=1) + dz = np.fromfile(f, dtype=np.float32, count=1) + ndz = np.fromfile(f, dtype=np.int32, count=1) + zmplt = np.fromfile(f, dtype=np.float32, count=1) + c0 = np.fromfile(f, dtype=np.float32, count=1) + nPade = np.fromfile(f, dtype=np.int32, count=1) + nStarter = np.fromfile(f, dtype=np.int32, count=1) + rs = np.fromfile(f, dtype=np.float32, count=1) + lz = np.fromfile(f, dtype=np.int32, count=1) + + rec_marker2 = np.fromfile(f, dtype=datatype, count=1) + + rr = np.arange(ndr*dr, rmax+ndr*dr, ndr*dr) + rd = ndz * dz * (np.arange(lz) + 1) + Nrr = len(rr) + Nrd = len(rd) + + tl = np.zeros([Nrr, Nrd],dtype=np.float32) + for it in range(Nrr): + marker1 = np.fromfile(f, dtype=datatype, count=1) + tl[it][:] = np.fromfile(f, dtype=np.float32, count=Nrd) + marker2 = np.fromfile(f, dtype=datatype, count=1) + + tl = np.transpose(tl) + + return [rd, rr, tl] + +def readTLline_RAM(filePath=None): + if filePath != None: + f = open(filePath + "/tl.line") + else: + f = open("tl.line") + + lines = f.readlines() + f.close() + rr = np.zeros(len(lines)) + tl = np.zeros(len(lines)) + for i in range(len(lines)): + data = lines[i].strip().split(' ') + rr[i] = float(data[0]) + tl[i] = float(data[-1]) + + return rr, tl + + +class PosRAM: + def __init__(self, zs, zr, rmax, zmax, dr, dz, zmplt): + self.zs = zs + self.zr = zr + self.rmax = rmax + self.zmax = zmax + self.dr = dr + self.dz = dz + self.zmplt = zmplt + + +class SSPRAM: + def __init__(self, rp, cw_p, cb_p, rhob_p, attn_p): + self.rp = rp + self.cw = cw_p + self.cb = cb_p + self.rhob = rhob_p + self.attn = attn_p + + +class PROF: + def __init__(self, depth, value): + self.z = depth + self.v = value + + + diff --git a/pyat/readwrite.py b/pyat/readwrite.py new file mode 100755 index 0000000..f1cc918 --- /dev/null +++ b/pyat/readwrite.py @@ -0,0 +1,1507 @@ +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 +import re +from matplotlib import pyplot as plt +from math import floor + +""" + translation of basic read functions from Mike Porter's + Matlab ReadWrite folder + Currently lacking bellhop support (just need to translate read_bell) + and bellhop3d (need to translate readsxrx) +""" + + +def equally_spaced( x ): + # routine to see if the vector x is composed of roughly equally spaced + # values + n = len( x ) + xtemp = np.linspace( x[ 0 ], x[ -1 ], n ) + delta = abs( x[:] - xtemp[:] ) # the colon converts both to column vectors + + if ( np.max( delta ) < 1e-9 ): + return 1 + else: + return 0 + +# Write a field-parameters file +def write_fieldflp(flpfil, Option, Pos, **kwargs): + if flpfil[-4:] != '.flp': + flpfil = flpfil + '.flp' + + if 'scooter' in kwargs.keys(): + with open(flpfil, 'w') as f: + f.write('\'' + '{:4}'.format(Option) + ' \'' + ' ! Option \r\n') + f.write('{:5d}'.format(len(Pos.r.range)) + ' \t \t \t \t ! NRR') + if ( len( Pos.r.range ) > 2) and equally_spaced( Pos.r.range ): + f.write('\r\n ' + '{:6f}'.format(Pos.r.range[0])) + f.write('\r\n ' + '{:6f}'.format(Pos.r.range[-1])) + elif len(Pos.r.range) == 1: + f.write('\r\n ' + '{:6f}'.format(Pos.r.range[0]) + ' ') + else: + for i in range(len(Pos.r.range)): + f.write('\r\n ' + '{:6f}'.format(Pos.r.range[i]) + ' ') + + + f.write( '/ \t ! RR(1) ... (km) \r\n' ) + return + with open(flpfil, 'w') as f: + f.write('/ ! Title \r\n' ) + f.write('\'' + '{:4}'.format(Option) + ' \'' + ' ! Option \r\n') + f.write('999999 ! Mlimit (number of modes to include) \r\n' ) + if 'NProf' in kwargs.keys(): + Nprof = kwargs['NProf'] + f.write(str(Nprof)+' ! NProf' ) + rProf = kwargs['rProf'] + if ( len( rProf ) > 2) and equally_spaced( rProf ): + f.write('\r\n ' + '{:6f}'.format(rProf[0])) + f.write('\r\n ' + '{:6f}'.format(rProf[-1])) + elif len(rProf) == 1: + f.write('\r\n ' + '{:6f}'.format(rProf[0]) + ' ') + else: + for i in range(len(rProf)): + f.write('\r\n ' + '{:6f}'.format(rProf[i]) + ' ') + f.write(' / ! rProf (km) \r\n') + else: + f.write('1 ! NProf \r\n' ) + f.write('0.0 / ! rProf (km) \r\n' ) + + # receiver ranges + f.write('{:5d}'.format(len(Pos.r.range)) + ' \t \t \t \t ! NRR') + + if ( len( Pos.r.range ) > 2) and equally_spaced( Pos.r.range ): + f.write('\r\n ' + '{:6f}'.format(Pos.r.range[0])) + f.write('\r\n ' + '{:6f}'.format(Pos.r.range[-1])) + elif len(Pos.r.range) == 1: + f.write('\r\n ' + '{:6f}'.format(Pos.r.range[0]) + ' ') + else: + for i in range(len(Pos.r.range)): + f.write('\r\n ' + '{:6f}'.format(Pos.r.range[i]) + ' ') + + + f.write( ' / \t ! R(1) ... (km) \r\n' ) + + # source depths + + f.write('{:5d}'.format(len(Pos.s.depth))+'\t \t \t \t ! NSD') + + if ( len( Pos.s.depth ) > 2 and equally_spaced( Pos.s.depth ) ): + f.write('\r\n '+'{:6f}'.format(Pos.s.depth[ 0 ])+' ') + f.write('\r\n '+'{:6f}'.format(Pos.s.depth[ -1 ])+' ') + elif len(Pos.s.depth) == 1: + f.write( '\r\n '+ '{:6f}'.format(Pos.s.depth[0])+' ') + else: + for i in range(len(Pos.s.depth)): + f.write( '\r\n '+ '{:6f}'.format(Pos.s.depth[i])+' ') + + f.write( '/ \t ! SD(1) ... (m) \r\n' ) + + # receiver depths + + f.write('{:5}'.format(len(Pos.r.depth)) + ' \t \t \t \t ! NRD') + + if ( len( Pos.r.depth ) > 2 and equally_spaced( Pos.r.depth ) ): + f.write( '\r\n '+'{:6f}'.format(Pos.r.depth[0])+' ') + f.write( '\r\n '+'{:6f}'.format(Pos.r.depth[-1])+' ') + elif (len(Pos.r.depth) == 1): + f.write( '\r\n '+'{:6f}'.format(Pos.r.depth[0])+' ') + else: + for i in range(len(Pos.r.depth)): + f.write( '\r\n '+'{:6f}'.format(Pos.r.depth[i])) + + f.write( '/ \t ! RD(1) ... (m) \r\n' ) + + # receiver range offsets + f.write( '{:5}'.format(len( Pos.r.depth ))+' \t \t ! NRR', ) + if not hasattr(Pos.r, 'offsets'): + f.write( '\r\n '+'{:6.2f}'.format(0.0) +' ') # just zeros + else: + if ( len( Pos.r.offsets) > 2 and equally_spaced( Pos.r.offsets ) ): + f.write( '\r\n '+'{:6f}'.format(Pos.r.offsets[1])+' ') + f.write( '\r\n '+'{:6f}'.format(Pos.r.offsets[-1])+' ') + elif (len(Pos.r.offsets) == 1): + f.write( '\r\n '+'{:6f}'.format(Pos.r.offsets[0])+' ') + else: + for i in range(len(Pos.r.offsets)): + f.write( '\r\n '+'{:6f}'.format(Pos.r.offsets[i])) + f.write( '/ \t \t ! RR(1) ... (m) \r\n' ) + return + +def fileparts(fname): + fpath = os.path.dirname(os.path.abspath(fname)) + if '.' not in fname: + ext = None + else: + ext_ind = [i for i in range(len(fname)) if fname[i]=='.'] + if len(ext_ind) > 1: + print('Warning: lots of periods...') + #raise ValueError("something fishy about a filename with two periods") + ext_ind = ext_ind[-1] + ext = fname[ext_ind:] + #print('ext', ext) + fname = fname[:ext_ind] + return fpath, fname, ext + +class Empty(): + def __init__(self): + return + +def write_sbp(fname, pattern,num_dec_pts=2): + """ + From BELLHOP docs, + ``The source beam pattern file has the format + NSBPPts + angle1 power1 + .... + with angle following the BELLHOP convention, i.e. declination angle in degrees (So that 90 degrees points to the bottom. The power is in dB. + '' + Input- + fname - string + basic file name to write. + ends with .sbp, if not it's appended + pattern - np 2d array of floats + first column is angles IN DEGREES + second column is power IN DB + """ + if fname[-4:] != '.sbp': + fname += '.sbp' + num_angles = pattern.shape[0] + np.savetxt(fname, pattern, fmt='%.'+str(num_dec_pts) + 'f', delimiter=' ', header=str(num_angles) + ' # Num pts',comments='') + return + + +def write_env( envfil, model, TitleEnv, freq, ssp, bdry, pos, beam, cint, RMax, NMESH=None,*varargin ): + ''' + write_env - writes an environment file given the environment parameters + + Parameters + ---------- + envfil : str + name of the environment file to save, should not include extensions + model : str + type of model to run. example ['kraken'] + TitleEnv : str + title of environment file (printed at top of .env file) + freq : float + frequency of environment + ssp : pyat.pyat.env.SSP + sound speed profile object + bdry : pyat.pyat.env.Bndry + boundary condition (top and bottom) object + pos : pyat.pyat.env.Pos + Position object? contains information about source and reciever depths and ranges + beam : (unknown) + unknown, but passing [] works + cint : an empty class? + unkown, but passing Empty where Empty is defined by (below) works + class Empty: + def __init__(self): + return + RMax : int + maximum range + NMESH : int + number of mesh points used for sound speed profile interpolation. Default is None, which + which uses the number of defined points in the sound speed profile object + + ''' + if (envfil[-4:] != '.env' ): + envfil = envfil + '.env' # append extension + + if ( len( varargin ) == 0 ): + f = open( envfil, 'w' ) # create new env file + else: + f = open( envfil, 'a' ) # append to existing envfil + + + model = model.upper() # convert to uppercase + + + f.write('\'' + TitleEnv + '\' ! Title \r\n') + f.write('{:8.2f}'.format(freq) +' \t \t \t ! Frequency (Hz) \r\n') + f.write('{:5d}'.format(ssp.NMedia)+ ' \t \t \t ! NMedia \r\n') + f.write('\'' + bdry.Top.Opt + '\''+ ' \t \t \t ! Top Option \r\n') + + if ( bdry.Top.Opt[0] == 'A' ): # analytic boundary + f.write(' {:6.2f}'.format(ssp.depth[0]) + \ + ' {:6.2f}'.format(bdry.Top.hs.alphaR) + \ + ' {:6.2f}'.format(bdry.Top.hs.betaR) + \ + ' {:6.2g}'.format(bdry.Top.hs.rho) + \ + ' {:6.2f}'.format(bdry.Top.hs.alphaI) + \ + ' {:6.2f}'.format(bdry.Top.hs.betaI) + \ + ' \t ! upper halfspace \r\n') + + # SSP + for medium in range(ssp.NMedia): + + # get NMESH + if NMESH is None: + NMESH = ssp.N[ medium ] + f.write('{:5d}'.format(NMESH) + \ + ' {:4.2f}'.format(ssp.sigma[ medium ]) + \ + ' {:6.2f}'.format(ssp.depth[ medium+1 ]) + ' \t ! N sigma depth \r\n') + for ii in range(len( ssp.raw[ medium].z )): + f.write('\t {:6.2f} '.format(ssp.raw[ medium ].z[ ii ]) + \ + '{:6.2f} '.format(ssp.raw[ medium ].alphaR[ ii ]) + \ + '{:6.2f} '.format(ssp.raw[ medium ].betaR[ ii ]) + \ + '{:6.2g}'.format(ssp.raw[ medium ].rho[ ii ] ) + \ + ' {:10.6f} '.format(ssp.raw[ medium ].alphaI[ ii ]) + \ + '{:6.2f} '.format(ssp.raw[ medium ].betaI[ ii ]) + \ + '/ \t ! z c cs rho \r\n') + + # lower halfspace + f.write('\''+bdry.Bot.Opt + '\'' + ' {:6.2f}'.format(ssp.sigma[1]) + ' \t \t ! Bottom Option, sigma\r\n') # ssp.sigma( 2 ) ) + + fmtr = {'all': lambda x: ' {:6.2f}'.format(float(x))}#, 'void': lambda x: ''} + a = np.array([]) + d = np.array([float(ssp.depth[ssp.NMedia])]) + vls = [d, bdry.Bot.hs.alphaR, bdry.Bot.hs.betaR, bdry.Bot.hs.rho, \ + bdry.Bot.hs.alphaI, bdry.Bot.hs.betaI] + strings = [np.array2string(x, formatter = fmtr) for x in vls] + s = ' {:6.2f} '.format(ssp.depth[ssp.NMedia]) + for string in strings[1:]: + if string == '': + string = ' 0.00' + s += string + if ( bdry.Bot.Opt[0] == 'A' ): +# print('{:6.2f}'.format(bdry.Bot.hs.betaI[0])) + f.write(' ' + s + ' \t / \t ! lower halfspace \r\n') + + if model in [ 'SCOOTER', 'KRAKEN', 'KRAKENC', 'SPARC' ]: + f.write('{:6.0f} '.format(cint.Low) + '{:6.0f} \t \t ! cLow cHigh (m/s) \r\n'.format(cint.High)) # phase speed limits + f.write('{:8.2f} \t \t \t ! RMax (km) \r\n'.format(RMax )) # maximum range + + # source depths + f.write('{:5d} \t \t \t \t ! NSD'.format(len( pos.s.depth ) )) + + if ( len( pos.s.depth ) >= 2) and equally_spaced( pos.s.depth ): + f.write('\r\n {:6f} '.format(pos.s.depth[0]) +' {:6f} '.format(pos.s.depth[-1])) + else: + f.write('\r\n ')#{:6f} '.format( pos.s.depth[0] )) + for tmp_depth in pos.s.depth: + f.write('{:6f} '.format(tmp_depth) ) + + f.write('/ \t ! SD(1) ... (m) \r\n' ) + + # receiver depths + + f.write('{:5d} \t \t \t \t ! NRD'.format(len( pos.r.depth ) )) + + if ( len( pos.r.depth ) >= 2) and equally_spaced( pos.r.depth ) : + f.write('\r\n {:6f} '.format(pos.r.depth[0]) + '{:6f} '.format( pos.r.depth[-1])) + else: + f.write('\r\n ') + for tmp_depth in pos.r.depth: + f.write('{:6f} '.format(tmp_depth) ) + + f.write('/ \t ! RD(1) ... (m) \r\n' ) + + # receiver ranges + if (model == 'BELLHOP' ) or (model == 'FirePE' ) : + # receiver ranges + f.write('{:5d} \t \t \t \t ! NRR'.format(len( pos.r.range ) )) + + if ( len( pos.r.range ) >= 2) and equally_spaced( pos.r.range ) : + f.write('\r\n {:6f} '.format(pos.r.range[0]) + '{:6f} '.format(pos.r.range[-1])) + else: + f.write('\r\n {:6f} '.format(pos.r.range[0] )) + f.write('/ \t ! RR(1) ... (km) \r\n' ) + write_bell(f, beam ) + + f.close() + +def write_ssp(sspfile, cw, r_arr): + """ + Input + sspfile - string + path of file to write + cw - numpy 2d array + ssp profiles, columns are the depth profiles for given range + r_arr - numpy 1d array + ranges (in km) for updating ssp + """ + if sspfile[-3:] != 'ssp': + sspfile += '.ssp' + with open(sspfile, 'w') as f: + f.write(str(r_arr.size)+'\r\n') + for val in r_arr: + f.write(str(val) + '\t') + f.write('\r\n') + for row in range(cw.shape[0]): + for val in cw[row,:]: + f.write(str(val) + '\t') + f.write('\r\n') + return + +def write_bathy(btyfile, bathm): + """ + Input + btyfile - string + path of file to write + range_depth_array - numpy 2d array + first column is ranges (in km), second column is depths + """ + if btyfile[-3:] != 'bty': + btyfile += '.bty' + with open(btyfile, 'w') as f: + f.write('CS\r\n') + f.write(str(len(bathm.d)) + '\r\n') + for row in range(len(bathm.r)): + f.write(' ' + str(bathm.r[row]) + ' ' + str(bathm.d[row]) + '\n') + return + +def write_bell(f, beam): + f.write('\''+beam.RunType+'\''+ ' \t \t \t \t ! Run Type \r\n'); + + if (beam.Ibeam != None): + f.write('{: 1) and ((beam.RunType[1] != 'G') and (beam.RunType[1] != 'B') and (beam.RunType[1] != 'S')) : + f.write('\''+beam.Type[0:1]+'\''+' {:f}'.format(beam.epmult)+' {:f}'.format(beam.rLoop )+' \t \t ! ''Min/Fill/Cer, Sin/Doub/Zero'' Epsmult RLoop (km) \r\n') + f.write('{:d}'.format(beam.Nimage)+' {:d}'.format(beam.Ibwin)+' \t \t \t \t ! Nimage Ibwin \r\n') + +def read_shd_bin(*varargin): + s = Source(0) + r = Dom(0,0) + pos = Pos(s,r) + ''' + Read TL surfaces from a binary Bellhop/Kraken .SHD file + without having to convert to ASCII first. + Useage: + ... = read_shd_bin( filename, xs, ys ) + where (xs, ys) is the source coordinate in km + (xs, ys) are optional + Output is a 4-D pressure field p( Ntheta, Nsd, Nrd, Nrr ) + ''' + if (len(varargin) < 1) or (len(varargin) > 3): + raise ValueError("Can only pass one to three arguments: filename; (filename, xs, ys); or (filename, freq)") + + filename = varargin[0] + + # optional frequency + if len(varargin) == 2: + freq = varargin[1] + else: + freq = np.NaN + + # optional source (x,y) coordinate + if len(varargin) >= 3: + xs = varargin[1] + ys = varargin[2] + else: + xs = np.NaN + ys = np.NaN + + ## + f = open( filename, 'rb' ) + + recl = unpack(' 3): + raise ValueError("Can only pass one to three arguments: filename; (filename, xs, ys); or (filename, freq)") + + filename = varargin[0] + + # optional frequency + if len(varargin) == 2: + freq = varargin[1] + else: + freq = np.NaN + + # optional source (x,y) coordinate + if len(varargin) >= 3: + xs = varargin[1] + ys = varargin[2] + else: + xs = np.NaN + ys = np.NaN + + PlotType = []; # in case this was not set + + [holder , name, ext ] = fileparts( filename ) + if ( ext == '.mat' ) : + [ holder, holder1, ext2 ] = fileparts( name ) + + if ext2 == '.shd': + FileType = 'shdmat' + elif ext2 == '.grn': + FileType = 'grnmat' + else: + pass + else: + if filename == 'ASCFIL': + FileType = 'asc' + elif filename == 'SHDFIL': + FileType = 'shd' + elif filename == 'tl.grid': + FileType = 'RAM' + else: + endchar = len( filename ) + if ( endchar >= 4 ): + FileType = filename[ endchar - 3 : endchar ].lower() + + ## + PlotTitle, PlotType, freqVec, atten, pos, pressure = [], [], [], [], [], [] + if FileType in ['shd', 'grn' ]: # binary format + if len(varargin) == 1: + [ PlotTitle, PlotType, freqVec, atten, pos, pressure ] = read_shd_bin( filename ) + if len(varargin) == 2: + [ PlotTitle, PlotType, freqVec, atten, pos, pressure ] = read_shd_bin( filename, freq ) + if len(varargin) == 3: + [ PlotTitle, PlotType, freqVec, atten, pos, pressure ] = read_shd_bin( filename, xs, ys ) + elif FileType == 'shdmat': # Shade function mat file + loadmat( filename ) + + # has a specific source xs, ys been given? + if not np.isnan( xs ): + xdiff = abs( pos.s.x - xs * 1000. ) + [ holder, idxX ] = min( xdiff ) + ydiff = abs( pos.s.y - ys * 1000. ) + [ holder, idxY ] = min( ydiff ) + + # extract the appropriate source index + pressureT = pressure[ idxX, idxY, :, :, :, : ] + pressure = np.reshape( pressureT, [ pos.Nsd, pos.Nrd, pos.Ntheta, pos.Nrr ] ) + + # has a specific frequency been given? + if not np.isnan( freq ): + freqdiff = abs( freqVec - freq ) + [ holder, ifreq ] = min( freqdiff ) + # extract the appropriate source index + pressure = pressure[ ifreq, 1, :, : ] + + elif FileType == 'asc': # ascii format + [ PlotTitle, PlotType, freqVec, atten, pos, pressure ] = read_shd_asc( filename ) + + elif FileType == 'grnmat': # Green's function mat file + loadmat( filename ) + pos.r.range = np.array(pos.r.range.T); # make it a column vector to match read_shd_bin + + elif FileType == 'RAM': + [ PlotTitle, PlotType, freqVec, atten, pos, pressure ] = read_ram_tlgrid( filename ) + else: + raise ValueError( 'Unrecognized file extension' ) + + # clean up PlotTitle by taking only the part up inside the quotes + # nchars = strfind( PlotTitle, '''' ); # find quotes + # PlotTitle = [ PlotTitle( nchars( 1 ) + 1 : nchars( 2 ) - 1 ) ] + return [ PlotTitle, PlotType, freqVec, atten, pos, pressure ] + +def crci( c, alpha, freq, AttenUnit ): + ''' + Convert real wave speed and attenuation to a single complex wave speed + 5 cases: + N for Nepers/meter + M for dB/meter (M for Meters) + F for dB/m-kHZ (F for frequency dependent) + W for dB/wavelength (W for Wavelength) + Q for Q + T for Thorpe + ''' + omega = 2.0 * np.pi * freq + + # *** Convert to Nepers/m *** + + alphaT = 0.0 + + if AttenUnit[0] == 'N' : # Nepers/m + alphaT = alpha + elif AttenUnit[0] == 'M' : # dB/meter + alphaT = alpha / 8.6858896 + elif AttenUnit[0] == 'F' : # dB/m-kHZ + alphaT = alpha * freq / 8685.8896 + elif AttenUnit[0] == 'W' : # dB/wavelength + if ( c != 0.0 ): + alphaT = alpha * freq / ( 8.6858896 * c ) + elif AttenUnit[0] == 'Q' : + if( c * alpha != 0.0 ): + alphaT = omega / ( 2.0 * c * alpha ) + else: + print('Unknown attenuation unit') + + # added volume attenuation + if AttenUnit[1] == ( 'T' ) : # Thorp + f2 = np.square(( freq / 1000.0 )) + # Original Thorp (1967) formula + # alphaT = 40.0 * f2 / ( 4100.0 + f2 ) + 0.1 * f2 / ( 1.0 + f2 ); + # alphaT = alphaT / 914.4; # dB / m + # alphaT = alphaT / 8.6858896; # Nepers / m + + # Updated formula from JKPS Eq. 1.34 + Thorp = 3.3e-3 + 0.11 * f2 / ( 1.0 + f2 ) + 44.0 * f2 / ( 4100.0 + f2 ) + 3e-4 * f2 # dB/km + Thorp = Thorp / 8685.8896 # Nepers / m + alphaT = alphaT + Thorp + # *** Convert Nepers/m to equivalent imaginary sound speed *** + alphaT = alphaT * c*c / omega + crci = complex(c, alphaT) + return crci + +def topbot( lines, line_ind, freq, BCType, AttenUnit ): + ''' Handles top and bottom boundary conditions + Input: + ENVFIL: Environmental file + freq: frequency # BCType: Boundary condition type + + Output: + Bdry.cp: P-wave speed in halfspace + Bdry.cs: S-wave speed in halfspace + Bdry.rho: density in halfspace + + BumDen: Bump density + eta: Principal radius 1 + xi: Principal radius 2 + ''' + global alphaR, betaR, rhoR, alphaI, betaI + + # *** Echo to PRTFIL user's choice of boundary condition *** + bctypes = ['S', 'H', 'T', 'I', 'V', 'R', 'A', 'F', 'W', 'P'] + if BCType not in bctypes: + raise ValueError( 'Fatal error: Unknown boundary condition type' ) + + # ****** Read in BC parameters depending on particular choice ****** + + cp = 0.0 + cs = 0.0 + rho = 0.0 + ice = Ice(0.0, 0.0, 0.0) + hs = None + + # *** Twersky ice model parameters *** + + if ( BCType == 'S' or BCType == 'H' or BCType == 'T' or BCType == 'I' ): + tmp = lines[line_ind] + ice.BumDen, ice.eta, ice.xi = float(tmp.split()[0]), float(tmp.split()[1]), float(tmp.split()[2]) + line_ind += 1 + + # *** Half-space properties *** + + if ( BCType == 'A' ): + tmp = lines[line_ind] + line_ind += 1 + tmp = tmp.split() + tmp = [x for x in tmp if x[0].isnumeric()] # filter out non numbers + num_vals = len(tmp) + if num_vals == 6: + ztmp, alphaR, betaR, rhoR, alphaI, betaI = [float(x) for x in tmp[0:6]] + elif num_vals == 5: + ztmp, alphaR, betaR, rhoR, alphaI = [float(x) for x in tmp[0:5]] + elif num_vals == 4: + ztmp, alphaR, betaR, rhoR = [float(x) for x in tmp[0:4]] + elif num_vals == 3: + ztmp, alphaR, betaR = [float(x) for x in tmp[0:3]] + elif num_vals == 2: + ztmp, alphaR = [float(x) for x in tmp[0:2]] + elif num_vals == 1: + ztmp = [float(x) for x in tmp][0] + else: # there were no vals to read in so defaults will be used + pass + + cp = crci( alphaR, alphaI, freq, AttenUnit ) + cs = crci( betaR, betaI, freq, AttenUnit ) + rho = rhoR + + # raw values in user units + hs = HS(alphaR, betaR, rho, alphaI, betaI) + return [ cp, cs, rho, ice, hs, line_ind] + +def read_env_core( envfil ): + global HV, NFirstAcoustic, NLastAcoustic + global alphaR, betaR, rhoR, alphaI, betaI + + alphaR = 1500 # defaults + betaR = 0 + rhoR = 1 + alphaI = 0 + betaI = 0 + + bdry = Bndry(None,None) + bdry.Top = TopBndry('dummy_opt') + bdry.Bot = BotBndry('dummy_opt', 10) + + NFirstAcoustic = 0 + bdry.Top.cp = 0.0 + bdry.Top.cs = 0.0 + bdry.Top.rho = 0.0 + bdry.Bot.cp = 2000.0 + bdry.Bot.cs = 0.0 + bdry.Bot.rho = 2.0 + + ssp = None + + with open(envfil, 'r') as f: + lines = f.readlines() + # first line is title + line_ind = 0 + TitleEnv = re.findall('\'(.*)\'',lines[line_ind] )[0] + # second line is freq + line_ind += 1 + freq = float((lines[line_ind].split())[0]) + line_ind += 1 + + # third line is number of media + tmp = lines[line_ind] + line_ind += 1 + NMedia = int(tmp.split()[0]) + # initialize depth array + ssp_depth = [0]*(NMedia+1) + + # fourth line is Top option + tmp = lines[line_ind] + line_ind += 1 + TopOpt = re.findall('\'(.*)\'', tmp)[0] + if len(TopOpt) < 5: + TopOpt = TopOpt + ' '*(4-len(TopOpt)) # whitespace pad to 4 characters + # convert the deprecated '*' option to '~' + bdry.Top.Opt = TopOpt.replace('*', '~') + + SSPType = bdry.Top.Opt[0] + bdry.Top.BC = bdry.Top.Opt[1] + # this should throw error, I may need to pad the end with spaces + AttenUnit = bdry.Top.Opt[2:4] + + # *** SSP approximation options *** + ssp_types = ['N', 'C', 'S', 'Q', 'H', 'A'] + if SSPType not in ssp_types: + raise ValueError('SSP must be in ' + str(ssp_types)) + + + # *** Attenuation options *** + atten_opts = ['N', 'F', 'M', 'W', 'Q'] + if AttenUnit[0] not in atten_opts: + raise ValueError( 'Unknown attenuation units, must be ' + str(atten_opts) ) + + [ bdry.Top.cp, bdry.Top.cs, bdry.Top.rho, holder, bdry.Top.hs, line_ind ] = topbot(lines, line_ind, freq, bdry.Top.BC, AttenUnit ) + + ssp_z = [] + ssp_c = [] + ssp_cs = [] + ssp_rho = [] + ssp_raw_list = [] + Loc = [0] * NMedia + ssp_Npts = [0]*NMedia + ssp_N = [0]*NMedia + ssp_sigma = [0]*(NMedia+1) + for medium in range(NMedia): + if ( medium == 1 ): + Loc[medium] = 0 + else: + Loc[ medium ] = Loc[ medium - 1 ] + ssp_Npts[ medium - 1 ] + + # next line is the num points in ssp, sigma, depth + tmp = lines[line_ind] + line_ind += 1 + ssp_N[medium], ssp_sigma[medium], ssp_depth[medium+1] = int(float(tmp.split()[0])), float(tmp.split()[1]), float(tmp.split()[2]) + + # read in the SSP + ssp_raw_list.append(SSPraw([], [], [], [], [], [])) + + for ii in range(1000): + tmp = lines[line_ind] + line_ind += 1 + tmp = tmp.split() + tmp = [x for x in tmp if x[0].isnumeric()] # filter out non numbers + num_vals = len(tmp) + if num_vals == 6: + ztmp, alphaR, betaR, rhoR, alphaI, betaI = [float(x) for x in tmp[0:6]] + elif num_vals == 5: + ztmp, alphaR, betaR, rhoR, alphaI = [float(x) for x in tmp[0:5]] + elif num_vals == 4: + ztmp, alphaR, betaR, rhoR = [float(x) for x in tmp[0:4]] + elif num_vals == 3: + ztmp, alphaR, betaR = [float(x) for x in tmp[0:3]] + elif num_vals == 2: + ztmp, alphaR = [float(x) for x in tmp[0:2]] + elif num_vals == 1: + ztmp = [float(x) for x in tmp[0]][0] + else: # there were no vals to read in so defaults will be used + pass + + + ## TRANSLATE CRCI, IN BELLHOP DIR + cp = crci( alphaR, alphaI, freq, AttenUnit ) + cs = crci( betaR, betaI, freq, AttenUnit ) + ssp_z.append(ztmp) # add in to existing vector + ssp_c.append(cp) + ssp_cs.append(cs) + ssp_rho.append(rhoR) + + ssp_raw_list[ medium].z.append(ztmp) # add in to existing vector + ssp_raw_list[ medium].alphaR.append(alphaR) + ssp_raw_list[ medium].alphaI.append(alphaI) + ssp_raw_list[ medium].betaR.append(betaR) + ssp_raw_list[ medium].betaI.append(betaI) + ssp_raw_list[ medium].rho.append(rhoR) + + # check for end of this layer + if ( ztmp == ssp_depth[ medium+1 ] ): + ssp_depth[0] = ssp_z[0] + # break + break + + # end of loop for that medium + if ( ssp_N[ medium ] == 0 ): # calculate mesh automatically + # choose a reference sound speed + C = alphaR + if ( betaR > 0.0 ): + C = betaR # shear? + deltaz = 0.05 * C / freq # default sampling: 20 points per wavelength + ssp_N[ medium ] = round( ( ssp_depth[ medium + 1 ] - ssp_depth[ medium ] ) / deltaz ) + ssp_N[ medium ] = np.max( [ssp_N[ medium ], 10] ) # require a minimum of 10 points + + # keep track of first and last acoustic medium + if ( not np.any( cs ) ): # shear anywhere? + if ( NFirstAcoustic == 0 ): + NFirstAcoustic = medium + NLastAcoustic = medium + + # stuff for Bellhop + if ( medium == 0 ): + HV = np.diff( np.array(ssp_z) ) # layer thicknesses + ssp_cz = np.diff( np.array(ssp_c )) / HV # gradient of ssp (centered-difference approximation) + ssp_Npts[ medium ] = ii + if ( medium == 0 ): + ssp_depth[0] = ssp_z[0] + + ## lower halfspace + + # read in next line + tmp = lines[line_ind] + line_ind += 1 + BotOpt = re.findall('\'(.*)\'', tmp)[0] + + if len(BotOpt) < 5: + BotOpt = BotOpt + ' '*(4-len(BotOpt)) # whitespace pad to 4 characters + # convert the deprecated '*' option to '~' + bdry.Bot.Opt = BotOpt.replace('*', '~') + bdry.Bot.BC = bdry.Bot.Opt[0] + + ## TRANSLATE TOPBOT IN BELLHOP DIR + [ bdry.Bot.cp, bdry.Bot.cs, bdry.Bot.rho, holder, bdry.Bot.hs, line_ind ] = topbot( lines, line_ind, freq, bdry.Bot.BC, AttenUnit ) + + + bdry.Top.depth = ssp_depth[0] + bdry.Bot.depth = ssp_depth[ NMedia ] + + # Get rho, c just INSide the boundary (used later for reflection + # coefficients) + I = NFirstAcoustic + bdry.Top.rhoIns = ssp_rho[ I ] + bdry.Top.cIns = ssp_c[ I ] + + I = Loc[ NLastAcoustic ] + ssp_Npts[ NLastAcoustic ] + bdry.Bot.rhoIns = ssp_rho[ I ] + bdry.Bot.cIns = ssp_c[ I ] + ssp = SSP(ssp_raw_list, ssp_depth, NMedia, N=ssp_N, sigma=ssp_sigma) + return [TitleEnv, freq, ssp, bdry, lines, line_ind] + +def readvector(lines, line_ind): + ''' + Reads a row-vector + This routine emulates a Fortran language capability that allows a '/' + to be used to terminate input + A user can type: + 5 ! Number of values + 0 1000 / + as a shorthand for + 5 + 0 250 500 750 1000 + Here the '/' terminates the reading and the code then creates 4 equally + spaced values between 0 and 1000. + However, one also has the option of using specific (non-equispaced + values) + ''' + line_ind = int(line_ind) + tmp = lines[line_ind] + line_ind += 1 + Nx = int(tmp.split()[0]) + tmp = lines[line_ind] + line_ind += 1 + + if '/' in tmp: + tmp = tmp.split() + tmp = [x for x in tmp if x[0].isnumeric()] + if len(tmp) == 1: + x = float(tmp[0]) + else: + x = [float(tmp[0]), float(tmp[1])] + if Nx > 2: + x = np.linspace( x[ 0 ], x[ 1 ], Nx ) # generate the vector + else: + x = [float(x) for x in tmp if x[0].isnumeric()] + return x, Nx, line_ind + +def readsdrd(lines, line_ind): + ''' + Read source depths and receiver depths + Variable 'Pos' is a structure: + Pos.r.depth = vector of receiver depths + Pos.Nrd = number of receiver depths + Pos.s.depth = vector of source depths + Pos.Nsd = number of source depths + ''' + [ depth, Nx, line_ind] = readvector(lines, line_ind); + source = Source(depth) + + [ depth, Nrd, line_ind ] = readvector(lines, line_ind); + dom = Dom(None, depth) + pos = Pos(source, dom) + return pos, line_ind + +def read_env( envfil, model ): + ''' Read an environmental file + + should check here that model is one of the known ones + ''' + if ( envfil != 'ENVFIL' ) and (envfil[-4:] != '.env' ): + envfil = envfil + '.env' # append extension + + model = model.upper() # convert to uppercase + + [ TitleEnv, freq, ssp, bdry, lines, line_ind] = read_env_core( envfil ) # read in the environmental file + + model_cat = ['SCOOTER', 'KRAKEN', 'KRAKENC', 'KRAKEL', 'SPARC' ] + + if model in model_cat: + tmp = lines[line_ind] + line_ind += 1 + tmp = [x for x in tmp.split() if x[0].isnumeric()] + low, high = [float(x) for x in tmp] + cint = cInt(low, high) + tmp = lines[line_ind] + line_ind += 1 + RMax = float(tmp.split()[0]) + + + else: # dummy values for BELLHOP + cint = cInt(1500, 1e9) + RMax = 100 + + # BELLHOP3D has x-y coordinates of sources as well + + + if ( model == 'BELLHOP3D' ) : + raise ValueError("Bellhop 3d not supported atm") + + # !!! check: does this delete Pos.s.x, etc. when executed + [pos, line_ind] = readsdrd(lines, line_ind) # read in the source and receiver depths + + if ( model == 'BELLHOP' ) : + raise ValueError("Bellhop not supported at the moment") +# Pos.r.range = readr( fid ) # read in receiver ranges +# Pos.Nrr = length( Pos.r.range ) +# Beam = read_bell( fid, Bdry, freq, Bdry.Bot.depth, Bdry.Top.depth, Pos.r.range( end ) ) +# elseif ( strcmp( model, 'BELLHOP3D' ) ) +# Pos.r.range = readr( fid ) # read in receiver ranges +# Pos.Nrr = length( Pos.r.range ) +# Pos.theta = readtheta( fid ) +# Pos.Ntheta = length( Pos.theta ) +# Beam = read_bell( fid, Bdry, freq, Bdry.Bot.depth, Bdry.Top.depth, Pos.r.range( end ) ) + else: # dummy value for models that don't use Beam parameters + beam = Beam('CG', Nbeams=0, alpha= [-15,15], box=Box(1.05*RMax, 1.05*np.max(pos.r.depth)), deltas = 0) + pos.r.range = np.linspace( 0, RMax, 501 ) # set up receiver range vector + return [ TitleEnv, freq, ssp, bdry, pos, beam, cint, RMax] + +def read_modes(**kwargs): + ''' + Read the modes produced by KRAKEN + usage: + keys are 'fname', 'freq', 'modes' + 'fname' and 'freq' are mandatory, 'modes' is if you only want a subset of modes + [ Modes ] = read_modes_bin( filename, modes ) + filename is without the extension, which is assumed to be '.moA' + freq is the frequency (involved for broadband runs) + (you can use freq=0 if there is only one frequency) + modes is an optional vector of mode indices + + derived from readKRAKEN.m Feb 12, 1996 Aaron Thode + + Translated to python by Hunter Akins 2019 + + Modes.M number of modes + Modes.k wavenumbers + Modes.z sample depths for modes + Modes.phi modes + + Modes.Top.bc + Modes.Top.cp + Modes.Top.cs + Modes.Top.rho + Modes.Top.depth + + Modes.Bot.bc + Modes.Bot.cp + Modes.Bot.cs + Modes.Bot.rho + Modes.Bot.depth + + Modes.N Number of depth points in each medium + Modes.Mater Material type of each medium (acoustic or elastic) + Modes.Nfreq Number of frequencies + Modes.Nmedia Number of media + Modes.depth depths of interfaces + Modes.rho densities in each medium + Modes.freqVec vector of frequencies for which the modes were calculated + ''' + + + filename = kwargs['fname'] + if 'freq' in kwargs.keys(): + freq = kwargs['freq'] + if 'modes' in kwargs.keys(): + modes = kwargs['modes'] + with open(filename, 'rb') as f: + iRecProfile = 1; # (first time only) + + lrecl = 4*unpack('= len(lines)-1: + break + else: + angle = float(lines[line_ind]) + line_ind += 1 + tmp = lines[line_ind].split(' ') + tmp = [x for x in tmp if x != ''] + nsteps = int(tmp[0]) + num_top_bnc= int(tmp[1]) + num_bot_bnc= int(tmp[2]) + line_ind += 1 + x,y = np.zeros(nsteps), np.zeros(nsteps) + counter = 0 + for line_ind in range(line_ind, line_ind + nsteps): + tmp = lines[line_ind].split(' ') + tmp = [x for x in tmp if x != ''] + x[counter],y[counter] = float(tmp[0]), float(tmp[1]) + counter +=1 + + if num_top_bnc == 0 and num_bot_bnc==0: + axis.plot(x,y, color='k') + elif num_top_bnc == 1 and num_bot_bnc == 0: + axis.plot(x,y, color='r', alpha=.85) + elif num_top_bnc == 0 and num_bot_bnc == 1: + axis.plot(x,y, color='b', alpha=.85) + else: + axis.plot(x,y, color='g',alpha=.7) + line_ind += 1 + + return fig, axis + + +def get_rays(fname): + with open(fname, 'r') as f: + lines = f.readlines() + TITLE = lines[0] + freq = float(lines[1]) + tmp = lines[2].split(' ') + tmp = [x for x in tmp if x != ''] + Nsx = int(tmp[0]) + Nsy = int(tmp[1]) + Nsz = int(tmp[2]) + tmp = lines[3].split(' ') + tmp = [x for x in tmp if x != ''] + Nalpha = int(tmp[0]) + Nbeta = int(tmp[1]) + deptht = float(lines[4]) + depthb = float(lines[5]) + + + + line_ind = 7 + ray_collections = [] + for i in range(Nsz): + source_collection=[] + for ibeam in range(Nalpha): + if line_ind >= len(lines)-1: + break + else: + angle = float(lines[line_ind]) + line_ind += 1 + tmp = lines[line_ind].split(' ') + tmp = [x for x in tmp if x != ''] + nsteps = int(tmp[0]) + num_top_bnc= int(tmp[1]) + num_bot_bnc= int(tmp[2]) + line_ind += 1 + xy = np.zeros((2,nsteps)) + counter = 0 + for line_ind in range(line_ind, line_ind + nsteps): + tmp = lines[line_ind].split(' ') + tmp = [x for x in tmp if x != ''] + xy[0,counter],xy[1,counter] = float(tmp[0]), float(tmp[1]) + counter +=1 + ray = Eigenray(angle, num_top_bnc, num_bot_bnc, xy) + source_collection.append(ray) + line_ind += 1 + ray_collections.append(source_collection) + + return ray_collections + + +def read_shd_asc(filename): + raise NotImplementedError + +def read_ram_tlgrid(filename): + raise NotImplementedError diff --git a/pyat/rlcore.py b/pyat/rlcore.py new file mode 100755 index 0000000..e3b6c1e --- /dev/null +++ b/pyat/rlcore.py @@ -0,0 +1,74 @@ +import numpy as np + +"""混响级计算""" +def call_RLs_Ss(v,SL,freq,tau,theta,l,r): #海面混响 + """ + v:风速，m/s + SL：声源级,dB + freq:声纳频率，KHz + tau: 发射脉冲宽度,ms + theta：海底掠射角,du + l: 阵列长度,m + r： 斜向距离,m""" + freq = freq*1000 + beta=158*np.power((v*np.power(freq,1.0/3)),-0.58) + Ss=3.3*beta*np.log10(theta/30)-42.4*np.log10(beta)+2.6#海面散射强度经验公式，Chapman,Harris.水声学原理，p237 + C=1500 #声速 + tau=tau/1000 #单位换算 + lam=C/freq#波长 + phi= np.power(10,np.log10(lam/2/l/np.pi)+0.92)#l>lambda水平线阵/边长>>lambda等效波束宽度，水声学原理p225 + TL=20*np.log10(r) #传播损失经验公式 + A = C*tau*phi*r/2 #海面散射面积 + rls = SL-2*TL+Ss+10*np.log10(A) + return rls +def call_RLs_Sb(SL,freq,tau,theta,l,r): #海底混响 + """ + SL：声源级,dB + freq:声纳频率，KHz + tau: 发射脉冲宽度,ms + theta：海底掠射角,du + l: 阵列长度,m + r： 斜向距离,m""" + Sb=-28+10*np.log10(np.sin(theta)*np.sin(theta))#海底散射强度经验公式，Mackenzie.水声学原理，p243 + C=1500 #声速 + tau=tau/1000 #单位换算 + lam=C/(freq*1000)#波长 + phi= np.power(10,np.log10(lam/2/l/np.pi)+0.92)#l>lambda水平线阵/边长>>lambda等效波束宽度，水声学原理p225 + TL=20*np.log10(r) #传播损失经验公式 + A = C*tau*phi*r/2 #海底散射面积 + rls = SL-2*TL+Sb+10*np.log10(A) + return rls + +def call_RLs_Sv(SL,Sv,freq,tau,l,r): #海面混响 + """ + SL：声源级,dB + SV:散射强度,用户输入 + freq:声纳频率，KHz + tau: 发射脉冲宽度,ms + theta：海底掠射角,du + l: 阵列长度,m + r： 斜向距离,m""" + freq = freq*1000 + C=1500 #声速 + tau=tau/1000 #单位换算 + lam=C/freq#波长 + phi= np.power(10,np.log10(lam/2/l/np.pi)+0.92)#l>lambda水平线阵/边长>>lambda等效波束宽度，水声学原理p225 + TL=20*np.log10(r) #传播损失经验公式 + A = C*tau*phi*r*r/2 #海面散射面积 + rls = SL-2*TL+Sv+10*np.log10(A) + return rls + +def call_rl(v,SL,Sv,freq,tau,theta,l,r): + """ + v:风速，m/s + SL：声源级,dB + SV:散射强度,用户输入 + freq:声纳频率，KHz + tau: 发射脉冲宽度,ms + theta：海底掠射角,du + l: 阵列长度,m + r： 斜向距离,m""" + rlss = call_RLs_Ss(v,SL,freq,tau,theta,l,r) + rlsb = call_RLs_Sb(SL,freq,tau,theta,l,r) + rlsv = call_RLs_Sv(SL,Sv,freq,tau,l,r) + return [rlss, rlsb, rlsv] diff --git a/pyat/setup.py b/pyat/setup.py new file mode 100644 index 0000000..c02d355 --- /dev/null +++ b/pyat/setup.py @@ -0,0 +1,7 @@ +from distutils.core import setup +from Cython.Build import cythonize +setup(ext_modules = cythonize(["bellhop.py", "config.py", "cscore.py", "czcore.py", "env.py", "krakenc.py", + "OANPd.py", "ProfileDataType.py", "ProfileRead.py", "ProfileTool.py", "ram.py", "ramRW_new.py", + "ramRW.py", "RANDIcore.py", "readwrite.py", "rlcore.py", "tlcore.py", "WGNPd.py" + ])) +# setup(ext_modules = cythonize(["cscore.py"])) \ No newline at end of file diff --git a/pyat/tlcore.py b/pyat/tlcore.py new file mode 100644 index 0000000..cf4e1ff --- /dev/null +++ b/pyat/tlcore.py @@ -0,0 +1,87 @@ +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 +import matplotlib.pyplot as plt + +def call_tl(ETOPO1Path, STPath, freq, sourceDepth, sourceLat, sourceLon, receiveDepth, + receiveLat=None, receiveLon=None, range=None, azi=None, timeType='Annual', modelType='kraken', isPlot=False): + """ Input: + freq - float, Frequency + sourceDepth - float, Source depth + sourceLat - float, Source lat + sourceLon - float, Source lon + receiveDepth - float, Receive depth + receiveLat - float, Receive lat + receiveLon - float, Receive lon + range - float, Range from source to receive + azi - float, azimuth angle from receive to source + timeType - char, #全年平均 'Annual', #1-12月'January''February''March''April''May''June''July''August''September' + October': '10''November''December', #北半球四季'Winter''Spring''Summer''Autumn' + modelType - char, 'bellhop'、'kraken'、'ram' + isPlot - plot or not + """ + """ Output: + ranges 距离 + depths 深度 + TL 传播损失 + bathym类 地形数据 bathym.r距离 bathym.d深度 + """ + # 定义时间段 + timeIndex = timeType + # 创建类 + data = DataBase(ETOPO1Path, STPath) + # 设置时间段 + data.set_timePeriod(timeIndex) + + bathym, ssp, sed, base = [], [], [], [] + if (receiveLat != None) and (receiveLon != None): + [bathym, ssp, sed, base] = GetEnv(data, 25, Coord(sourceLon, sourceLat), coordR=Coord(receiveLon, receiveLat)) + elif (range != None) and (azi != None): + [bathym, ssp, sed, base] = GetEnv(data, 25, Coord(sourceLon, sourceLat), R=range, azi=azi) + else: + print('no receive info!') + 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) + # Plot + if isPlot: + plt.pcolor(Pos.r.range, Pos.r.depth, TL[0,0,:,:], vmin=40, vmax=130, cmap='jet_r', shading='auto') + plt.xlabel("Range/m") + plt.ylabel("Depth/m") + plt.gca().invert_yaxis() + plt.colorbar() + plt.show() + 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) + # Plot + if isPlot: + plt.pcolor(Pos.r.range, Pos.r.depth, TL[0,0,:,:], vmin=40, vmax=130, cmap='jet_r', shading='auto') + plt.xlabel("Range/m") + plt.ylabel("Depth/m") + plt.gca().invert_yaxis() + plt.colorbar() + plt.show() + 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) + # Plot + if isPlot: + plt.pcolormesh(rr, rd, TL, vmin=60, vmax=130, cmap='jet_r', shading='auto') + plt.xlabel("Range/m") + plt.ylabel("Depth/m") + plt.gca().invert_yaxis() + plt.colorbar() + plt.show() + return rr, rd, TL, bathym + else: + return [], [], [], bathym + \ No newline at end of file diff --git a/rl.py b/rl.py new file mode 100755 index 0000000..c8d8d78 --- /dev/null +++ b/rl.py @@ -0,0 +1,58 @@ +# -*- coding: utf-8 -*- +import warnings +warnings.filterwarnings('ignore') +import argparse, json + +def _argparse(): + parser = argparse.ArgumentParser(description='混响级计算') + 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: + paramData = json.load(jsonFile) + if ("windSpeed" not in paramData) or ("soundLevel" not in paramData) or\ + ("SV" not in paramData) or ("freq" not in paramData) or\ + ("tau" not in paramData) or ("theta" not in paramData) or\ + ("length" not in paramData) or ("range" not in paramData) or\ + ("outParam" not in paramData): + print(f'{argfile}中的参数有误，请核实！') + return + v = float(paramData["windSpeed"]) + SL = float(paramData["soundLevel"]) + Sv = float(paramData["SV"]) + freq = float(paramData["freq"]) + tau = float(paramData["tau"]) + theta = float(paramData["theta"]) + l = float(paramData["length"]) + r = float(paramData["range"]) + outParam = paramData["outParam"] + + from pyat.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) + + #输出json文件 + outData = {} + outData["param"] = paramData + outData["海面混响"] = f'{round(rlss,2)}' + outData["海底混响"] = f'{round(rlsb,2)}' + outData["体积混响"] = f'{round(rlsv,2)}' + with open(("" if outParam["outPath"]=="" else outParam["outPath"] + "/")+"rlOut.json", 'w') as jsonFile: + json.dump(outData, jsonFile, indent=2, ensure_ascii=False) + + #海面混响 + print("海面混响计算结果:", round(rlss,2)) + #海底混响 + print("海底混响计算结果:", round(rlsb,2)) + #立体混响 + print("体积混响计算结果:", round(rlsv,2)) + + +if __name__ == '__main__': + main() diff --git a/tl.py b/tl.py new file mode 100644 index 0000000..aedb543 --- /dev/null +++ b/tl.py @@ -0,0 +1,79 @@ +# -*- coding: utf-8 -*- +import warnings +warnings.filterwarnings('ignore') +import argparse, json +from common import NoIndent, MyEncoder + +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: + 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): + print(f'{argfile}中的参数有误，请核实！') + return + ETOPO1Path = paramData["ETOPO1Path"] + STPath = paramData["STPath"] + freq = float(paramData["freq"]) + sourceDepth = float(paramData["sourceDepth"]) + receiveDepth = float(paramData["receiveDepth"]) + sLat = float(paramData["sLat"]) + sLon = float(paramData["sLon"]) + rLat = float(paramData["rLat"]) + rLon = float(paramData["rLon"]) + timeType = paramData["timeType"] + modelType = paramData["modelType"] + png = paramData["png"] + outParam = paramData["outParam"] + + from pyat.tlcore import call_tl + import numpy as np + [ranges, depths, TL, bathym] = call_tl(ETOPO1Path, STPath, freq=freq, sourceDepth=sourceDepth, sourceLat=sLat, sourceLon=sLon, receiveDepth=receiveDepth, + receiveLat=rLat, receiveLon=rLon, timeType=timeType, modelType=modelType, isPlot=png) + + # write data to file + # 保存地形数据 + bathym.r = np.around(bathym.r, decimals=0) + bathym.d = np.around(bathym.d, decimals=0) + bathymData = { + "param": paramData, + "count": len(bathym.d), + "descript": "距离、深度", + "data": NoIndent(list(zip(bathym.r, bathym.d))) + } + 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) + + # 保存传播损失数据 + tlpoint = [] + 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())) + tlData = { + "param": paramData, + "count": len(ranges)*len(depths), + "descript": "距离、深度、损失", + "data": NoIndent(tlpoint) + } + 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() \ No newline at end of file diff --git a/tl_realtime.py b/tl_realtime.py new file mode 100644 index 0000000..8b6156b --- /dev/null +++ b/tl_realtime.py @@ -0,0 +1,207 @@ +# -*- coding: utf-8 -*- +import warnings +warnings.filterwarnings('ignore') +import argparse, json +from common import NoIndent, MyEncoder + +def _argparse(): + parser = argparse.ArgumentParser(description='tl realtime 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: + paramData = json.load(jsonFile) + if ("freq" not in paramData) or ("sourceDepth" not in paramData) or\ + ("receiveDepth" not in paramData) or ("modelType" not in paramData) or\ + ("bathyfile" not in paramData) or ("sspfile" not in paramData) or\ + ("z_min" not in paramData) or ("z_max" not in paramData) or\ + ("cp_min" not in paramData) or ("cp_max" not in paramData) or\ + ("cs_min" not in paramData) or ("cs_max" not in paramData) or\ + ("rho_min" not in paramData) or ("rho_max" not in paramData) or\ + ("ap_min" not in paramData) or ("ap_max" not in paramData) or\ + ("as_min" not in paramData) or ("as_max" not in paramData) or\ + ("Cb_p" not in paramData) or ("Cb_s" not in paramData) or\ + ("Rhob" not in paramData) or ("Ab_p" not in paramData) or\ + ("Ab_s" not in paramData) or ("png" not in paramData) or\ + ("outParam" not in paramData): + print(f'{argfile}中的参数有误，请核实！') + return + freq = float(paramData["freq"]) + sourceDepth = float(paramData["sourceDepth"]) + receiveDepth = float(paramData["receiveDepth"]) + modelType = paramData["modelType"] + btyfile = paramData["bathyfile"] + sspfile = paramData["sspfile"] + z_min = float(paramData["z_min"]) + z_max = float(paramData["z_max"]) + cp_min = float(paramData["cp_min"]) + cp_max = float(paramData["cp_max"]) + cs_min = float(paramData["cs_min"]) + cs_max = float(paramData["cs_max"]) + rho_min = float(paramData["rho_min"]) + rho_max = float(paramData["rho_max"]) + ap_min = float(paramData["ap_min"]) + ap_max = float(paramData["ap_max"]) + as_min = float(paramData["as_min"]) + as_max = float(paramData["as_max"]) + Cb_p = float(paramData["Cb_p"]) + Cb_s = float(paramData["Cb_s"]) + Rhob = float(paramData["Rhob"]) + Ab_p = float(paramData["Ab_p"]) + Ab_s = float(paramData["Ab_s"]) + png = paramData["png"] + outParam = paramData["outParam"] + + 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 + import matplotlib.pyplot as plt + + bathym, ssp, sed, base = [], [], [], [] + with open(btyfile, 'r') as f: + lines = f.readlines() + if len(lines) < 1: + print(f'{btyfile}中的参数有误，请核实！') + return + # 读取数据 + data = np.loadtxt(btyfile, skiprows=1, usecols=(0, 1)) + bathym = BathyM(data[:, 0], data[:, 1]) + + profile_num = 0 + with open(sspfile, 'r') as f: + try: + line = f.readline() + profile_num = int((line.split())[0]) + while line: + line = f.readline() + if '距离' in line: + ranges, depths, vels = [], [], [] + line = f.readline() + while line and ('分隔符' not in line): + ranges.append(1) + depths.append(float((line.split())[0])) + vels.append(float((line.split())[1])) + line = f.readline() + ssp.append(SSP_(ranges, depths, vels)) + except Exception as err: + print(f'{sspfile}中的参数有误，请核实！') + print(err) + return + + """akima 0829""" + z_ssp = np.array([]) + for i in range(len(ssp)): + z_ssp = np.append(z_ssp, ssp[i].z) + z_ssp = np.unique(z_ssp) + z_ssp = np.sort(z_ssp) + + for i in range(len(ssp)): + dep = ssp[i].z + vel = ssp[i].c + # 对剖面进行垂向插值处理 + depInterp = z_ssp[z_ssp<=ssp[i].z[-1]] + # 应用Akima插值算法，对ssp进行插值 + func = Akima1DInterpolator(dep, vel) + velInterp = func(depInterp) + # 剖面赋值 + ssp[i].z = depInterp + ssp[i].c = velInterp + len_real = len(depInterp)-1 + #for it in range(len_real,int(deep_max)): + # ssp[i].c = np.append(ssp[i].c,0) + + # Sediment + z = np.array([z_min, z_max]) + cp = np.array([cp_min, cp_max]) + cs = np.array([cs_min, cs_max]) + rho = np.array([rho_min, rho_max]) + a_p = np.array([ap_min, ap_max]) + a_s = np.array([as_min, as_max]) + for it in range(profile_num): + sed.append(Sed(z, cp, cs, rho, a_p, a_s)) + + # Base + base = [] + for it in range(profile_num): + base.append(Base(Cb_p, Cb_s, Rhob, Ab_p, Ab_s)) + + ranges, depths, TL = [], [], [] + 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) + ranges, depths, TL = Pos.r.range, Pos.r.depth, TL[0,0,:,:] + # Plot + if png: + plt.pcolor(Pos.r.range, Pos.r.depth, TL[0, 0, :, :], vmin=40, vmax=130, cmap='jet_r', shading='auto') + plt.xlabel("Range/m") + plt.ylabel("Depth/m") + plt.gca().invert_yaxis() + plt.colorbar() + plt.show() + print(TL[0, 0, :, :]) + 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) + ranges, depths, TL = Pos.r.range, Pos.r.depth, TL[0,0,:,:] + # Plot + if png: + plt.pcolor(Pos.r.range, Pos.r.depth, TL[0, 0, :, :], vmin=40, vmax=130, cmap='jet_r', shading='auto') + plt.xlabel("Range/m") + plt.ylabel("Depth/m") + plt.gca().invert_yaxis() + plt.colorbar() + plt.show() + elif modelType == 'ram': + # 调用ram模块 + [rd, rr, TL] = call_ram(freq, sourceDepth, receiveDepth, bathym, ssp, sed) + ranges, depths, TL = rr, rd, TL + # Plot + if png: + plt.pcolormesh(rr, rd, TL, vmin=60, vmax=130, cmap='jet_r', shading='auto') + plt.xlabel("Range/m") + plt.ylabel("Depth/m") + plt.gca().invert_yaxis() + plt.colorbar() + plt.show() + + # 保存地形数据 + bathym.r = np.around(bathym.r, decimals=0) + bathym.d = np.around(bathym.d, decimals=0) + bathymData = { + "param": paramData, + "count": len(bathym.d), + "descript": "距离、深度", + "data": NoIndent(list(zip(bathym.r, bathym.d))) + } + with open(("" if outParam["outPath"]=="" else outParam["outPath"] + "/")+"tlRealtimeOut_bathym.json", 'w') as jsonFile: + json_data = json.dumps(bathymData, cls=MyEncoder, indent=2, ensure_ascii=False) + jsonFile.write(json_data) + + # 保存传播损失数据 + tlpoint = [] + 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())) + tlData = { + "param": paramData, + "count": len(ranges)*len(depths), + "descript": "距离、深度、损失", + "data": NoIndent(tlpoint) + } + with open(("" if outParam["outPath"]=="" else outParam["outPath"] + "/")+"tlRealtimeOut_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() \ No newline at end of file diff --git a/wgnp.py b/wgnp.py new file mode 100644 index 0000000..4d4c169 --- /dev/null +++ b/wgnp.py @@ -0,0 +1,60 @@ +# -*- coding: utf-8 -*- +import warnings +warnings.filterwarnings('ignore') +import argparse, json + +def _argparse(): + parser = argparse.ArgumentParser(description='Wind Generated Noise Prediction') + 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: + paramData = json.load(jsonFile) + if ("dataFile" not in paramData) or ("freq" not in paramData) or\ + ("depth" not in paramData) or ("windSpeed" not in paramData) or\ + ("month" not in paramData) or\ + ("png" not in paramData) or \ + ("outParam" not in paramData): + print(f'{argfile}中的参数有误，请核实！') + return + dataFile = paramData["dataFile"] + freq = float(paramData["freq"]) + depth = float(paramData["depth"]) + windSpeed = float(paramData["windSpeed"]) + month = int(paramData["month"]) + png = paramData["png"] + outParam = paramData["outParam"] + + from pyat.WGNPd import WGNPd + from pyat.DataBase import DataBase + import matplotlib.pyplot as plt + import numpy as np + + freq = 50 * np.power(2, np.arange(freq)/3) + frequency, spl = WGNPd(dataFile, depth, windSpeed, month, freq) + + #输出json文件 + outData = {} + outData["param"] = paramData + outData["descript"] = "frequency and noise level" + frequency = np.around(frequency, decimals=0) + spl = np.around(spl, decimals=2) + outData["data"] = list(zip(frequency, spl)) + with open(("" if outParam["outPath"]=="" else outParam["outPath"] + "/")+"wgnpOut.json", 'w') as jsonFile: + json.dump(outData, jsonFile, indent=2, ensure_ascii=False) + + if png: + plt.semilogx(freq, spl) + plt.xlabel('Frequency / Hz') + plt.ylabel('Noise Level / dB') + plt.show() + +if __name__ == '__main__': + main() \ No newline at end of file