Coverage for python/bellhop/bellhop.py: 84%

1##############################################################################

6# This file was originally part of arlpy, released under Simplified BSD License.

7# It has been relicensed in this repository to be compatible with the Bellhop licence (GPL).

9##############################################################################

11"""Underwater acoustic propagation modeling toolbox.

13This toolbox uses the Bellhop acoustic propagation model. For this model

14to work, the conplete bellhop-hub package must be built and installed

15and `bellhop.exe` should be in your PATH.

16"""

18import os as _os

19import re as _re

20import subprocess as _proc

22from tempfile import mkstemp as _mkstemp

23from struct import unpack as _unpack

24from sys import float_info as _fi

26import numpy as _np

27from scipy import interpolate as _interp

28import pandas as _pd

30import matplotlib.pyplot as _pyplt

31import matplotlib.cm as _cm

32import bokeh as _bokeh

34from bellhop.constants import _Strings, _Maps

35import bellhop.environment

36import bellhop.plotutils as _plt

38# this format to explicitly mark the functions as public:

39from bellhop.readers import read_env2d as read_env2d

40from bellhop.readers import read_ssp as read_ssp

41from bellhop.readers import read_bty as read_bty

42from bellhop.readers import read_refl_coeff as read_refl_coeff

44# models (in order of preference)

45_models = []

47def create_env2d(**kv):

48 """Create a new 2D underwater environment.

50 Parameters

51 ----------

52 **kv : dict

53 Keyword arguments for environment configuration.

55 Returns

56 -------

57 env : dict

58 A new 2D underwater environment dictionary.

60 Example

61 -------

63 To see all the parameters available and their default values:

65 >>> import bellhop as bh

66 >>> env = bh.create_env2d()

67 >>> bh.print_env(env)

69 The environment parameters may be changed by passing keyword arguments

70 or modified later using dictionary notation:

72 >>> import bellhop as bh

73 >>> env = bh.create_env2d(depth=40, soundspeed=1540)

74 >>> bh.print_env(env)

75 >>> env['depth'] = 25

76 >>> env['bottom_soundspeed'] = 1800

77 >>> bh.print_env(env)

79 The default environment has a constant sound speed.

80 A depth dependent sound speed profile be provided as a Nx2 array of (depth, sound speed):

82 >>> import bellhop as bh

83 >>> env = bh.create_env2d(depth=20,

84 >>>. soundspeed=[[0,1540], [5,1535], [10,1535], [20,1530]])

86 A range-and-depth dependent sound speed profile can be provided as a Pandas frame:

88 >>> import bellhop as bh

89 >>> import pandas as pd

90 >>> ssp2 = pd.DataFrame({

91 0: [1540, 1530, 1532, 1533], # profile at 0 m range

92 100: [1540, 1535, 1530, 1533], # profile at 100 m range

93 200: [1530, 1520, 1522, 1525] }, # profile at 200 m range

94 index=[0, 10, 20, 30]) # depths of the profile entries in m

95 >>> env = bh.create_env2d(depth=20, soundspeed=ssp2)

97 The default environment has a constant water depth. A range dependent bathymetry

98 can be provided as a Nx2 array of (range, water depth):

100 >>> import bellhop as bh

101 >>> env = bh.create_env2d(depth=[[0,20], [300,10], [500,18], [1000,15]])

102 """

103 env = bellhop.environment.new()

104 for k, v in kv.items():

105 if k not in env.keys():

106 raise KeyError('Unknown key: '+k)

107 env[k] = _np.asarray(v, dtype=_np.float64) if not isinstance(v, _pd.DataFrame) and _np.size(v) > 1 else v

108 env['depth_max'] = env['depth_max'] or _np.max(env['depth'])

109 env = check_env2d(env)

110 return env

111

112

113

114def check_env2d(env):

115 """Check the validity of a 2D underwater environment definition.

116

117 :param env: environment definition

118

119 Exceptions are thrown with appropriate error messages if the environment is invalid.

120

121 >>> import bellhop as bh

122 >>> env = bh.create_env2d()

123 >>> check_env2d(env)

124 """

125 try:

126 assert env['type'] == '2D', 'Not a 2D environment'

127 max_range = _np.max(env['rx_range'])

128 if env['surface'] is not None:

129 assert _np.size(env['surface']) > 1, 'surface must be an Nx2 array'

130 assert env['surface'].ndim == 2, 'surface must be a scalar or an Nx2 array'

131 assert env['surface'].shape[1] == 2, 'surface must be a scalar or an Nx2 array'

132 assert env['surface'][0,0] <= 0, 'First range in surface array must be 0 m'

133 assert env['surface'][-1,0] >= max_range, 'Last range in surface array must be beyond maximum range: '+str(max_range)+' m'

134 assert _np.all(_np.diff(env['surface'][:,0]) > 0), 'surface array must be strictly monotonic in range'

135 assert env['surface_interp'] == _Strings.curvilinear or env['surface_interp'] == _Strings.linear, 'Invalid interpolation type: '+str(env['surface_interp'])

136 if _np.size(env['depth']) > 1:

137 assert env['depth'].ndim == 2, 'depth must be a scalar or an Nx2 array'

138 assert env['depth'].shape[1] == 2, 'depth must be a scalar or an Nx2 array'

139 assert env['depth'][0,0] <= 0, 'First range in depth array must be 0 m'

140 assert env['depth'][-1,0] >= max_range, 'Last range in depth array must be beyond maximum range: '+str(max_range)+' m'

141 assert _np.all(_np.diff(env['depth'][:,0]) > 0), 'Depth array must be strictly monotonic in range'

142 assert env['depth_interp'] == _Strings.curvilinear or env['depth_interp'] == _Strings.linear, 'Invalid interpolation type: '+str(env['depth_interp'])

143 max_depth = _np.max(env['depth'][:,1])

144 env["_bottom_bathymetry"] = "from-file"

145 else:

146 max_depth = env['depth']

147 if isinstance(env['soundspeed'], _pd.DataFrame):

148 # For DataFrames, apply the same minimum point requirements as numpy arrays

149 if env['soundspeed_interp'] == 'spline':

150 assert env['soundspeed'].shape[0] > 3, 'soundspeed profile must have at least 4 points for spline interpolation'

151 else:

152 assert env['soundspeed'].shape[0] > 1, 'soundspeed profile must have at least 2 points'

153 assert env['soundspeed'].index[0] <= 0, 'First depth in soundspeed array must be 0 m'

154 assert env['soundspeed'].index[-1] >= max_depth, 'Last depth in soundspeed array must be beyond water depth: '+str(max_depth)+' m'

155 assert _np.all(_np.diff(env['soundspeed'].index) > 0), 'Soundspeed array must be strictly monotonic in depth'

156 elif _np.size(env['soundspeed']) > 1:

157 assert env['soundspeed'].ndim == 2, 'soundspeed must be a scalar or an Nx2 array'

158 assert env['soundspeed'].shape[1] == 2, 'soundspeed must be a scalar or an Nx2 array'

159 # Minimum points depend on interpolation type

160 if env['soundspeed_interp'] == 'spline':

161 assert env['soundspeed'].shape[0] > 3, 'soundspeed profile must have at least 4 points for spline interpolation'

162 else:

163 assert env['soundspeed'].shape[0] > 1, 'soundspeed profile must have at least 2 points'

164 assert env['soundspeed'][0,0] <= 0, 'First depth in soundspeed array must be 0 m'

