USML
wave_queue Class Reference
Collaboration diagram for wave_queue:

Detailed Description

Wavefront propagator for the WaveQ3D model.

Implemented using a third order Adams-Bashforth algorithm which estimates the new position and direction from the previous three iterations.

\[ \vec{r}_{n+1} = \vec{r}_{n} + \Delta t \left( \frac{23}{12} \frac{ d \vec{r} }{ dt }_n - \frac{16}{12} \frac{ d \vec{r} }{ dt }_{n-1} + \frac{ 5}{12} \frac{ d \vec{r} }{ dt }_{n-2} \right) \]

\[ \vec{\xi}_{n+1} = \vec{\xi}_{n} + \Delta t \left( \frac{23}{12} \frac{ d \vec{\xi} }{ dt }_n - \frac{16}{12} \frac{ d \vec{\xi} }{ dt }_{n-1} + \frac{ 5}{12} \frac{ d \vec{\xi} }{ dt }_{n-2} \right) \]

where:

  • $ \vec{r} $ = position in spherical earth coordinates
  • $ \vec{\xi} $ = normalized direction in spherical earth coordinates = direction divided by speed of sound
  • $ \frac{d}{dt}_n$ = time derivative at step "n"
  • $ \Delta t $ = marching algorithm time step

The actual work of computing these derivatives and the effects of the ocean environment is done by the wave_front class. The wave_front entries are stored in a circular queue so that environmental parameters computations and wavefront derivatives can be reused by subsequent time steps.

References:
S.M. Reilly, G. Potty, Sonar Propagation Modeling using Hybrid Gaussian Beams in Spherical/Time Coordinates, January 2012.

Member Enumeration Documentation

Type of spreading model to be used.

Enumerator
CLASSIC_RAY 
HYBRID_GAUSSIAN 

Constructor & Destructor Documentation

wave_queue ( ocean_model ocean,
const seq_vector freq,
const wposition1 pos,
const seq_vector de,
const seq_vector az,
double  time_step,
const wposition targets = NULL,
spreading_type  type = HYBRID_GAUSSIAN 
)

Initialize a propagation scenario.

Parameters
oceanReference to the environmental parameters.
freqFrequencies over which to compute propagation (Hz).
posLocation of the wavefront source in spherical earth coordinates.
deInitial depression/elevation angles at the source location (degrees, positive is up). Requires a minimum of 3 angles if eigenrays are being computed.
azInitial azimuthal angle at the source location (degrees, clockwise from true north). Requires a minimum of 3 angles if eigenrays are being computed. Ray fans that wrap around all azimuths should include rays for both 0 and 360 degrees.
time_stepPropagation step size (seconds).
targetsList of acoustic targets.
typeType of spreading model to use: CLASSIC_RAY or HYBRID_GAUSSIAN.
~wave_queue ( )
virtual

Destroy all temporary memory.

Member Function Documentation

bool addProplossListener ( proplossListener pListener)

Add a proplossListener to the _proplossListenerVec vector.

void build_eigenray ( unsigned  t1,
unsigned  t2,
unsigned  de,
unsigned  az,
double  distance2[3][3][3] 
)
private

Used by detect_eigenrays() to compute eigneray parameters and add a new eigenray entry to the current target.

Used by detect_eigenrays() to compute eigenray parameters and add a new eigenray entry to the current target.

Inverts the 3-D Taylor series to find offsets in time, source D/E, and source AZ that minimize the distance to the target.

\[ d^2( \vec{\rho} ) = c + \vec{b} \cdot \vec{\rho} + \frac{1}{2} \vec{\rho} \cdot A \cdot \vec{\rho} \]

\[ \frac{ \delta d^2 }{ \delta \vec{\rho} } = \vec{b} + A \cdot \vec{\rho} = 0 \]

\[ A \cdot \vec{\rho} = - \vec{b} \]

\[ \vec{\rho} = - A^{-1} \vec{b} \]

