EIGENRAY: FORTRAN computer code for computing long-range, deep-ocean acoustic propagation by ray tracing

Ray Tracing for Ocean Acoustic Tomography Technical Memorandum APL-UW TM 3-98 December 1998 Applied Physics Laboratory University of Washington 1013 NE 40th Street Seattle, Washington 98105-6698 Approved for public release; distribution is unlimited. DARPA Grant MDA 972-93-1-003 ONR Grant N00014-97-1-0259 by Brian D. Dushaw and John A. Colosi _______________UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY_______________ ACKNOWLEDGEMENTS This work is conducted as part of the Acoustic Thermometry of Ocean Climate (ATOC) and North Pacific Acoustic Laboratory (NPAL) projects supported by DARPA (Grant MDA 972-93-1-0003) and ONR (Grant N00014-97-1-0259), respectively. John Colosi is grateful for a Young Investigator Award from the Office of Naval Research. Matthew Dzieciuch made several constructive comments on this technical report and on the development of the code. Bruce Cornuelle urged the inclusion of the effects of ocean currents in the ray calculations. Bob Odom provided most of the discussion of finite-frequency effects on rays passing near the ocean surface. Questions concerning this report and its associated FORTRAN code may be addressed to: Brian Dushaw Applied Physics Laboratory College of Ocean and Fisheries Sciences University of Washington 1013 N.E. 40th Street Seattle, WA 98105-6698 (206) 543-1300 dusha[email protected] John Colosi Applied Ocean Physics and Engineering Department Woods Hole Oceanographic Institution MS #11 Woods Hole, MA 02543 (508) 289-2317 [email protected] ii TM 3-98 _______________UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY_______________ ABSTRACT This report describes a new, flexible computer code in the FORTRAN computer language to make ray calculations for ocean acoustic tomography. The Numerical Recipes software package provided the basis for much of this computer code. The ray equations are reviewed, and ray equations that include the effects of ocean current are derived. Methods are derived for rapidly integrating those equations to obtain time front and eigenray information for long-range, deep-water acoustic transmissions. These methods include a look-up table for sound speed, sound speed gradient, second derivative of sound speed, and range-dependent information. Cubic spline methods are used to interpolate sound speed with depth and to obtain the derivatives of sound speed. The choice of the step size increments used to integrate the equations is a critical aspect of the integration, affecting both the accuracy of the prediction and the speed of computation. A predetermined, user-specified step size appears to allow more efficient calculations than "adaptive step" methods. "Adaptive step" methods adjust the step size automatically to maintain a given accuracy in the integration of the ray equations, while user-specified step sizes allow one to use prior knowledge of the integration problem to achieve the desired accuracy with much less computational overhead. Several integration methods were explored, but the classical 4th order Runge-Kutta method appears to be the most efficient and best method for this integration problem. Appendices describe detailed aspects of the computer code, as well as the methods used for deriving eigenray information and for parallelizing the ray calculations. The computer code is designed to be unstable so that the user can easily modify it to his or her own porpoises. TM 3-98 iii _______________UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY_______________ TABLE OF CONTENTS Page Motivation .................................................................................................................. 1 Ray Equations, or Equations of Motion..................................................................... 2 Look-Up Tables and Sound Speed Interpolation....................................................... 4 Integration of the Differential Equations................................................................... 6 Surface and Bottom Reflections .............................................................................. 8 Benchmarking vs Accuracy....................................................................................... 9 Conclusions ................................................................................................................ 10 Appendix A. A Technical Summary and a Flow Chart of the Computer Code........ 12 Appendix B. Calculation of Eigenrays ..................................................................... 14 Appendix C. Modifying the Code for Computations in Parallel .............................. 17 Appendix D. Input and Output Files and Other Operational Information................ 18 Appendix E. The Ray Equations in Terms of Sound Slowness ................................ 21 Appendix F. The Ray Equations with Current.......................................................... 22 References .................................................................................................................. 26 iv TM 3-98 _______________ _______________ LIST OF FIGURES Page Figure 1. Step sizes determined by an adaptive ray trace and the predetermined step size presently implemented in the code............................................ 28 Figure 2. Time front predictions associated with various step size scalings .......... 29 Figure 3. Flow chart of ray trace code.................................................................... 30 Figure 4. Ray path increments and sound speed..................................................... 31 TM 3-98 v UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY _______________UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY_______________ vi TM 3-98 This is a blank page. MOTIVATION Although numerous ray tracing codes are available, none satisfy all of the present requirements of long-range ocean acoustic tomography. Ocean acoustic tomography is described by Munk, Worcester, and Wunsch (1995). What is required for long-range tomography is a fast, accurate, and flexible code. The requirement of flexibility necessitates coding that can be easily modified by the user; thus FORTRAN is the preferred computer language. Flexibility allows the user to easily implement his or her own algorithms, such as a better search for eigenrays or an alternative (perhaps faster or more accurate) integration routine. Issues of numerical accuracy in ray predictions are discussed in the section on integrating the ray equations. Finally, the requirement of speed necessitates coding that can sometimes become opaque. This report therefore also describes the methods employed to achieve computer code that is highly efficient. The present and immediate goal of the authors of this report is to achieve fast, accurate wav efront and eigenray travel time predictions at basin scale ranges (3−5 Mm) in the North Pacific Ocean as part of the North Pacific Acoustic Laboratory and Acoustic Thermometry of Ocean Climate projects (Dushaw 1999; Dushaw et al. 1999; Worcester et al. 1998, Colosi et al. 1999). These projects require eigenray predictions for a half dozen or so time series of sound speed sections derived from acoustic data obtained at 3−5 Mm range. For a single sound speed section, calculation of the time series of ray travel times takes about 30 hours on a 200-MHz Pentium Pro computer, so the need for the fastest possible code is evident. Many of the ideas used in implementing the code to be discussed here originated in the RAY code and associated technical report by Bowlin et al. (1992). Alas, this code is in C. The new computer code described here uses the FORTRAN cubic spline and integration routines from Numerical Recipes (Press et al. 1992; Numerical Recipes hereinafter); thus use of this code requires a license to use the Numerical Recipes software ($40). CREDO: The code described here is meant to be easily modified by the user, and so it will never be a stable ray code. The code is meant to be transparent and fast. Metaphorically, if the Bowlin RAY code is a Fiat 2000 with a Fiat engine, the code described here is a '67 Chevy Impala with JAT O (Jet Assisted Take-Off) propulsion.* *Urban Legend No. 37 TM 3-98 1 _______________ _______________ UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY _______________UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY_______________ RAY EQUATIONS, OR EQUATIONS OF MOTION According to Bowlin et al. (1992): The equations of motion for a ray traveling through the ocean can be cast in Cartesian coordinates as follows: d θ dr =∂rc ctan θ −∂zc c(1a) dz dr =tan θ (1b) dt dr =sec θ c(1c) where θ is the angle of the ray with respect to the horizontal raxis, and zis the vertical coordinate [positive upward]. Bowlin et al. (1992) continue: For long range ocean acoustics, the curvature of the Earth’s surface makes non-Cartesian coordinates more suitable for ray tracing. Let new z’ axes lie along radii passing through the center of the Earth with z′=0at sea level and z′=R eat the Earth’s center, where Reis the radius of the Earth, and let the new r′be the range measured along a circular arc at sea level. ..... The new equations of motion which include the geometrical effects due to a spherical Earth are d θ dr′=fe∂z′c c−∂r′c ctan θ −1 Re(2a) dz dr′=fetan θ (2b) dt dr′=fesec θ c(2c) where fe=dr dr′=(Re−z′) Re(3) 2 TM 3-98 _______________UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY_______________ These are the equations that RAY integrates.... [These latter equations have z positive downward.] However, the RAY code actually integrates not θ , but cos θ and sin θ , thus avoiding the calculation of the transcendental functions. An alternate way to integrate these equations is to use the equations in Cartesian coordinates (the equations above with Re→∞and fe→1) but apply the well-known Earth flattening transformation (Aki and Richards 1980). This transformation is applied once to the initial sound speeds and associated depths, and the subsequent integration in the Cartesian coordinates is then mathematically equivalent to the above equations. Physically, this transformation is a stretching of depth and sound speed equivalent to the curvature of the Earth’s surface. If ε =z/Re(z positive downward) then the Earth flattening transformation is z=z*(1+ ε /2 + ε * ε /3) and c=c*(1+ ε + ε * ε ). Since some computation can be saved if this transformation is applied once before integrating the differential equations, this transformation together with the Cartesian differential equations is preferred. Ray predictions using the two methods agree to within 1 ms at 3-Mm range, and the integration of the Cartesian equations appears to be roughly 10% faster. If cs =cos θ ,sn =sin θ and we apply the flat Earth transformation to zand c, the equations above can be re-written as dcs dr′=−sn   ∂z′c c−∂r′c c sn cs (4a) dsn dr′=cs   ∂z′c c−∂r′c c sn cs (4b) dz dr′=sn cs (4c) dt dr′=1 ccs (4d) No trigonometric functions need to be calculated while integrating these equations. Since integration of these equations is most sensitive to the angle integration, the redundant equations for angle (the equations for cs and sn) are actually helpful for the stability of the integration. The above equations can be cast in only three equations, where the first equation is for d(tan θ )/dr, but this formulation involves a term sqrt(1+tan θ *tan θ )inthe integration for travel time (see Appendix B). This square root is computationally expensive (for one test case it increased computation time by about 15%), perhaps more expensive than integrating the additional differential equation. Similarly, it is important to code TM 3-98 3 _______________UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY_______________ (M. Dzieciuch personal communication 1999, Worcester et al. 1998). However, it is difficult to model the near-surface interaction with rays in a general way, particularly for broad-band acoustic transmissions. For the time being, those concerned about this effect may make an empirical assessment of its magnitude by comparing ray and mode or parabolic equation predictions. BENCHMARKING VS ACCURACY The speed of a ray prediction depends on the accuracy required by the user. Thus, it is virtually impossible to compare the speeds of various ray codes, unless one has a means of ensuring that the various predictions are to the same accuracy. The table below shows computation times for 5000 rays when using an early version of the ray code. The calculations are for a 3-Mm range, range-dependent sound speed section. For these calculations, the Step Size Scale is a scaling factor that is applied to the step size function shown in Figure 1. The table shows that the computation time depends linearly on the step size. Step Size Scale 0.5 1.0 2.0 4.0 Computation Time (min.) 32.0 16.5 8.75 4.5 Figure 2 shows that the time front predictions associated with these various step size scalings deteriorate considerably as step size is increased. By comparison, an adaptive step size prediction took about 1 hour, the Bowlin RAY prediction (Figure 2) took 34.0 minutes, and the Colosi prediction (Figure 2) (paranoid for accuracy) took over 5 hours. While for all of these predictions the location of the individual ray arrivals may vary, all of the time fronts agree within milliseconds. The shape and absolute travel times of time fronts appear to be fairly robust. The lesson here is that the user of a ray code needs to balance his or her desire for accuracy against the speed of computation. The astute reader will have noticed in Figure 2 occasional gaps in the time fronts, particularly in the earliest part of the arrival pattern. These gaps apparently result from the ray approximation; similar predictions using the parabolic equation do not show these gaps. In the ray approximation, rays which pass near a surface layer will either travel into that surface layer or miss it altogether, even though the launch angle in the two cases may differ only infinitesimally. This property of the ray approximation is a likely cause of the gaps in the time front prediction. When such gaps occur near the depth of a receiver, it can be difficult to obtain an associated eigenray prediction. Indeed, difficulty in obtaining the eigenray prediction for a ray known to exist probably results from such an unphysical gap in the time front. 10 TM 3-98 _______________UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY_______________ CONCLUSIONS We hav e described a method for making fast, accurate computations of acoustic rays as a tool for ocean acoustic tomography. We described ray equations, including a derivation of ray equations that include the effects of weak ocean currents (Appendix F). Methods using cubic spline interpolation and a look-up table allow sound speed and sound speed gradient to be calculated rapidly and accurately at arbitrary range and depth. The choice of the step size used in integrating the differential equations is critical, affecting both the computation time and the accuracy of the ray predictions. The best method for integrating the ray equations appears to be ordinary, classical 4th order Runge-Kutta integration. Presently known integration methods using the adaptive-step techniques to maintain a user-specified numerical accuracy appear to carry too much computational cost to be competitive. A user-specified step size greatly increases the efficiency of the computation, but a more efficient adaptive-stepping method may yet be devised. Integration accuracy can be initially checked for a particular problem by repeating the ray predictions with several trial step-size functions. For 5-Mm range acoustic transmissions across the North Pacific, eigenray travel times calculated with the code described here agreed with the travel times of the Bowlin code to within 20 ms. Similar agreement was found for other calculations at similar multimegameter range. The upper and lower turning depths of eigenray paths calculated using the code here and the Bowlin code agreed to within a few meters. Sound speeds calculated using the annual mean Levitus ’94 ocean atlas were used for these calculations. When calculating time fronts to similar accuracy, the present code appears to be comparable in speed to the Bowlin code. However, when calculating eigenrays (Appendix B), the present code appears to be 2−3 times faster than the Bowlin code because most of the time it integrates only two differential equations. In addition, the present eigenray code appears to be a little more effective at finding the eigenrays than the Bowlin code, so fewer rays need to be calculated to define the initial time front used to find the eigenrays. The code described here is inherently unstable, yet it is highly flexible. The computer code is easily modified for particular problems. One application of this code may be to calculate the forward problem matrices used for travel time inversion while calculating the eigenrays. This calculation would result in more accurate matrices and it will not require saving the ray paths to one’s hard disk. As described in some of the appendices, the code can calculate eigenrays, can be modified to run on a parallel computer for rapidly obtaining ray predictions, or can be modified to include the effects of ocean current. The present suite of software consists of about 1800 lines (including lots of commentary) of FORTRAN code. TM 3-98 11 _______________UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY_______________ APPENDIX A A TECHNICAL SUMMARY AND A FLOW CHART OF THE COMPUTER CODE A flow chart sketching the computer code is shown in Figure 3. The main program ray.f first sets up the sound speed look-up tables (subroutine speed.f) before integrating the rays (subroutine dodeint.f) with the desired ray launch angles. The code is in double precision throughout. Note that for the code to be universally portable the common tables must have the double precision arrays and variables listed before the integer variables. The sound speed look-up tables (ctab) are calculated using the Numerical Recipes cubic spline routines dspline.f and dsplint.f. The subroutine dspline.f calculates the second derivative of sound speed, from which the cubic splined sound speeds are obtained. The subroutine dsplint.f has been modified to obtain the subroutine dsplint_both.f which returns the values of both sound speed and sound speed gradient. While reading in the values for sound speed, the subroutine also obtains the table of profile ranges (range). A step size look-up table (steps) is also determined in the subroutine speed.f.As mentioned earlier, this predetermined step size is linear from the surface to a user-specified depth and then has a tanh functional form below that depth. The user is asked to specify the transition depth from linear to tanh forms and the step size values at the surface, at the transition depth, and at 5500-m depth. The bathymetry look-up table (btab) is determined in the subroutine bathy.f. This subroutine reads in the maximum allowable bottom bounces (MBONK), the tolerance for the ray to miss the ocean bottom (btol), and the bathymetry data. The bathymetry lookup table consists of the five variables: range, depth, slope, twice the slope angle, and a value of an intercept of the line segment. Once the sound speed, range, step size and bathymetry look-up tables have been defined, the rays are integrated using the subroutine dodeint.f. Presently the code loops over a set of rays at equally incremented launch angles. The subroutine dodeint.f has been heavily modified from the Numerical Recipes subroutine example of the same name. The role of this subroutine is as a driver for taking the Runge-Kutta steps (step size h) with the subroutine drk4.f; the user wishing to change the method of integration from Runge-Kutta may merely substitute an alternate subroutine to drk4.f. The subroutine drk4.f has been only slightly modified from the Numerical Recipes subroutine of the same name. The subroutine drk4.f calls the subroutine derivs.f which calculates the coupled differential equations (the ray equations) at arbitrary depth and ray angle. The subroutine dodeint.f saves the ray paths derived during the integration in variables xp and yp and returns the cosine and sine ray angle at the receiver range, the depth of the ray at the receiver range, and the ray travel time for writing to a file. Range dependence is implemented in the derivs.f subroutine using the variable deltaX;deltaX is the horizontal range between a step in the integration and the range of the sound speed profile most recently exceeded by the ray. Since range dependence is implemented using a constant sound speed gradient between profiles, the sound speed at 12 TM 3-98 _______________UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY_______________ any step in the integration can be found by c(r,z)=c(ri,z)+deltaX dc/dr. Here, c(ri,z) is the sound speed profile at range ri. The value for the vertical gradient of sound speed can be found similarly. In general, the step size h, O(5−200) m, is much less than deltaX, O(10 km). TM 3-98 13 _______________UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY_______________ APPENDIX B CALCULATION OF EIGENRAYS There are several ways in which eigenrays can be computed accurately, yet with far fewer calculations than are required for a complete time front prediction. In all cases, eigenrays must be found by first tracing a fan of rays and then determining the ray launch angles that result in rays arriving at the receiver depth. However, there are a number of shortcuts that can be implemented. In a ray fan equally spaced in launch angle, most of the rays concentrate in the cusps of the time front and so frequently do little to resolve the time front in the depth region of interest. Thus, one of the simplest ways to speed up eigenray calculations is to use an initial prediction with a small number of rays in order to define the range of ray angles that arrive near the depth of the receiver. Using these selected angle ranges, a second prediction with many rays can be done that is far more efficient at defining the time front near the depth of the reciever. A second way to gain efficiency is to omit the calculation of travel time in the initial fan of rays. All that is required to derive eigenrays is the angle that produces a ray that arrives at the receiver depth, so that travel time is not necessary. The present incarnation of the eigenray code implements this using only two coupled differential equations (for tangent of ray angle and ray depth; see Eq. B1 below), rather than the four required for obtaining travel time as well. Once the launch angles resulting in rays that arrive at the receiver depth are found, a small number of integrations with all four differential equations can then be performed to obtain all the information about the eigenrays. Eigenray predictions are probably more efficient if fewer, more accurate rays are used to calculate the time front rather than many, inaccurate rays. The more linear the relation between launch angle and arrival depth, the easier it will be to obtain the eigenrays. After the two ray arrivals that span the receiver depth are determined, eigenrays are probably best found by a sequence of bisections of ray launch angle giving a sequence of rays that converge on the receiver depth. An initial interpolation to solve for the eigenray launch angle might improve efficiency, but in general the ray arrivals are not linear in launch angle and arrival depth, so the interpolation will not be accurate. The code presently uses ordinary bisection on rays that span the receiver depth. Because of nonlinearities in the relation between launch angle and receiver depth, the sequence of bisections sometimes misses the receiver depth. Various contingencies have been built into the code to account for this, but it may sometimes happen that the eigenray cannot be found at all. It appears that at 3-Mm range, usually less than five bisections are required to find a ray that arrives within 10 m of the receiver depth. If one wants only the travel time of the eigenrays, one could interpolate the predicted time front to get the travel times associated with the depth of the receiver. Travel times determined in this way will (for benign oceans) be accurate to within a few milliseconds. 14 TM 3-98 _______________UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY_______________ Thus, the travel times of the eigenrays can be determined without having to actually solve for the eigenray path. This technique has not been implemented. The code here is flexible enough that the user can easily implement ideas for finding eigenrays. If someone develops code that works well for determining eigenrays, he or she is requested to forward those methods to the authors of this report. The code for finding eigenrays is similar to that described in Appendix A, with the addition of subroutines dodeint_short.f, drk4_short.f, and derivs_short.f. These subroutines use only two differential equations in tangent of ray angle and ray depth. If tn =tan( θ ), then these equations are dtn dr =(1+tn2)∂z′c c(B1a) dz dr =tn (B1b)    dt dr =√  1+tn2 c   (B1c) The subroutine used to find the eigenray ray launch angles is find_arr.f. This routine uses the time front information from the initial ray traces to select pairs of ray angles with depths at the receiver range that span the receiver depth. A series of bisections on ray launch angle is used to converge on the desired eigenray ray launch angle. The rays at the selected angles are then recalculated using all four differential equations to obtain all the information about the ray. Since this entire procedure relies on these two different integrations giving nearly identical ray paths, the procedure is actually one test of integration accuracy. If the eigenrays in the second integration miss the receiver depth, one problem may be that the step sizes are too big. However, the integrations are so sensitive that the series of step sizes during the two integrations must be fairly well synchronized. If even slightly different step sizes are used in the second integration, the rays can miss the the receiver depth by O(100 m). Once the eigenrays have been calculated, frequently the next step for the tomographer is to calculate the forward problem matrix, Gij = Γi ∫ η j(x)ds c(x)2(B2) TM 3-98 15 _______________UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY_______________ Here Γiis the ith ray path, η jis one of a set of functions used to model the ocean, ds is an element of ray path length (equal to dr/cos( θ )), and cis sound speed. Traditionally, to facilitate this calculation the parameter ds/c2is calculated for each saved step in the ray trace. With this parameter, the above integral can be calculated with relative ease. Howev er, ds/c2is not entirely easy to calculate with rigorous accuracy. The problem of calculating ds/c2is illustrated in Figure 4. While a ray path is calculated using very small horizontal step sizes, O(5−200 m), the ray path is saved at a much large range increment, O(1000 m). What is required, then, is a value for ds/c2 appropriate for this much larger horizontal range increment. It is not obvious how this parameter is to be calculated rigorously. It is easy, of course, to approximate ds by summing the small dsithat occur during the integration, but what of 1/c2? Should the values of c,1/c,or1/c 2be averaged over the several small steps? Sound speed varies so little over a path increment that it probably does not matter too much. The present code calculates dsi/c2at each integration step and sums these values over the several integration steps that comprise the saved range increment. The value of ds/c2is saved together with the ray path range and depth pairs. Possible errors may arise at depths where the sound speed or mode functions are varying rapidly, howev er. With the present set of subroutines, it would be relatively easy to merge the ray tracing routines with the calculation of the forward problem. In this way, the forward problem matrices can be calculated at the very small step sizes of the integration on the fly, as it were, without having to fill one’s hard disk up with the stored ray paths. This is the ideal way to calculate the forward problem matrices. The Bowlin et al. (1992) RAY code appears to address the problem of ds/c2by integrating a variable q=s/c2. With this variable, ds/c2may presumably be calculated at the nth step by qn−qn−1. This procedure is not correct, since d(s/c2)=ds/c2−2sdc/c 3. The second term has a dangerous dependence on the path length itself. The second term oscillates in sign because dc changes sign from upward-going to downward-going rays. Calculation of this variable requires the integration of an additional differential equation. For small values of s, the values of ds/c2derived from present code and from the Bowlin code appear agree to 1 part in 10,000 for the North Pacific environment modeled by the annual mean Levitus atlas. 16 TM 3-98 _______________UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY_______________ APPENDIX C MODIFYING THE CODE FOR COMPUTATIONS IN PARALLEL Calculating ray predictions for long-range transmissions is a natural application for parallel computation. Each ray trace is independent of all the others, so several computer processors can calculate subsets of rays separately. A parallel version of the ray code described here has been developed using the MPI suite of subroutines (Gropp et al. 1994). In order for several processors to be able to trace rays, all that is required is that they hav e access to the sound speed look-up table. Thus, the parallel version of the code broadcasts this look-up table, together with a few other variables, to all available processorss. With this information, the available processors can trace subsets of the desired rays. When they hav e completed these calculations, they then communicate the results (travel times, ray paths, etc.) back to the master process for printing or other use. The MPI routines perform these communications tasks (broadcasting to all processors or sending and receiving information between processors) relatively easily. The speed of the ray calculations scales well with the number of available processors. With a homogeneous set of processors, the complete set of rays to be calculated can be broken up into equally sized subsets. Each processor then gets a subset of rays to calculate. The present parallelized version of the ray code adopts this strategy. With an inhomogeneous set of processors, "load balancing" becomes an issue. Fewer rays should be sent to slower processors. It may be better to assign rays singly to processors and assign additional rays to a particular processor as its ray calculations are completed. MPI includes "non-blocking" send and receive communication subroutines which may be used for this purpose. That is, the code will not wait for a particular ray trace to be completed before assigning rays to other processors. In this way, a processor that is four times faster than another processor will complete four times as many rays without requiring the user to preallocate unequally sized subsets of rays to the processors. This strategy has not been implemented, but it would not be difficult to do so. TM 3-98 17 _______________UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY_______________ APPENDIX D INPUT AND OUTPUT FILES AND OTHER OPERATIONAL INFORMATION As noted often earlier, this code is meant to be modified by the user for his or her own purposes. This includes the formats of the sound speed data files and the means to input other information such as receiver depth, range, and other parameters of the ray trace as well as the output formats. The present setup of the code is described here as an example of a way to implement the code. The code is presently run by executing the command time eigenray < in.ray "time", of course, reports the time it takes to execute the raytrace; "eigenray" is the executable, and the file "in.ray" contains all of the relevant parameters that "eigenray" asks for. An example of a file for input to eigenray: 8. 12. Starting and ending values for ray angle (degrees) 100 1500. # rays to sketch out time front, and depth below which to omit rays in the second raytrace. The second ray trace is not performed for time front calculations. 100 # rays in second ray trace. 1000 1000 Source and receiver depths (meters) 0.0 3000000 Start and end ranges (meters) 1000.0 The range increment to save the ray paths (meters) 10.0 Tolerance for missing receiver depth (meters) 1 Sav e eigenray paths? 0=no, 1=yes 0 Calculate either eigenrays (=0) or a time front (=1) 1 Enable ocean bottom? 0=no, 1=yes test.ssp The sound speed filename 1500. Transistion depth for step size, linear to tanh (meters) 5. 100. 500. Values of step size at the surface, the transition depth, and at the bottom (all in meters) 1.0 Step size scaling e.g., 0.5,1.0,2.0; for testing accuracy 30 The maximum number of bottom bounces before dropping a ray 5.0 The tolerance for missing the ocean bottom (meters) test.bth The bathymetry file name. These last three lines are not required if the flag disabling bathymetry is set. 18 TM 3-98 _______________UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY_______________ The file of sound speeds has the same format as that of the Bowlin RAY code, with "-1 range" marking the start of a profile, followed by depth and sound speed. An example sound speed file: -1 0.000 0 1500.597 10 1500.465 20 1499.793 30 1498.795 50 1495.540 75 1491.698 100 1489.109 125 1487.366 150 1486.174 200 1484.450 250 1483.021 300 1481.752 400 1479.939 500 1479.103 600 1479.168 700 1479.679 800 1479.737 900 1480.225 1000 1480.835 1100 1481.586 1200 1482.276 1300 1483.210 1400 1484.137 1500 1485.137 1750 1487.678 2000 1491.187 2500 1498.126 3000 1506.161 3500 1514.534 4000 1523.252 4500 1532.259 5000 1541.404 5500 1550.529 -1 50.963 A new sound speed profile at 50.963 km range. 0 1500.751 10 1500.612 20 1499.982 30 1499.038 (etc., repeated as necessary) TM 3-98 19 _______________UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY_______________ REFERENCES Aki, K., and P. Richards, Quantitative Seismology, Theory and Methods, 2 vols., Freeman, San Francisco, 1980, 932 pp. Bowlin, J. B., J. L. Spiesberger, T. F. Duda, and L. F. Freitag, Ocean Acoustical RayTracing Software RAY, Tech. Rep. WHOI-93-10, Woods Hole Oceanographic Institution, Woods Hole, Mass., 1992, 49 pp. Colosi, J. A., E. K. Scheer, S. M. Flatte´, B. D. Cornuelle, M. A. Dzieciuch, W. Munk, P. F. Worcester, and A. B. Baggeroer, "Comparisons of measured and predicted acoustic fluctuations for a 3250-km propagation experiment in the eastern North Pacific," J. Acoust. Soc. Am., in press, 1999. Dushaw, B. D., Inversion of multimegameter range acoustic data for ocean temperature, IEEE J. Ocean. Eng., in press, 1999. Dushaw, B. D., B. M. Howe, J. A. Mercer, R. C. Spindel, and the ATOC Group, Multimegameter range acoustic data obtained by bottom-mounted hydrophone arrays for measurement of ocean temperature, IEEE J. Ocean. Eng., in press, 1999. Franchi, E. R., and M. J. Jacobson, Ray propagation in a channel with depth-variable sound speed and current, J. Acoust. Soc. Am., 52, 316−331, 1972. Gropp, W., E. Lusk, and A. Skjellum, Using MPI: Portable Parallel Programming with the Message-Passing Interface, The MIT Press, Cambridge, Massachusetts, 1994, 307 pp. Moler, C. B., and L. P. Solomon, Use of splines and numerical integration in geometrical acoustics, J. Acoust. Soc. Am., 48, 739−744, 1970. Munk, W., P. Worcester, and C. Wunsch, Ocean Acoustic Tomography, Cambridge University Press, New York, 1995, 433 pp. Murphy, E. L., and J. A. Davis, Modified ray theory for bounded media, J. Acoust. Soc. Am.,56, 1747−1760, 1974. Odom, R., Near Surface Turning Rays, Technical Memorandum, Applied Physics Laboratory, University of Washington, in preparation, 1998. Parker, R., Geophysical Inverse Theory, Princeton University Press, Princeton, New Jersey, 1994, 386 pp. 26 TM 3-98 _______________UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY_______________ Pierce, A. D., Acoustics: An Introduction to Its Physical Principles and Applications, Acoustical Society of America, Woodbury, New York, 1989, 678 pp. Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in FORTRAN: The Art of Scientific Computing, 2nd Ed., Cambridge Univ. Press., 1992, 963 pp. Worcester, P. W., B. D. Cornuelle, M. A. Dzieciuch, W. H. Munk, B. M. Howe, J. A. Mercer, R. C. Spindel, J. A. Colosi, K. Metzger, and T. G. Birdsall, A test of basinscale thermometry using a large-aperture vertical array at 3250-km range in the eastern North Pacific, J. Acoust. Soc. Am., in press, 1998. TM 3-98 27 0 500 1000 1500 3000 2500 2000 1500 1000 500 0 0 500 1000 1500 Depth (m) Step Size (m) Step Size (m) Figure 1. Step sizes determined by an adaptive ray trace and the predetermined step size presently implemented in the code. The left panel shows the step sizes as a function of depthderived by an adaptive step size algorithm used to calculate a single ray with a fairly small error tolerance. The right panel shows an analytic step size function that produces ray predictions as accurate as the adaptive step sizing but that require much less computation time. 28 TM 3-98 _______________UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY_______________ 3 2 1 0 Receiver Depth (km) 3 2 1 0 Receiver Depth (km) 3 2 1 0 Receiver Depth (km) 2021 2023 2025 3 2 1 0 Travel Time (s) Receiver Depth (km) S=4.0 S=2.0 S=1.0 S=0.5 Bowlin Colosi 2021 2023 2025 Travel Time (s) Figure 2. Time front predictions associated with various step size scalings. The left panels show the deterioration of the time front as the step size is increased. For very small step sizes the ray arrivals are evenly spaced. For comparison, predictions using the Bowlin and Colosi codes are shown in the right panels. In all cases the shape and travel time of the time front are nearly identical. TM 3-98 29 Program EIGENRAY Subroutine SPEED step size table - steps(#depths) (read in steps parameters) sound speed table - ctab(#depths,#profiles,6) (read in sound speed profiles) Subroutine SPLINE (2nd Deriv.) Subroutine SPLINT_BOTH (c(z), dc/dz) Information about source, receiver depth, range, etc. Subroutine ODEINT integrate over a path (save information at regular intervals if desired) Subroutine RK4 Runge-Kutta step Subroutine DERIVS Subroutine DERIVS (ray eqns.) Loop over ray launch angles Stop Write out ray information A multitude of steps One call to initialize Figure 3. A flow chart of one implementation of the ray trace code. This code is designed to be easily modified for the user?s purposes. 30 TM 3-98 _______________UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY_______________ Saved range increment hi h1h2h3h4hN-1 s1 s2 s3 s4 s5 sN ... ... ds1 ds2 ds3 ds4 dsN-1 ... c(z) or z) Depth Range i Figure 4. A schematic figure showing the relation of increments of ray path, s, to the sound speed profile or mode function. The ray path is calculated at small range increments hi (e.g., 5-200 m) but saved at a much larger increments (e.g., 1000 m). A value of ds/c2 associated with the saved path increment is required to calculate the matrices used in the inversion of travel times. It is probably better to calculate the matrices while performing a ray trace; better accuracy may be obtained by using the small step sizes. TM 3-98 31 _______________UNIVERSITY OF WASHINGTON • APPLIED PHYSICS LABORATORY_______________ REPORT DOCUMENTATION PAGE Form Approved OPM No. 0704-0188 Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the data needed, and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302, and to the Office of Information and Regulatory Affairs, Office of Management and Budget, Washington, DC 20503. 1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED 5. FUNDING NUMBERS 4. TITLE AND SUBTITLE 6. AUTHOR(S) 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION REPORT NUMBER 9. SPONSORING / MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSORING / MONITORING AGENCY REPORT NUMBER 11. SUPPLEMENTARY NOTES 12a. DISTRIBUTION / AVAILABILITY STATEMENT12a. DISTRIBUTION / AVAILABILITY STATEMENT 12b. DISTRIBUTION CODE 13. ABSTRACT (Maximum 200 words) 14. SUBJECT TERMS 17. SECURITY CLASSIFICATION OF REPORT 18. SECURITY CLASSIFICATION OF THIS PAGE 19. SECURITY CLASSIFICATION OF ABSTRACT 15. NUMBER OF PAGES 16. PRICE CODE 20. LIMITATION OF ABSTRACT NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89) Prescribed by ANSI Std. 239-18 299-01 This report describes a new, flexible computer code in the FORTRAN computer language to make ray calculations for ocean acoustic tomography. The Numerical Recipes software package provided the basis for much of this computer code. The ray equations are reviewed, and ray equations that include the effects of ocean current are derived. Methods are derived for rapidly integrating those equations to obtain time front and eigenray information for long-range, deep-water acoustic transmissions. These methods include a look-up table for sound speed, sound speed gradient, second derivative of sound speed, and range-dependent information. Cubic spline methods are used to interpolate sound speed with depth and to obtain the derivatives of sound speed. The choice of the step size increments used to integrate the equations is a critical aspect of the integration, affecting both the accuracy of the prediction and the speed of computation. A predetermined, user-specified step size appears to allow more efficient calculations than“adaptive step” methods. “Adaptive step” methods adjust the step size automatically to maintain a given accuracy in the integration of the ray equations, while user-specified step sizes allow one to use prior knowledge of the integration problem to achieve the desired accuracy with much less computational overhead. Several integration methods were explored, but the classical 4th order Runge-Kutta method appears to be the most efficient and best method for this integration problem. Appendices describe detailed aspects of the computer code, as well as the methods used for deriving eigenray information and for parallelizing the ray calculations. The computer code is designed to be unstable so that the user can easily modify it to his or her own purposes. 36 Unclassified Unclassified Unclassified SAR APL-UW TM 3-98 December 1998 Technical Approved for public release; distribution is unlimited. Ray tracing, ocean acoustics, geometrical optics Office of Naval Research Defence Advanced Research Projects Agency Ballston Tower 1 3701 N. Fairfax Drive 800 N. Quincy Street Arlington, VA 22203 Arlington, VA 22217 Applied Physics Laboratory University of Washington 1013 NE 40th Street Seattle, WA 98105-6698 Ray Tracing for Ocean Acoustic Tomography DARPA Grant MDA 972-93-1-003 ONR Grant N00014-97-1-0259 Brian D. Dushaw and John A. Colosi