165 assert env['soundspeed'][-1,0] >= max_depth, 'Last depth in soundspeed array must be beyond water depth: '+str(max_depth)+' m'

166 assert _np.all(_np.diff(env['soundspeed'][:,0]) > 0), 'Soundspeed array must be strictly monotonic in depth'

167 assert env['soundspeed_interp'] in _Maps.interp_rev, 'Invalid interpolation type: '+str(env['soundspeed_interp'])

168 if max_depth not in env['soundspeed'][:,0]:

169 indlarger = _np.argwhere(env['soundspeed'][:,0]>max_depth)[0][0]

170 if env['soundspeed_interp'] == _Strings.spline:

171 tck = _interp.splrep(env['soundspeed'][:,0], env['soundspeed'][:,1], s=0)

172 insert_ss_val = _interp.splev(max_depth, tck, der=0)

173 else:

174 insert_ss_val = _np.interp(max_depth, env['soundspeed'][:,0], env['soundspeed'][:,1])

175 env['soundspeed'] = _np.insert(env['soundspeed'],indlarger,[max_depth,insert_ss_val],axis = 0)

176 env['soundspeed'] = env['soundspeed'][:indlarger+1,:]

177 assert _np.max(env['tx_depth']) <= max_depth, 'tx_depth cannot exceed water depth: '+str(max_depth)+' m'

178 assert _np.max(env['rx_depth']) <= max_depth, 'rx_depth cannot exceed water depth: '+str(max_depth)+' m'

179 assert env['min_angle'] > -180 and env['min_angle'] < 180, 'min_angle must be in range (-180, 180)'

180 assert env['max_angle'] > -180 and env['max_angle'] < 180, 'max_angle must be in range (-180, 180)'

181 if env["bottom_reflection_coefficient"] is not None:

182 env["bottom_boundary_condition"] = "from-file"

183 if env['tx_directionality'] is not None:

184 assert _np.size(env['tx_directionality']) > 1, 'tx_directionality must be an Nx2 array'

185 assert env['tx_directionality'].ndim == 2, 'tx_directionality must be an Nx2 array'

186 assert env['tx_directionality'].shape[1] == 2, 'tx_directionality must be an Nx2 array'

187 assert _np.all(env['tx_directionality'][:,0] >= -180) and _np.all(env['tx_directionality'][:,0] <= 180), 'tx_directionality angles must be in [-90, 90]'

188 return env

189 except AssertionError as e:

190 raise ValueError(e.args)

191

192def print_env(env):

193 """Display the environment in a human readable form.

194

195 :param env: environment definition

196

197 >>> import bellhop as bh

198 >>> env = bh.create_env2d(depth=40, soundspeed=1540)

199 >>> bh.print_env(env)

200 """

201 env = check_env2d(env)

202 keys = ['name'] + sorted(list(env.keys()-['name']))

203 for k in keys:

204 v = str(env[k])

205 if '\n' in v:

206 v = v.split('\n')

207 print('%20s : '%(k) + v[0])

208 for v1 in v[1:]:

209 print('%20s '%('') + v1)

210 else:

211 print('%20s : '%(k) + v)

212

213def plot_env(env, surface_color='dodgerblue', bottom_color='peru', tx_color='orangered', rx_color='midnightblue', rx_plot=None, **kwargs):

214 """Plots a visual representation of the environment.

215

216 :param env: environment description

217 :param surface_color: color of the surface (see `Bokeh colors <https://bokeh.pydata.org/en/latest/docs/reference/colors.html>`_)

218 :param bottom_color: color of the bottom (see `Bokeh colors <https://bokeh.pydata.org/en/latest/docs/reference/colors.html>`_)

219 :param tx_color: color of transmitters (see `Bokeh colors <https://bokeh.pydata.org/en/latest/docs/reference/colors.html>`_)

220 :param rx_color: color of receviers (see `Bokeh colors <https://bokeh.pydata.org/en/latest/docs/reference/colors.html>`_)

221 :param rx_plot: True to plot all receivers, False to not plot any receivers, None to automatically decide

222

223 Other keyword arguments applicable for `bellhop.plot.plot()` are also supported.

224

225 The surface, bottom, transmitters (marker: '*') and receivers (marker: 'o')

226 are plotted in the environment. If `rx_plot` is set to None and there are

227 more than 2000 receivers, they are not plotted.

228

229 >>> import bellhop as bh

230 >>> env = bh.create_env2d(depth=[[0, 40], [100, 30], [500, 35], [700, 20], [1000,45]])

231 >>> bh.plot_env(env)

232 """

233

234 min_x = 0

235 max_x = _np.max(env['rx_range'])

236 if max_x-min_x > 10000:

237 divisor = 1000

238 min_x /= divisor

239 max_x /= divisor

240 xlabel = 'Range (km)'

241 else:

242 divisor = 1

243 xlabel = 'Range (m)'

244 if env['surface'] is None:

245 min_y = 0

246 else:

247 min_y = _np.min(env['surface'][:,1])

248 if _np.size(env['depth']) > 1:

249 max_y = _np.max(env['depth'][:,1])

250 else:

251 max_y = env['depth']

252 mgn_x = 0.01*(max_x-min_x)

253 mgn_y = 0.1*(max_y-min_y)

254 oh = _plt.hold()

255 if env['surface'] is None:

256 _plt.plot([min_x, max_x], [0, 0], xlabel=xlabel, ylabel='Depth (m)', xlim=(min_x-mgn_x, max_x+mgn_x), ylim=(-max_y-mgn_y, -min_y+mgn_y), color=surface_color, **kwargs)

257 else:

258 # linear and curvilinear options use the same altimetry, just with different normals

259 s = env['surface']

260 _plt.plot(s[:,0]/divisor, -s[:,1], xlabel=xlabel, ylabel='Depth (m)', xlim=(min_x-mgn_x, max_x+mgn_x), ylim=(-max_y-mgn_y, -min_y+mgn_y), color=surface_color, **kwargs)

261 if _np.size(env['depth']) == 1:

262 _plt.plot([min_x, max_x], [-env['depth'], -env['depth']], color=bottom_color)

263 else:

264 # linear and curvilinear options use the same bathymetry, just with different normals

265 s = env['depth']

266 _plt.plot(s[:,0]/divisor, -s[:,1], color=bottom_color)

267 txd = env['tx_depth']

268 _plt.plot([0]*_np.size(txd), -txd, marker='*', style=None, color=tx_color)

269 if rx_plot is None:

270 rx_plot = _np.size(env['rx_depth'])*_np.size(env['rx_range']) < 2000

271 if rx_plot:

272 rxr = env['rx_range']

273 if _np.size(rxr) == 1:

274 rxr = [rxr]

275 for r in _np.array(rxr):

276 rxd = env['rx_depth']

277 _plt.plot([r/divisor]*_np.size(rxd), -rxd, marker='o', style=None, color=rx_color)

278 _plt.hold(oh)

279

280def plot_ssp(env, **kwargs):

281 """Plots the sound speed profile.

282

283 :param env: environment description

284

285 Other keyword arguments applicable for `bellhop.plot.plot()` are also supported.

286

287 If the sound speed profile is range-dependent, this function only plots the first profile.

288

289 >>> import bellhop as bh

290 >>> env = bh.create_env2d(soundspeed=[[ 0, 1540], [10, 1530], [20, 1532], [25, 1533], [30, 1535]])

291 >>> bh.plot_ssp(env)

292 """