where:

  • $ \vec{ \rho } $ = offsets in time, D/E, AZ
  • $ A $ = Hessian matrix (2nd derivative)
  • $ \vec{b} $ = gradient vector (1st derivative)
  • $ c $ = constant

This routine uses these offsets and a Taylor series for direction to compute the target D/E and AZ. It avoids extrapolating source and target angles too far into shadow zone by limiting the D/E and AZ offsets to the width of one beam.

Intensity is computed using either a classic ray theory or a hybrid Gaussian beam summation across the wavefront. Attenuation is incorporated into the ray intensity using an interpolation in time along the CPA ray. Phase is copied from the CPA as are the counts for surface, bottom, and caustics.

There is an interesting degenerate case that happens for targets near the source location. The CPA calculation from detect_eigenrays() creates non-physical eigenrays around the point at which the 2nd wavefront transitions from occurring before surface reflection to including just after it. This routines examines the offsets in time and ensonified area to eliminate these non-physical eigenrays.

Parameters
t1Row number of the current target.
t2Column number of the current target.
deD/E angle index number. Can not equal a value at the edge of the ray fan.
azAZ angle index number. Can not equal a value at the edge of the ray fan.
distance2Distance squared to each of the 27 neighboring points. The first index is time, the second is D/E and the third is AZ. Warning: this array is modified during the computation.

In order to speed up the code, it required splitting the code to run faster. This then only checks the _az_boundary condition once per function call.

In order to speed up the code, it required splitting the code to run faster. This then only checks the _az_boundary condition once per function call.

void close_netcdf ( )

Close netCDF wavefront log.

void compute_offsets ( const double  distance2[3][3][3],
const c_vector< double, 3 > &  delta,
c_vector< double, 3 > &  offset,
c_vector< double, 3 > &  distance,
bool &  unstable 
)
staticprivate

Find relative offsets and true distances in time, D/E, and AZ.

Find relative offsets and true distances in time, D/E, and azimuth.

Uses the analytic solution for the inverse of a symmetric 3x3 matrix to solve

             H x = g
             x = inv(H) g

where

  • x is the relative offset,
  • g is the gradient (1st derviative) of the distances around CPA, and
  • H is the Hessian (2nd derviative) of the distances around CPA.

This implementation reverts to a simplier calculation using just the diagonals of the Hessian if the inverse can not be computed.

Compute distances from offsets, for each coordinate, by assuming that the other two offsets are zero.

             d(n)^2 = - g(n) x(n) - 0.5*H(n,n) x(n)^2

The DE distance calculations can be unstable if the target is far from CPA. This routine deals with this by computing the DE distance from the time and AZ distances.

             d(DE)^2 = d(total)^2 - d(time)^2 - d(AZ)^2
Parameters
distance2Distance squared to each of the 27 neighboring points. The first index is time, the second is D/E and the third is AZ.
deltaAxis step size in each dimension. The first index is time, the second is D/E and the third is AZ.
offsetDistance to the current target in ray coordinate units. The first index is time offset, the second is launch angle D/E offset, and the third is launch angle AZ offset. Maximum limited to +/- 1 beam. (output)
distanceDistance to the current target in world coordinate units (meters). Give the distance the same sign as the relative offset. (output)
unstableTrue if full inverse is expected to be unstable. May be updated by this routine if target is far outside of ray fan. (input/output)
const wave_front* curr ( )
inline

Return current element in the wavefront.

void detect_caustics ( )
private

Detects and processes all of the logic necessary to determine points along the wavefronts that have folded over and mark them as caustics.

Detects and processes the caustics along the next wavefront.

This logic determines if any two points have crossed over each other when going from current wavefront to the next.

void detect_eigenrays ( )
private

Detect and process wavefront closest point of approach (CPA) with target.

Requires a minimum of three rays in the D/E and AZ directions. Targets beyond the edge of the wavefront are matched to the next ray inside the fan.

void detect_reflections ( )
private

Detect and process boundary reflections and caustics.

Loops through all of the "next" wavefront elements to see if any are on the wrong side of a boundary.

