***********************************
Program package ANRAY, version 4.01
***********************************

Dirk Gajewski *),  Ivan Psencik **)

*)  Inst.of Geophysics, University of Hamburg,
    Bundesstrasse 55, 20146 Hamburg, FRG

**) Geophysical Institute, Czechosl.Acad.Sci.,
    Bocni II, 141 31 Praha 4, Czechoslovakia

***************************************
Working version, not yet fully debugged
***************************************


Main differences from the version 3.02
--------------------------------------
There are several new options offerred in this version. The main
news are the following:

(a) The distribution of elastic parameters in individual layers can
be specified in grid points of a 3-D rectangular grid. Each of the
21 elastic parameters is specified in the grid individually. In case
of an isotropic layer, only specification of P and S wave velocities
is sufficient. If larger grid than specified in the present version is
to be used, changes in the dimensions in the routines of the block
MODBS.FOR must be made, see descritpion of input data for the B-spline
approximation.

(b) Graphical output in all the programs of the package is produced in
the form of postscript files. Calcomp-PostScript interface was developed by
L.Klimes.

(c) In order to increase the precision of the determination of points
of intersection of rays and dynamic ray tracing with interfaces and
boreholes, and thus to improve the convergence of the two-point ray
tracing, procedures proposed in [4] are adopted in this package.

(d) Tests, making possible a check of the precission of the ray tracing
and dynamic ray tracing, are offered. The tests can be used both
passively and actively, see the new switch IPREC in the description of
input data for the program ANRAY, sub 10. Further investigation of the
precission of the ray and dynamic ray tracing and its control is under
way.

(e) Many errors and problems have been removed and other small scale
modifications have been made. Most errors were detected and removed in
the routines for transformation of dynamic ray tracing across interfaces.
Most of them were spotted by P.Bakker (Shell Int., Rijswijk). I am very
grateful to Peter for it and I shall be grateful to anybody else for
indication of other problems since with the new options, new bugs were
surely introduced.


Description of the package
--------------------------
     Program package ANRAY can be used for computation of rays,
travel times, ray amplitudes and ray synthetic seismograms in 3-D
laterally varying structures containing isotropic and/or anisotropic
layers. Synthetic seismograms can be constructed at receivers
distributed regularly or irregularly along surface, interfaces or
vertical profiles.

    Program package consists of nine programs, basic program ANRAY
and programs ANRAYPL, SYNTAN, SEISPLOT, POLARPLOT, FRESAN, SYNFAN,
BPLOT and VELPL.

    Program ANRAY is an updated version of program with the same
name which has been used in packages ANRAY86 and ANRAY89 written by
Gajewski & Psencik [1], [2]. It is designed for ray, travel time and
ray amplitude computations. Two ways of approximation of distribution
of elastic parameters are available. In the first, "B-spline
approximation", elastic parameters in an arbitrary point of a layer
are determined by B-spline approximation from values of parameters
specified at grid points of a 3-D rectangular network. B-splines with
tension by A.K.Cline [3] with additions by L.Klimes, also used in
[4], are used. In the second way, "isosurface interpolation", elastic
parameters are determined by vertical linear interpolation of values
specified at interfaces, which represent surfaces of constant values
of elastic parameters. Rays can be computed in two modes. In the first
one, rays are specified by the point source location and the initial
orientation of the slowness vector at the source: initial-value ray
tracing. In the second mode, rays are specified by the point source
location and a system of regularly or irregularly distributed receivers
situated on surface or on an interface or on a vertical profile:
two-point ray tracing. The point source can be situated at any point
of the model. Polarization vectors, geometrical spreading and reflection,
transmission and conversion coefficients may be evaluated along the rays.
Program ANRAY can produce three files, one for plotting ray diagrams,
travel times and ray amplitudes, second for computation of synthetic
seismograms and third for plotting plane sections of slowness, phase
velocity or group velocity surfaces. Program ANRAY consists of 7 blocks
of routines, ANRAY.FOR, A2.FOR, A3.FOR, A4.FOR, A5.FOR,
CALCOPS.FOR and MODBS.FOR or MODIS.FOR, which must be linked together.
If B-spline approximation is to be used, MODBS.FOR must be used, if
isosurface interpolation is to be used, MODIS.FOR must be used.
Compile source code file anrayis.for
for the isosurface interpolation, and anraybs.for
for the B-spline interpolation.
Refer to the description of program ANRAY.

    The program ANRAYPL can be used for plotting of horizontal and
vertical ray diagrams, time-distance and amplitude-distance curves
of individual elementary waves computed in the program ANRAY.
Several types of point sources, namely explosive source, single
force, double couple source can be considered. Program ANRAYPL is
a modification of the program RAYPLOT from the package SEIS83 [5].
Program ANRAYPL consists of routines: ANRAYPL.FOR, BORDER.FOR,
SOURCE.FOR and CALCOPS.FOR.
Refer to the description of program ANRAYPL.

    The program SYNTAN can be used for computation of ray synthetic
seismograms from the results computed in the program ANRAY.
Gabor wavelet is used as a source-time function. Hilbert transform
of this wavelet is evaluated by approximate formula, see [6].
Several types of point sources, namely explosive source, single
force, double couple source, can be considered. Program SYNTAN is
a modification of the program SYNTPL from the program package SEIS83
[5]. Program SYNTAN consists of routines: SYNTAN.FOR,
SOURCE.FOR and CALCOPS.FOR.
Refer to the description of program SYNTAN.

    Programs FRESAN and SYNFAN can be used for frequency domain