293

294 svp = env['soundspeed']

295 if isinstance(svp, _pd.DataFrame):

296 svp = _np.hstack((_np.array([svp.index]).T, _np.asarray(svp)))

297 if _np.size(svp) == 1:

298 if _np.size(env['depth']) > 1:

299 max_y = _np.max(env['depth'][:,1])

300 else:

301 max_y = env['depth']

302 _plt.plot([svp, svp], [0, -max_y], xlabel='Soundspeed (m/s)', ylabel='Depth (m)', **kwargs)

303 elif env['soundspeed_interp'] == _Strings.spline:

304 ynew = _np.linspace(_np.min(svp[:,0]), _np.max(svp[:,0]), 100)

305 tck = _interp.splrep(svp[:,0], svp[:,1], s=0)

306 xnew = _interp.splev(ynew, tck, der=0)

307 _plt.plot(xnew, -ynew, xlabel='Soundspeed (m/s)', ylabel='Depth (m)', hold=True, **kwargs)

308 _plt.plot(svp[:,1], -svp[:,0], marker='.', style=None, **kwargs)

309 else:

310 _plt.plot(svp[:,1], -svp[:,0], xlabel='Soundspeed (m/s)', ylabel='Depth (m)', **kwargs)

311

312def compute_arrivals(env, model=None, debug=False, fname_base=None):

313 """Compute arrivals between each transmitter and receiver.

314

315 :param env: environment definition

316 :param model: propagation model to use (None to auto-select)

317 :param debug: generate debug information for propagation model

318 :param fname_base: base file name for Bellhop working files, default (None), creates a temporary file

319 :returns: arrival times and coefficients for all transmitter-receiver combinations

320

321 >>> import bellhop as bh

322 >>> env = bh.create_env2d()

323 >>> arrivals = bh.compute_arrivals(env)

324 >>> bh.plot_arrivals(arrivals)

325 """

326 env = check_env2d(env)

327 (model_name, model) = _select_model(env, _Strings.arrivals, model)

328 if debug:

329 print('[DEBUG] Model: '+model_name)

330 return model.run(env, _Strings.arrivals, debug, fname_base)

331

332def compute_eigenrays(env, tx_depth_ndx=0, rx_depth_ndx=0, rx_range_ndx=0, model=None, debug=False, fname_base=None):

333 """Compute eigenrays between a given transmitter and receiver.

334

335 :param env: environment definition

336 :param tx_depth_ndx: transmitter depth index

337 :param rx_depth_ndx: receiver depth index

338 :param rx_range_ndx: receiver range index

339 :param model: propagation model to use (None to auto-select)

340 :param debug: generate debug information for propagation model

341 :param fname_base: base file name for Bellhop working files, default (None), creates a temporary file

342 :returns: eigenrays paths

343

344 >>> import bellhop as bh

345 >>> env = bh.create_env2d()

346 >>> rays = bh.compute_eigenrays(env)

347 >>> bh.plot_rays(rays, width=1000)

348 """

349 env = check_env2d(env)

350 env = env.copy()

351 if _np.size(env['tx_depth']) > 1:

352 env['tx_depth'] = env['tx_depth'][tx_depth_ndx]

353 if _np.size(env['rx_depth']) > 1:

354 env['rx_depth'] = env['rx_depth'][rx_depth_ndx]

355 if _np.size(env['rx_range']) > 1:

356 env['rx_range'] = env['rx_range'][rx_range_ndx]

357 (model_name, model) = _select_model(env, _Strings.eigenrays, model)

358 if debug:

359 print('[DEBUG] Model: '+model_name)

360 return model.run(env, _Strings.eigenrays, debug, fname_base)

361

362def compute_rays(env, tx_depth_ndx=0, model=None, debug=False, fname_base=None):

363 """Compute rays from a given transmitter.

364

365 :param env: environment definition

366 :param tx_depth_ndx: transmitter depth index

367 :param model: propagation model to use (None to auto-select)

368 :param debug: generate debug information for propagation model

369 :param fname_base: base file name for Bellhop working files, default (None), creates a temporary file

370 :returns: ray paths

371

372 >>> import bellhop as bh

373 >>> env = bh.create_env2d()

374 >>> rays = bh.compute_rays(env)

375 >>> bh.plot_rays(rays, width=1000)

376 """

377 env = check_env2d(env)

378 if _np.size(env['tx_depth']) > 1:

379 env = env.copy()

380 env['tx_depth'] = env['tx_depth'][tx_depth_ndx]

381 (model_name, model) = _select_model(env, _Strings.rays, model)

382 if debug:

383 print('[DEBUG] Model: '+model_name)

384 return model.run(env, _Strings.rays, debug, fname_base)

385

386def compute_transmission_loss(env, tx_depth_ndx=0, mode=_Strings.coherent, model=None, tx_type="default", debug=False, fname_base=None):

387 """Compute transmission loss from a given transmitter to all receviers.

388

389 :param env: environment definition

390 :param tx_depth_ndx: transmitter depth index

391 :param mode: coherent, incoherent or semicoherent

392 :param model: propagation model to use (None to auto-select)

393 :param tx_type: point or line

394 :param debug: generate debug information for propagation model

395 :param fname_base: base file name for Bellhop working files, default (None), creates a temporary file

396 :returns: complex transmission loss at each receiver depth and range

397

398 >>> import bellhop as bh

399 >>> env = bh.create_env2d()

400 >>> tloss = bh.compute_transmission_loss(env, mode=bh.incoherent)

401 >>> bh.plot_transmission_loss(tloss, width=1000)

402 """

403 env = check_env2d(env)

404 if mode not in [_Strings.coherent, _Strings.incoherent, _Strings.semicoherent]:

405 raise ValueError('Unknown transmission loss mode: '+mode)

406 if tx_type not in _Maps.source_rev:

407 raise ValueError(f'Unknown source type: {tx_type!r}')

408 if env['tx_type'] == 'default':

409 env['tx_type'] = tx_type

410 else:

411 if not(tx_type == 'default') and not(env['tx_type'] == tx_type):

412 raise ValueError('ENV file defines source type "'+env['tx_type']+'" inconsistent with Python argument tx_type="'+tx_type+'"')

413 if _np.size(env['tx_depth']) > 1:

414 env = env.copy()

415 env['tx_depth'] = env['tx_depth'][tx_depth_ndx]

416 (model_name, model) = _select_model(env, mode, model)

417 if debug:

418 print('[DEBUG] Model: '+model_name)

419 return model.run(env, mode, debug, fname_base)

420

421def arrivals_to_impulse_response(arrivals, fs, abs_time=False):

422 """Convert arrival times and coefficients to an impulse response.

423

424 :param arrivals: arrivals times (s) and coefficients

425 :param fs: sampling rate (Hz)

426 :param abs_time: absolute time (True) or relative time (False)

427 :returns: impulse response

428

429 If `abs_time` is set to True, the impulse response is placed such that

430 the zero time corresponds to the time of transmission of signal.

431

432 >>> import bellhop as bh

433 >>> env = bh.create_env2d()

434 >>> arrivals = bh.compute_arrivals(env)

435 >>> ir = bh.arrivals_to_impulse_response(arrivals, fs=192000)

436 """

437 t0 = 0 if abs_time else min(arrivals.time_of_arrival)

438 irlen = int(_np.ceil((max(arrivals.time_of_arrival)-t0)*fs))+1