Relies on detect_reflections_surface() and detect_reflections_bottom() to do the actual work of detecting and processing reflections. These routines work recursively with their opposite so that multiple reflections can take place in a single time step. This is critical in very shallow water where the reflected position may already be beyond the opposing boundary.

At the end of this process, the wave_front::find_edges() routine is used to break the wavefront down into ray families. A ray family is defined by a set of rays that have the same surface, bottom, or caustic count.

bool detect_reflections_bottom ( unsigned  de,
unsigned  az 
)
private

Detect and process reflection for a single (DE,AZ) combination.

The attenuation and phase of reflection loss are added to the values currently being stored in the next wave element. Works recursively with detect_reflections_surface() so that multiple reflections can take place in a single time step.

Parameters
deD/E angle index number.
azAZ angle index number.
Returns
True if first recursion reflects from bottom.
bool detect_reflections_surface ( unsigned  de,
unsigned  az 
)
private

Detect and process surface reflection for a single (DE,AZ) combination.

Detect and process reflection for a single (DE,AZ) combination.

The attenuation and phase of reflection loss are added to the values currently being stored in the next wave element. Works recursively with detect_reflections_bottom() so that multiple reflections can take place in a single time step.

Parameters
deD/E angle index number.
azAZ angle index number.
Returns
True if first recursion reflects from surface.
double getIntensityThreshold ( )
inline

getIntensityThreshold

Returns
Returns current value of the intensity threshold in dB
void init_netcdf ( const char *  filename,
const char *  long_name = NULL 
)

Initialize recording to netCDF wavefront log.

Opens the file and records initial conditions. The file structure is illustrated by the netCDF sample below:

  netcdf sample_test {
  dimensions:
         frequency = 1 ;
         source_de = 25 ;
         source_az = 9 ;
         travel_time = UNLIMITED ; // (### currently)
  variables:
         double frequency(frequency) ;
                 frequency:units = "hertz" ;
         double source_de(source_de) ;
                 source_de:units = "degrees" ;
                 source_de:positive = "up" ;
         double source_az(source_az) ;
                 source_az:units = "degrees_true" ;
                 source_az:positive = "clockwise" ;
         double travel_time(travel_time) ;
                 travel_time:units = "seconds" ;
         double latitude(travel_time, source_de, source_az) ;
                 latitude:units = "degrees_north" ;
         double longitude(travel_time, source_de, source_az) ;
                 longitude:units = "degrees_east" ;
         double altitude(travel_time, source_de, source_az) ;
                 altitude:units = "meters" ;
                 altitude:positive = "up" ;
         short surface(travel_time, source_de, source_az) ;
                 surface:units = "count" ;
         short bottom(travel_time, source_de, source_az) ;
                 bottom:units = "count" ;
         short caustic(travel_time, source_de, source_az) ;
                 caustic:units = "count" ;
  // global attributes:
                 :Conventions = "COARDS" ;
  data:
     frequency = 2000 ;
     source_de = 0, 1, 2, ...
     source_az = -1, 0, 1 ;
     travel_time = 0, 0.1, 0.2, ...
     latitude = 45, 45, 45, ...
     longitude = -45, -45, -45, ...
     altitude = -75, -75, -75, ...
  }
Parameters
filenameName of the file to write to disk.
long_nameOptional global attribute for identifying data-set.
void init_wavefronts ( )
private

Initialize wavefronts at the start of propagation using a 3rd order Runge-Kutta algorithm.

The Runge-Kutta algorithm is much more computationally expensive than the Adams-Bashforth algorithm used during propagation. But Runge-Kutta is self starting, only happens at initialization, and it avoids introducing the start-up errors that would be present with cheaper methods.

Assumes that all of the elements of the _curr wavefront have been initialized prior to this initialization. When this method is complete, the wavefront elements for _past, _prev, _and _next will all have valid position, direction, and closest point of approach data. However, the _next element will not have been checked for interface collisions or eigenray collisions with targets.

