Calling from Python

Installation

While UnderwaterAcoustics.jl is written in Julia, it can be called from Python using JuliaCall. For this, we need to install the juliacall package. For some of the examples below, we will also need to numpy and matplotlib, so let’s install them as well:

pip install juliacall numpy matplotlib


Next, we will need to install UnderwaterAcoustics.jl and optionally Plots.jl. To do this, start a Python interpreter and run the following commands:

import juliapkg
juliapkg.add("UnderwaterAcoustics", "0efb1f7a-1ce7-46d2-9f48-546a4c8fbb99")
juliapkg.add("Plots", "91a5bcdd-55d7-5caf-9e0b-520d859cae80")
Tip

While Plots.jl is not necessary if you plan to use matplotlib for plotting, we install it here to demonstrate its use from Python. The plot recipes for Plots.jl that come with UnderwaterAcoustics.jl are very convenient and work nicely if you using Plots.jl from Python. If you use matplotlib, you will need to write your own plotting code based on the data returned from UnderwaterAcoustics.jl. We will show examples of both approaches below.

Getting started

Now we are ready to use UnderwaterAcoustics.jl from Python. To load the package from Python:

from juliacall import Main as jl
jl.seval("using UnderwaterAcoustics")


Let’s start with a simple example similar to the one shown in Quickstart:

env = jl.UnderwaterEnvironment(bathymetry=20, seabed=jl.SandyClay)
pm = jl.PekerisRayTracer(env)
tx = jl.AcousticSource(0.0, -5.0, 1000.0)
rx = jl.AcousticReceiver(100.0, -10.0)
rays = jl.arrivals(pm, tx, rx)


We can display the ray arrivals:

rays
7 element Vector{RayArrival}:
∠ -2.9°  0⤒  0⤓ ∠  2.9° |  65.06 ms | -40.0 dB ϕ   0.0°↝
∠  8.5°  1⤒  0⤓ ∠  8.5° |  65.71 ms | -40.1 dB ϕ 180.0°↝
∠-14.0°  0⤒  1⤓ ∠-14.0° |  66.98 ms | -62.6 dB ϕ 170.9°↝
∠ 19.3°  1⤒  1⤓ ∠-19.3° |  68.85 ms | -74.3 dB ϕ -23.4°↝
∠-24.2°  1⤒  1⤓ ∠ 24.2° |  71.26 ms | -80.4 dB ϕ-145.6°↝
∠ 28.8°  2⤒  1⤓ ∠ 28.8° |  74.16 ms | -73.5 dB ϕ  10.8°↝
∠-33.0°  1⤒  2⤓ ∠-33.0° |  77.50 ms | -100.7 dB ϕ-167.2°↝

or plot them using Plots.jl:

jl.seval("using Plots")               # python doesn't support "!", so we create
jl.seval("plot_add = plot!")          #  an alias plot_add for plot! in Julia

p = jl.plot(env, xlims=(-10,110))
jl.plot_add(tx)
jl.plot_add(rx)
jl.plot_add(rays)
jl.display(p)

We can also get the acoustic field or the transmission loss:

jl.acoustic_field(pm, tx, rx)
0.011324625487986332 + 0.013441467920062827im
jl.transmission_loss(pm, tx, rx)
35.101536894459066

We can also create a grid of receivers and compute the acoustic field over the grid:

rxs = jl.AcousticReceiverGrid2D(jl.range(1.0, 100, step=0.1), jl.range(-20, 0, step=0.1))
x = jl.transmission_loss(pm, tx, rxs)


We can display the transmission loss grid or plot it using Plots.jl:

x
991×201 Matrix{Float64}:
 19.1283  19.2328  19.5428  19.8492  …   8.82646  11.1099  16.4209  87.1227
 19.1329  19.2373  19.5471  19.8536      8.87496  11.1689  16.4855  87.123
 19.1381  19.2424  19.552   19.8586      8.92702  11.2323  16.5549  87.1233
 19.1442  19.2483  19.5577  19.8643      8.9824   11.2997  16.6288  87.1237
 19.151   19.255   19.5641  19.8709      9.04085  11.3708  16.7068  87.124
 19.1587  19.2626  19.5715  19.8784  …   9.10213  11.4454  16.7886  87.1244
 19.1674  19.2712  19.5799  19.8869      9.16604  11.5232  16.8739  87.1248
 19.1772  19.2809  19.5893  19.8965      9.23241  11.604   16.9625  87.1253
 19.1882  19.2918  19.6     19.9073      9.30109  11.6875  17.0541  87.1258
 19.2005  19.3039  19.6119  19.9194      9.37202  11.7737  17.1485  87.1263
  ⋮                                  ⋱                               ⋮
 37.3745  37.4387  37.4769  37.4921     57.8844   61.3578  67.2495  98.1624
 37.4079  37.4741  37.5144  37.5317     57.9056   61.3788  67.2688  98.1837
 37.4422  37.5104  37.5527  37.5722     57.9268   61.4001  67.2893  98.205
 37.4774  37.5476  37.592   37.6136  …  57.9481   61.4217  67.3109  98.2263
 37.5134  37.5857  37.6321  37.6558     57.9693   61.4435  67.3336  98.2477
 37.5504  37.6247  37.6732  37.6989     57.9903   61.4655  67.3573  98.269
 37.5882  37.6646  37.7151  37.7428     58.0112   61.4877  67.3819  98.2904
 37.6269  37.7053  37.7578  37.7876     58.0319   61.5099  67.4075  98.3118
 37.6665  37.7469  37.8014  37.8332  …  58.0524   61.5322  67.4338  98.3333