439 ir = _np.zeros(irlen, dtype=_np.complex128)

440 for _, row in arrivals.iterrows():

441 ndx = int(_np.round((row.time_of_arrival.real-t0)*fs))

442 ir[ndx] = row.arrival_amplitude

443 return ir

444

445def plot_arrivals(arrivals, dB=False, color='blue', **kwargs):

446 """Plots the arrival times and amplitudes.

447

448 :param arrivals: arrivals times (s) and coefficients

449 :param dB: True to plot in dB, False for linear scale

450 :param color: line color (see `Bokeh colors <https://bokeh.pydata.org/en/latest/docs/reference/colors.html>`_)

451

452 Other keyword arguments applicable for `bellhop.plot.plot()` are also supported.

453

454 >>> import bellhop as bh

455 >>> env = bh.create_env2d()

456 >>> arrivals = bh.compute_arrivals(env)

457 >>> bh.plot_arrivals(arrivals)

458 """

459 t0 = min(arrivals.time_of_arrival)

460 t1 = max(arrivals.time_of_arrival)

461 oh = _plt.hold()

462 if dB:

463 min_y = 20*_np.log10(_np.max(_np.abs(arrivals.arrival_amplitude)))-60

464 ylabel = 'Amplitude (dB)'

465 else:

466 ylabel = 'Amplitude'

467 _plt.plot([t0, t1], [0, 0], xlabel='Arrival time (s)', ylabel=ylabel, color=color, **kwargs)

468 min_y = 0

469 for _, row in arrivals.iterrows():

470 t = row.time_of_arrival.real

471 y = _np.abs(row.arrival_amplitude)

472 if dB:

473 y = max(20*_np.log10(_fi.epsilon+y), min_y)

474 _plt.plot([t, t], [min_y, y], xlabel='Arrival time (s)', ylabel=ylabel, ylim=[min_y, min_y+70], color=color, **kwargs)

475 _plt.hold(oh)

476

477def plot_rays(rays, env=None, invert_colors=False, **kwargs):

478 """Plots ray paths.

479

480 :param rays: ray paths

481 :param env: environment definition

482 :param invert_colors: False to use black for high intensity rays, True to use white

483

484 If environment definition is provided, it is overlayed over this plot using default

485 parameters for `bellhop.plot_env()`.

486

487 Other keyword arguments applicable for `bellhop.plot.plot()` are also supported.

488

489 >>> import bellhop as bh

490 >>> env = bh.create_env2d()

491 >>> rays = bh.compute_eigenrays(env)

492 >>> bh.plot_rays(rays, width=1000)

493 """

494 rays = rays.sort_values('bottom_bounces', ascending=False)

495 max_amp = _np.max(_np.abs(rays.bottom_bounces)) if len(rays.bottom_bounces) > 0 else 0

496 if max_amp <= 0:

497 max_amp = 1

498 divisor = 1

499 xlabel = 'Range (m)'

500 r = []

501 for _, row in rays.iterrows():

502 r += list(row.ray[:,0])

503 if max(r)-min(r) > 10000:

504 divisor = 1000

505 xlabel = 'Range (km)'

506 oh = _plt.hold()

507 for _, row in rays.iterrows():

508 c = int(255*_np.abs(row.bottom_bounces)/max_amp)

509 if invert_colors:

510 c = 255-c

511 c = _bokeh.colors.RGB(c, c, c)

512 _plt.plot(row.ray[:,0]/divisor, -row.ray[:,1], color=c, xlabel=xlabel, ylabel='Depth (m)', **kwargs)

513 if env is not None:

514 plot_env(env)

515 _plt.hold(oh)

516

517def plot_transmission_loss(tloss, env=None, **kwargs):

518 """Plots transmission loss.

519

520 :param tloss: complex transmission loss

521 :param env: environment definition

522

523 If environment definition is provided, it is overlayed over this plot using default

524 parameters for `bellhop.plot_env()`.

525

526 Other keyword arguments applicable for `bellhop.plot.image()` are also supported.

527

528 >>> import bellhop as bh

529 >>> import numpy as np

530 >>> env = bh.create_env2d(

531 rx_depth=np.arange(0, 25),

532 rx_range=np.arange(0, 1000),

533 min_angle=-45,

534 max_angle=45

535 )

536 >>> tloss = bh.compute_transmission_loss(env)

537 >>> bh.plot_transmission_loss(tloss, width=1000)

538 """

539 xr = (min(tloss.columns), max(tloss.columns))

540 yr = (-max(tloss.index), -min(tloss.index))

541 xlabel = 'Range (m)'

542 if xr[1]-xr[0] > 10000:

543 xr = (min(tloss.columns)/1000, max(tloss.columns)/1000)

544 xlabel = 'Range (km)'

545 oh = _plt.hold()

546 _plt.image(20*_np.log10(_fi.epsilon+_np.abs(_np.flipud(_np.array(tloss)))), x=xr, y=yr, xlabel=xlabel, ylabel='Depth (m)', xlim=xr, ylim=yr, **kwargs)

547 if env is not None:

548 plot_env(env, rx_plot=False)

549 _plt.hold(oh)

550

551def pyplot_env(env, surface_color='dodgerblue', bottom_color='peru', tx_color='orangered', rx_color='midnightblue',

552 rx_plot=None, **kwargs):

553 """Plots a visual representation of the environment with matplotlib.

554

555 :param env: environment description

556 :param surface_color: color of the surface (see `Bokeh colors <https://bokeh.pydata.org/en/latest/docs/reference/colors.html>`_)

557 :param bottom_color: color of the bottom (see `Bokeh colors <https://bokeh.pydata.org/en/latest/docs/reference/colors.html>`_)

558 :param tx_color: color of transmitters (see `Bokeh colors <https://bokeh.pydata.org/en/latest/docs/reference/colors.html>`_)

559 :param rx_color: color of receviers (see `Bokeh colors <https://bokeh.pydata.org/en/latest/docs/reference/colors.html>`_)

560 :param rx_plot: True to plot all receivers, False to not plot any receivers, None to automatically decide

561

562 Other keyword arguments applicable for `bellhop.plot.plot()` are also supported.

563

564 The surface, bottom, transmitters (marker: '*') and receivers (marker: 'o')

565 are plotted in the environment. If `rx_plot` is set to None and there are

566 more than 2000 receivers, they are not plotted.

567

568 >>> import bellhop as bh

569 >>> env = bh.create_env2d(depth=[[0, 40], [100, 30], [500, 35], [700, 20], [1000,45]])

570 >>> bh.plot_env(env)

571 """

572

573 if _np.array(env['rx_range']).size > 1:

574 min_x = _np.min(env['rx_range'])

575 else:

576 min_x = 0

577 max_x = _np.max(env['rx_range'])

578 if max_x - min_x > 10000:

579 divisor = 1000

580 min_x /= divisor

581 max_x /= divisor

582 xlabel = 'Range (km)'

583 else:

584 divisor = 1

585 xlabel = 'Range (m)'

586 if env['surface'] is None:

587 min_y = 0

588 else:

589 min_y = _np.min(env['surface'][:, 1])

590 if _np.size(env['depth']) > 1:

591 max_y = _np.max(env['depth'][:, 1])

592 else:

593 max_y = env['depth']

594 mgn_x = 0.01 * (max_x - min_x)

595 mgn_y = 0.1 * (max_y - min_y)