bool is_closest_ray ( unsigned  t1,
unsigned  t2,
unsigned  de,
unsigned  az,
const double &  center,
double  distance2[3][3][3] 
)
private

Used by detect_eigenrays() to discover if the current ray is the closest point of approach (CPA) to the current target.

Used by detect_eigenrays() to discover if the current ray is the closest point of approach to the current target.

Computes the distance to each of the 27 neighboring wavefront points for later use in the calculation of eigenray interpolation products.

Exits early if the central ray is on the edge of a ray family. Exits early if the distance to any of the surrounding points is smaller than the distance to the central point, unless that point is on the edge of the ray family, to which the search continues. Ties are awarded to the higher time, higher D/E, and higher AZ. Extrapolates outside of ray families by not testing to see if points on the edge of the family are closer to the target than the central ray.

Assumes that the wavefront has three rays in the D/E and AZ directions. Also assumes that the calling routine does not search the rays on the edges of the wavefront. This ensures that each of the tested rays has valid rays to either side of it.

Parameters
t1Row number of the current target.
t2Column number of the current target.
deD/E angle index number.
azAZ angle index number.
centerReference to the center of the distance2 cube.
distance2Distance squared to each of the 27 neighboring points. The first index is time, the second is D/E and the third is AZ (output).
Returns
True if central point is closest point of approach.

In order to speed up the code, it required splitting the code to run faster. This then only checks the _az_boundary condition once per function call.

void make_taylor_coeff ( const double  value[3][3][3],
const c_vector< double, 3 > &  delta,
double &  center,
c_vector< double, 3 > &  gradient,
c_matrix< double, 3, 3 > &  hessian,
bool  diagonal_only = false 
)
staticprivate

Computes the Taylor series coefficients used to compute eigenrays.

Computes the Taylor coefficients used to compute eigenrays.

The vector Taylor series uses the first derivative (gradient) and second derivative (Hessian) to estimate eigenray products as a function of time, D/E, and AZ from the closest point of approach (CPA).

Parameters
valueValue to interpolate around the CPA.
deltaAxis step size in each dimension. The first index is time, the second is D/E and the third is AZ.
centerValue at the center of the grid (output).
gradientFirst derivative in 3 dimensions (output).
hessianSecond derivative in 3 dimensions (output).
diagonal_onlyZeros out Hessian cross terms when true.
References:
Weisstein, Eric W. "Hessian." From MathWorld–A Wolfram Web Resource. http://mathworld.wolfram.com/Hessian.html
const wave_front* next ( )
inline

Return next element in the wavefront.

bool notifyProplossListeners ( unsigned  targetRow,
unsigned  targetCol,
eigenray  pEigenray 
)

For each proplossListener in the _proplossListenerVec vector call the addEigenray method to provide eigenrays.

unsigned num_az ( ) const
inline

Number of AZ angles in the ray fan.

unsigned num_de ( ) const
inline

Number of D/E angles in the ray fan.

const wave_front* past ( )
inline

Return past element in the wavefront.

const wave_front* prev ( )
inline

Return previous element in the wavefront.

bool removeProplossListener ( proplossListener pListener)

Remove a proplossListener from the _proplossListenerVec vector.

void save_netcdf ( )

Write current record to netCDF wavefront log.

Records travel time, latitude, longtiude, altitude for the current wavefront.

void set_bottom_reverb ( reverb_model model)

Register a bottom reverberation model.

void set_surface_reverb ( reverb_model model)

Register a surface reverberation model.

void setIntensityThreshold ( double  dThreshold)
inline

setIntensityThreshold

Parameters
dThresholdThe new value of the intensity threshold in dB.
double source_az ( unsigned  az)
inline

Initial azimuthal angle at the source location.

Parameters
azIndex of the element to access.
Returns
Depression/elevation angle. (degrees, clockwise from true north)
double source_de ( unsigned  de)
inline

Initial depression/elevation angle at the source location.

