	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