596 if env['surface'] is None:

597 _pyplt.plot([min_x, max_x], [0, 0], color=surface_color, **kwargs)

598 _pyplt.xlabel(xlabel)

599 _pyplt.ylabel('Depth (m)')

600 print(min_x, mgn_x, max_x, mgn_x)

601 _pyplt.xlim([min_x - mgn_x, max_x + mgn_x])

602 _pyplt.ylim([-max_y - mgn_y, -min_y + mgn_y])

603 else:

604 # linear and curvilinear options use the same altimetry, just with different normals

605 s = env['surface']

606 _pyplt.plot(s[:, 0] / divisor, -s[:, 1], color=surface_color, **kwargs)

607 _pyplt.xlabel(xlabel)

608 _pyplt.ylabel('Depth (m)')

609 _pyplt.xlim([min_x - mgn_x, max_x + mgn_x])

610 _pyplt.ylim([-max_y - mgn_y, -min_y + mgn_y])

611 if _np.size(env['depth']) == 1:

612 _pyplt.plot([min_x, max_x], [-env['depth'], -env['depth']], color=bottom_color, **kwargs)

613 else:

614 # linear and curvilinear options use the same bathymetry, just with different normals

615 s = env['depth']

616 _pyplt.plot(s[:, 0] / divisor, -s[:, 1], color=bottom_color, **kwargs)

617 txd = env['tx_depth']

618 # print(txd, [0]*_np.size(txd))

619 _pyplt.plot([0] * _np.size(txd), -txd, marker='*', markersize=6, color=tx_color, **kwargs)

620 if rx_plot is None:

621 rx_plot = _np.size(env['rx_depth']) * _np.size(env['rx_range']) < 2000

622 if rx_plot:

623 rxr = env['rx_range']

624 if _np.size(rxr) == 1:

625 rxr = [rxr]

626 for r in _np.array(rxr):

627 rxd = env['rx_depth']

628 _pyplt.plot([r / divisor] * _np.size(rxd), -rxd, marker='o', color=rx_color, **kwargs)

629

630def pyplot_ssp(env, **kwargs):

631 """Plots the sound speed profile with matplotlib.

632

633 :param env: environment description

634

635 Other keyword arguments applicable for `bellhop.plot.plot()` are also supported.

636

637 If the sound speed profile is range-dependent, this function only plots the first profile.

638

639 >>> import bellhop as bh

640 >>> env = bh.create_env2d(soundspeed=[[ 0, 1540], [10, 1530], [20, 1532], [25, 1533], [30, 1535]])

641 >>> bh.plot_ssp(env)

642 """

643

644 svp = env['soundspeed']

645 if isinstance(svp, _pd.DataFrame):

646 svp = _np.hstack((_np.array([svp.index]).T, _np.asarray(svp)))

647 if _np.size(svp) == 1:

648 if _np.size(env['depth']) > 1:

649 max_y = _np.max(env['depth'][:, 1])

650 else:

651 max_y = env['depth']

652 _pyplt.plot([svp, svp], [0, -max_y], **kwargs)

653 _pyplt.xlabel('Soundspeed (m/s)')

654 _pyplt.ylabel('Depth (m)')

655 elif env['soundspeed_interp'] == _Strings.spline:

656 ynew = _np.linspace(_np.min(svp[:, 0]), _np.max(svp[:, 0]), 100)

657 tck = _interp.splrep(svp[:, 0], svp[:, 1], s=0)

658 xnew = _interp.splev(ynew, tck, der=0)

659 _pyplt.plot(xnew, -ynew, **kwargs)

660 _pyplt.xlabel('Soundspeed (m/s)')

661 _pyplt.ylabel('Depth (m)')

662 _pyplt.plot(svp[:, 1], -svp[:, 0], marker='.', **kwargs)

663 else:

664 _pyplt.plot(svp[:, 1], -svp[:, 0], **kwargs)

665 _pyplt.xlabel('Soundspeed (m/s)')

666 _pyplt.ylabel('Depth (m)')

667

668def pyplot_arrivals(arrivals, dB=False, color='blue', **kwargs):

669 """Plots the arrival times and amplitudes with matplotlib.

670

671 :param arrivals: arrivals times (s) and coefficients

672 :param dB: True to plot in dB, False for linear scale

673 :param color: line color (see `Bokeh colors <https://bokeh.pydata.org/en/latest/docs/reference/colors.html>`_)

674

675 Other keyword arguments applicable for `bellhop.plot.plot()` are also supported.

676

677 >>> import bellhop as bh

678 >>> env = bh.create_env2d()

679 >>> arrivals = bh.compute_arrivals(env)

680 >>> bh.plot_arrivals(arrivals)

681 """

682 t0 = min(arrivals.time_of_arrival)

683 t1 = max(arrivals.time_of_arrival)

684 if dB:

685 min_y = 20 * _np.log10(_np.max(_np.abs(arrivals.arrival_amplitude))) - 60

686 ylabel = 'Amplitude (dB)'

687 else:

688 ylabel = 'Amplitude'

689 _pyplt.plot([t0, t1], [0, 0], color=color, **kwargs)

690 _pyplt.xlabel('Arrival time (s)')

691 _pyplt.ylabel(ylabel)

692 min_y = 0

693 for _, row in arrivals.iterrows():

694 t = row.time_of_arrival.real

695 y = _np.abs(row.arrival_amplitude)

696 if dB:

697 y = max(20 * _np.log10(_fi.epsilon + y), min_y)

698 _pyplt.plot([t, t], [min_y, y], color=color, **kwargs)

699 _pyplt.xlabel('Arrival time (s)')

700 _pyplt.ylabel(ylabel)

701

702def pyplot_rays(rays, env=None, invert_colors=False, **kwargs):

703 """Plots ray paths with matplotlib

704

705 :param rays: ray paths

706 :param env: environment definition

707 :param invert_colors: False to use black for high intensity rays, True to use white

708

709 If environment definition is provided, it is overlayed over this plot using default

710 parameters for `bellhop.plot_env()`.

711

712 Other keyword arguments applicable for `bellhop.plot.plot()` are also supported.

713

714 >>> import bellhop as bh

715 >>> env = bh.create_env2d()

716 >>> rays = bh.compute_eigenrays(env)

717 >>> bh.plot_rays(rays, width=1000)

718 """

719 rays = rays.sort_values('bottom_bounces', ascending=False)

720 max_amp = _np.max(_np.abs(rays.bottom_bounces)) if len(rays.bottom_bounces) > 0 else 0

721 if max_amp <= 0:

722 max_amp = 1

723 divisor = 1

724 xlabel = 'Range (m)'

725 r = []

726 for _, row in rays.iterrows():

727 r += list(row.ray[:, 0])

728 if max(r) - min(r) > 10000:

729 divisor = 1000

730 xlabel = 'Range (km)'

731 for _, row in rays.iterrows():

732 c = _np.abs(row.bottom_bounces) / max_amp

733 if invert_colors:

734 c = 1.0 - c

735 c = _cm.gray(c)

736 if "color" in kwargs.keys():

737 _pyplt.plot(row.ray[:, 0] / divisor, -row.ray[:, 1], **kwargs)

738 else:

739 _pyplt.plot(row.ray[:, 0] / divisor, -row.ray[:, 1], color=c, **kwargs)

740 _pyplt.xlabel(xlabel)

741 _pyplt.ylabel('Depth (m)')

742 if env is not None:

