TRIMAIN-83: A Variable-Bottom, Multiprofile Raytace Program.
1. Introduction
TRIMAIN-83 is a FORTRAN-IV program in use on the VAX-11/780
computer in the Large Aperture Acoustics Branch at NRL. Its purpose
is to apply acoustic ray theory to an ocean environment in which both
sound speed structure and bottom depth vary with horizontal range.
The program can compute transmission loss by direct "eigenray" addition
or by probabilistic distribution of ray arrivals, and can be used to
compute intensity level, travel time, and source and receiver angles
for individual eigenrays.
The original program TRIMAIN was written in FORTRAN-63 for the
CDC-3800 computer [1] and documented in an NRL Report [2]. The purpose
of this note is to document the conversion to FORTRAN-IV, with some
modifications. A brief outline of the computational methods is given
here, along with instructions for preparation of input data, and the
input data and resulting output for a sample run. For further details
about individual subroutines, the reader may refer to the original
report.
2. Computations in TRIMAIN-83
Single-profile raytrace computations begin with sound speed data
in the form of a profile {c , z }; i.e., a table of sound speed vs depth
i i
at discrete depths from ocean surface to bottom. An interpolating
function is used to define c(z) for all depths; the form of the function
used and the method of fitting the data points define the raytracing
algorithm. For some interpolating functions the differential equation
for the ray trajectories can be solved in closed form [3]. In
particular, if the function
-1/2
c(z) = ( a + bz ) (1)
is used to define c(z) for z < z < z , solutions for the range,
i i+1
travel time, and spreading loss in the layer can be written out in
closed form [4]. Furthermore, if we start with a profile defined by
the range-dependent interpolating formula
-1/2
c(z,r) = ( a + b z + b r ) (2)
1 2
we can apply the rotation transformation
r = ( b r' + b z' )/b'
1 2
z = ( b z' - b r' )/b'
1 2
in order to transform the problem into the range-independent case
-1/2
c(z') = ( a + b'z' ) ,
2 2 1/2
where b' = ( b + b ) ,
1 2
so that the closed-form solutions mentioned above can be used.
The scheme used in TRIMAIN-83 is to divide the range-depth plane
into a network of right triangles, within each of which an interpolating
function of the form (2) applies. The bottom is modeled as piecewise-
linear, so that each bottom segment serves a a side of a triangle at the
deepest layer (See Fig. 1). With this interpolating scheme, the sound
speed is made continuous along the triangle edges, but its derivative
cannot be. The vertical triangle boundaries are located at the ranges
of input and interpolated sound speed profiles, an interpolated profile
being inserted at the range of each break point of the piecewise-linear
bottom.
The TRIMAIN ray trace starts with a set of rays at the source,
with initial angles specified by the user. The program advances the
whole set of rays, triangle by triangle, from each receiver or profile
range to the next. When the rays are at a receiver range, eigenrays are
found by linear interpolation, in the ray arrival set, to the receiver
depths. The eigenrays are described fully in terms of intensity level,
travel time, and source and receiver angles. Intensity at the receiver
may be computed by addition (either coherent or incoherent) of the
eigenrays. Another useful output is the ray history listing, which
is a tabulation, at each specified range, of the location, ray angle,
travel time, ray intensity level, and turnaround history of each of the
rays actually traced.
An alternative method of intensity calculation, suitable for
long ranges, is the 'depth-smoothing' computation, which distributes the
energy of each ray traced in a quasi-gaussian about its computed depth;
this approach is recommended for long-range problems not only because
it seems more physically plausible than the deterministic addition of
a discontinuous set of eigenrays, but also because it can be done with
a much sparser set of rays. A "range/depth-smoothing' method, using
convergence-zone averaging, may be used also if a 'mean' transmission
loss curve at long ranges is desired.
The recommendation of the depth-smoothed method for computing
transmission loss does not mean that eigenray calculations at long
range are useless; they can still be valuable for comparing travel
times and relative signal levels of multipaths. It should be borne in
mind that the density of the initial ray angle distribution required
to guarantee finding (almost) all of the eigenrays increases greatly
with receiver range.
3. Input Data
Note: The decimal place portion of the F-formats below
has been omitted, to make this section more readable.
Record 1. Any alphanumeric title in columns 1-72.
Record 2. Source and Environmental Parameters:
SD, FHZ, SORLEV, DATD, SLDB, XATT, XBOT
FORMAT(10F8)
SD Source depth (m).
FHZ Frequency (Hz).
SORLEV Source level (dB); default=0.
DATD Source beam pattern steering angle(deg).
(+ up, - down, default=0 = none).
Read pattern in as a table if non-zero.
SLDB Surface loss (dB); constant for all angles.
XATT Volume attenuation indicator:
=0 for no volume attenuation.
>0 for Thorp-type attenuation formula:
a(dB/km) = .10936 Y/(1+Y) + 50 Y/(5000 + Y),
where Y = square of frequency in KHz.
XBOT Bottom loss indicator:
<0 for perfectly absorbing bottom.
=0 for perfectly reflecting bottom.
1-9 for MGS or user-defined bottom loss class.
Record 3. Output Option Flags:
ICOH, IRAN, IT2, IT3, IGNRAY, IRD, IRF, IPL, NSRS, NBRS, LIMDB
FORMAT(11I5)
ICOH Transmission loss by coherent ray addition.
IRAN Transmission loss by random ray addition.
IT2 Transmission loss by depth smoothing.
IT3 Transmission loss by range/depth smoothing.
IGNRAY Eigenray output on print file TRIMAIN.IGN.
IRD Ray histories on print file TRIMAIN.DST.
IRF Ray histories on unformatted file TRIMAIN.RAY.
(This file is used by ray plot program.)
IPL Plot(s) of transmission loss vs. range.
NSRS = Surface reflections allowed before killing ray.
NBRS = Bottom reflections allowed before killing ray.
LIMDB= Loss in dB allowed before killing ray.
Record 4. Receiver Ranges and Depths:
R1KM, DRKM, R2KM, (RCD(I),I=1,NRCD) FORMAT(10F8)
R1KM, DRKM, R2KM are starting, increment, and final
ranges in km for equally-spaced receivers. If no
transmission loss calculation is requested, they are
still required for ray history or eigenray output.
RCD(I) are receiver depths in m. Up to 27 depths
may be entered; terminate input with a zero or negative
depth unless all 27 are used. These depths are used
only for transmission loss or eigenray calculations.
Record 5a. Ray Initialization (Fan) Parameters:
GAMDD, DGAMD, GAMUD, SL, PH, XCONT FORMAT(10F8)
GAMDD = Lower angle limit of this fan in degrees.
DGAMD = Angular increment of fan in degrees.
GAMUD = Upper angle limit of this fan in degrees.
SL = Source gain (dB) for this fan (default=0).
PH = Initial phase (radians) for this fan (default=0).
XCONT = Continuation flag; =0 if no more fans follow.
(Note that + angles are upgoing, - downgoing and
fans must be entered in order of increasing angle.
They need not be contiguous but may not overlap.)
Record 5b. Source Beam Pattern (if DATD is nonzero.)
(VBP(I),I=1,NBP) FORMAT(16F5)
(Read in the decrements (dB) for a symmetric beam
pattern, starting with the axial value, at 1-degree
spacing; last value read is extended to 90 degrees).
Record 6a. Bottom Loss Class Specification (if XBOT > 0).
RBLKM, ICLASS FORMAT(F8,I2)
RBLKM = Maximum range (km) for this bottom class.
(A non-positive range means to end of run).
ICLASS = 0 for perfectly reflecting bottom
1-5 are MGS classes, 6-9 are user-defined.
Note that this input is not used unless XBOT > 0.
Record 6b. Bottom Loss Table (if 5 < ICLASS < 10).
(RB(I),I=1,N) FORMAT(16F5)
RB = Values of bottom loss in dB, tabulated in 1-degree
steps in grazing angle, starting at 0 degrees. The
last non-zero value will be extended to 90 degrees.
Record 6c. Bottom Phase-Shift Table (for ICOH non-zero)
(RP(I),I=1,N) FORMAT(16F5)
RP = Values of bottom phase-shift in radians, tabulated
in the same way as RB. This input is used only with
user-defined bottom-loss classes, only when coherent-
phase transmission loss calculation is specified.
Each user-defined bottom-loss table (with phase-shift
table, if used) follows immediately the first Record
6a that specifies it. Once defined, The same class may
be specified for a subsequent range interval without
reading in the table again.
Record 7. Bathymetry Profile:
(RBKM(I),ZB(I),I=1,NRB) FORMAT(10F8)
RBKM and ZB are ranges in km, depths in m.
Record 8a. Sound speed Profile Title Card:
NCUR, RKM, TITLE FORMAT(I8, F8, 10A4)
NCUR = 0 means curved-earth correction is applied
to the profile [5].
NCUR > 0 means omit curved-earth correction.
NCUR < 0 means end of data (use after last profile).
RKM = Range of this profile in km.
TITLE = Any 40-character title.
Record 8b. Sound speed Profile:
(Z(I),C(I),I=1,NZ) FORMAT(10F8)
Z(I) and C(I) are depths in m, sound speeds in m/sec.
Notes: Depths must be in increasing order, starting
with 0. A maximum of 50 input points may be read in.
Terminate the profile with a zero or negative depth
(add an extra record if necessary), unless all 50
points are used. If an input profile does not extend
to the bottom, the program will extrapolate linearly
from the last two points; so it is important to be
sure that the input data goes deep enough. It
improves program efficiency to use input profiles
tabulated at common depths. In a single-profile
problem the profile need be entered only once.
REFERENCES:
1. E. L. Wright: "Ray-tracing with Horizontal and Vertical
Gradients", J. Acoust. Soc. Am. 48, 92(A) (1970).
2. B. G. Roberts Jr.: "Horizontal-Gradient Acoustical Ray-
Trace Program TRIMAIN", Nav. Res. Lab. Report 7827(1974).
3. M. A. Pederson: "Ray Theory Applied to a Wide Class of
Velocity Functions", J. Acoust. Soc. Am. 43,619-634(1968).
4. M. A. Pederson and D. F. Gordon: "Normal-Mode Theory
Applied to Short-Range Propagation in an Underwater
Acoustic Surface Duct", J. Acoust. Soc. Am. 37,105-118
(1965).
5. A. G. D. Watson: "The Effect of the Earth's Sphericity
on the Propagation of Sound in the Sea", Admiralty
Research Lab. Report ARL/N 29/L (1958).
Ocean Surface
------------------------------------------------------------------------
| . | . | . | . | . |
| . | . | . | . | . |
| . | . | . | . |. |
|---------------|-----------|-------------------|---------------|------|
~ ~ ~ ~ ~ ~
(About 30 layers in deep-water profile)
~ ~ ~ ~ ~ ~
|---------------|-----------|-------------------|---------------|------|
| . | | | . | . |
| . | . | . | . | |
| . | | | . | |
| . | . | . | . | . |
| . | | | . | |
| . | . | . | . | |
| . | | | . |. |
|---------------|-----------|-------------------.---------------|------|
| . | . | .......... | . |
| . | . | ................... | . |
| . | . | ............................ |. |
|---------------.................................................------|
| ................ ............................ |
| .................... BOTTOM .............................. |
| ........................ ................................ |
........................................................................
Figure 1
Partition of the range-depth plane into right triangles. The grid
depths are the set of sound speed profile depths and bottom profile
depths. The grid ranges are the set of sound speed profile ranges and
bottom profile ranges. (The bottom is piecewise linear.)
------------------------------------------------------------
Changes (22/11/88) MBP
I've found it necessary to make a few changes to the TRIM80 version
of TRIMAIN which came from NRL:
Changes in plotting routines to accomodate MINDIS
Changes in the unit numbers
Also, it's now set up with a separate command file and environmental file.
Thus, one creates a MUNK.ENV file and runs TRIMAIN via:
$ @TRIMAIN MUNK
Then you can obtain a ray plot via:
$ @TRIPLT MUNK