The underwater acoustic propagation modeling toolbox (uwapm) in arlpy is integrated with the popular Bellhop ray tracer distributed as part of the acoustics toolbox. In this notebook, we see how to use arlpy.uwapm to simplify the use of Bellhop for modeling.
Start off with checking that everything is working correctly:
import arlpy.uwapm as pm
import arlpy.plot as plt
import numpy as np
pm.models()
The bellhop model should be listed in the list of models above, if everything is good. If it isn't listed, it means that bellhop.exe is not available on the PATH, or it cannot be correctly executed. Ensure that bellhop.exe from the acoustics toolbox installation is on your PATH (updated .profile or equivalent, if necessary, to add it in).
From here on we assume that the bellhop model is available, and proceed...
We next create an underwater 2D environment (with default settings) to model:
env = pm.create_env2d()
pm.print_env(env)
We can see the default values above. Most numbers are in SI units. The only ones that aren't are:
bottom_absoption is in dB/wavelengthmin_angle and max_angle are in degreesThe default environment is an isovelocity Pekeris waveguide with a water depth of 25 m, a omnidirectional transmitter at 5 m depth, and a receiver at 10 m depth 1 km away. We can visualize the environment (not to scale):
pm.plot_env(env, width=900)
Let's simulate it and see what the eigenrays between the transmitter and receiver look like:
rays = pm.compute_eigenrays(env)
pm.plot_rays(rays, env=env, width=900)
We can also compute the arrival structure at the receiver:
arrivals = pm.compute_arrivals(env)
pm.plot_arrivals(arrivals, width=900)
The arrivals are returned in pandas data frame that we can query, if we like. For example, we can look up the time of arrival, angle of arrival, and the number of surface/bottom bounces for the first 10 entries in the arrivals list (not necessarily sorted by time):
arrivals[arrivals.arrival_number < 10][['time_of_arrival', 'angle_of_arrival', 'surface_bounces', 'bottom_bounces']]
We can also convert to a impulse response time series, if we want to use it for further signal processing:
ir = pm.arrivals_to_impulse_response(arrivals, fs=96000)
plt.plot(np.abs(ir), fs=96000, width=900)
So far, we have a simple isovelicty Pekeris waveguide. Let us next have something more interesting - an environment with some bathymetric structure and a more complicated soundspeed profile.
Let's first start off by defining our bathymetry, a steep up-slope for the first 300 m, and then a gentle downslope:
bathy = [
[0, 30], # 30 m water depth at the transmitter
[300, 20], # 20 m water depth 300 m away
[1000, 25] # 25 m water depth at 1 km
]
and then our soundspeed profile:
ssp = [
[ 0, 1540], # 1540 m/s at the surface
[10, 1530], # 1530 m/s at 10 m depth
[20, 1532], # 1532 m/s at 20 m depth
[25, 1533], # 1533 m/s at 25 m depth
[30, 1535] # 1535 m/s at the seabed
]
Now we can create our environment with a muddy bottom, and a transmitter that is at 20 m depth:
env = pm.create_env2d(
depth=bathy,
soundspeed=ssp,
bottom_soundspeed=1450,
bottom_density=1200,
bottom_absorption=1.0,
tx_depth=15
)
pm.print_env(env)
pm.plot_env(env, width=900)
We can also plot the soundspeed profile used:
pm.plot_ssp(env)
Looks more interesting! Let's see what the eigenrays look like, and also the arrival structure:
rays = pm.compute_eigenrays(env)
pm.plot_rays(rays, env=env, width=900)
We could also ignore the receiver, and plot rays launched at various angles:
rays = pm.compute_rays(env)
pm.plot_rays(rays, env=env, width=900)
or place lots of receivers in a grid to visualize the acoustic pressure field (or equivalently transmission loss). We can modify the environment (env) without having to recreate it, as it is simply a Python dictionary object:
env['rx_range'] = np.linspace(0, 1000, 1001)
env['rx_depth'] = np.linspace(0, 30, 301)
tloss = pm.compute_transmission_loss(env)
pm.plot_transmission_loss(tloss, env=env, clim=[-60,-30], width=900)
We see a complicated interference pattern, but an interesting focusing at 800 m at a 15 m depth. The detailed interference pattern is of course sensitive to small changes in the environment. A less sensitive, but more averaged out, transmission loss estimate can be obtained using the incoherent mode:
tloss = pm.compute_transmission_loss(env, mode='incoherent')
pm.plot_transmission_loss(tloss, env=env, clim=[-60,-30], width=900)
Now, let's use a directional transmitter instead of an omni-directional one:
beampattern = np.array([
[-180, 10], [-170, -10], [-160, 0], [-150, -20], [-140, -10], [-130, -30],
[-120, -20], [-110, -40], [-100, -30], [-90 , -50], [-80 , -30], [-70 , -40],
[-60 , -20], [-50 , -30], [-40 , -10], [-30 , -20], [-20 , 0], [-10 , -10],
[ 0 , 10], [ 10 , -10], [ 20 , 0], [ 30 , -20], [ 40 , -10], [ 50 , -30],
[ 60 , -20], [ 70 , -40], [ 80 , -30], [ 90 , -50], [100 , -30], [110 , -40],
[120 , -20], [130 , -30], [140 , -10], [150 , -20], [160 , 0], [170 , -10],
[180 , 10]
])
env['tx_directionality'] = beampattern
tloss = pm.compute_transmission_loss(env)
pm.plot_transmission_loss(tloss, env=env, clim=[-60,-30], width=900)
Now you can see the directionality and the sidelobe structure of the transmitter.
Finally, let's try adding a long wavelength swell on the water surface:
surface = np.array([[r, 0.5+0.5*np.sin(2*np.pi*0.005*r)] for r in np.linspace(0,1000,1001)])
env['surface'] = surface
tloss = pm.compute_transmission_loss(env)
pm.plot_transmission_loss(tloss, env=env, clim=[-60,-30], width=900)
Now, if I placed a receiver at 800 m, and 15 m depth, roughly where we see some focusing, what would the eigenrays and arrival structure look like?
env['rx_range'] = 800
env['rx_depth'] = 15
rays = pm.compute_eigenrays(env)
pm.plot_rays(rays, env=env, width=900)
arrivals = pm.compute_arrivals(env)
pm.plot_arrivals(arrivals, dB=True, width=900)
We plotted the amplitudes in dB, as the later arrivals are much weaker than the first one, and better visualized in a logarithmic scale.