743 pyplot_env(env)

744

745def pyplot_transmission_loss(tloss, env=None, **kwargs):

746 """Plots transmission loss with matplotlib.

747

748 :param tloss: complex transmission loss

749 :param env: environment definition

750

751 If environment definition is provided, it is overlayed over this plot using default

752 parameters for `bellhop.plot_env()`.

753

754 Other keyword arguments applicable for `bellhop.plot.image()` are also supported.

755

756 >>> import bellhop as bh

757 >>> import numpy as np

758 >>> env = bh.create_env2d(

759 rx_depth=np.arange(0, 25),

760 rx_range=np.arange(0, 1000),

761 min_angle=-45,

762 max_angle=45

763 )

764 >>> tloss = bh.compute_transmission_loss(env)

765 >>> bh.plot_transmission_loss(tloss, width=1000)

766 """

767 xr = (min(tloss.columns), max(tloss.columns))

768 yr = (-max(tloss.index), -min(tloss.index))

769 xlabel = 'Range (m)'

770 if xr[1] - xr[0] > 10000:

771 xr = (min(tloss.columns) / 1000, max(tloss.columns) / 1000)

772 xlabel = 'Range (km)'

773 trans_loss = 20 * _np.log10(_fi.epsilon + _np.abs(_np.flipud(_np.array(tloss))))

774 x_mesh, ymesh = _np.meshgrid(_np.linspace(xr[0], xr[1], trans_loss.shape[1]),

775 _np.linspace(yr[0], yr[1], trans_loss.shape[0]))

776 trans_loss = trans_loss.reshape(-1)

777 # print(trans_loss.shape)

778 if "vmin" in kwargs.keys():

779 trans_loss[trans_loss < kwargs["vmin"]] = kwargs["vmin"]

780 if "vmax" in kwargs.keys():

781 trans_loss[trans_loss > kwargs["vmax"]] = kwargs["vmax"]

782 trans_loss = trans_loss.reshape((x_mesh.shape[0], -1))

783 _pyplt.contourf(x_mesh, ymesh, trans_loss, cmap="jet", **kwargs)

784 _pyplt.xlabel(xlabel)

785 _pyplt.ylabel('Depth (m)')

786 _pyplt.colorbar(label="Transmission Loss(dB)")

787 if env is not None:

788 pyplot_env(env, rx_plot=False)

789

790

791def _quoted_opt(*args: str) -> str:

792 """Concatenate N input _Strings. strip whitespace, surround with single quotes

793 """

794 combined = "".join(args).strip()

795 return f"'{combined}'"

796

797

798def models(env=None, task=None):

799 """List available models.

800

801 :param env: environment to model

802 :param task: arrivals/eigenrays/rays/coherent/incoherent/semicoherent

803 :returns: list of models that can be used

804

805 >>> import bellhop as bh

806 >>> bh.models()

807 ['bellhop']

808 >>> env = bh.create_env2d()

809 >>> bh.models(env, task=coherent)

810 ['bellhop']

811 """

812 if env is not None:

813 env = check_env2d(env)

814 if (env is None and task is not None) or (env is not None and task is None):

815 raise ValueError('env and task should be both specified together')

816 rv = []

817 for m in _models:

818 if m[1]().supports(env, task):

819 rv.append(m[0])

820 return rv

821

822def _select_model(env, task, model):

823 if model is not None:

824 for m in _models:

825 if m[0] == model:

826 return (m[0], m[1]())

827 raise ValueError('Unknown model: '+model)

828 for m in _models:

829 mm = m[1]()

830 if mm.supports(env, task):

831 return (m[0], mm)

832 raise ValueError('No suitable propagation model available')

833

834def load_shd(fname_base):

835 with open(fname_base+'.shd', 'rb') as f:

836 recl, = _unpack('i', f.read(4))

837 # _title = str(f.read(80))

838 f.seek(4*recl, 0)

839 ptype = f.read(10).decode('utf8').strip()

840 assert ptype == 'rectilin', 'Invalid file format (expecting ptype == "rectilin")'

841 f.seek(8*recl, 0)

842 nfreq, ntheta, nsx, nsy, nsd, nrd, nrr, atten = _unpack('iiiiiiif', f.read(32))

843 assert nfreq == 1, 'Invalid file format (expecting nfreq == 1)'

844 assert ntheta == 1, 'Invalid file format (expecting ntheta == 1)'

845 assert nsd == 1, 'Invalid file format (expecting nsd == 1)'

846 f.seek(32*recl, 0)

847 pos_r_depth = _unpack('f'*nrd, f.read(4*nrd))

848 f.seek(36*recl, 0)

849 pos_r_range = _unpack('f'*nrr, f.read(4*nrr))

850 pressure = _np.zeros((nrd, nrr), dtype=_np.complex128)

851 for ird in range(nrd):

852 recnum = 10 + ird

853 f.seek(recnum*4*recl, 0)

854 temp = _np.array(_unpack('f'*2*nrr, f.read(2*nrr*4)))

855 pressure[ird,:] = temp[::2] + 1j*temp[1::2]

856 return _pd.DataFrame(pressure, index=pos_r_depth, columns=pos_r_range)

857

858### Bellhop propagation model ###

859

860class _Bellhop:

861

862 def __init__(self):

863 pass

864

865 def supports(self, env=None, task=None):

866 if env is not None and env['type'] != '2D':

867 return False

868 fh, fname = _mkstemp(suffix='.env')

869 _os.close(fh)

870 fname_base = fname[:-4]

871 self._unlink(fname_base+'.env')

872 rv = self._bellhop(fname_base)

873 self._unlink(fname_base+'.prt')

874 self._unlink(fname_base+'.log')

875 return rv

876

877 def run(self, env, task, debug=False, fname_base=None):

878 taskmap = {

879 _Strings.arrivals: ['A', self._load_arrivals],

880 _Strings.eigenrays: ['E', self._load_rays],

881 _Strings.rays: ['R', self._load_rays],

882 _Strings.coherent: ['C', self._load_shd],

883 _Strings.incoherent: ['I', self._load_shd],

884 _Strings.semicoherent: ['S', self._load_shd]

885 }

886 fname_flag=False

887 if fname_base is not None:

888 fname_flag = True

889

890 fname_base = self._create_env_file(env, taskmap[task][0], fname_base)

891

892 results = None

893 if self._bellhop(fname_base):

894 err = self._check_error(fname_base)

895 if err is not None:

896 print(err)

897 else:

898 try:

899 results = taskmap[task][1](fname_base)

900 except FileNotFoundError:

901 print('[WARN] Bellhop did not generate expected output file')

902 if debug:

903 print('[DEBUG] Bellhop working files: '+fname_base+'.*')

904 elif fname_flag:

905 print('[CUSTOM FILES] Bellhop working files: '+fname_base+'.*')

906 else:

907 self._unlink(fname_base+'.env')

908 self._unlink(fname_base+'.bty')

909 self._unlink(fname_base+'.ssp')

910 self._unlink(fname_base+'.ati')

911 self._unlink(fname_base+'.sbp')

912 self._unlink(fname_base+'.prt')

913 self._unlink(fname_base+'.log')

914 self._unlink(fname_base+'.arr')

915 self._unlink(fname_base+'.ray')

916 self._unlink(fname_base+'.shd')

917 return results

918

919 def _bellhop(self,*args):

920 try:

921 runcmd = f'bellhop.exe {" ".join(list(args))}'

