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