Parabolic Equation

This module contains methods for numerically solving the parabolic wave equation of Thomson and Chapman.

class kadlu.sound.parabolic_equation.Grid(dr, r_max, dz, z_max, dq=None, q_max=6.283185307179586, q=None)[source]

Bases: object

Regular grid in cylindrical coordinates r,q,z, where

r: radial coordinate (distance) in meters

z: vertical coordinate (depth) in meters

q: azimuthal coordinate (angle) in radians

Note: The sea surface is at z = 0. We use positive z values below the sea surface and negative values above.

Args:

dr: float: Radial grid spacing in meters
r_max: float: Radial range in meters
dz: float: Vertical grid spacing in meters
z_max: float: Vertical range in meters
dq: float: Angular grid spacing in radians
q_max: float: Angular range in radians
q: float, list or numpy.array: Angular coordinates in radians. If specified, dq and q_max are ignored.

Attrs:

dr: float: Radial grid spacing in meters
r: numpy.array: Radial coordinates in meters
nr: int: Number of radial grid points
dz: float: Vertical grid spacing in meters
z: numpy.array: Vertical coordinates in meters
nz: int: Number of vertical grid points
dq: float: Angular grid spacing in radians
q: numpy.array: Angular coordinates in radians
nq: int: Number of angular grid points
kz: numpy.array: Vertical wavenumber coordinates in inverse meters
below: numpy.array: Indices of z axis corresponding to values below (z>=0) the surface
q_qz, z_qz: numpy.array: Coordinates of q-z meshgrid
below_qz: numpy.array: Indices of q-z meshgrid corresponding to values below (z>=0) the surface

mirror(x, z_axis=0)[source]

Replace values above sea surface with mirrored values from below surface.

Args:

x: numpy.array: Array containing the data to be mirrored
z_axis: int: z axis of array

Returns:

x: numpy.array: Mirrored array

class kadlu.sound.parabolic_equation.TransmissionLoss(freq, bathy_func, bathy_deriv_func, sound_speed_func, bottom, propagation_range=50, angular_bin=10, c0=1500, **kwargs)[source]

Bases: object

Compute the transmission loss by solving the parabolic wave equation.

The seafloor acoustic properties are specified in a dictionary with the following keys,

sound_speed (bottom sound speed in m/s)

density (bottom density in g/cm3)

attenuation (bottom attenuation in dB/lambda)

Keyword arguments can be used to specify the computational grid kadlu.sound.parabolic_equation.Grid.

Args:

freq: float: Frequency in Hz
bathy_func: function: Bathymetry interpolation function in variables x,y
bathy_deriv_func: function: Bathymetry derivative interpolation function in variables x,y
sound_speed_func: function: Sound speed interpolation function in variables x,y,z
bottom: dict: Seafloor acoutic properties. Must contain the keys ‘sound_speed’, ‘density’, and ‘attenuation’.
propagation_range: float: Propagation range in km. Default is 50 km.
angular_bin: float: Angular bin size in degrees. Default is 10 degrees.
c0: float: Reference sound speed in m/s
source_depth: array-like: Source depths in meters.

Attributes:

Example:

calc(source_depth=None, rec_depth=[0.1], vertical=False, aperture=86, progress_bar=True, nz_max=250, nr_max=250, size_limit=True, return_field=False)[source]

Calculate the transmission loss in the horizontal plane at the specified depth(s).

Args:

source_depth: array-like: Source depths in meters
rec_depth: array-like: Receiver depths in meters
vertical: array-like: Compute the transmission loss in the vertical plane. Note: This will slow down the computation.
aperture: float: Half-beamwidth of the starter field in degrees.
progress_bar: bool: Show progress bar. Default is True
nz_max: int: Maximum number of vertical bins for output arrays
nr_max: int: Maximum number of radial bins for output arrays
size_limit: bool: Limit of the size of the output arrays to the dimensions specified in nz_max and nr_max. If False, nz_max and nr_max have no effect
return_field: bool: Instead of returning the transmission loss, TL, return the complex pressure field, f, where TL = 20 * np.log10(np.abs(f)). Default is False.

Returns:

tl_h: numpy.array: Transmission loss in dB in the horizontal plane; has shape (len(source_depth), len(rec_depth), len(q), len(r)).
ax_h: dict: Axes of the horizontal transmission loss array, (source_depth, rec_depth, q, r). q and r are the azimuthal and radial coordinate axes.
tl_v: numpy.array: Transmission loss in dB in the vertical plane; has shape (len(source_depth), len(z), len(r), len(q)). Only returned if vertical is True.
ax_v: dict: Axes of the vertical transmission loss array, (source_depth, z, r, q). z and r are the vertical and radial coordinate axes. q is the azimuthal axes. Only returned if vertical is True.

plot_horiz(source_depth_idx=0, rec_depth_idx=0)[source]

Plot the transmission loss on a horizontal plane in polar coordinates.

Args:

source_depth_idx: int: Source depth index
rec_depth_idx: int: Receive depth index

Returns:

fig: matplotlib.figure.Figure: A figure object.

plot_vert(angle=0, source_depth_idx=0, show_bathy=True, max_depth=None, show_ssp=False)[source]

Plot the transmission loss on a vertical plane in carthesian coordinates.

Args:

angle: float: Angle in degrees
source_depth_idx: int: Source depth index
show_bathy: bool: Superimpose bathymetry on plot
show_ssp: bool: Superimpose sound speed profile on plot

Returns:

fig: matplotlib.figure.Figure: A figure object.

kadlu.sound.parabolic_equation.artificial_absorp(z, z_max, D, alpha)[source]

