Acoustic Channel Simulator (January 2014)
(P. Qarabaqi and M. Stojanovic, Northeastern University)
(see also http://millitsa.coe.neu.edu/?q=projects)
Uses arrival time information from BELLHOP or a simplified multipath calculator for generating large-scale channel realizations.
arlpy, a Python interface to BELLHOP (New September 2018)
(Mandar Chitre, National University of Singapore)
A Julia interface to BELLHOP (New December 2020)
(Mandar Chitre, National University of Singapore)
https://org-arl.github.io/UnderwaterAcoustics.jl/stable/pm_bellhop.html
BELLHOP gaussian beam/finite element beam code
(M. Porter, Heat, Light, and Sound Research, Inc.)
- Fortran95 source in Acoustics Toolbox. (updated 6 July 2018)
- BELLHOP User's Guide/Manual (rough draft) (31 January 2011)
- C++ version The Jaffe Lab for Underwater Imaging, Scripps Institution of Oceanography UC San Diego
- Description of BELLHOP algorithm
- General Description (A brief user's guide) (O. Rodríguez, University of Algarve)
- World Ocean Simulation System (WOSS). A GUI for BELLHOP applied to underwater communications and networking (F. Guerra, P. Casari, M. Zorzi, U. Padova, ITALIA).
- Python version (This is a very preliminary version done as a demo by an HLS summer student. It is not for production use.)
- Python Scripts for reading/writing input/output files
(Orlando Rodriguez/Univ. of Algarve)
VirTEX: Virtual Timeseries EXperiment code (contributions from J. Peterson, M. Porter, and M. Siderius, Orig. HLS Research)
- ASA presentation
- Readme file
- Complete package 1Mb zip file
February 2013
When a source is placed on an SSP point, which is also a local
minimum then a horizontally launched ray is a pathological
situation. If its launch angle or position were slightly perturbed
then it would generate a ray that cycles repeatedly in tiny loops.
The logic for dynamically adjusting the step size along a ray has
always been confused by this and Bellhop would take
infinitessimally small steps until running out of storage space
for the ray. The new version detects this pathological case and
uses the standard stepsize. Hopefully this will not have other bad
side effects.
March 2013
A term involving the rz derivative of the SSP had been neglected
when the Quad option was used to interpolate a 2D SSP. Similarly,
curvature changes in crossing range-segments in the SSP had been
ignored. Fixed ... Logic was also added so that the ray stepsize
would be automatically adjusted to land on the ranges defining the
SSP. Previously we had required the user to control that through
the bathymetry or altimetry file.
July 2013
Changes made to make it easier to call Bellhop as a subroutine
from Matlab. In particular all the variables can easily be set
from the main routine or by reading them from an envfil.
Fixed an error in the formula that handles changes in beam
curvature during boundary reflection. The earlier formula
implicitly assumed the sound speed gradient was perpendicular to
the boundary. This was also fixed in Bellhop3D for both Nx2D and
3D options. The results are clearly improved when there is a
complicated bathymetry together with ocean SSP gradients.
February 2014
There is a further curvature correction to beams when they reflect
from a curved boundary. (Boundaries in Bellhop can either be
curved or piecewise linear.) The curvature correction had the
wrong sign for the curvature.
May 2014
The arrivals file had inaccurate receive angles in the case when
Gaussian beams were used. Bellhop was storing the receiver angle
associated with the last beam contributing to a given receiver
location. This was not very accurate. Bellhop has been modified to
use the average angle of all the contributing beams (in a group of
beams associated with a certain eigenray). Bellhop3D has also been
modified accordingly.
April 2015
The logic in the routine that calculates the beam influence for
Geometric Gaussian beams had been rearranged for greater
efficiency. In the process a factor of 2 was dropped in
calculating the intensity with the incoherent field option. Thus
incoherent calculations with Geometric Gaussian beams were coming
out 3 dB too low. This error affected the 2015 release of Bellhop.
June 2015
A user wanted to see ray traces for multiple source depths in the
same plot. A minor modification was done to permit this.
November 2015
Bellhop had used an approximate treatment of ocean volume
attenuation based directly on the cylindrical range to the
receiver. The volume attenuation should really be based on the
path length of the ray. This has been a standard approximation in
ray models; however, it does not always work well. For instance,
for HF acoustic modems in deep water being used over a vertical
path, this approximation yields 0 loss (since the cylindrical
range to the receiver is 0), thus ignoring significant loss in the
depth direction. Furthermore, Bellhop had no provision for a
depth-dependent volume attenuation. Both Bellhop and Bellhop3D
have been modified to allow a depth-dependent attenuation and to
incorporate losses in terms of the length of the ray path. The
changes to do this were simple but numerous. All the routines for
interpolating the soundspeed profile were modified to return a
complex sound speed. The option to write an arrivals file was
modified to write a complex delay where the imaginary part
represents the attenuation along the eigenray path. Matlab
routines that read the arrivals file were also modified. All the
beam influence routines were modifed to incorporate the complex
delay. A new test case ('VolAtt') was added to exercise the
new capability. With so many changes, I assume some bugs will show
up, despite the extensive testing.
April 2016
The incoherent calculations for Geometric Gaussian beams were
still about 1.6 dB too low. I had assumed loosely that two
adjacent beams were contributing. I modified the formulat to use
the sum of an infinite number of beams (with increasing distance
form the receiver). This was actually the process that had always
been used for the *coherent* sum, but I had been less careful with
the incoherent sum. The updated formula now increases to a
fracyion of a dB with what had been in the Geometric Hat beams.
Further refinements are possible for refracted media.
July 2016
Fixed an error in Bellhop3D using Geometric Hat beams in Cartesian
coordinates: the old version would generally have a beam make a
contribution to a radial and its continuation in the other half
plane. For instance, if the beam made a contribution to the
45-degree radial, it would also show an erroneous contribution to
the 225-degree radial.
August 2016
Added an option for geometric hat beams in ray-centered
coordinates (rather than Cartesian coordinates). The geometric
beams in Cartesian coordinates have been a standard. Only the
Cerveny-style beams had been implemented in ray-centered
coordinates. The interpolation of the beams works better with
ray-centered coordinates, especially for beams that are nearly
vertical. However, the Cartesian beams require less computation.
September 2016
Fixed a bug in Bellhop3D in which the 2D option was not clearing
the pressure field for a radial. As a result, the field on new
radials had the field in storage from the previous radial in areas
that would have been shadows.
Modified the Gaussian beams in ray-centered coordinates to use the
same algorithm as that used for hat beams in ray-centered
coordinates. The new code is faster by about 4x. (That factor
depends on lots of particulars of the problem and the set up.)
Bellhop reads the SSP data to a user-specified depth. That depth
was being converted to single precision and then in a rare case
did not produce a match at the last SSP point due to the lower
precision. Fixed.
Added (but did not test) curvilinear altimetry option based on the
analogous bathymetry option.
For Geometric, (Gaussian and hat beams) in ray-centered
coordinates the 'Arrivals' option was incorrectly using a ray
matrix from the 2D case (ray2D vs. ray3D). Fixed.
October 2016
The Bellhop3D altimetry option had probably never been tested. I
created a directory 'ParaBot' with a parabolic altimetry and
bathymetry file based on the same case in the Bellhop tests. It
turned out this had never worked properly. The logic is, of
course, almost identical to the bathymetry case but it wasn't
quite right. Based on the ParaBot test it looks like it is now
working correctly. The pressure field is symmetric about the
horizontal axis, indicating the altimetry and bathymetry files are
being treated the same way. the match with the normal Bellhop is
not quite perfect and deserves another look. However, it is
awfully close.
Bellhop and Bellhop3D have always ignored the density in the
volume assuming it was seawater and could be treated as a
constant. That constant was taken as 1.0 gm/cm^3, which is very
slightly off. In running an aero-acoustic problem for the air that
density is about 1000x too high which led to an incorrect boundary
reflection coefficient. Both models have been modified to use the
density as specified in the environmental file including its depth
variation. This effects about 6 different types of interpolation
options in both models; however, I have just tested the c-linear
option and only in Bellhop.
July 2017
The logic in Bellhop3D had always been sloppy in terms of how it
would handle a ray that exited the domain where the bathymetry or
altimetry was defined. This could result in accesses to undefined
variables. (Bellhop handles this by extending the domain as a
constant.) Bellhop3D has been modified to terminate the ray trace
for the particular take-off angle involved when this happens. The
ray trace then continues to the next launch angle.
Often problems in the ray trace occur because of user errors in
defining the bathymetry. For instance, NaNs are often present.
Bellhop3D now checks for such NaNs when they're first read and
generates a warning message. It also tries to detect when a ray
has entered a patch where the bathymetry is improperly defined
(because of the Nan) and then terminates the trace of that
particular ray.
A problem was also present in that Bellhop3D did not always check
for each new ray which bathymetry/altimetry segment the source was
over/under. This has been fixed; however, I need to review this
since Bellhop uses a different process where it *always* searches
for the active segment. This is less efficient but safer. In any
case, Bellhop and Bellhop3D should likely be set up the same way.
The Bellhop3D ASCII arrivals file has been modified so that it is
stylistically consistent with the Bellhop version. (Previously the
delay information was written as a complex number by Bellhop3D vs.
two reals in Bellhop.)
December 2017
The arrivals option for the case of an Nx2D run had not been
implemented. This feature was added so that an arrivals
calculation can be done with either Nx2D or 3D options.
May 2018
The geometric Gaussian beam options in Bellhop and Bellhop3D use a
'stent' which sets a lower limit on how narrow a beam can be. That
lower limit is based on arc-length from the source so that it
allows narrow beams close to the source. Otherwise the beams add
up to too high a level. Bellhop and Bellhop3D had been using just
a coarse estimate of arc-length based on the number of steps times
the step-size. However, that estimate was way off in cases where
the default step-size is chosen to be really large. That is
normally OK because the codes adjust the step-size down
automatically as the ray steps along. However, the number of steps
times the step-size then overestimates the arc-length in a big
way. Both codes have been modified to estimate the arc-length
using the travel time delay, which is tracked as part of the ray
trace and incorporates the dynamic changes in step-length.
The semi-coherent option in Bellhop was over-riding any source
beam pattern if that was provided. Fixed.
The Gaussian beam option in Bellhop3D (in 3D mode) now can be done
in the original ray-centered coordinates or in Cartesian
coordinates. The latter is vastly faster.
June 2018
The treatment of volume attenuation in the Matlab version of
Bellhop has been upgraded to be consistent with the Fortran
version. (See note dated November 2015.)
The fix of Feb. 2013 was also incorporated in BellhopM
July 2018
The Matlab version of Bellhop (BellhopM) has the wrong sign on the
phase changes that happen when a ray crosses a caustic (because I
used ' instead of .' on the associated vector, inadvertently
conjugating it). Fixed. The arc-length correction from May 2018
was also implemented in BellhopM.
A problem showed up in Bellhop3D when rays were traced over a
bathymetry given in Eastings and Northings. Those values are in
the millions, but described (in this case) a small patch. The ray
steps are around 1 m and when written in single precision to the
ray file, the ray coordinates in Easting and Northings don't have
enough precision. The rays then look jittery. The fix was simply
to write the ray file in double precision. This was changed in the
regulard 2D Bellhop as well, although it's less likely one would
create a Bellhop file with that kind of problem.
There were lines in WriteArrivals for the 3D case that factored in
1 / sqrt( r ) for a point source and something else for a line
source. Neither one makes sense in a 3D run so they were removed.
The 3D code does the equivalent of a sqrt( r) through the actual
spreading of the beams in each direction; a line source would
require a different treatment because of the 3D.
The format for ARRFiles generated by Bellhop3D has been modified
to write the number of bearings (as well as receiver depths and
ranges). Azimuthal angles at the source and receiver were also
added to the ARRFile.
When Bellhop3D was run in the Nx2D mode it was neglecting to make
the normals to the boundaries have unit length. (They were
projections of the 3D normal vectors onto the 2D propagation
plane; that process changed the length of the normal vectors.) As
a result the angles of reflection were slightly off producing
errors in the results.
August 2018
The previous version created pyramidal beams with a rectangular
base. Thus rotations of the q_tilde, q_hat coordinates were
ignored. The new version takes these into account yielding
pyramidal beams with a parallelogram base. The resulting beams
then produce what could be called a tiling, with essentially no
gaps or overlaps between adjacent beams. A lot of work was also
done to handle boundary curvature. When the curvilinear boundary
option is selected it now automatically calculates the boundary
curvature in each direction and applies jumps to the q functions
(related to beam curvature) when the beams reflect off those
boundaries. The previous version only incorporated curvature in
calculating the angle of the reflected rays.
November 2018
The option for a bottom type that varied with range was not
working properly. A confusion in the order of the real and
imaginary parts of the P-wave and S-wave speeds was causing the
sound speed attenuation to be ignored. An additional bug had been
introduced when the broadband option was introduced into the
acoustic toolbox-- the routine for converting attenuation to a
complex sound speed was missing a newly introduced variable. Fixed
...
December 2018
The density interpolation with the hexahedral option was not
interpolating density correctly. Since a varying density is rarely
used in ocean acoustics this problem would likely not have shown
up in any practical case. A modification has been made so that it
is now interpolated ln linearly in depth. However, the
interpolation is piecewise constant in x and y.
January 2019
PCHIP, Spline, and N2-linear SSP options added to the Matlab
version of BELLHOP.
February 2019
2D or 3D soundspeeds are read in from an auxillary SSPFILE. In
that case, a depth-dependent volume attenuation can be input in
the envfil. However, the attenuation was being interpolated
incorrectly due to a confusion between depth into a layer and the
proportional depth in the layer. This has been fixed but not
tested. Note that the depths for the SSP values must match between
the ENVFILE and the SSPFILE. There is currently no check to
enforce that.
The François-Garrison formulas for volume attenuation have been
added as an alternative to the older Thorp formula. In the process
of doing this I found that I had volume attenuation incorporated
twice in the Matlab version of BELLHOP. This error was introduced
in June 2018 (see above) when I switched from handling attenuation
by a simple formula based on cylindrical range to a more accurate
one based on attenuation along the arc of the ray path. I failed
to remove the older range formula when I added the new one.
May 2019
The incoherent option with paraxial beams (Cerveny) was producing
levels that were way too low. The code was summing the absolute
values of the pressures, rather than the intensities (pressure
squared). With that corrected for both beams in Cartesian
coordinates and ray-centered beams, the levels are better, but
they still seem a bit off. However, the idea of incoherently
summing beams that are in the same family does not really make
sense.
In the case of a full 360 degree sweep of rays, BELLHOP and
BELLHOP3D were supposed to automatically remove a duplicate beam
at, e.g. -180 and 180 degrees. The coding was done lazily, testing
for exact equality and with round-off that test was not always
satisfied. The code has been modified to test for approximate
equality. This modification is of little importance, but it does
eliminate an artifact in at/tests/BeamPattern/omni.
BELLHOP3D was not including the volume attenuation when the
incoherent/semi-coherent option was selected.
June 2019
Elasticity in the halfspaces has been included. Note that BOUNCE
can be used to generate reflection coefficients for elastic
halfspaces and even more complicated layered sub-bottoms. Then the
resulting reflection coefficient can be used in BELLHOP. However,
this modification simplifies the process for halfspaces.
July 2020
There is a hidden option in BELLHOP to include a beam shift on
boundary reflection. (It is hidden or undocumented because it is
not considered to be reliable.) That option was being ignored when
a ray trace run was selected. As a result, the ray plots did not
reflect the true trajectories that were using in a TL calculation.
Fixed ... This beamshift option is available in both BELLHOP and
BELLHOP3D, but will not produce reliable results with complicated
bottoms. There is no logic in the software to deal with the fact
that a displaced ray/beam may be moved to the next bathymetric
domain. However, this change at least ensures the ray plots will
represent the underlying behavior for the TL.
August 2020
In the test case at/tests/BeamPattern/omni.env there are
homogeneous lossless halfspaces for the top and bottom.
Calculating the boundary impedance uses a square root to get a
vertical wavenumber. In this special case it it taking the square
root of a negative real. It turns out the gfortran and Intel
compilers do different things with those negative reals (which lie
on the the branch cut of the square root). They both give a purely
imaginary value, but with different signs. The Fortran standard
does not lock this down and the behavior can also be a function of
the machine architecture. The effect was to generate a divide by
zero with the Intel compiler. This was trapped with all the traps
enabled. Logic has been added to force the desired sign and fix
this inconsistency.
This error could occur with elastic layers and/or with other codes
in the Acoustics Toolbox. However, one would rarely use a lossless
halfspace in those models. Also, SCOOTER, for instance always
works a little off the real axis and wouldn't experience this.
This also might be an issue in the Matlab version of BELLHOP
depending on the architecture.
September 2020
An error was showing up in the Fortran BELLHOP under gfortran 9.2,
but not earlier versions of gfortran and not the Intel Fortran.
The error caused the ray trace to take a large number of small
steps. I found that the stepsize, h, had been declared as INTENT(
OUT ), when it should have been INTENT( INOUT ) since a maximum
step was being passed, that was supposed to be reduced by
ReduceStep. So the step size is both an input value and an output
value. Fixed. This error also was present in BELLHOP3D although
there had been no reports of problems. It has been fixed there
also.
October 2020
When the number of beams selected was 1, BELLHOP and BELLHOP3D
interpreted this as a full 360 degree sweep (since the difference
between the first and last beam MOD( 360 ) is 0. In such cases it
decrements the number of beams resulting in 0 beams. Using just
one declination angle would be a rare choice in BELLHOP; however,
in BELLHOP3D it can make sense to trace just one bearing angle for
an Nx2D run. In addition, a small change was made to ensure that
the vectors for these angles are always allocated to a minimum
length of 3. This is necessary because of the coding to
automatically generate linearly interpolated values, which is
flagged in the third element.
November 2020
A user was finding that BELLHOP3D was failing for a specific
azimuthal take-off angle for the beams. The error message was that
there was not enough storage for the ray. The cause of this was
that the ray was in a valley with the source precisely on the line
of the bathymetry that defined the minimum of the valley, and
launching along the same angle as that line. As a result, the ray
tracer sends the ray ping-ponging back and forth across that line
segment, always taking micro-steps since it tries to ensure that a
ray steps onto all interfaces.
BELLHOP and BELLHOP3D have now been modified so that the
'microsteps' are larger. The above described behavior still
occurs, but fewer steps are consumed in bouncing around the line
of the valley minimum. As a result, there is much less likelihood
of running out of storage. The same problem can occur if a ray is
launched along line defining a minimum in the ocean SSP. This can
happen in BELLHOP and BELLHOP3D. In cases where the number of
steps is a problem, this can be resolved by offsetting the source.
April 2021
We were trying to model an air-over-water case by putting an air
layer into the same medium as the water. (Recall that BELLHOP only
allows one 'medium'; in KRAKEN and SCOOTER we would handle this with
separate media.) Other people have done this successfully for ocean
over sediment problems. BELLHOP is not really designed for this and
it will not produce both transmitted and reflected rays at the
air/water interface. Instead, it gets a sort of reflection from the
gradient between air and water for rays below the critical angle; it
gets a sort of transmission between air and water for rays above the
critical angle.
If you do this, you generally want to spread that transition between
air and water over a meter or so. The reason is that the numerical
integrator is very inaccurate with sharp gradients in the SSP. You
can also address that by putting a lot of SSP points in the
transition layer. That forces BELLHOP to use a step size that lands
on each SSP point, i.e. the numerical integrator uses small steps
through the transition layer. So this is not a bug fix, but an alert
about properly handling SSPs with extremely strong gradients.
Regardless of how you sample the SSP, you must not have duplicate
depths. In other words there must be a transition layer of nonzero
thickness between air and water.
June 2021
When BELLHOP3D was run in the Nx2D mode for an arrivals calculation,
it was not writing the source and receiver azimuthal angles (which
are equal to each other in an Nx2D run) to the arrivals file This
causes problems in post-processing if the user is interpolating the
field between bearing lines.
July 2021
When BELLHOP3D was run in the Nx2D mode for an eigenray calculation,
it was writing the ray file in the 2D format used by the standard
BELLHOP. As a result, the rays could not be plotted using
plotray3d.m. They also could not be plotted by plotray.m used by
BELLHOP because the header for the ray file said it was a 3D ray
file.
When BELLHOP3D was run in the Nx2D mode for an arrivals calculation
it was omitting a necessary sqrt( r ) factor to account for
cylindrical spreading. As a result the Nx2D and 3D results were
typically very similar, but with different levels.
If a ray hits the boundary very close to a node were the bathymetry
(or altimetry) is sampled, that ray can escape (slip through the
boundary). This is a detail of the algorithm, but it always steps
just barely over the boundary and then applies the change in
direction due to the reflection. However, if that happens right near
a node then that ray hits other logic that forces it to land on
interfaces in range. Then with two adjacent ray points outside the
boundary, the algorithm kills the ray, deciding it has managed
to escape. This happened in at/tests/parabot for ray 118 (but
interestingly only with gfortran and without the runtime checks on).
It has been fixed by reducing the amount of the (small) overstep
that BELLHOP does on reflections. However, this is not a universal
fix, so it may (rarely) show up in other cases.
September 2021
In a harbor scenario that was 20 m x 20 m x 200 m we found that a
ray could step way past the bounding box. The particular case has a
50 Hz source leading to a 2-m default stepsize. The bounding box
only kills off a ray when it has stepped outside the box, which is a
fairly rough test. Both BELLHOP and BELLHOP3D have been modified to
treat the bounding box like other interfaces and make careful
adjustments to make sure the ray steps very precisely on the
bounding box (actually a very small fraction of the default stepsize
past that). This allows more precise treatment for confined cases
like this.
December 2021
A user created a BELLHOP3D bathymetry file with the
y-coordinates decreasing monotonically. The code assumes x
and y coordinates increase monotonically and fails if they
don't. A test has been added to flag this error, which otherwise
would cause the code to terminate mysteriously.
April 2023
A user found that with the Intel compiler an ASCII arrivals file was
wrapped at 80 characters causing the output info to be split between
two lines. This in turn affected their reader for the arrival file.
It turns out that the Fortran standard does not regulate this and
wrapping is not done by gfortran. A switch was added to the Makefile
-no-wrap-margin to change this. However, people should be aware that
this problem may occur with other compilers. It can also be fixed by
adding a logical record length in the OPEN statement for the file or
by using an explicit format statement instead of the list-directed
output. None of these solutions is ideal.
A user converted BELLHOP/BELLHOP3D to C++ and did meticulous checks
to make sure ray paths and other outputs were 'identical' with
different architectures and between the Fortran and the C++
versions. This was a very nice cross-check. Most of the differences
had no practical import or they addressed parts of the code that are
no longer used or supported. Trivial changes are not reporter here;
however, three changes were more significant bugs:
1) the ray normal calculation in BELLHOP3D was using the rotation
angle ray2%phi (which was uninitialized). It was supposed to use
ray1%phi,
2) the ray phase in the 2D Gaussian Cartesian beam was using the
(undefined) phase at step iS instead of the correct phase at step
iS-1 used for all the other beams,
3) the interpolation of the phase (phaseInt) in a zone that had just
crossed a caustic was using phaseInt = phase + pi / 2 instead of the
correct formula phaseInt = phaseInt + pi / 2.
These errors have been corrected and may have affected the
interference patterns in the coherent TL plots for different cases.
October 2024
The capability to treat a bottom type that varies with lat./long.
has been added. This is incorporated by adding additional lines to
the *.bty file that specify a matrix of bottom-type indices. This is
followed by a table that defines the geoacoustic bottom properties
associated with every index.
December 2024
As noted above on June 2019 a capability for elastic halfspaces was
added in BELLHOP. That same capability has now been added to
BELLHOP3D in both the Nx2D and 3D options. The test case
at/tests/Bellhop3d/provinces exercises this capability.
Previously range-dependent bottom types were hardwired to have
attenuation units in dB/wavelength. Now in BELLHOP and BELLHOP3D the
attenuation units are inherited from those selected in the envfil.
April 2025
Some incompatibilities (single vs. double) were introduced in the
binary arrivals format written by BELLHOP3D vs. what was expected in
the Matlab read routines. Fixed.
The plotsssp.m routine could misrepresent the SSP by under-sampling
it in the case of a spline fit. The fix was to simply increase the
sampling to that the spline shape was always visible.
August 2025
A 0/0 divide was causing a problem for an axial ray occurring
precisely at a cusp in the SSP (where the bi-linear profile had a
local minimum). This is a pathological situation where the ray
oscillates with an infinitesimal loop length. This affected arrivals
runs only, for both BELLHOP and BELLHOP3D. A test has been added to
ArrMod.f90 to bypass this situation.
November 2025
BELLHOP3D now allows bottom grain size to vary as a function of
lat./long. Previously as noted under October 2024 BELLHOP3D was
modified to allow the bottom geo-acoustic parameters to vary as a
function of lat./long. Also the grain size formula calculates the
sound speed in the bottom using a ratio based on the local sound
speed at the bottom of the water column. Previously the ratio was
based on a hardwired value of 1500 m/s relying on the user to update
the hardwired value by hand. This issue of the sound speed ratio has
also been fixed in BELLHOP.
A problem was fixed in the calculation of receiver depth indices
that received beam contributions in BELLHOP3D. This caused the first
received to always receive zero arrivals. It only affected the case
when ray-centered beams were selected.
December 2025
P-chip interpolation (similar to cubic spline, but fewer artifacts)
turned on in BELLHOP3D. This had originally been implemented in
BELLHOP but never implemented in BELLHOP3D. Implementing it in
BELLHOP3D turned out to be trivial. We are still gaining experience
with this option, but it is likely to be the best interpolation
choice in general cases. It has only been implemented in BELLHOP and
BELLHOP3D for stratified sound-speed profiles; it is not implemented
for range-dependent or 3D sound-speed profiles.
Fixed a slightly off curvature correction formula in BELLHOP3D.
After many surface interactions in a surface duct case it caused the
beams to become far too wide.
February 2026
In the December 2025 release I inadvertently introduced a bug in the
influence.f90 routine. This affected the Geometric Gaussian beams in
Cartesian coordinates. Since this routine is used by both BELLHOP
and BELLHOP3D it affected both. Fixed ...
June 2026
A check has been added to the section that reads in grain size to
ensure the user units are within the allowed bounds.
An error in BELLHOP3D was fixed in which the sound speed ratio and
attenuation were not being set.
On running the test battery I noticed a small error (a few dB) in
the near field for test case at/tests/Calib. This error occurred in
the Gaussian beam case in Cartesian coordinates. It was due to a
factor of lambda being dropped in the application of the stent,
which limits the minimum allowable beam size. This was fixed;
however, in checking also BELLHOP3D I note that stents are not
implemented at all in the 3D influence routine, which will degrade
the results. In the Nx2D case BELLHOP3D uses the same influence
routine as BELLHOP so this change corrects the error in both places.
August 2026
BELLHOP3D previously was not able to handle the case of a single
receiver using beams (hat or Gaussian) in ray-centered coordinates
(because it calculated the receiver index based on the spacing
between adjacent receivers). Logic has been added to handle this
case.
The various boundary reflection formulas in BELLHOP3D (both Nx2D and
3D) were computing the sound speed and gradient one step in advance
of the true step. However, one step in advance meant that we were
actually using a step from the previous beam since the ray
trajectory for the new beam had not yet been calculated. This
problem was particularly acute in a case with an extremely strong
near-surface gradient. Corrected ...
October 2026
I don't generally maintain the Matlab version of BELLHOP (called bellhopM) and it is not current with bug fixes in recent years. However, at a user's request I modifed bellhopM to output ray files in the standard format, rather than plotting them directly in bellhopM. The resulting ray files can be plotted in the usual way in Matlab using the plotray.m script.
Various bug fixes in BELLHOP3D: The beam influence routine for GeoHat beams in ray-centered coordinates (InfluenceGeoHatRayCen.f90) had a possible error in tracking caustic phase changes. The error occured for very coarse field sampling.
The ApplyContribution routine that adds in contributions to arrivals in an arrivals run was incorrectly over-writing a ray-tangent vector. This affected arrivals runs that used beams in Cartesian coordinates.
A 2 pi factor was missing in the ApplyContribution routine in the case of Gaussian beams in ray-centered coordinates. (It was included in the case of Gaussian beams in Cartesian coordinates.)