Program ANRAY - Version 4.01
****************************

   Program ANRAY is a modified version  of packages ANRAY86 and
ANRAY89 by Gajewski and Psencik [1],[2]. It is designed for ray,
travel time, ray amplitude, Green's function and ray synthetic
seismograms computations. 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.
This mode is called 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 along profiles
on surface or an interface, or along vertical line profiles. This
mode is called two-point ray tracing. Polarization vectors,
geometrical spreading and reflection, transmission and conversion
coefficients may be evaluated along the rays. Basic program of the
package ANRAY also called ANRAY can produce three files: the first for
plotting ray diagrams, travel times and ray amplitudes, the second for
computation of synthetic seismograms and the third for plotting plane
sections of slowness, phase velocity or group velocity surfaces.

   Since the program is based on the ray method it does not work
properly in the so-called singular regions like vicinity of caustics,
of critical points, of transitions from shadow to illuminated region,
in anisotropic media in directions in which the two qS waves propagate
with nearly the same phase velocity. The existence of the latter
singularities does not affect only the amplitude computations but
may also cause problems with convergence of iterations to two-point
rays.

Most important characteristics of the package ANRAY:

  - the model is 3-D, laterally inhomogeneous, with curved layers
    of nonzero thickness;
  - elastic parameters inside layers can be determined either by
    a vertical linear interpolation between interfaces, on which
    constant values of elastic parameters are specified or by a
    3-D interpolation by B-splines with tension from the values
    of elastic parameters specified in grid points of a 3-D grid;
  - a point source can be situated anywhere in an isotropic or
    anisotropic layer;
  - receivers can be distributed not only along surface or interface
    profiles, but also along vertical profiles;
  - surface and interface profiles do not need to pass through the
    vertical projection of the point source on the surface or the
    corresponding interface;
  - 3-D two-point ray tracing allowing to search for rays arriving
    from a source to a prescribed vicinity of a receiver specified
    on a surface, interface or vertical profile, based on paraxial
    ray approximation;
  - the output of the program ANRAY can be processed by modified
    programs from the modified 2-D package SEIS83 [3] for time-domain
    or frequency-domain computations of ray synthetics.

Some useful information

  - when using ANRAY in the two-point ray-tracing mode, it is
    desirable that the box containing the model is constructed
    so that the source is always situated inside the model, away
    from the vertical boundaries of the model; this makes possible
    an effective use of the shooting procedure. It is also
    recommended to avoid receivers situated in the corners of
    the model box;
  - the bottom boundary of the model cannot be used as a reflector
    when ray amplitudes are to be evaluated; calculation of
    reflection coefficients requires knowledge of parameters
    below the reflector.

Coordinate systems used in the program
--------------------------------------

    Several coordinate systems are used in the program. The most
important ones from the user's point of view are the following:

A) "Model" Cartesian coordinate system

    This coordinate system is used for the description of the model,
i.e., for the delineation of interfaces and the elastic parameter
distribution. It is defined as follows. The x and y coordinate
axes are situated in a horizontal plane, with the vertical (z) axis
chosen positive downwards. The positive x and y axes are chosen such
that the coordinate system is a right handed one.

B) "Receiver" cylindrical coordinate system

    In the two-point ray-tracing mode, i.e. in the mode, in which
rays arriving to a priori specified receivers are sought, receivers
must be specified. Receivers must be always situated on profiles.
These profiles are always situated in a vertical plane. In case of
receivers situated along a surface profile, receivers are located
along the intersection of the vertical plane with the surface. The
same holds for a profile of receivers distributed along an interface.
In case of receivers distributed along a vertical profile, the
profile is straight vertical line specified by its x- and
y-coordinates. The coordinates of receivers are specified in a
cylindrical coordinate system with its vertical axis passing through
a point specified by the user. In this way, receivers along a surface
or interface profile are fully specified by the azimuth of the profile
and distances of receivers from the vertical axis. In case of a vertical
profile, receivers are fully specified by x- and y-coordinates of the
vertical profile and by "model" z-coordinates of receivers.

C) "Elastic tensor" Cartesian coordinate system

    A convenient and frequently used coordinate system for the
description of the anisotropic properties of a material is the system
whose axes coincide with the principal axes of the elasticity tensor.
Such a system is very practical since in it many of the 21 elastic
parameters are zero. In this program the data related to the elastic
parameters may be specified in such a coordinate system. Additional
information is necessary specifying how the "elasticity tensor"
coordinate system is to be rotated to coincide with the "model"
coordinate system, see description of input data for isosurface
interpolation sub 7a.

Description of the model and the source
---------------------------------------

    The model is situated in a "box" bounded by plane vertical boundaries
parallel to the (x,z) and (y,z) planes, and by, possibly curved, top
and bottom boundaries formed by interfaces. The model consists of a
system of layers separated by curved non-intersecting interfaces.
The layers as well as interfaces are numbered sequentially from the
top to the bottom. All interfaces must be defined on the whole (x,y)
cross-section of the model.

    The interfaces are specified by points on a rectangular grid
completely covering the (x,y) model cross-section. Thus, the first
and last lines of this grid must be situated on vertical boundaries
of the model. Bicubic spline interpolation is used to approximate
the interfaces between the grid points. The interfaces are therefore
smooth without corner points.

    The medium inside layers may be either isotropic or anisotropic.
All 21 elastic parameters and density may vary spatially within the
anisotropic layers and the P and S wave velocities and density may
vary spatially in the isotropic layers. The elastic parameters are
taken to have the dimensions of velocity squared, which corresponds
to the elastic parameters normally associated with an anisotropic
medium divided by the density of the medium. The P and S wave
velocities in an isotropic layer are, for simplicity, also referred
to as elastic parameters. Two methods of the determination of the
elastic parameters in individual layers may be used.

    (a) Linear interpolation of elastic parameters along vertical lines
between interfaces representing isosurfaces of elastic parameters.
In this approach, the distribution of elastic parameters is laterally
varying within the layers if their boundaries are curved or inclined.
The vertical gradient of the elastic parameters within a layer is
large when the interfaces bounding the layer are close to one another
and vice versa.

    (b) B-spline interpolation of elastic parameters from their values
specified in grid points of a 3-D grid. B-splines with tension by
A.K.Cline [4].

    The length and velocity units km and km/s, m and m/s or m and m/ms
must be used consistently. The density should be always specified in
g/cm**3 (or 1000 kg/m**3).

    A point source may be situated at any point of the model (for
the effective functioning of the two-point ray-tracing procedure
it is recommended that it is situated inside the model box, away
from its vertical boundaries). If it is situated at an interface,
it is automatically considered to be situated in the underlying
layer. A unit isotropic radiation function is incorporated in the
program ANRAY. Radiation functions of single force, double couple
and explosive (implosive) point sources may be introduced during
further processing of the data obtained from the program ANRAY in
the programs ANRAYPL, SYNTAN or FRESAN.

    The polarization vectors at the point source situated in an
isotropic layer are specified in the following way. The initial
polarization vector of a P wave is specified by a unit vector tangent
to the ray at the source. The initial polarization vectors of an S
wave are specified by unit vectors situated in the plane perpendicular
to the ray. The first of these vectors is situated along the
intersection of a vertical plane with the plane perpendicular to
the ray. It is orientated so that it has always a negative component
into the z axis in the "model" coordinates The other unit vector is
chosen such that the triplet formed by the first and second unit
vectors in the plane perpendicular to the ray and the vector tangent
to the ray at the source form a right-handed system. For this choice,
the second unit vector in the plane perpendicular to the ray is
horizontal.