Compute complex term, to be added to the square of the index of refraction in order to simulate the effect of an artificial absorption layer at the bottom of the physical domain.

Args:

z: numpy.array: Depth values in meters
z_max: float: Maximum depth in meters
D: float: Characteristic width of the absorption layer in meters
alpha: float: Attenuation coefficient

Returns:

: numpy.array: Absorption term, to be added to the square of the physical index of refraction.

kadlu.sound.parabolic_equation.complex_sound_speed(c, alpha)[source]

Computes complex sound speed representing the effect of volume attenuation.

Args:

c: array-like: Sound speed in m/s
alpha: float: Attenuation coefficient in dB/lambda

Returns:

ci: array-like: Complex sound speed term

kadlu.sound.parabolic_equation.eff_index_refr_sq(z, zb, dzb2, k0, L, n2, n2b, r, rb, return_density=False)[source]

Compute the effective index of refraction squared.

Args:

z: numpy.array: Depth values in meters
zb: numpy.array: Seafloor depth in meters
dzb2: numpy.array: Gradient of the seafloor depth squared
k0: float: Reference wavenumber in inverse meters
L: float: Characteristic smoothing length in meters
n2: numpy.array: Water refractive index squared
n2b: numpy.array: Bottom refractive index squared
r: numpy.array: Water density in g/cm^3
rb: numpy.array: Bottom density in g/cm^3
return_density: bool: Return the smoothened density array. Default is False.

Returns:

n2e: numpy.array: Effective index of refraction squared
rs: numpy.array: Smoothened density. Only returned if return_density is True.

kadlu.sound.parabolic_equation.index_refr_sq(c0, c, alpha=0)[source]

Compute index of refraction.

If alpha > 0, the index of refraction receives a complex term, which represents the effect of volume attenuation.

Args:

c0: float: Reference sound speed in m/s
c: array-like: Sound speed in m/s
alpha: float: Attenuation coefficient in dB/lambda

Returns:

n2: array-like: Index of refraction squared

kadlu.sound.parabolic_equation.prop_defr(x, k0, kz, nq)[source]

Compute the defractive propagation matrix, U_D.

Args:

x: float: Step size in meters
k0: float: Reference wavenumber in inverse meters
kz: array-like: Vertical wavenumber in inverse meters
nq: int: Number of angular bins

Returns:

U: numpy.array: Defractive propagation matrix, has shape (1,nz,nq)

kadlu.sound.parabolic_equation.prop_refr(x, k0, n)[source]

Compute the refractive propagation matrix, U_R.

Args:

x: float: Step size in meters
k0: float: Reference wavenumber in inverse meters
n: numpy.array: Index of refraction, has shape (nz,nq)

Returns:

U: numpy.array: Refractive propagation matrix, has shape (1,nz,nq)

kadlu.sound.parabolic_equation.smooth_density(z, zb, dzb2, L, r, rb)[source]

Smoothen the density transition from the water column to the bottom layer and compute first and second derivatives.

Args:

z: numpy.array: Depth values in meters
zb: numpy.array: Seafloor depth in meters
dzb2: numpy.array: Gradient of the seafloor depth squared
L: float: Characteristic smoothing length in meters
r: numpy.array: Water density in g/cm^3
rb: numpy.array: Bottom density in g/cm^3

Returns:

rs: numpy.array: Density in g/cm^3, with smooth transition from water to bottom
dr2: numpy.array: Square of density gradient in (g/cm^3/m)^2
d2r: numpy.array: Divergence of density gradient (g/cm^3/m^2)

kadlu.sound.parabolic_equation.smooth_index_refr(z, zb, n2, n2b, L)[source]

Smoothen the index of refraction from the water column to the bottom layer.

Args:

z: numpy.array: Depth values in meters
zb: numpy.array: Seafloor depth in meters
n2: numpy.array: Index of refraction squared in water
n2b: numpy.array: Index of refraction squared in bottom
L: float: Characteristic smoothing length in meters

Returns:

: numpy.array: Index of refraction squared, with smooth transition from water to bottom

kadlu.sound.parabolic_equation.thomson_starter(k0, kz, dz, zs, theta1)[source]

Compute Thomson starter field \psi (0, k_z) as defined in Jensen Sec. 6.4.2.3.

Args:

k0: float: Reference wavenumber in inverse meters
kz: numpy.array: Vertical wavenumber in inverse meters
dz: float: Vertical bin size in meters
zs: array-like: Source depth(s) in meters
theta1: float: Half-beamwidth in radians

Returns:

psi: numpy array: Starter field, \psi (0, k_z). Has shape (len(zs), len(kz), 1)

Example:

>>> import numpy as np
>>> from kadlu.sound.parabolic_equation import Grid, thomson_starter
>>>
>>> # Create a grid with depth of 500 m, range of 1 km,
>>> # and grid spacings of 100 meters and 10 degrees
>>> grid = Grid(100., 1000., 100., 500., 10.*np.pi/180.)
>>>
>>> # compute Thomson starter field with an aperture of 86 degrees,
>>> # for a reference wavenumber of 0.04 m^-1 and a source depth of 9 meters,
>>> psi = thomson_starter(k0=0.04, kz=grid.kz, dz=grid.dz, zs=9, theta1=86*np.pi/180.)
>>> psi = np.round(psi, 4) # round to 4 decimals
>>> print(psi[0])
[[ 0.    +0.j    ]
 [ 0.0101-0.0101j]
 [ 0.0205-0.0205j]
 [ 0.0319-0.0319j]
 [ 0.0451-0.0451j]
 [ 0.    +0.j    ]
 [-0.0451+0.0451j]
 [-0.0319+0.0319j]
 [-0.0205+0.0205j]
 [-0.0101+0.0101j]]