February 2008 - Brian Dushaw Version 3.0 of the eigenray code has had a major revision to its cubic splining routines. The Numerical Recipes routines have been abandoned in favor of routines from netlib.org. The instigating reason for the change was that the netlib.org routine implemented the "not-a-knot" boundary condition, which seems more appropriate for sound speed profiles. (Matlab's spline uses the same, by default.) The effect, I believe, is a slightly better spline over the lowest several points of the sound speed profile. The input files should all be unchanged from the previous code. As I looked at the new routines, it became apparent that there was a much better way of going about meshing the tabulated spline to the Runge-Kutta integration using the look-up table. I've implemented that, and the calculations work much better now. Speeds significantly faster than the previous version, and results that are better behaved (depending on your value for epsilon, of course). The computation time to complete the test case went from 1 min. using the previous code to 20 s using the new code, with results that looked far more accurate! The new code is elegant and slick and offers a number of advantages: (a) No Numerical Recipes routines (sort of; at least the splining is not used.) (b) Choose the depth interval for interpolation in the look-up tables, as you like. No longer locked to 3-m intervals. (c) MUCH cleaner source code. Previous versions of the code used a look-up table that was locked-in at 3-m sampling (or at least it was not straightforward to change it). Now if you want to change this increment, just edit the appropriate values in include/table.h (Ndepmax and depth_inc; Ndepmax points separated by depth_inc meters). 1-m intervals and 5500 m depth seems to use up 250MB of RAM in the test case. The code is now inherently more accurate, so that the values of epsilon can be relaxed somewhat from what was appropriate for the older code. Use epsilon=1e-4, say, for some quick answers, and 1e-7 for more accuracy. The main drawback with smaller epsilon is longer computation time; too small and the code may just thrash around trying to get accuracy the the Runge-Kutta-Merson integration can't achieve. (If you really want accuracy, you might also reduce the step sizes (10. 400. 700.) in the input file; that will save some of the thrashing.) *** Standard recommended values for both speed and accuracy: epsilon=1e-5 and step sizes (10 400 700). You may have to experiment with these numbers to get good results for your own situation. The integration seems to be inherently fairly accurate, however, so that excessively small values for the parameters may needlessly use up the computer cycles. The routine src/speed.f has been completely rewritten. There is an option "iprintflag" you can set at the top that will write out the sound speed profile ("c_out.dat") as interpolated by the code. The option turns off the flat-earth correction as well. The mex file src/meig.c compiled and ran on my machine (Suse linux 10.3, intel ifort 10.1, matlab 2007b). I'll leave the binary here - perhaps you can run it; this is really tempermental, however. (a) execute make to make the object files. (b) check the mexopts.sh for your own situation. (c) execute mex src/meig.c to make the mex file. (d) run the test case "runeig.m" in matlab (and watch matlab crash...) Incidentally, the file "ray.info.accurate" stores the results of an highly accurate run of the test case (for the older version of the code) which can be used to test the effect of epsilon, etc. END