Description of elementary waves and their codes
-----------------------------------------------

    Various refracted, multiply reflected and converted elementary
waves may be considered. Each of these waves is described by a
unique code which can be manually generated.

    The codes are defined as follows. The whole ray is divided
into segments, each of which lies between two successive points
at which the ray strikes an interface. If the endpoints of a
segment lie on different interfaces, the segmnent is called a simple
segment. If the end points lie on the same interface, the segment
is called compound segment. Any compound segment is regarded as
to be formed of two simple segments. The part of the ray between
the source and the point at which the ray first strikes an interface
is considered as the first segment of the ray. If the source is
situated on an interface and the ray strikes first this interface
from below, the first segment is regarded as a compound segment.

    The code of a wave consists of chain of doublets corresponding
to individual simple segments of its rays, successively, from
the source towards the termination points. The first number in the
doublet gives the number of the layer, in which the corresponding
segment is situated. The second number specifies the type of the
wave along the segment. The numbers 1 and 2 correspond to quasishear
waves (in an isotropic layer it is sufficient to specify the number 1
for the shear wave, the number 2 gives the same wave) the number 3
corresponds to the quasicompressional wave (or compressional wave
in an isotropic layer).

Description of the ray computations
-----------------------------------

    Individual rays are specified by the source position and by two
angles specifying the orientation of the slowness vector at the
source. For simplicity, one of the angles is called here azimuth,
and the other declination. Both angles are specified in "model"
coordinate system. Let us note that in anisotropic layers, the above
angles generally differ from the take-off angles of rays. In
isotropic layers, the angles specifying the slowness vector coincide
with the take-off angles of corresponding rays.

    The azimuth is measured in the (x,y) plane, positively clockwise
from the positive x axis. The declination is defined as an angle
between the slowness vector and its projection into the horizontal
plane. It is measured positively clockwise, zero declination
corresponding to a ray leaving the source with horizontal slowness
vector.

    The program ANRAY may be run in two different modes:
(a) Initial-value ray-tracing mode,
(b) Two-point ray-tracing mode.

    In the initial-value ray-tracing mode, the rays are computed for
a system of specified initial angles, without a priori knowledge of
where they will terminate.

    In the two-point ray-tracing mode, the rays which terminate in a
prescribed vicinity of receivers regularly or irregularly distributed
along a surface, interface or vertical profile are computed.
Two-point ray tracing is based on the shooting method. Modified
versions of the routine TRAMP from the program package SEIS83 [2] are
used for this purpose.

    In both modes, a system of initial angles must be specified,
both azimuths and declinations. In the initial-value ray tracing
mode, ray paths and related quantities corresponding directly to
these initial angles are computed. In the two-point ray tracing mode,
the system of initial angles serves only as a basic system, within
which an iterative search for rays terminating at the specified
receivers is performed.

    The procedure for finding rays terminating in a vicinity of
prescribed receivers works as follows. First, search is made within
an iterative procedure (routine PROFIL) for a ray terminating on the
profile containing the receiver. During this search declination is
held constant. If a ray terminating on the profile is found, it
enters into the second iterative procedure (routine RECEIV), which
deals only with rays terminating on the profile. The second iterative
procedure searches for a ray terminating at the receiver. The
iterative procedure is based on the paraxial ray approximation.

                                                     
Description of input and output files
-------------------------------------

    Main input data are read from the standard input by list-directed
input (free format) and consist of a single line containing following
data:
    'LIN' 'LOU' 'LU1' 'LU2' 'LU3'/
Here:
    'LIN' is the name of the input data file LIN.
    'LOU' is the name of the output log file LOU.
    'LU1' is the name of the output data file LU1.
    'LU2' is the name of the output data file LU2.
    'LU3' is the name of the output data file LU3.
    / is a slash recomended in batch and script files to enable future
        extensions.
Defaults:
    'LIN'='anray.dat'
    'LOU'='anray.out'
    'LU1'='lu1.out'
    'LU2'='lu2.out'
    'LU3'='lu3.out'
Examples of the main input data for anrayis.for:
    'anrayis.sch' /
Examples of the main input data for anraybs.for:
    'anraybs.grd' /

    The data for the program ANRAY are read in from the file LIN,
specified by the user. Output data describing the computations are
stored in the file LOU. Output data for plotting ray diagrams, travel
time curves and amplitude-distance curves in the program ANRAYPL are
stored in the file LU1. Output data for computing ray synthetic
seismograms in the program SYNTAN are stored in the file LU2.
Output data for plotting sections of velocity surfaces are stored
in the file LU3.

                                                      
Input data in the file LIN
--------------------------

1) MTEXT                                          FORMAT(20A4)

   MTEXT...    description of the problem under consideration,
         (alphanumeric text).

2) INULL,ISURF                                    FORMAT(17I5)

   INULL...    in some tests (e.g. on zero elements of matrices),
         the numbers less than or equal RNULL=10**(-INULL)
         are taken as zero. Default value (INULL=4).
   ISURF...    ISURF=0: no calculation of sections of velocity
               surfaces.
               ISURF=1: calculation of sections of velocity
               surfaces.

3) MPRINT,NINT                                    FORMAT(14I5)

   MPRINT...  controls storage of the description of the model
         in the file LOU.
         MPRINT=0: Only a copy of input data, no other description
         of the model is stored.
         MPRINT=1: Writes the most important input data plus printer
         plots of the form of interfaces in the model (see ZMIN,
         ZMAX in line 5).
         MPRINT=2: The same as for MPRINT=1 plus description of the
         unrotated and rotated compressed matrices of the elastic
         parameters.

   NINT...  number of interfaces in the model (including first
         interface representing the model surface and the last
         interface representing the model bottom).

4) The lines 4a-4d repeat NINT times, each for one interface, from
   the top to the bottom of the model.

4a)MX,MY                                          FORMAT(14I5)

   MX,MY... number of grid lines x=const and y=const for
         approximation of the interfaces.

4b)SX(1),...,SX(MX)                               FORMAT(8F10.5)

   SX(I)...  x-coordinates of grid lines for approximation of
         the interfaces.

4c)SY(1),...SY(MY)                                FORMAT(8F10.5)

   SY(I)...  y-coordinates of grid lines for approximation of
         the interfaces.

4d)Z(1),...,Z(MX*MY)                              FORMAT(8F10.5)

   Z(I)...  z-coordinates (depths) of the interface at the grid
         points, successively for x=SX(1) from y=SY(1) to SY(MY),
         then on a new line for x=SX(2) from y=SY(1) to SY(MY),
         a.s.o.

5) This line repeats NINT times, each one for one interface, from
   the top to the bottom of the model.

   ZMIN,ZMAX                                      FORMAT(8F10.5)

   ZMIN,ZMAX... the depth range at which the interfaces are displayed
         in the printer plots. The printer plots consist of systems
         of digits from 0 to 9 and possibly asterisks covering the
         whole horizontal extent of the model with 101 points in the
         x direction and 51 points in the y direction. The digits are
         determined from the computed z coordinates of interfaces on
         the 101x51 grid using the relation
                  I=IFIX(10.*(Z-ZMIN)/(ZMAX-ZMIN).
         Outside the range (ZMIN,ZMAX), the asterisks are printed.
         These lines should be included (possibly as blank lines)
         even when MPRINT=0, i.e. when the printer plots are not
         required.

6) Cards controlling the approximation of elastic parameters and
   density distribution