922 #print(f"RUNNING {runcmd}")

923 result = _proc.run(runcmd,

924 stderr=_proc.STDOUT, stdout=_proc.PIPE,

925 shell=True)

926 #print(f"RETURN CODE: {result.returncode}")

927 if result.returncode == 127:

928 return False

929 except OSError:

930 return False

931 return True

932

933 def _unlink(self, f: str) -> None:

934 try:

935 _os.unlink(f)

936 except FileNotFoundError:

937 pass

938

939 def _print(self, fh, s, newline=True):

940 _os.write(fh, (s+'\n' if newline else s).encode())

941

942 def _print_array(self, fh, a, label="", nn=None):

943 if nn is None:

944 nn = _np.size(a)

945 if nn == 1:

946 self._print(fh, "1")

947 self._print(fh, f"{a} / ! {label} (single value)")

948 else:

949 self._print(fh, f"{nn}")

950 for j in a:

951 self._print(fh, f"{j} ", newline=False)

952 self._print(fh, f"/ ! {label} ({nn} values)")

953

954 def _create_env_file(self, env, taskcode, fname_base=None):

955

956 # get env file name

957 if fname_base is not None:

958 fname = fname_base+'.env'

959 fh = _os.open(_os.path.abspath(fname), _os.O_WRONLY | _os.O_CREAT | _os.O_TRUNC)

960 else:

961 fh, fname = _mkstemp(suffix='.env')

962 fname_base = fname[:-4]

963

964 self._print(fh, "'"+env['name']+"'")

965 self._print(fh, f"{env['frequency']} ! FREQ (Hz)")

966 self._print(fh, "1 ! NMedia=1 always for Bellhop")

967

968 svp = env['soundspeed']

969 svp_depth = 0.0

970 svp_interp = _Maps.interp_rev[env['soundspeed_interp']]

971 svp_boundcond = _Maps.boundcond_rev[env['top_boundary_condition']]

972 svp_attunits = _Maps.attunits_rev[env['attenuation_units']]

973 svp_volatt = _Maps.volatt_rev[env['volume_attenuation']]

974 if isinstance(svp, _pd.DataFrame):

975 svp_depth = svp.index[-1]

976 if len(svp.columns) > 1:

977 assert svp_interp == 'Q', "SVP DataFrame with multiple columns implies quadrilateral interpolation."

978 else:

979 svp = _np.hstack((_np.array([svp.index]).T, _np.asarray(svp)))

980 if env['surface'] is None:

981 comment = "SSP parameters: Interp / Top Boundary Cond / Attenuation Units / Volume Attenuation)"

982 self._print(fh, f"'{svp_interp}{svp_boundcond}{svp_attunits}{svp_volatt}' ! {comment}")

983 else:

984 comment = "SSP parameters: Interp / Top Boundary Cond / Attenuation Units / Volume Attenuation)"

985 self._print(fh, f"'{svp_interp}VWT*' ! {comment}")

986 self._create_bty_ati_file(fname_base+'.ati', env['surface'], env['surface_interp'])

987

988 max_depth = env["depth_max"] if env["depth_max"] else env['depth'] if _np.size(env['depth']) == 1 else max(_np.max(env['depth'][:,1]), svp_depth)

989

990 if env['volume_attenuation'] == _Strings.francois_garrison:

991 comment = "Francois-Garrison volume attenuation parameters (sal, temp, pH, depth)"

992 self._print(fh,f"{env['fg_salinity']} {env['fg_temperature']} {env['fg_pH']} {env['fg_depth']} ! {comment}")

993 if _np.size(svp) > 1: # why is this necessary? Included to match bellhop examples but seems erroneous/misplaced/redundant

994 self._print(fh, f"{svp[0,0]} {svp[0,1]} / ! MAXDEPTH SSP")

995

996 # max depth should be the depth of the acoustic domain, which can be deeper than the max depth bathymetry

997 comment = "DEPTH_Npts DEPTH_SigmaZ DEPTH_Max"

998 self._print(fh, f"{env['depth_npts']} {env['depth_sigmaz']} {env['depth_max']} ! {comment}")

999

1000 if _np.size(svp) == 1:

1001 self._print(fh, f"0.0 {svp} / ! '0.0' SSP_Const")

1002 self._print(fh, f"{max_depth} {svp} / ! MAXDEPTH SSP_Const")

1003 elif svp_interp == 'Q':

1004 sspenv = env['ssp_env']

1005 # if the SSP data was provided in the ENV file, use that:

1006 if sspenv is not None:

1007 for j in range(sspenv.shape[0]):

1008 self._print(fh, f"{sspenv[j,0]} {sspenv[j,1]} / ! ssp_{j}")

1009 # otherwise use the SSP data specified in the dataframe:

1010 else:

1011 for j in range(svp.shape[0]):

1012 self._print(fh, f"{svp.index[j]} {svp.iloc[j,0]} / ! ssp_{j}")

1013 self._create_ssp_file(fname_base+'.ssp', svp)

1014 else:

1015 for j in range(svp.shape[0]):

1016 self._print(fh, f"{svp[j,0]} {svp[j,1]} / ! ssp_{j}")

1017

1018 depth = env['depth']

1019 bot_bc = _Maps.boundcond_rev[env['bottom_boundary_condition']]

1020 dp_flag = _Maps.bottom_rev[env['_bottom_bathymetry']]

1021 comment = "BOT_Boundary_cond / BOT_Roughness"

1022 self._print(fh, f"{_quoted_opt(bot_bc,dp_flag)} {env['bottom_roughness']} ! {comment}")

1023

1024 comment = "DEPTH_Max BOT_SoundSpeed BOT_SoundSpeed_Shear BOT_Density [ BOT_Absorp [ BOT_Absorp_Shear ] ]"

1025 if _np.size(depth) > 1:

1026 self._create_bty_ati_file(fname_base+'.bty', depth, env['depth_interp'])

1027

1028 if env['bottom_boundary_condition'] == "acousto-elastic":

1029 if env['bottom_absorption'] is None:

1030 self._print(fh, f"{env['depth_max']} {env['bottom_soundspeed']} {env['bottom_soundspeed_shear']} {env['bottom_density']/1000} / ! {comment}")

1031 elif env['bottom_absorption_shear'] is None:

1032 self._print(fh, "%0.6f %0.6f %0.6f %0.6f %0.6f /" % (env['depth_max'], env['bottom_soundspeed'], env['bottom_soundspeed_shear'], env['bottom_density']/1000, env['bottom_absorption']))

1033 else:

1034 self._print(fh, "%0.6f %0.6f %0.6f %0.6f %0.6f %0.6f /" % (env['depth_max'], env['bottom_soundspeed'], env['bottom_soundspeed_shear'], env['bottom_density']/1000, env['bottom_absorption'], env['bottom_absorption_shear']))

1035

1036 if env['bottom_boundary_condition'] == "from-file":

1037 self._create_refl_coeff_file(fname_base+".brc", env['bottom_reflection_coefficient'])

1038

1039 self._print_array(fh, env['tx_depth'], nn=env['tx_ndepth'], label="TX_DEPTH")

1040 self._print_array(fh, env['rx_depth'], nn=env['rx_ndepth'], label="RX_DEPTH")

1041 self._print_array(fh, env['rx_range']/1000, nn=env['rx_nrange'], label="RX_RANGE")

1042

1043 beamtype = _Maps.beam_rev[env['beam_type']]

