The following methods should be defined for MyPropagationModel:
acoustic_field
The acoustic field is represented by complex numbers with amplitude that is related to the source level (spl) and transmission loss, and angle that is related to the acoustic phase at the source frequency.
TipAdditional options
The acoustic_field(), transmission_loss(), arrivals() and impulse_response() methods may support propagation model specific keyword arguments (options) to control finer details of the propagation model.
A transmission_loss() method may be optionally defined for MyPropagationModel:
transmission_loss(pm, tx, rxs)
Compute the transmission loss from the source tx to the receivers rxs using propagation model pm. If rxs denotes a single receiver, the result is a scalar. If rxs is an AbstractArray, the result is an array of transmission losses (in dB) with the same shape as rxs.
If it is not defined, the transmission loss is automatically computed as 20 * log10(acoustic_field(...)).
A propagation model should also typically define:
arrivals(pm, tx, rx; paths=true)
Compute the arrivals at the receiver rx due to the source tx using propagation model pm. Returns an array of arrivals.
For ray models, eigenray paths are typically included in the arrivals. However, if they are not needed, one may set paths=false to allow the propagation model to avoid computing them.
The returned arrivals should be an array of arrivals that extend:
AbstractAcousticArrival
Superclass for all acoustic arrivals.
The information held in an arrival is propagation model dependent. For example, ray models may return arrivals that contain ray information such as time of arrival, angle of arrival, amplitude and phase of arrival, eigenpath, etc. On the other hand, models based on normal modes may return arrivals containing mode information such as mode number, horizontal and vertical wavenumber, etc. For ray and mode arrivals, the following concrete subtypes should be used when possible:
RayArrival
Type representing a single acoustic ray arrival.
Properties:
t / time: arrival time (s)
ϕ / phasor: complex amplitude
ns / surface_bounces: number of surface bounces
nb / bottom_bounces: number of bottom bounces
θₛ / launch_angle: launch angle at source (rad)
θᵣ / arrival_angle: arrival angle at receiver (rad)
path: ray path (optional, vector of 3-tuples or missing)
The properties are accessible with the short names for brevity, and longer more descriptive names where readability is desired or unicode symbols are undesired.
ModeArrival
Type representing a single acoustic mode arrival.
Properties:
m / mode: mode number
kᵣ / hwavenumber: horizontal wavenumber (rad/m)
ψ(z) / mode_function: mode function
v / group_velocity: group velocity (m/s)
vₚ / phase_velocity: phase velocity (m/s)
The properties are accessible with the short names for brevity, and longer more descriptive names where readability is desired or unicode symbols are undesired.
Some propagation models may be able to generate an impulse response. If so, they should define:
Compute the impulse response at the receiver rx due to the source tx using propagation model pm at the given sampling frequency fs. If abstime is true, the result is in absolute time from the start of transmission. Otherwise, the result is relative to the earliest arrival time of the signal at the receiver (possibly with some guard period to accommodate acausal response). ntaps specifies the number of taps in the impulse response. If not specified, the number of taps is chosen automatically based on the arrival times.
Compute the impulse response at the receiver rx due to the source tx using propagation model pm at the given sampling frequency fs.
Several kwargs may be specified:
If abstime is true (default: false), the result is in absolute time from the start of transmission. Otherwise, the result is relative to the earliest arrival time of the signal at the receiver (with some guard period to accommodate acausal response).
ntaps (default: nothing for automatic) specifies the number of taps in the impulse response.
fmin and fmax specifies the bandwidth of interest (default: 0 to fs/2). If the impulse response is used to convolve with bandlimited signals, it is recommended that the impulse response bandwidth be reduced to match the signal bandwidth to reduce computational load and improve numerical stability.
nmodes (default: nothing for no limit) is the maximum of modes used in estimating the duration of the impulse response. Reducing the number of modes reduces the duration of the impulse response, but may cause aliasing.
threshold (default: -60 dB) is used to drop attenuated modes to manage computational load.
acausal (default: 0.02 s) controls the computation of acausal impulse response (before the estimated time of arrival of first mode).
taper (default: 0.1) applies a tukey window to the impulse response to limit it to the estimated delay spread.
The impulse response is computed only for positive frequencies. Such an impulse response is suitable for convolution with passband complex analytic signals. If convolved with real signals, the resulting signal is approximately equivalent to converting the real signal to a complex analytic form and then convolving it with the impulse response.
If defined, the impulse response may be used to generate a channel model automatically (by calling channel()).
Channel model
If a propagation model can estimate a received signal from a transmit signal without having to compute an impulse response and convolve it, it may wish to implement the channel modeling API directly:
Compute a channel model from the sources txs to the receivers rxs using propagation model pm. The result is a channel model with the same number of input channels as the number of sources and output channels as the number of receivers. The channel model accepts signals sampled at rate fs and returns signals sampled at the same rate.
An additive noise model may be optionally specified as noise. If specified, it is used to corrupt the received signals.
Propagation model specific keyword arguments kwargs supported by impulse_response() can be passed through when generating a channel.
Simulate the transmission of passband signal x through the channel model ch. If txs is specified, it specifies the indices of the sources active in the simulation. The number of sources must match the number of channels in the input signal. If rxs is specified, it specifies the indices of the receivers active in the simulation. Returns the received signal at the specified (or all) receivers.
fs specifies the sampling rate of the input signal. The output signal is sampled at the same rate. If fs is not specified but x is a SampledSignal, the sampling rate of x is used. Otherwise, the signal is assumed to be sampled at the channel’s sampling rate.
If abstime is true, the returned signals begin at the start of transmission. Otherwise, the result is relative to the earliest arrival time of the signal at any receiver. If noisy is true and the channel has a noise model associated with it, the received signal is corrupted by additive noise.
In some cases (e.g. channel replay techniques), a channel model may be defined without the need to derive it from a propagation model. In such a case, one may extend UnderwaterAcoustics.jl by directly defining the channel model.
Scatterers and shapes
Environments may contain scatterers — objects in the water column described by a geometric shape and an acoustic boundary condition:
Create a scatterer with the given shape and acoustic boundary condition. A scatterer is an object in the water column, described by its geometry (an AbstractShape) and the acoustic property of its surface (an AbstractAcousticBoundary, e.g. RigidBoundary for a sound-hard object or PressureReleaseBoundary for a sound-soft object such as a bubble cloud). How the boundary condition is applied — including whether the object is treated as impenetrable — is determined by the propagation model.
Propagation models that support scatterers are expected to apply the boundary condition under a local tangent-plane (Kirchhoff) approximation: the reflection coefficient of the boundary is applied using the local surface normal at each interaction point. This approximation is valid when the local radius of curvature of the scatterer is much larger than the acoustic wavelength.
Propagation models that support scattering should query has_scatterers(env) and consume the scatterer geometry through the shape API described below. Models that do not support scattering should reject environments with scatterers with a clear error message (all built-in models do so).
The parametric boundary map
The fundamental geometric primitive of a shape is a parametric map of its closed boundary over a unit parameter domain. A new shape type (let’s call it MyShape) extends:
AbstractShape
Superclass for all scatterer shapes.
A shape describes the geometry of a closed object through a parametric map of its boundary over a unit parameter domain. The boundary of a 2-D shape is a curve parametrized by u ∈ [0,1], and the boundary of a 3-D shape is a surface parametrized by (u, v) ∈ [0,1]². A concrete shape must define:
Base.ndims(shape) — spatial dimension of the shape (2 or 3)
boundary_point(shape, u...) — parametric boundary position map
Optionally, a shape may provide analytic normal() and curvature() methods; if absent, they are automatically derived from boundary_point() using automatic differentiation. A shape may also override the derived queries distance(), boundary_projection(), intersect_ray(), is_inside(), boundary_points(), bounding_box() and location() with analytic implementations for speed and accuracy.
and must define just two methods:
Base.ndims(shape::MyShape) =2# spatial dimension (2 or 3)boundary_point(shape::MyShape, u) # parametric boundary position map
boundary_point(shape, u)
boundary_point(shape, u, v)
Get the position on the boundary of shape at the given boundary parameters. The boundary of a 2-D shape is parametrized by a single parameter u ∈ [0,1], and the boundary of a 3-D shape by two parameters (u, v) ∈ [0,1]². The parametric map must trace the full closed boundary as the parameters sweep the unit domain, with u = 0 and u = 1 mapping to the same point. The position is returned as an XYZ named tuple with coordinates in meters.
2-D shapes lie in the x–z plane (y = 0), with z negative downward. The boundary of a 2-D shape must be traversed counter-clockwise (with x rightward and z upward) as u increases, so that automatically derived normals point outward. Similarly, the parametrization of a 3-D shape must be right-handed (∂p/∂u × ∂p/∂v pointing outward).
Everything else is derived automatically. In particular, boundary normals and curvatures are computed from boundary_point() using automatic differentiation if the shape does not provide them. A shape may provide analytic implementations for accuracy and speed:
normal(shape, u)
normal(shape, u, v)
Get the unit outward normal to the boundary of shape at the given boundary parameters, as an XYZ named tuple. If a shape does not provide an analytic normal, it is derived from boundary_point() using automatic differentiation.
curvature(shape, u)
curvature(shape, u, v)
Get the signed curvature (in 1/m) of the boundary of a 2-D shape at boundary parameter u. The curvature is positive where the boundary curves toward the interior of the shape (convex), and negative where it curves away (concave). If a shape does not provide an analytic curvature, it is derived from boundary_point() using automatic differentiation.
The return convention for the curvature of 3-D shapes is not yet settled, and missing is returned for 3-D shapes that do not provide a curvature.
For simple shapes, one may avoid defining a new type altogether by wrapping a parametric function:
ParametricCurve(f)
Wrap a function f(u) as a 2-D shape. The function must map a boundary parameter u ∈ [0,1] to a position on the boundary of the shape, tracing a closed curve counter-clockwise in the x–z plane (with x rightward and z upward) as u increases, so that automatically derived normals point outward. The position may be given in any form accepted by xyz() (e.g. an (x, z) tuple or an XYZ named tuple). The function must be differentiable if normals, curvatures or the generic derived queries are used.
Examples
julia> circle = ParametricCurve(u -> (10 * cospi(2u), -40 + 10 * sinpi(2u)));
julia> boundary_point(circle, 0.25)
(x = 0.0, y = 0.0, z = -30.0)
ParametricSurface(f)
Wrap a function f(u, v) as a 3-D shape. The function must map boundary parameters (u, v) ∈ [0,1]² to a position on the boundary of the shape, covering the full closed surface with a right-handed parametrization (∂p/∂u × ∂p/∂v pointing outward), so that automatically derived normals point outward. The position may be given in any form accepted by xyz(). The function must be differentiable if normals are used.
Derived geometry queries
The framework provides a set of derived queries with generic fallbacks based on boundary_point(). These are properties of the shape (not of any particular solver), so any propagation model gets them for free. A shape may override any of them with analytic implementations (as Ellipse does):
distance(p1, p2)
Compute distance between two positions.
distance(shape, pos)
Get the signed distance (in meters) from position pos to the boundary of shape. The distance is negative if pos is inside the shape, positive if outside, and zero on the boundary. The magnitude is the Euclidean distance to the nearest boundary point. A shape may return a conservative lower bound on the magnitude far away from the boundary, but the distance must be exact (with the correct sign and zero level set) near the boundary.
For 2-D shapes, the distance is computed in the x–z plane, ignoring the y coordinate of pos.
boundary_projection(shape, pos)
Project position pos onto the boundary of shape, i.e., find the boundary point nearest to pos. Returns a named tuple with fields:
u: boundary parameter(s) of the nearest boundary point
point: position of the nearest boundary point
normal: unit outward normal at that point
curvature: curvature at that point
distance: signed distance from pos to the boundary (negative inside)
For 2-D shapes, the projection is computed in the x–z plane, ignoring the y coordinate of pos.
intersect_ray(shape, origin, dir)
Find the nearest intersection of the ray origin + t * dir (for t ≥ 0) with the boundary of shape. Returns nothing if the ray does not intersect the boundary, and otherwise a named tuple with fields:
t: ray parameter at the intersection (distance along the ray in units of the length of dir)
u: boundary parameter(s) of the intersection point
point: position of the intersection point
normal: unit outward normal at the intersection point
curvature: curvature at the intersection point
For 2-D shapes, the intersection is computed in the x–z plane, ignoring the y coordinates of origin and dir.
is_inside(shape, pos)
Return true if position pos is strictly inside shape, and false otherwise. For 2-D shapes, the test is performed in the x–z plane, ignoring the y coordinate of pos.
boundary_points(shape; n=128)
Sample points on the boundary of shape, e.g. for plotting. For 2-D shapes, n+1 points tracing the closed boundary are returned (with the first and last point coinciding). For 3-D shapes, points sampled on an (n+1) × (n+1) grid over the parameter domain are returned.
bounding_box(shape)
Get an axis-aligned bounding box of shape, as a named tuple (min, max) of positions. The generic implementation estimates the box from sampled boundary points, and may slightly underestimate the extent of the shape.
WarningThe distance() contract
Solvers that trace curved rays are expected to find boundary interactions as root-crossings of the signed distance() along the ray. A custom shape that overrides distance() may return a conservative lower bound on the magnitude far from the boundary, but near the boundary the distance must be sign-exact with the correct zero level set — a globally loose bound would make a solver fire the boundary event at the wrong place, or never detect it. The generic fallback and the analytic Ellipse implementation satisfy this automatically.
The generic fallbacks converge their inner minimizations to true stationary points, so ForwardDiff duals propagate correctly through them (by the envelope theorem) — custom shapes defined by a differentiable boundary_point() are fully differentiable with respect to both query positions and shape parameters.
Example: a custom shape from one function
A closed curve defined by a single (differentiable) function is a fully functional shape:
usingUnderwaterAcoustics# a lobed blob centered at (500, -40), with radius varying between 16 and 24 mblob =ParametricCurve() do u r =20+4*cospi(6u) (500+ r *cospi(2u), -40+ r *sinpi(2u))end@infois_inside(blob, (x=520, z=-40))@infodistance(blob, (x=500, z=-70))@infointersect_ray(blob, (x=0, z=-40), (x=1, z=0))
[ Info: true
[ Info: 9.095787053632662
[ Info: (t = 484.0, u = 0.5, point = (x = 484.0, y = 0.0, z = -40.0), normal = (x = -1.0, y = 0.0, z = -0.0), curvature = -0.078125)
Normals, curvatures, distances and ray intersections are all computed automatically from the parametric function.