6a) ISQRT,IRHO                                    FORMAT(14I5)

   ISQRT... controls the approximation of elastic parameters.
         ISQRT=0:  Approximation in elastic parameters.
         ISQRT=1:  Approximation in square roots of elastic
         parameters (velocity approximation).
         Note: In anisotropic layers (IANI.NE.0), ISQRT=0
         automatically.

   IRHO...     controls the density distribution in the model.
         IRHO=0:   The density RHO at any point of the model is
         determined from the relation RHO=1.7+0.2*VP, where VP
         is the square root of the elastic parameter A11 (which
         in an isotropic layer corresponds to the P-wave velocity).
         IRHO=1:   The density is constant throughout each layer
         with value specified in line 6b.

6b) This line is read only if IRHO.NE.0.
    RHO(1),RHO(2),...,RHO(NINT-1)                 FORMAT(8F10.5)

   RHO(I)... density in the I-th layer. It has meaning only if
         IRHO=1. If RHO(I) is not specified, RHO(I)=1.0 by default.

7) Description of distribution elastic parameters in individual
   layers. The input data differ for isosurface interpolation and
   B-spline approximation.

**********************
B-spline approximation
**********************

   The elastic parameters are specified at grid points of a rectangular
   network. The network may be different in different layers and it
   must always cover the whole layer. The input data 7a-7d repeat
   (NINT-1)times, one set for each layer, from the top to the bottom of
   the model.

7a)IANI,SIGMA                                     FORMAT(I10,F10.5)

   IANI...     specifies properties of the layer.
         IANI=0:   The layer is isotropic.
         IANI=1:   The layer is anisotropic.

   SIGMA...    the tension factor, see [4]. This value indicates
         the curviness desired. If ABS(SIGMA) is nearly zero
         (e. g. .001) the resulting surface is approximately the
         tensor product of cubic splines. If ABS(SIGMA) is large
         (e. g. 50.) the resulting surface is approximately tri-
         linear. If SIGMA equals zero tensor products of cubic
         splines result. A standard value for SIGMA is approximately
         1. In absolute value.

7b)NX,NY,NZ                                      FORMAT(26I3)

   NX... number of grid lines in the x-direction (NX.GE.1). For
         NX=1, the medium is homogeneous in the direction of the
         x-axis.
   NY... number of grid lines in the y-direction (NY.GE.1). For
         NY=1, the medium is homogeneous in the direction of the
         y-axis.
   NZ... number of grid lines in the z-direction (NZ.GE.1). For
         NZ=1, the medium is homogeneous in the direction of the
         z-axis.

7c)X1(1),X1(2),...,X1(NX)                        FORMAT(8F10.5)

   X1(I)...  coordinate of the I-th grid line in the x-direction.

7c)Y1(1),Y1(2),...,Y1(NY)                        FORMAT(8F10.5)

   Y1(I)...  coordinate of the I-th grid line in the y-direction.

7c)Z1(1),Z1(2),...,Z1(NZ)                        FORMAT(8F10.5)

   Z1(I)...  coordinate of the I-th grid line in the z-direction.

7d) Values of elastic parameters at the grid points of the 3-D grid.
   In an isotropic layer, when IANI=0, values of the squares of the P
   and S wave velocities are read in successively. In case of an
   anisotropic layer, individual elastic parameters are read in
   successively. The elastic parameters are ordered in the following
   way: A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33,A34,A35,A36,A44,
   A45,A46,A55,A56,A66. Each elastic parameter is read in in the
   following way

   (((W(I,J,K),K=1,NZ),J=1,NY),I=1,NX)           FORMAT(8F10.5)

   W(I,J,K) contains the value of II-th elastic parameter at the
         grid point (X(I),Y(J),Z(K)).

   Note: The maximum numbers of grid lines to specify the distribution
   of elastic parameters in a single layer cannot exceed, in the present
   version, 10. In order to change this number, the dimensions of arrays
   X1,Y1,Z1,W,WG1,WG2,C,VX,VY,VZ and TEMP must be changed correspondingly.
   The dimension of the array C must be of at least NX*NY*NZ and of the
   array TEMP of at least 3*MAX(NX, NY, NZ). The values of variables NW1
   and NW2 in the beginning of the routine MODEL (program block MODBS.FOR)
   must be also changed to the new value. The total number of grid points
   in the whole model cannot exceed, in the present version, 100. In order
   to change this number, the dimensions 100 of the arrays in the common
   block /EPAR/ must be changed correspondingly.

************************
Isosurface interpolation
************************

   The elastic parameters are specified at interfaces representing
   their isosurfaces. The input data 7a-7c repeat (NINT-1)times,
   one set for each layer, from the top to the bottom of the model.

7a)IANI,ANGU(1),ANGU(2),ANGU(3),ANGL(1),ANGL(2),ANGL(3)
                                                  FORMAT(I10,6F10.5)

   IANI...     specifies properties of the layer.
         IANI=0:   The layer is isotropic.
         IANI=1:   The layer is anisotropic.

   ANGU(1),...,ANGU(3)... the rotation angles, by which the
         "elasticity tensor" coordinate system is to be rotated into
         the "model" coordinate system, specified at the upper (U)
         interface bounding the layer. The angles are specified in
         degrees. ANGU(1) is an angle between positive direction of
         the "model" x-axis and horizontal projection of positive
         part of "crystal" z-axis. ANGU(1) is positive if measured
         from positive direction of "model" x-axis towards positive
         "model" y-axis. ANGU(2) is an angle between the positive
         directions of the "model" z-axis and "crystal" z-axis.
         ANGU(3) is a rotation angle in the plane perpendicular to
         the "crystal" z-axis. For ANGU(3)=0., "crystal" y-axis is
         horizontal, x-axis is perpendicular to it and points down.
         For ANGU(3) non-zero, the "crystal" x- and y-axes are
         rotated positvely from their positions for ANGU(3)=0.

   ANGL(1),...,ANGL(3)... the same as for ANGU, but at the lower (L)
         interface bounding the layer.

         Note 1: The values of rotation angles inside of a layer are
         determined in the same manner as the elastic parameters.
         They are obtained by a linear interpolation along vertical
         lines between the corresponding values of ANGU and ANGL.

         Note2: In the case of an isotropic layer the parameters ANGU
         and ANGL have no meaning.

7b)This line differs for isotropic and anisotropic layer.

   Anisotropic layer - IANI.NE.0.

   (A66U(J,K),J=1,6),K=1,6)                       FORMAT(6F10.5)

   A66U... elastic parameters (i.e. actual parameters divided by
         density) in a compressed 6x6 form at the upper (U) interface
         bounding the layer. It is sufficient to specify only right
         upper corner of the matrix A66U.

   Isotropic layer - IANI.EQ.0.

   A66U(1,1),A66U(4,4)                            FORMAT(6F10.5)

   A66U(1,1),A66U(4,4)... squares of P- and S-wave velocities at the
         upper (U) interface bounding the layer.

7c)This line differs for isotropic and anisotropic layer.

   Anisotropic layer - IANI.NE.0.

   (A66L(J,K),J=1,6),K=1,6)                       FORMAT(6F10.5)

   A66L... elastic parameters (i.e. actual parameters divided by
         density) in a compressed 6x6 form at the lower (L) interface
         bounding the layer. It is sufficient to specify only right
         upper corner of the matrix A66U.

   Isotropic layer - IANI.EQ.0.

   A66L(1,1),A66L(4,4)                            FORMAT(6F10.5)

   A66L(1,1),A66L(4,4)... squares of P- and S-wave velocities at the
         lower (L) interface bounding the layer.