1044 beampattern = " "

1045 txtype = _Maps.source_rev[env['tx_type']]

1046 gridtype = _Maps.grid_rev[env['grid']]

1047 if env['tx_directionality'] is not None:

1048 beampattern = "*"

1049 self._create_sbp_file(fname_base+'.sbp', env['tx_directionality'])

1050 runtype_str = taskcode + beamtype + beampattern + txtype + gridtype

1051 self._print(fh, f"'{runtype_str.rstrip()}' ! RUN TYPE")

1052 self._print(fh, f"{env['nbeams']} ! NBeams")

1053 self._print(fh, "%0.6f %0.6f / ! ALPHA1,2 (degrees)" % (env['min_angle'], env['max_angle']))

1054 step_size = env["step_size"] or 0.0

1055 box_depth = env["box_depth"] or 1.01*max_depth

1056 box_range = env["box_range"] or 1.01*_np.max(env['rx_range'])

1057 self._print(fh, f"{step_size} {box_depth} {box_range/1000} ! STEP (m), ZBOX (m), RBOX (km)")

1058 _os.close(fh)

1059 return fname_base

1060

1061 def _create_bty_ati_file(self, filename, depth, interp):

1062 with open(filename, 'wt') as f:

1063 f.write("'%c'\n" % ('C' if interp == _Strings.curvilinear else 'L'))

1064 f.write(str(depth.shape[0])+"\n")

1065 for j in range(depth.shape[0]):

1066 f.write("%0.6f %0.6f\n" % (depth[j,0]/1000, depth[j,1]))

1067

1068 def _create_sbp_file(self, filename, dir):

1069 with open(filename, 'wt') as f:

1070 f.write(str(dir.shape[0])+"\n")

1071 for j in range(dir.shape[0]):

1072 f.write("%0.6f %0.6f\n" % (dir[j,0], dir[j,1]))

1073

1074 def _create_refl_coeff_file(self, filename, rc):

1075 with open(filename, 'wt') as f:

1076 f.write(str(rc.shape[0])+"\n")

1077 for j in range(rc.shape[0]):

1078 f.write(f"{rc[j,0]} {rc[j,1]} {rc[j,2]}\n")

1079

1080 def _create_ssp_file(self, filename, svp):

1081 with open(filename, 'wt') as f:

1082 f.write(str(svp.shape[1])+"\n")

1083 for j in range(svp.shape[1]):

1084 f.write("%0.6f%c" % (svp.columns[j]/1000, '\n' if j == svp.shape[1]-1 else ' '))

1085 for k in range(svp.shape[0]):

1086 for j in range(svp.shape[1]):

1087 f.write("%0.6f%c" % (svp.iloc[k,j], '\n' if j == svp.shape[1]-1 else ' '))

1088

1089 def _readf(self, f, types, dtype=str):

1090 if type(f) is str:

1091 p = _re.split(r' +', f.strip())

1092 else:

1093 p = _re.split(r' +', f.readline().strip())

1094 for j in range(len(p)):

1095 if len(types) > j:

1096 p[j] = types[j](p[j])

1097 else:

1098 p[j] = dtype(p[j])

1099 return tuple(p)

1100

1101 def _check_error(self, fname_base: str):

1102 err = None

1103 try:

1104 with open(fname_base+'.prt', 'rt') as f:

1105 for lno, s in enumerate(f):

1106 if err is not None:

1107 err += '[BELLHOP] ' + s

1108 elif '*** FATAL ERROR ***' in s:

1109 err = '\n[BELLHOP] ' + s

1110 except FileNotFoundError:

1111 pass

1112 return err

1113

1114 def _load_arrivals(self, fname_base):

1115 with open(fname_base+'.arr', 'rt') as f:

1116 hdr = f.readline()

1117 if hdr.find('2D') >= 0:

1118 freq = self._readf(f, (float,))

1119 tx_depth_info = self._readf(f, (int,), float)

1120 tx_depth_count = tx_depth_info[0]

1121 tx_depth = tx_depth_info[1:]

1122 assert tx_depth_count == len(tx_depth)

1123 rx_depth_info = self._readf(f, (int,), float)

1124 rx_depth_count = rx_depth_info[0]

1125 rx_depth = rx_depth_info[1:]

1126 assert rx_depth_count == len(rx_depth)

1127 rx_range_info = self._readf(f, (int,), float)

1128 rx_range_count = rx_range_info[0]

1129 rx_range = rx_range_info[1:]

1130 assert rx_range_count == len(rx_range)

1131 else:

1132 freq, tx_depth_count, rx_depth_count, rx_range_count = self._readf(hdr, (float, int, int, int))

1133 tx_depth = self._readf(f, (float,)*tx_depth_count)

1134 rx_depth = self._readf(f, (float,)*rx_depth_count)

1135 rx_range = self._readf(f, (float,)*rx_range_count)

1136 arrivals = []

1137 for j in range(tx_depth_count):

1138 f.readline()

1139 for k in range(rx_depth_count):

1140 for m in range(rx_range_count):

1141 count = int(f.readline())

1142 for n in range(count):

1143 data = self._readf(f, (float, float, float, float, float, float, int, int))

1144 arrivals.append(_pd.DataFrame({

1145 'tx_depth_ndx': [j],

1146 'rx_depth_ndx': [k],

1147 'rx_range_ndx': [m],

1148 'tx_depth': [tx_depth[j]],

1149 'rx_depth': [rx_depth[k]],

1150 'rx_range': [rx_range[m]],

1151 'arrival_number': [n],

1152 # 'arrival_amplitude': [data[0]*_np.exp(1j * data[1]* _np.pi/180)],

1153 'arrival_amplitude': [data[0] * _np.exp( -1j * (_np.deg2rad(data[1]) + freq[0] * 2 * _np.pi * (data[3] * 1j + data[2])))],

1154 'time_of_arrival': [data[2]],

1155 'complex_time_of_arrival': [data[2] + 1j*data[3]],

1156 'angle_of_departure': [data[4]],

1157 'angle_of_arrival': [data[5]],

1158 'surface_bounces': [data[6]],

1159 'bottom_bounces': [data[7]]

1160 }, index=[len(arrivals)+1]))

1161 return _pd.concat(arrivals)

1162

1163 def _load_shd(self, fname_base):

1164 return load_shd(fname_base)

1165

1166 def _load_rays(self, fname_base):

1167 with open(fname_base+'.ray', 'rt') as f:

1168 f.readline()

1169 f.readline()

1170 f.readline()

1171 f.readline()

1172 f.readline()

1173 f.readline()

1174 f.readline()

1175 rays = []

1176 while True:

1177 s = f.readline()

1178 if s is None or len(s.strip()) == 0:

1179 break

1180 a = float(s)

1181 pts, sb, bb = self._readf(f, (int, int, int))

1182 ray = _np.empty((pts, 2))

1183 for k in range(pts):

1184 ray[k,:] = self._readf(f, (float, float))

1185 rays.append(_pd.DataFrame({

1186 'angle_of_departure': [a],

1187 'surface_bounces': [sb],

1188 'bottom_bounces': [bb],

1189 'ray': [ray]

1190 }))

1191 return _pd.concat(rays)

1192

1193

1194

1195_models.append(('bellhop', _Bellhop))

1196

1197__all__ = [

1198 name

1199 for name in globals()

1200 if not name.startswith("_") # ignore private names

1201]