Parameters
deIndex of the element to access.
Returns
Depression/elevation angle. (degrees, positive is up)
const wposition1& source_pos ( )
inline

Location of the wavefront source in spherical earth coordinates.

Returns
Common origin of all points on the wavefront
void step ( )

Marches to the next integration step in the acoustic propagation.

Uses the third order Adams-Bashforth algorithm to estimate the position and direction of each point on the next wavefront from the previous three iterations. Automatically rotates the queue and updates the environmental parameters on the new wavefront.

Accumulated non-spreading losses are computed by combining the individual values in the next wavefront with prior losses in the current wavefront.

If targets have been specified, this function calls detect_eigenrays() at the end of each step to search for wavefront collisions with those targets.

At the end of each step, the next iteration may extend beyond one of the boundaries. This allows the eigenray calculation to accurately portray targets near the interface. Reflections are computed at the beginning of the next iteration to ensure that the next wave elements are alway inside of the water column.

double time ( )
inline

Elapsed time for the current element in the wavefront.

Friends And Related Function Documentation

friend class reflection_model
friend
friend class spreading_hybrid_gaussian
friend
friend class spreading_ray
friend

Member Data Documentation

bool _az_boundary
private

Create an Azimuthal boundary loop condition upon initialization.

This condition will prevent the production of multiple eigenrays for instances where the first azimuthal angle is equivalent to the last azimuthal angle in the AZ vector that is passed.

wave_front * _curr
private
bool _de_branch
private

Treat all targets that are slightly away from directly above the source as special cases.

const seq_vector* _frequencies
private

Frequencies over which to compute propagation loss (Hz).

Defined as a pointer to support virtual methods in seq_vector class.

double _intensity_threshold
private

The value of the intensity threshold in dB Any eigenray intensity values that are weaker than this threshold are not sent the proplossListner(s); Defaults to -300 dB.

NcVar * _nc_altitude
private
NcVar * _nc_bottom
private
NcVar * _nc_caustic
private
NcFile* _nc_file
private

The netCDF file used to record the wavefront log.

Reset to NULL when not initialized.

NcVar * _nc_latitude
private
NcVar * _nc_longitude
private
NcVar * _nc_on_edge
private
int _nc_rec
private

Current record number in netDCF file.

NcVar * _nc_surface
private
NcVar* _nc_time
private

The netCDF variables used to record the wavefront log.

wave_front * _next
private
ocean_model& _ocean
private

Reference to the environmental parameters.

Assumes that the storage for this data is managed by calling routine.

wave_front* _past
private

Circular queue of wavefront elements needed by the third order Adams-Bashforth algorithm.

  • _past is for iteration n-2
  • _prev is for iteration n-1
  • _curr is for iteration n (the current wavefront)
  • _next is for iteration n+1
wave_front * _prev
private
std::vector<proplossListener *> _proplossListenerVec
private

Vector containing the references of objects that will be used to update classes that require eigenrays as they are built.

These classes must implement addEigenray method.

reflection_model* _reflection_model
private

Reference to the reflection loss model component.

const seq_vector* _source_az
private

Initial azimuthal angle (AZ) at the source location (degrees, clockwise from true north).

Defined as a pointer to support virtual methods in seq_vector class.

const seq_vector* _source_de
private

Initial depression/elevation angle (D/E) at the source location (degrees, positive is up).

Defined as a pointer to support virtual methods in seq_vector class.

const wposition1 _source_pos
private

Location of the wavefront source in spherical earth coordinates.

spreading_model* _spreading_model
private

Reference to the spreading loss model component.

Supports either classic ray theory or Hybrid Gaussian Beams.

const wposition* _targets
private

List of acoustic targets.

matrix<double> _targets_sin_theta
private

Intermediate term: sin of colatitude for targets.

By caching this value here, we avoid re-calculating it each time the that wave_front::compute_target_distance() needs to compute the distance squared from each target to each point on the wavefront.

double _time
private

Time for current entry in the wave_front circular queue (seconds).

double _time_step
private

Propagation step size (seconds).