8) This line is read only if ISURF=1 (generation of data for plots of
   velocity surfaces).

   NPAR                                          FORMAT(16I5)

   NPAR... switch indicating a section which wave surface is to
         be plotted.
         NPAR=1: slowness surface.
         NPAR=2: phase velocity surface.
         NPAR=3: group velocity surface.

9) This line is read only if  ISURF=1 (generation of data for plotsof
   velocity surfaces).

   X0,Y0,Z0,DDELTA                               FORMAT(4F10.5)

   X0,Y0,Z0... coordinates of the point for which a wave surface
         is to be plotted. The point should not be situated closer
         to the border of the model or interface than two time steps
         along the ray, DT, see input line 13.

   DDELTA... step in the angle of the phase normal (in radians) for
         plotting the section of a selected surface. In case of the
         slowness surface and phase velocity surface, it is the step
         with which points of the corresponding section are plotted.
         In case of group velocity, the density of plotted points
         changes according to the intensity of the energy flux. Where
         energy flux is denser, the points will be denser and vice
         versa. Default value, DDELTA=0.05 rad.

10) ICONT,MEP,MOUT,MDIM,METHOD,MREG,ITMAX,IPOL,IPREC,IRAYPL,
   IPRINT,IAMP,MTRNS,ICOEF,IRT,ILOC,MCOD,MORI    FORMAT(24I3)

   ICONT...  controls the continuation of the computations.
         ICONT=0: Indicates the termination of the computations.
           The line 10 is the last line in the input data set.
         ICONT=1:  The computations continue, the next line
           to be read in is the line 10.
         Note: For ICONT=0, the rest of the line 10 can be blank.

   MEP...  specifies whether initial-value ray tracing or
         two-point ray tracing is to be performed.
         MEP.EQ.0: Initial-value ray tracing is performed.
         MEP.NE.0: Two-point ray tracing is performed.
           ABS(MEP): The number of receivers, at which the rays
           should terminate.
           MEP=1: Computation of a ray terminating at a single
           receiver on a surface of the model or an internal
           interrface.
           MEP.LT.0 or MEP.GT.1: The receivers can be distributed
           along  a surface profile, a vertical profile or a
           profile along an interface, see parameter ILOC in
           this line; the orientation of the profile is specified
           by parameter PROF in line 11.
           MEP.GT.0: The receivers are distributed regularly
           along the profile (see line 11).
           MEP.LT.0: The receivers are distributed irregularly
           along the profile (see line 11).

   MOUT... controls printing of the results to the file LOU.
         MOUT=0:   Only the input data are reproduced and
           a list of codes is printed.
         MOUT=1:   Only the results for the rays terminating
           at specified receivers, and boundary rays are printed.
         MOUT=2:   The results for all rays terminating on a
           specified profile are printed.
         MOUT=3:   The results for all rays computed.
         Note: In case of initial-value ray tracing, MOUT=1
         already gives sufficient information about the results
         of computation.

   MDIM... controls the level of the computations.
         MDIM=0:   Only rays and travel times along them are
           computed.
         MDIM=1:   Spreading-free amplitudes, i.e. products of
           plane wave reflection/transmission/conversion
           coefficients are evaluated along the rays in addition
           to the quantities obtained for MDIM=0.
         MDIM=2:   Ray amplitudes including geometrical spreading

   METHOD...  specifies the method used for the iterative
         determination of rays in the two-point ray tracing
         procedure.
         METHOD=0: Paraxial ray approximation from the rays
            specified by the basic system of initial angles,
            see lines 14 and 15 or 16.
         METHOD=1: Paraxial ray approximation from the previously
            found two-point rays (fast but not always reliable).
         Note 1: In case of initial-value ray tracing (MEP=0),
         the parameter METHOD has no meaning.
         Note 2: Paraxial ray approximation can be used only if
         results of dynamic ray tracing are available. Therefore,
         for MEP.NE.0, MDIM is set automatically equal to 1,
         MDIM=1.

   MREG...  controls the evaluation of amplitudes at the termination
         points of the rays.
         MREG=0:   The conversion coefficients are applied if the
           termination point is situated on the surface.
         MREG=1:   The conversion coefficients are not applied.
         Note: MREG is set to be equal to 1 automatically
         for VSP computations and for receivers situated along
         internal interfaces (when ILOC.gt.0). In these cases
         conversion coefficients are not considered.

   ITMAX... number of iterations permitted in the two-point ray
         tracing procedure. The default value is ITMAX=10.
         Note: In case of initial-value ray tracing (MEP=0),
         the parameter ITMAX has no meaning.

   IPOL...  controls printing of polarization vectors in the routine
         POLAR.
         IPOL=0: Polarization vectors are not printed at all.
         IPOL=1: Polarization vectors are printed only at
         the points of incidence of a ray at an interface.
         IPOL=2: Polarization vectors are printed at all computed
         points along the ray (this option makes the computations
         more time consuming).

   IPREC...  switch controlling precission of the ray tracing and
         dynamic ray tracing calculations.
         IPREC=0: Precission is not tested.
         IPREC=1: Precission is tested and where inaccuracy exceeds
         the value RNULL (see description of the parameter INULL in
         input data 2), ray tracing and dynamic ray tracing are
         modified so that they satisfy the precission tests.
         IPREC=2: Precission is tested but no modification of results
         is performed. When inaccuracy exceeds the value RNULL, the
         values of unsatisfactory tests are printed at any point of
         any ray where the inaccuracy occurs.
         Note: As a test of precission the following relations are
         tested. (a) Scalar product of the slowness vector and the group
         velocity vector; its value should be p_i*v_i=1. (b) Derivative
         of the eikonal equation with respect to the two ray parameters
         r_1 and r_2 (take-off angles of the considered ray}; both
         derivatives should be zero. (c) Scalar products of the slowness
         vector and the vectors Q_{i1}=dx_i/dr_1 and Q_{i2}=dx_i/dr_2
         at constant travel time; their values should be p_i*Q_{iJ}=0.

   IRAYPL...   dummy parameter.

   IPRINT,IAMP,MTRNS,ICOEF,IRT... switches controlling storage of
         auxiliary quantities in the output file LOU. IPRINT controls
         output of quantities computed along the rays in routine
         and routine INV; IAMP controls the printing of
         the auxiliary quantities concerning amplitudes in the
         routine AMPL; MTRNS - transformation of the slowness vector
         at the interfaces in the routine TRANSL; ICOEF - reflection/
         transmission/conversion coefficient in the routine COEF;
         IRT - iterative determination of the point of intersection
         of a ray with an interface in the routine OUT.
         It is recommended to specify all these parameters zero,
         which means that no auxiliary quantities are printed.

   ILOC... controls the orientation and position of the profile, on
         which receivers are situated.
         ILOC=0:  Receivers along a surface profile.
         ILOC=1:  Receivers along a vertical profile.
         ILOC.GT.1:  Receivers along the ILOC-th interface.

   MCOD... controls specification of initial angles of rays of
         individual elementary waves.
         MCOD=0:  The system of initial angles of rays is determined
           once in the beginning of computations. It is the same for
           all elementary waves.
         MCOD=1:  The system of initial angles of rays is determined
           for each elementary wave separately.
   MORI... controls specification of ray parameters: angles specifying
         a ray at the source. The two angles are chosen as polar angles
         related to the "model" axis z and y. The reason for the two
         options is to avoid singular situations as that which occurs
         for vertical rays if polar angles are related to the z-axis.
         See specification of the initial slowness vector by ray
         angles specified in lines 14 and 15 or 16.
         MORI=0:  The ray angles are related to the "model" z-axis.
         MORI=1:  The ray angles are related to the "model" y-axis.

