From 888a80e69cc1e9c553c01f1078decb59f20271d8 Mon Sep 17 00:00:00 2001 From: Rose0428 <827933803@qq.com> Date: Fri, 7 Apr 2023 14:00:00 +0800 Subject: [PATCH] =?UTF-8?q?=E6=9B=B4=E6=94=B9=E4=BA=86=E5=A3=B0=E7=BA=BF?= =?UTF-8?q?=E8=AE=A1=E7=AE=97=E6=96=B9=E5=BC=8F=EF=BC=8C=E5=90=88=E7=90=86?= =?UTF-8?q?=E5=88=86=E9=85=8D=E4=BA=86=E8=A7=92=E5=BA=A6=E5=92=8C=E5=A3=B0?= =?UTF-8?q?=E7=BA=BF=E8=BF=9B=E8=A1=8C=E5=B9=B6=E8=A1=8C=E8=BF=90=E7=AE=97?= =?UTF-8?q?=E3=80=82?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- pyat/bellhop.py | 118 +++++++++++++++++++++++++++++++++--------------- pyat/krakenc.py | 10 ++-- pyat/ram.py | 18 ++++++-- pyat/tlcore.py | 32 ++++++++----- 4 files changed, 120 insertions(+), 58 deletions(-) diff --git a/pyat/bellhop.py b/pyat/bellhop.py index 7de2d15..d2628ee 100644 --- a/pyat/bellhop.py +++ b/pyat/bellhop.py @@ -3,7 +3,9 @@ from env import Pos, Source, Dom, cInt, SSPraw, SSP, HS, BotBndry, TopBndry, Bnd from os import system import numpy as np from config import AtBinPath - +from multiprocessing import Pool +from scipy.stats import norm +import math def write_ssp(sspfile, ssp, bathm, NZmax): if sspfile[-3:] != 'ssp': @@ -33,7 +35,9 @@ def write_bathy(btyfile, bathm): f.write(' ' + str(bathm.r[row]) + ' ' + str(bathm.d[row]) + '\r\n') return - +def call_Bellhop_p(filename): + system(AtBinPath + "/bellhop " + filename) + def calZmax(ssp): d = np.array([], dtype=int) for it in range(len(ssp)): @@ -96,20 +100,31 @@ def call_Bellhop(freq, sd, rd, bathm, ssp, sed, base): bdy = Bndry(top, bottom) # Beam params run_type = 'C' # incoherently sum beams, see AT docs for more info - nbeams = 3000 - alpha = np.linspace(-45, 45, 2) # min and max launch angle -20 degrees to 20 degrees + NAlphaRange = 12 + Alpha = alphadiv(NAlphaRange,Rmax) # min and max launch angle -20 degrees to 20 degrees + totalBeams = beamsnumber(freq, Rmax, max(bathm.r)) box = Box(Zmax, max(bathm.r)) # bound the region you let the beams go, depth in meters and range in km deltas = 0 # length step of ray trace, 0 means automatically choose - beam = Beam(RunType=run_type, Nbeams=nbeams, alpha=alpha, box=box, deltas=deltas) # package + Filenames = [] + for iAlphaRange in range(len(Alpha) - 1): + alpha = np.array([Alpha[iAlphaRange], Alpha[iAlphaRange + 1]]) + nbeams = int(totalBeams * (alpha[1] - alpha[0]) / 180) - # 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) + beam = Beam(RunType=run_type, Nbeams=nbeams, alpha=alpha, box=box, deltas=deltas) # package + filenameI = filename + str(iAlphaRange) + write_env(filenameI + '.env', 'BELLHOP', 'Pekeris profile', freq, sspB, bdy, pos, beam, cInt, Rmax) + write_ssp(filenameI, ssp, bathm, NZmax) + write_bathy(filenameI, bathm) + Filenames.append(filenameI) + + pool = Pool(NAlphaRange) + pool.map(call_Bellhop_p, Filenames) - system(AtBinPath + "/bellhop " + filename) # Read sound field - [x, x, x, x, Pos1, pressure] = read_shd(filename + '.shd') + [x, x, x, x, Pos1, pressure] = read_shd(filename + '0.shd') + for iAlphaRange in range(NAlphaRange - 1): + [x, x, x, x, Pos1, pressure1] = read_shd(filename + str(iAlphaRange) + '.shd') + pressure = pressure + pressure1 TL = -10 * np.log10(np.square(np.abs(pressure))) # Transmission Loss return Pos1, TL @@ -197,32 +212,61 @@ def find_cvgcRays(rays_total): Rays.append(ray) return Rays +def alphadiv(NalphaRange, Rmax): + angle = [-90, 90] + if Rmax > 100 and Rmax <= 1000: + sigma = 25 + elif Rmax > 1000: + sigma = 15 + else: + sigma = 40 + for i in range(1, NalphaRange): + angle.append(norm.ppf(i / NalphaRange, loc=0, scale=sigma)) + + angle.sort() + return angle + +def beamsnumber(freq, Rmax, depth): + if Rmax > 700: + nbeams = max(int(0.01 * Rmax * 1000 * freq / 1500), 300) + elif Rmax > 300 and Rmax <= 700: + nbeams = max(int(0.05 * Rmax * 1000 * freq / 1500), 300) + elif Rmax > 100 and Rmax <= 300: + nbeams = max(int(0.1 * Rmax * 1000 * freq / 1500), 300) + else: + nbeams = max(int(0.3 * Rmax * 1000 * freq / 1500), 300) + + d_theta_recommended = math.atan(depth / (10.0 * Rmax * 1000)) + nbeams = max(int(math.pi / d_theta_recommended), nbeams) + + return nbeams + # 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) +''' ********************************************** ''' +''' ******************** 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) +# [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 +# #[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/krakenc.py b/pyat/krakenc.py index 93e04e8..d1ee796 100755 --- a/pyat/krakenc.py +++ b/pyat/krakenc.py @@ -1,8 +1,8 @@ from os import system import numpy as np -from readwrite import write_env, write_fieldflp, read_shd -from env import Pos, Source, Dom, cInt, SSPraw, SSP, HS, BotBndry, TopBndry, Bndry -from config import AtBinPath +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([]) @@ -15,7 +15,7 @@ def calZmax(ssp, sed): return Zmax -def call_krakenc(freq, sd, rd, bathm, ssp, sed, base): +def call_krakenc(freq, sd, bathm, ssp, sed, base): filename = "envK" envFilename = "envK.env" model = "KRAKEN" @@ -98,7 +98,7 @@ def call_krakenc(freq, sd, rd, bathm, ssp, sed, base): system(AtBinPath + "/krakenc " + filename) # Write *.flp file Opt_flp = 'RC' - pos.r.depth = rd + pos.r.depth = RD write_fieldflp(filename, Opt_flp, pos, NProf=len(ssp), rProf=bathm.r) # Call field.exe system(AtBinPath + "/field " + filename) diff --git a/pyat/ram.py b/pyat/ram.py index dfefe10..1c992fb 100755 --- a/pyat/ram.py +++ b/pyat/ram.py @@ -1,9 +1,17 @@ -from ramRW import PosRAM, PROF, SSPRAM, write_env_RAM, readTL_RAM +from pyat.ramRW import PosRAM, PROF, SSPRAM, write_env_RAM, readTL_RAM from os import system import numpy as np -from config import AtBinPath - -def call_ram(freq, sd, rd, bathm, ssp, sed): +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_ram(freq, sd, bathm, ssp, sed): ''' ********************************************************************************* ''' ''' ************************** Sound Field Calculation ****************************** ''' ''' ********************************************************************************* ''' @@ -15,6 +23,8 @@ def call_ram(freq, sd, rd, bathm, ssp, sed): 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 + Zmax = calZmax(ssp, sed) + rd = Zmax pos = PosRAM(sd, rd, rmax, zmax, dr, dz, zmplt) # Sound speed profile for RAM diff --git a/pyat/tlcore.py b/pyat/tlcore.py index d698707..a03d494 100644 --- a/pyat/tlcore.py +++ b/pyat/tlcore.py @@ -1,12 +1,12 @@ import numpy as np -from AcousticPara import GetEnv, Coord -from DataBase import DataBase -from bellhop import call_Bellhop -from krakenc import call_krakenc -from ram import call_ram +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(database, freq, sourceDepth, sourceLat, sourceLon, receiveDepth, +def call_tl(database, freq, sourceDepth, sourceLat, sourceLon, receiveLat=None, receiveLon=None, range=None, azi=None, modelType='kraken', isPlot=False): """ Input: freq - float, Frequency @@ -40,8 +40,7 @@ def call_tl(database, freq, sourceDepth, sourceLat, sourceLon, receiveDepth, return [], [], [], bathym if modelType == 'bellhop': sd = np.array([sourceDepth]) # Source depth - rd = np.linspace(0, receiveDepth, int(receiveDepth/10)) # Receiving depth - [Pos, TL] = call_Bellhop(freq, sd, rd, bathym, ssp, sed, base) + [Pos, TL] = call_Bellhop(freq, sd, bathym, ssp, sed, base) # Plot if isPlot: plt.pcolor(Pos.r.range, Pos.r.depth, TL[0,0,:,:], vmin=40, vmax=130, cmap='jet_r', shading='auto') @@ -53,8 +52,7 @@ def call_tl(database, freq, sourceDepth, sourceLat, sourceLon, receiveDepth, return Pos.r.range, Pos.r.depth, TL[0,0,:,:], bathym elif modelType == 'kraken': sd = np.array([sourceDepth]) # Source depth - rd = np.linspace(0, receiveDepth, int(receiveDepth/10)) # Receiving depth - [Pos, TL] = call_krakenc(freq, sd, rd, bathym, ssp, sed, base) + [Pos, TL] = call_krakenc(freq, sd, bathym, ssp, sed, base) # Plot if isPlot: plt.pcolor(Pos.r.range, Pos.r.depth, TL[0,0,:,:], vmin=40, vmax=130, cmap='jet_r', shading='auto') @@ -66,7 +64,7 @@ def call_tl(database, freq, sourceDepth, sourceLat, sourceLon, receiveDepth, return Pos.r.range, Pos.r.depth, TL[0,0,:,:], bathym elif modelType == 'ram': # 调用ram模块 - [rd, rr, TL] = call_ram(freq, sourceDepth, receiveDepth, bathym, ssp, sed) + [rd, rr, TL] = call_ram(freq, sourceDepth, bathym, ssp, sed) # Plot if isPlot: plt.pcolormesh(rr, rd, TL, vmin=60, vmax=130, cmap='jet_r', shading='auto') @@ -78,4 +76,14 @@ def call_tl(database, freq, sourceDepth, sourceLat, sourceLon, receiveDepth, return rr, rd, TL, bathym else: return [], [], [], bathym - \ No newline at end of file + + + +# # 创建类 +# database = DataBase() +# # 设置时间段 +# database.set_timePeriod('06') +# # 设置时间段 +# #database.set_timePeriod(timeIndex) +# range, depth, tl, bathym = call_tl(database=database,freq=200, sourceDepth=60, sourceLat=19, sourceLon=127, +# receiveLat=17.5, receiveLon=127.1, modelType='ram', isPlot=True)