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