11) This line is read only if MEP.NE.0, i.e., in the two-point
   ray-tracing mode. Specification of a profile with receivers.
   Specification of receiver positions along this profile. Receivers
   are specified in the "receiver" cylindrical coordinate system
   with axis through the point [XREF,YREF].

   For MEP.LT.0, NDST=-MEP.

   PROF,DST(1),DST(2),...,DST(NDST),XPRF,YPRF     FORMAT(8F10.5)

   PROF...  azimuthal angle specifying a surface or interface profile,
         along which the receivers are situated - an angle between the
         vertical plane containing the profile and the positive
         "model" x axis. The profile passes through the vertical axis
         specified by the coordinates XPRF and YPRF. PROF has no meaning
         in the VSP mode, when receivers are distributed along a vertical
         profile (ILOC=1).
         PROF is measured in radians from the positive x axis of the
         "model" coordinates towards the positive y axis. For example,
         PROF=0.0 corresponds to the profile situated in the vertical
         plane containing positive "model" x axis, PROF=1.5707 to the
         profile situated in the vertical plane containing positive
         "model" y axis.

   DST(I)...   distances of receivers from the origin of the profile
         (XPRF,YPRF), see below, measured along the profile. In case
         of a surface or interface profile, DST(I) is the distance of
         the I-th receiver from the vertical axis specified by
         coordinates XPRF and YPRF. In case of a vertical profile,
         DST(I) is the "model" z-coordinate of the receiver.

   XPRF,YPRF... x- and y- "model" coordinates of the vertical axis,
         which contains an origin with respect to which receivers
         situated along surface or interface profiles are specified.
         The origin should be chosen so that the receivers are always
         to one side from it. Default values XPRF=XSOUR, YPRF=YSOUR.
         For XSOUR and YSOUR, see input data line 13.

   For MEP=1, NDST=1.

   XREC,YREC                                      FORMAT(8F10.5)

   XREC,YREC... x- and y-coordinates of the receiver situated on
         the surface of the model or an internal interface (the
         interface is specified by the wave code).

   For MEP.GT.1, NDST=MEP.

   PROF,RMIN,RSTEP,XPRF,YPRF                      FORMAT(8F10.5)

   PROF...  azimuthal angle specifying a surface or interface profile,
         along which the receivers are situated - an angle between the
         vertical plane containing the profile and the positive
         "model" x axis. The profile passes through the vertical axis
         specified by the coordinates XPRF and YPRF. PROF has no meaning
         in the VSP mode, when receivers are distributed along a vertical
         profile (ILOC=1).
         PROF is measured in radians from the positive x axis of the
         "model" coordinates towards the positive y axis. For example,
         PROF=0.0 corresponds to the profile situated in the vertical
         plane containing positive "model" x axis, PROF=1.5707 to the
         profile situated in the vertical plane containing positive
         "model" y axis.

   RMIN... distance of the first receiver from the origin of the profile
         (XPRF,YPRF), see below. In case of a surface or interface
         profile, RMIN is the distance of the first receiver from the
         vertical axis specified by coordinates XPRF and YPRF. In case
         of a vertical profile, RMIN is "model" z-coordinate of the first
         receiver.

   RSTEP... the step in distances of neighbouring receivers from the
         origin (XPRF,YPRF). The distance of the I-th receiver is
         determined from the relation DST(I)=RMIN+RSTEP*(I-1).

   XPRF,YPRF... x- and y- "model" coordinates of the vertical axis,
         which contains an origin with respect to which receivers
         situated along surface or interface profiles are specified.
         The origin should be chosen so that the receivers are always
         to one side from it. Default values XPRF=XSOUR, YPRF=YSOUR.
         For XSOUR and YSOUR, see input data line 13.

12) This line is read only when ILOC.eq.1 and MEP.ne.0.

   XVSP,YVSP                                      FORMAT(8F10.5)

   XVSP,YVSP... x and y coordinates of the vertical profile in case
         of VSP computations.

13) XSOUR,YSOUR,ZSOUR,TSOUR,DT,AC,REPS,PREPS       FORMAT(8F10.5)

   XSOUR,YSOUR,ZSOUR... (x,y,z) coordinates of the source.

   TSOUR... initial time.

   DT...  time step for the P wave in the integration of the ray
         tracing system. For S waves the time step is automatically
         multiplied by 1.7 to make the speed of computations comparable
         with P waves. DT should be positive. The default value is DT=1.

   AC...  the required accuracy of the integration of the ray
         tracing system. Recommended values of this parameter are in
         the range (0.0001, 0.001). The default value is AC=0.0001.

   REPS... the vicinity of a receiver on the profile. Iterative
         procedure terminates if a ray is found, which has its
         termination point within this vicinity.
         Note: In case of initial-value ray tracing (MEP=0), the
         parameter REPS has no meaning. The default value is
         REPS=0.05.

   PREPS... specifies maximum permitted deviation of the termination
         point of a ray from the profile with receivers. The deviation
         is measured as the length of the arc of the circle passing
         through the termination point, with the centre at epicentre,
         between the profile and the termination point.
         Note: In case of initial-value ray tracing (MEP=0), the
         parameter PREPS has no meaning. The default value is
         PREPS=0.05.

14) This line is read only when MCOD.EQ.0

   AMIN,ASTEP,AMAX                                FORMAT(8F10.5)

15) This line is read only when MCOD.EQ.0

   BMIN,BSTEP,BMAX                                FORMAT(8F10.5)

   AMIN,ASTEP,AMAX,BMIN,BSTEP,BMAX... define the basic system of
         initial angles (in radians) for ray tracing. The angles
         have different meaning in the modes MORI.EQ.0 and
         MORI.NE.0.
           For MORI.EQ.0, the angles A and B specify a slowness
         vector p at the source in the "model" coordinates in the
         following way:
                p=c**(-1)*(cosBcosA, sinBcosA, sinA).
           For MORI.NE.0, the angles A and B specify a slowness
         vector p at the source in the "model" coordinates in the
         following way:
                p=c**(-1)*(cosAcosB, sinB, sinAcosB).
         The same values of AMIN,ASTEP,AMAX and BMIN,BSTEP,BMAX are
         used for all considered waves.

16) KC,KREF,CODE(1,1),CODE(1,2),CODE(2,1),CODE(2,2),...
   ...CODE(KREF,1),CODE(KREF,2)                   FORMAT(26I3)

   KC... KC.EQ.0:  Ray of direct unconverted wave with no reflections
         and no conversions at interfaces.
         KC.NE.0:  Ray of multiply reflected, possibly converted wave.
         KC=1:  After leaving the source, the ray is to strike first
         the interface below the source.
         KC=-1: after leaving the source, the ray is to strike first
         the interface above the source.

   KREF... the number of simple segments of the ray. In case of the
         direct ray (KC=0), specify KREF=1.
         KREF=0 indicates that no elementary wave is to be computed.

   CODE(I,1)... the number of the layer in which the I-th simple
         segment of the ray is situated.

   CODE(I,2)... the wave type along the I-th simple segment of the
         ray.
         CODE(I,2)=1 or 2: quasishear waves (qS1,qS2).
         CODE(I,2)=3: quasicompressional wave (qP).
         Note 1: In an isotropic layer, CODE(I,2)=1 and 2 give the
         same shear wave.
         Note 2: CODE(I,J) is an integer.