computation of ray synthetic seismograms. Frequency response is
computed in the program FRESAN. In the program SYNFAN it is multiplied
by the spectrum of the considered source-time function. Inverse
Fourier transform is then used to evaluate synthetic seismograms.
Programs FRESAN and SYNTAN are modifications of the programs GB
and SYNTGB from the package BEAM87 by V.Cerveny. Program FRESAN
consists of routines FRESAN.FOR, BORDER.FOR, SOURCE.FOR,
and CALCOPS.FOR. Program SYNFAN consists of routines SYNFAN.FOR
and CALCOPS.FOR.
Refer to the description of program FRESAN.
Refer to the description of program SYNFAN.

    Programs SEISPLOT, BPLOT and POLARPLOT can be used for plotting
results computed in progams SYNTAN and SYNFAN. Program SEISPLOT can be
used for plotting synthetic seismograms generated by the program SYNTAN,
program BPLOT can be used for plotting synthetic seismograms generated
by the program SYNFAN. The program POLARPLOT serves for plotting particle
motion diagrams and can be used in both cases. Programs SEISPLOT and
POLARPLOT are modifications of the programs with the same name from
the package SEIS83, program BPLOT is a modification of the program
with the same name from the package BEAM87. Program SEISPLOT consists
of routines SEISPLOT.FOR, BORDER.FOR and CALCOPS.FOR.
Program BPLOT consists of routines BPLOT.FOR, BORDER.FOR and
CALCOPS.FOR. Program POLARPLOT consists of routines POLARPLOT.FOR,
BORDER.FOR and CALCOPS.FOR.
Refer to the description of program SEISPLOT.
Refer to the description of program BPLOT.
Refer to the description of program POLARPLOT.

    Program VELPLOT can be used for plotting plane sections of
slowness, phase velocity and group velocity surfaces from the file
generated in the program ANRAY. Program VELPLOT consists of two
routines: VELPL.FOR and CALCOPS.FOR.
Refer to the Description of program VELPLOT.

    All files with the main programs (anray.for, anraypl.for, syntan.for,
fresan.for, synfan.for, seispl.for, bplot.for, polar.for and velpl.for)
contain, at their ends, Fortran 90 INCLUDE commands to include all
subroutine files required.  It is thus no need to compile subroutine
files separately and to submit them to a linker.

    FORTRAN77 is used throughout the program package ANRAY.
Standard Calcomp routines, namely PLOTS,
PLOT, SYMBOL and NUMBER, are used in plotting programs. Calcomp-PostScript
interface containing the above listed routines and developed by L.Klimes,
is a part of the package. This makes possible to generate the graphical
output in the form of ASCII files coded in the PostScript level 2 language
- encapsulated postscript file format version 3.0. The Calcomp-PostScript
interface is a part of the block of routines CALCOPS.FOR. This block
must be linked with each program of the package.

    The package is supplemented by a complete set of test input data for
an "anisotropic" model, in which isosurface approximation is used. At
present version the test input data for a model, in which B-splines are
used, are given only for the program ANRAY. A simple three-layered model
consisting of two homogeneous isotropic layers bounding from above and
below an anisotropic layer with horizontal variation of elastic parameters
in x- and y-direction is considered. At the receivers situated in a
borehole both in isotropic and anisotropic layer all direct waves and
one reflected qS wave are calculated.

    The program package is NOT yet fully debugged. Problems can be
encountered especially in singular regions in which two slowness surfaces
of the shear waves are close to each other. The ray method gives in such
regions distorted results and the two-point ray tracing procedure can fail
there. Introduction of criteria checking precission of the ray tracing and
dynamic ray tracing indicates that present formulation of these procedures
is very sensitive to the above mentioned singular regions. Alternative
formulations have to be, therefore, studied. Search for alternative ways
of solving the polynomial equation of the sixth order, necessary for the
determination of slowness vectors of waves generated by reflection or
transmission at an interface, is also desirable. Any report on other
errors, problems or inconsistencies, which may be encountered by users,
is welcome.

    In case of publication of results obtained with the program
package ANRAY, the users are asked kindly to refer to this report.

Good luck!


Praha, September 1997


REFERENCES
                                                        
[1] Gajewski,D., Psencik,I., 1986. Numerical modelling of seismic
    wave fields in 3-D laterally varying layered anisotropic
    structures - Program ANRAY86, Internal Report Inst.of Earth and
    Planet.Phys., University of Alberta, Edmonton.
                                                        
[2] Gajewski,D.,  Psencik,I.,1989. Ray  synthetic seismograms  in 3-D
    laterally inhomogeneous anisotropic structures - Program ANRAY89,
    Internal  Report   Centre  for  Computational   Seismology,  LBL,
    Berkeley.
                                                        
[3] Cline,A.K.,1981. FITPACK, Dept. of Computer Sciences, Univ. of
    Texas at Austin.
                                                        
[4] Cerveny,V., Klimes,L., Psencik,I.,1988. Complete seismic ray
    tracing in three-dimensional structures. In: Seismological
    Algorithms, D.J.Doornbos edit., Academic Press, New York, 89-168.
                                                        
[5] Cerveny,V., Psencik,I.,1984. Documentation of earthquake
    algorithms. SEIS83 - Numerical modeling of seismic wave fields
    in 2-D laterally varying layered structures by the ray method.
    E.R.Engdahl, edit., Report SE-35, Boulder, 36-40.
                                                        
[6] Cerveny,V., Molotkov,I.A., Psencik,I.,1977. Ray method in seismology,
    Univerzita Karlova, Praha.