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