17) The following two lines are read only when MCOD.NE.0. They are
    read in separately for each considered wave.

   AMIN,ASTEP,AMAX                                FORMAT(8F10.5)

   BMIN,BSTEP,BMAX                                FORMAT(8F10.5)

         The meaning of the above parameters is the same as in lines 12
    and 13.

Example of data LIN with isosurface interpolation
Example of data LIN with B-spline interpolation


Termination of computations
---------------------------

    The sequence of lines 16-17 can be repeated any number times,
separately for each considered elementary wave. To terminate the
computations of individual elementary waves put KREF=0 on the line
16 (or simply insert a blank line in the position of the line 16).
    The sequence of lines 10-17 can be repeated any number times
to perform different computations for the same model. Each system
of lines 10-17 produces file LU1 (if LU1.NE.0.AND.MEP.NE.0) for
plotting in program ANRAYPL, and file LU2 (if LU2.NE.0) for
program SYNTAN, in which ray synthetic seismograms are computed.
To terminate completely the computations, specify ICONT=0 (or
simply insert a blank line in the position of the line 10).


                                                      
Output to the file LU1
----------------------

    For LU1.NE.0 and MEP.NE.0, i.e. in the two-point ray-tracing mode,
certain important data concerning rays, travel times and amplitudes
are stored in the file LU1 for possible further processing. The file
LU1 can be used in the program ANRAYPL, included in this package, for
plotting ray diagrams, time-distance and amplitude-distance curves.
The data are stored in LU1 in the following formatted form:

(A)Information about the source and receivers.

1) ICONT,NDST,ILOC                               FORMAT(26I3)

      ICONT... indicative index.

         ICONT=0: Indication of the end of the data stored in the file
           LU1. The line 1) is the last data line in the file LU1.
         ICONT=1: Indication that a data set corresponding to one
           elementary wave follows.

      NDST... number of receivers along the considered profile.

      ILOC... controls the orientation and position of the profile,
         on which receivers are situated.
         ILOC=0:  Receivers along a surface profile.
         ILOC=1:  Receivers along a vertical profile.
         ILOC.GT.1:  Receivers along the ILOC-th interface.

2) RO                                            FORMAT(8F10.5)

      VP,VS,RO... P-wave velocity, S-wave velocity and density at the
         source.

3) NPN,NPN,NPN                                   FORMAT(26I3)

      Dummy data, NPN=2 by default.

4) APN,APN,APN,APN,APN                           FORMAT(5E15.5)

      Dummy data, APN=0.0 by default. There are two lines of these
      dummy data.

5) XPRF,YPRF,0.,PROF                        FORMAT(8F10.5)

      XPRF,YPRF... x- and y- "model" coordinates of the vertical axis,
         with respect to which receivers situated along surface or
         interface profiles are specified, see input data line 11.

      PROF... azimuthal angle specifying a surface or interface
         profile, along which the receivers are situated - an angle
         between the vertical plane containing the profile and the
         positive "model" x axis. PROF has no meaning in the VSP mode,
         when receivers are distributed along a vertical profile
         (ILOC=1). PROF is measured in radians from the positive x
         axis of the "model" coordinates towards the positive y axis.
         For example, PROF=0.0 corresponds to the profile situated
         in the vertical plane containing positive x axis,
         PROF=1.5707 to the profile situated in the vertical plane
         containing positive y axis.

6) (DST(I),I=1,NDST)                             FORMAT(8F10.5)

      DST(I)... coordinates of receivers along the profile. In case
         of a surface or interface profile, DST(I) is the radius of
         the receiver measured from the vertical axis passing through
         the source, see the description of the "receiver" coordinate
         system. In case of a vertical profile DST(I) is z-coordinate
         of the receiver in the "model" coordinate system.

(B)Information about rays. For each ray with a termination point at
   the specified receiver, the quantities given in the following
   lines 7 and 8 are stored. The data given in 7 and 8 are stored
   successively for all rays terminating at specified receivers.
   After the last ray, the parameter N (see 7) is set to zero.
   Then, the quantities in section C are stored.

7) N,II                                          FORMAT(26I3)

      N.NE.0... the number of points along the ray.
      N.EQ.0... indication of the end of the information about rays.
      II... dummy parameter.

      Note: Maximum value of N is 400.

8) (AY(2,I),AY(3,I),AY(4,I),I=1,N)               FORMAT(3E15.5)

      The (x,y,z) coordinates of the points along the rays, starting
      from the source (I=1) and ending with the termination point
      (I=N).

(C)Information about travel times and amplitudes.

9) INDEX                                         FORMAT(26I3)

      ABS(INDEX)... the total number of rays of the considered elementary
      wave, for which travel times and ray amplitudes are stored.
      INDEX.LT.0 indicates a shear wave in an isotropic medium.

      Note: Maximum value of ABS(INDEX) is 500.

10)(INDI(I),XCOOR(I),ZCOOR(I),TIME(I),AMPX(1,I),AMPY(1,I),AMPZ(1,I)
                                                 FORMAT(I5,9E15.5)

      INDI(I)... specifies the successive number of the receiver at
        which the I-th ray terminates.

      XCOOR(I)... radius of the termination point of the ray in the
        "receiver" cylindrical coordinate system, measured from the
        vertical axis through the source.

      ZCOOR(I)... depth of the termination point in the "model"
        coordinates.

      TIME(I)... travel time at the termination point of the I-th ray.

      AMPX(1,I),AMPY(1,I),AMPZ(1,I)...  complex amplitudes of the x,
         y and z components of the polarization vector in the "model"
         coordinates; for an S wave in an isotropic medium, complex
         amplitudes corresponding to the first polarization vector in
         the plane perpendicular to the ray, see the description of the
         model and the source in the introduction.

11) This line is stored only if INDEX.LT.0.

   AMPX(2,I),AMPY(2,I),AMPZ(2,I)                 FORMAT(6E15.5)

         AMPX(2,I),AMPY(2,I),AMPZ(2,I)... complex amplitudes of the
         components of the second polarization vector of an S wave,
         see the description of the model and the source in the
         introduction.

12)(P(I),I=1,3)                                  FORMAT(6E15.5)

         P(I)... components of the slowness vector of the considered
         wave at the source.

13)(POL(I),I=1,3)                                FORMAT(6E15.5)

         POL(I)... components of the polarization vector of the
         considered wave at the source. In case of an S wave in
         an isotropic medium, components of its first polarization
         vector in the plane perpendicular to the ray, see the
         description of the model and the source in the introduction.

14)This line is stored only if INDEX.LT.0.
   (POL1(I),I=1,3)                               FORMAT(6E15.5)

         POL1(I)... components of the second polarization vector of
         an S wave in an isotropic medium, see the description of the
         model and the source in the introduction.

   The lines 10-14 are repeated ABS(INDEX)-times.

   The sequence of data 1-14 repeats for each considered elementary
wave. Indication of the end of information in LU1 is ICONT=0 in 1).


                                                      
Output to the file LU2
----------------------

    For LU2.NE.0 and MEP.NE.0, i.e., in the two-point ray-tracing
mode, data necessary for construction of synthetic seismograms are
stored in the file LU2. Program SYNTAN, included in this package,
can be used for this purpose. The data are stored in LU2 in the
following formatted form:

1) MTEXT                                         FORMAT(17A4)

      MTEXT... alphanumeric test describing the computations.