p = jl.plot(env, xlims=(0,100))
jl.plot_add(rxs, x)
jl.display(p)

Plotting with matplotlib

If you prefer to use matplotlib for plotting, you can use the data returned from UnderwaterAcoustics.jl to create your own plots:

import numpy as np
import matplotlib.pyplot as plt

xloss = np.flipud(np.array(x).transpose())
nz, nx = xloss.shape
x_axis = np.linspace(0, 100, nx)
z_axis = np.linspace(0, -20, nz)
plt.figure()
mesh = plt.pcolormesh(x_axis, z_axis, xloss, shading='auto', cmap='viridis_r', vmin=0, vmax=42)
cbar = plt.colorbar(mesh)
cbar.set_label('Transmission loss (dB)')
plt.xlabel('x (m)')
plt.ylabel('z (m)')
plt.tight_layout()
plt.show()

More propagation models

Using other acoustic propagation models such as Bellhop, Kraken, etc is a breeze. During installation, include the relevant Julia packages. For example, to use Bellhop or Kraken, include AcousticToolbox.jl:

import juliapkg
juliapkg.add("AcousticsToolbox", "268a15bc-5756-47d6-9bea-fa5dc21c97f8")
Note

To use models like Bellhop, Kraken, etc, you will NOT need to install the FORTRAN binaries separately. They are automatically installed when you install the relevant Julia packages.

To use the propagation models from AcousticToolbox.jl:

from juliacall import Main as jl
jl.seval("using UnderwaterAcoustics")
jl.seval("using AcousticsToolbox")

bathy = jl.SampledField(np.array([200, 150]), x=np.array([0, 1000]))
ssp = jl.SampledField(np.array([1500, 1480, 1495, 1510, 1520]), z=jl.range(0, -200, step=-50), interp=jl.CubicSpline())
env = jl.UnderwaterEnvironment(bathymetry=bathy, soundspeed=ssp, seabed=jl.SandyClay)
pm = jl.Bellhop(env)
tx = jl.AcousticSource(0.0, -5.0, 1000.0)
rx = jl.AcousticReceiver(100.0, -10.0)
rays = jl.arrivals(pm, tx, rx)
6 element Vector{RayArrival}:
∠ -9.0°  0⤒  0⤓ ∠ -3.8° | 672.12 ms | -60.5 dB ϕ   0.0°↝
∠-13.2°  0⤒  1⤓ ∠-14.3° | 677.43 ms | -87.3 dB ϕ 138.1°↝
∠ 10.3°  1⤒  0⤓ ∠  6.3° | 680.38 ms | -60.3 dB ϕ 180.0°↝
∠ 16.2°  1⤒  1⤓ ∠-18.9° | 693.31 ms | -89.9 dB ϕ-149.1°↝
∠-20.9°  1⤒  1⤓ ∠ 24.4° | 724.26 ms | -89.7 dB ϕ-165.2°↝
∠ 24.9°  2⤒  1⤓ ∠ 29.0° | 749.86 ms | -87.1 dB ϕ   7.5°↝
jl.seval("using Plots")
jl.seval("plot_add = plot!")

p = jl.plot(env, xlims=(-10,110))
jl.plot_add(tx)
jl.plot_add(rx)
jl.plot_add(rays)
jl.display(p)

Channel modeling

We can also access the channel modeling API from Python. For example, we can obtain a channel from a propagation model:

env = jl.UnderwaterEnvironment(bathymetry=20, seabed=jl.SandyClay)
pm = jl.PekerisRayTracer(env)
tx = jl.AcousticSource(0.0, -5.0, 1000.0, spl=170)
rx = jl.AcousticReceiver(100.0, -10.0)
fs = 192000
ch = jl.channel(pm, tx, rx, fs, noise=jl.RedGaussianNoise(0.5e6))


We can then create a signal and pass it through the channel:

import numpy as np

x = np.sin(2 * np.pi * 10000 * np.arange(0, 0.001, 1/fs))
y = jl.transmit(ch, x, fs=fs)  # returned signal is a Julia SampledSignal object
y = np.array(y)                # so we convert it to a numpy array


and plot it using matplotlib:

import matplotlib.pyplot as plt

fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
ax1.plot(x)
ax2.plot(y[:1000])
ax2.set_xlabel("Samples")
plt.tight_layout()
plt.show()