2) NDST,KSH,ILOC...                              FORMAT(26I3)

      NDST... number of receivers along the profile.

      KSH... dummy parameter, automatically KSH=2.

      ILOC... controls the orientation and position of the profile,
         on which receivers are situated.
         ILOC=0:  Receivers along a surface profile.
         ILOC=1:  Receivers along a vertical profile.
         ILOC.GT.1:  Receivers along the ILOC-th interface.

3) XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,ROS             FORMAT(8F10.5)

      XSOUR,YSOUR,ZSOUR... coordinates of the source in the "model"
         coordinate system.

      TSOUR... initial time.

      RSTEP... an average distance between neighbouring receivers.

      ROS... density at the source.

4) DST(I),I=1,NDST                               FORMAT(8F10.5)

      DST(I)... coordinates of receivers along the profile. In case
         of a surface or interface profile, DST(I) is the radius of
         the receiver measured from the vertical axis passing through
         the source, see the description of the "receiver" coordinate
         system. In case of a vertical profile, DST(I) is z-coordinate
         of the receiver in the "model" coordinate system.

7) NC,II,T,AMPX1,AMPY1,AMPZ1,TAST     FORMAT(2I3,F10.5,12E12.6,F10.6)

      IABS(NC)... number of the considered elementary wave. NC.LT.0
         indicates that the wave starts as an S wave from the source
         situated in an isotropic layer.

      II... successive number of the receiver position.

      T... arrival time.

      AMPX1,AMPY1,AMPZ1... complex amplitudes of the x, y and z
         components of the polarization vector in the "model"
         coordinates; for an S wave in an isotropic medium, complex
         amplitudes corresponding to the first polarization vector in
         the plane perpendicular to the ray, see the description of the
         model and the source in the introduction.

      TAST... the global absorption factor t*; in this version,
         automatically t*=0.

8) This line is stored only if NC.LT.0.

   AMPX2,AMPY2,AMPZ2                             FORMAT(6E12.6)

         AMPX2,AMPY2,AMPZ2... complex amplitudes of the components of
         the second polarization vector of an S wave, see the description
         of the model and the source in the introduction.

9) (P(I),I=1,3)                                  FORMAT(3F10.5)

         P(I)... components of the slowness vector of the considered
         wave at the source.

10)(POL(I),I=1,3)                                FORMAT(3F10.5)

         POL(I)... components of the polarization vector of the
         considered wave at the source. In case of an S wave in
         an isotropic medium, components of its first polarization
         vector in the plane perpendicular to the ray, see the
         description of the model and the source in the introduction.

11)This line is stored only if NC.LT.0.
   (POL1(I),I=1,3)                               FORMAT(3F10.5)

         POL1(I)... components of the second polarization vector of
         an S wave in an isotropic medium, see the description of the
         model and the source in the introduction.

   The lines 7-11 are repeated ABS(NC)-times.

                                                      
Output to the file LU3
----------------------

    For ISURF=1, data necessary for construction of sections of
slowness surfaces, phase velocity and group velocity surfaces are
stored in the file LU3. Program VELPL, included in this package,
can be used for plotting the sections. The data are stored in LU3
in the following formatted form:

1) ITYPE,NPAR,INUM                               FORMAT(16I5)

      ITYPE...  ITYPE=1: qS1 wave.
                ITYPE=2: qS2 wave.
                ITYPE=3: qP wave.

      NPAR...   NPAR=1: plot of section of slowness surface.
                NPAR=2: plot of section of phase velocity surface.
                NPAR=1: plot of section of group velocity surface.

2) DDELTA,XX(I),YY(I),ZZ(I)                       FORMAT(4E15.5)

      DDELTA... angle (in radians), for which the vector (XX,YY,ZZ)
                is calculated.

      XX(I),YY(I),ZZ(I)... x,y and z components of the slowness vector
                (if NPAR=1), phase velocity vector (if NPAR=2) or
                group velocity vector (if NPAR=3) for the angle
                DDELTA. Maximum points: 300.

   The line 2 is repeated for all the angles of all the three waves.


                                                      
Output to the file LOU
----------------------

    Storage of results in the file LOU is controlled by the switches
MPRINT and MOUT. The switch MPRINT controls storage of parameters
describing the model. The switch MOUT controls the storage of
quantities during the initial-value or two-point ray tracing.

MPRINT:
   MPRINT=0: Only a copy of input data, no other description of the
         model is presented.
   MPRINT=1: Writes the most important input data plus printer plots
         of the form of interfaces in the model (see ZMIN,ZMAX in
         line 5).
   MPRINT=2: The same as for MPRINT=1 plus description of the
         unrotated and rotated compressed matrices of the elastic
         parameters.

MOUT:
   MOUT=0: Only input data are reproduced and list of the codes of
         the computed waves.
   MOUT=1: For each code, results describing all rays in the case of
         initial-value ray tracing and the rays arriving at specified
         receivers in the case of two-point ray tracing are stored.
         In the latter case, in addition, all boundary and critical
         rays are stored. All the results are stored in the routine
         RECEIV.

      The results for a ray terminating at a prescribed receiver are
      stored in two lines as follows:

         DD,R,XX,YY,ZZ,T,DEC,AZIM,IND,IND1,ITER,II,INDEX
         UU,(APX(I),APY(I),APZ(I),I=1,2),SPR,KMAH

         DD: Coordinate of the receiver (in the "receiver" coordinate
             system), at which the ray should arrive. For initial-
             value ray tracing, DD has no meaning, therefore, DD=0.0
             by default.

         R:  In the two-point ray tracing to receivers situated along
             the surface or an interface, R is the radius of the
             termination point in the "receiver" coordinates. In the
             two-point ray tracing to receivers distributed along a
             vertical profile, R is z-coordinate in "model"
             coordinates. R should differ from DD by less then REPS,
             see line 13 in input data file LIN.

         XX,YY,ZZ: "Model" coordinates of the termination point of
             the ray.

         T:   Corresponding arrival time.

         AZIM,DEC: Azimuth and declination specifying the slowness
             vector of the considered ray at the source (in radians).

         IND: Parameter describing the history of the ray. The most
             important values of IND are:
             (1 or 2) - the ray intersects the vertical (y,z)
                       boundary of the model with lower (1) or larger
                       (2) value of the coordinate x.
             (3) -     termination of the ray on the model surface.
             (4) -     termination of the ray on the lower boundary
                       of the model.
             (5) -     in the Runge-Kutta method, the basic time step
                       DT is halved more than 10 times, without reaching
                       the required accuracy AC, see line 13. There
                       is probably something wrong with your velocity
                       distribution in the model.
             (6) -     in line 13, DT=0 (not allowed). The computation
                       of the ray diagram does not start.
             (7) -     in line 13, DT.lt.0 (not allowed). The computation
                       of the ray diagram does not start.
             (8) -     termination of the ray at an interface because
                       its ray path does not correspond to the prescribed
                       code.
             (9) -     overcritical incidence of the ray. Termination of
                       calculation of the ray.
             (10) -    coinciding values of the phase velocities of
                       quasishear waves for the given direction.
                       Termination of calculation of ray.
             (11) -    matrix for the evaluation of R/T coefficients
                       singular. Termination of calculation of ray.
             (12) -    travel time along the ray exceeded the value TMAX.
                       This value is set TMAX=10000. Termination of
                       calculation of the ray.
             (13) -    the ray is specified by more than 400 points.
                       Its computation is terminated.
             (14) -    specified code does not correspond to the
                       position of the source in the model.
             (15) -    prescribed reflection from the bottom boundary
                       of the model while MDIM.NE.0. Termination of
                       calculation of the ray.
             (20) -    two neighbourhood interfaces in the model intersect
                       each other. The computation terminates.
             (30 or 31) - the ray intersects the vertical (x,z)
                       boundary of the model with lower (30) or larger
                       (31) value of the coordinate y.
             (39) -    number of iterations necessary for the determination
                       of the point of intersection of a ray with an
                       interface exceeded 10. Termination of calculation
                       of the ray. Proable cause: too large integration
                       step along a ray, DT, see input data line 13.
             (43) -    the ray terminates at the vertical plane perpendicular
                       to the profile source-borehole in the VSP mode.
             (50) -    the source is situated outside the model, the
                       computation of the ray diagram does not start.
             (IND.GT.100) - termination of the ray on the (IND-100)th
                       interface.

         IND1: ABS(IND1) - number of the interface last hit by the
                 ray. Negative sign of IND1 indicates at least one
                 compound segment of the ray, see description of
                 waves and wave codes.

         ITER: Number of iterations necessary to find the ray arriving
                 to the prescribed receiver.

         II: Dummy parameter.

         INDEX: Current number of the ray in the system of rays for
             the given code.

         UU: Square root of the product of the ratios of impedances
             of the reflected/transmitted to the incident waves at
             individual points of incidence along the ray path,
             including the ratios of cosines of angles of reflection/
             transmission to the cosines of angles of incidence at
             the same points.

         APX(J,I),APY(J,I),APZ(J,I)...  for termination point in an
             anisotropic medium or for P wave arrival to the
             termination point situated in an isotropic medium,
             APX(1,I),APY(1,I),APZ(1,I) are complex amplitudes of
             the x,y and z components of the polarization vector of
             the considered wave in the "model" coordinates; for
             an S wave arrival to termination point situated in an
             isotropic medium, APX(J,I),APY(J,I),APZ(J,I) are complex
             amplitudes of the x,y and z components corresponding to
             the two (J=1,2) unit vectors in the plane perpendicular
             to the ray, describing the polarization of S waves, see
             the description of the model and the source in the
             introduction.

         SPR... geometrical spreading.

         KMAH... index specifying how many times the considered ray
         touched a caustic.

      The information about boundary or critical rays determined in
      iterations along the profile is stored in one line of the form:

         R,ZZ,T,DEC,IND,IND1,IRES   'BOUNDARY (CRITICAL) RAY'

         The meaning of the individual parameters is the same as
         above.

         IRES: Indicates whether the boundary ray terminating on the
             considered profile was found (IRES=1) or not (IRES=0).

    MOUT=2: This option has meaning only in the two-point ray-tracing
         mode. For each code the results concerning all of the rays
         ariving to the considered profile are stored. All the
         results are stored in the routine RECEIV.

      The results for the rays specified by declinations from the
      basic system of declinations (see line 14 in input data
      file LIN) are stored in one line as follows:

         IND,IND1,KMAH,R,T,DEC.

      The results for the rays in the iteration to a specified
      receiver are stored in one line as follows:

         'ITERATIONS'  IND,IND1,ITER,KMAH,DD,R,T,DEC,AZIM.

      The results for the rays iterating to the boundary or critical
      rays are stored in one line as follows:

         'BOUNDARY (CRITICAL) RAY'  IND,IND1,ITER,R,T,DEC.

      The results for the rays arriving to prescribed receivers, for
      boundary and critical rays are stored in the same form as for
      MOUT=1. In addition to that information the results of the
      dynamic ray tracing are stored in the form of two 3x3 matrices

         XDYN

         YDYN

      The first matrix contains derivatives of the Cartesian coordinates
      with respect to the ray coordinates, the second one the derivatives
      of the components of the slowness vector with respect to the ray
      coordinates.

      The meaning of other symbols in the above results is the same as
      explained for the parameter MOUT=1.

    MOUT=3: Similar information as above is stored for all rays in the
         iterative procedure. The results stored in the routine PROFIL
         are marked by asterisk * in the beginning of line to
         distinguish them from results stored in the routine RECEIV.
         The use of the option MOUT=3 is not recommended.

Description of COMMON block/RAY/
--------------------------------

  COMMON/RAY/ AY(28,400),DS(20,20),KINT(20),HHH(3,3),TMAX,PS(3,7,20),
              IS(8,20),DINC,DTOLD,N,IREF,IND,IND1

  AY(I,N)...  values of quantities along the ray.
              N: successive number of a point on the ray.
              I=1: travel time
              I=2-4: x,y and z coordinates of the point of the ray.
              I=5-7: x,y and z components of the slowness vector.
   IANI.NE.0  I=8-28: elastic parameters of an anisotropic medium
   IANI=0     I=8-11: P wave velocity and its first derivatives with
                      respect to x, y and z coordinates.
              I=12-15: S wave velocity and its first derivatives with
                      respect to x, y and z coordinates;
              I=16-28 unspecified.

  DS(I,IREF)...  values of quantities at the point of incidence at
              an interface.
              IREF: successive number of the point of incidence.
              I=1-3: x,y and z components of the unit normal to the
                interface.
              I=4-6: x,y and z components of the projection of the
                slowness vector of the incident wave into the plane
                tangent to the interface.
              I=7: positive value of the projection of the group
                velocity vector of the incident wave to the normal
                to the interface.
              I=8: value of density on the side of incident wave.
              I=9: no meaning.
              I=10: positive value of the projection of the group
                velocity vector of the generated wave to the normal
                to the interface.
              I=11: value of density on the side of generated wave.
              I=12: no meaning.

   KINT(N)... successive numbers of points along the ray at points
              of incidence.
              N: successive number of the element of the ray.
              For compound element of the ray KINT=0.

   HHH(I,J)... I-th component of the J-th vector of the ray-centered
              vector basis e1, e2, t.

   TMAX...    maximum value of travel time; if it is exceeded, the
              calculation of the ray terminates.

   PS(I,J,K)... slowness vectors of the incident and generated
              waves at points of incidence at interfaces;
              I: index specifying the x-, y- and z-component
              of the slowness vector (I=1,2,3).
              J: index specifying the type of the considered wave;
              J=1-3: reflected wave, J=4-6: transmitted wave;
              J=7: incident wave.
              K: successive number of the point of incidence.

   IS(I,J)... values of quantities related to the J-th point of incidence
              at an interface.
              J: successive number of the point of incidence.
              IS(1,J): number of the layer, in which generated wave
                propagates.
              IS(2,J): specifies behaviour of the wave at an interface,
                specified in TRANSL, see variable ITRANS
                IS(2,J)=0: reflection,
                IS(2,J)=1: transmission,
                IS(2,J)=999: surface conversion.
              IS(3,J): number of the layer in which incident wave
                propagates.
              IS(4,J): dummy parameter
              IS(5,J): number of inhomogeneous waves generated in the
                layer where reflected wave propagates.
              IS(6,J): number of inhomogeneous waves generated in the
                layer where transmitted wave propagates.
              IS(7,J): number of the layer in which reflected wave
                would propagate.
              IS(8,J): number of the layer in which transmitted wave
                would propagate.

   N...       curent number of a point along the ray.

   IREF...    curent number of the point of incidence.

   IND...     parameter describing the history of the ray, see
              the description of the output to the file LOU.

   IND1...    ABS(IND1) - number of the last interface hit by the
              ray on its way to the receiver, see the description
              of the output to the file LOU.


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] 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.

[4] Cline,A.K.,1981. FITPACK, Dept. of Computer Sciences, Univ. of
    Texas at Austin.