Program ANRAY - Version 4.70
****************************

   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;
  - density-normalized 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 density-normalized elastic parameters
    specified in grid points of a 3-D grid;
  - in the B-spline mode, a water layer can be considered;
  - 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;
  - receivers can record either displacement vectors or pressure
    (if pressure is recorded, it is stored in the x-component of
    the amplitude vector AMPX);
  - conversion coefficients can be considered at receivers situated
    on the surface as well as at receivers situated at internal
    interfaces;
  - 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.
  - the package allows to perform calculations of the qS wave
    synthetics in the quasi-isotropic approximation.
  - the results of the computations can be displayed in the form
    of 2-D plots generated in the PostScript format or in the form
    of 3-D plots using the graphics system of the CRT package.

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 density-normalized
elastic parameters 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

    Tensor of density-normalized elastic parameters may be specified in
a coordinate system whose axes coincide with the principal axes of the
elasticity tensor. In such a case, in anisotropic media with higher
symmetry many of the 21 density-normalized elastic parameters are zero.
Since the axes of the "elastic tensor" coordinate system may differ from
the "model" coordinate system an additional information specifying how
the "elasticity tensor" coordinate system is to be rotated within the
model must be supplied, see description of input data for isosurface
interpolation sub 7a. 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
right handed.


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 density-normalized 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 density-
normalized elastic parameters have the dimensions of velocity squared.
The P and S wave velocities in an isotropic layer are, for simplicity,
also referred to as density-normalized elastic parameters. Two methods
of the determination of the density-normalized elastic parameters in
individual layers may be used.

    (a) Linear interpolation of density-normalized elastic parameters
along vertical lines between interfaces representing isosurfaces of
density-normalized elastic parameters. In this approach, the distribution
of density-normalized elastic parameters is laterally varying within
the layers if their boundaries are curved or inclined. The vertical
gradient of the density-normalized 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 density-normalized 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' 'LUBRD' 'LUGRD' 'LUIND' 'LURAY'/
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.
    'LUBRD' is the name of the output data file used in generating a VRML file.
    'LUGRD' is the name of the output data file used in generating a VRML file.
    'LUIND' is the name of the output data file used in generating a VRML file.
    'LURAY' is the name of the output data file used in generating a VRML file.
    / is a slash recomended in batch and script files to enable future
        extensions.
Defaults:
    'LIN'='anray.dat'
    'LOU'='anray.out'
    'LU1'='lu1.dat'
    'LU2'='lu2.dat'
Example of the main input data:
    'anray.an' / (B-spline interpolation; modbs.for)
or
    'anray.qi' / (B-spline interpolation; modbs.for)
or
    'anray.chs' / (linear interpolation; modis.for)
or
    'anray.vel' / (B-spline interpolation; modbs.for)
or
    'anray.rfr' / (linear interpolation; modis.for)

    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. If the file LU1 is specified, output data for
plotting ray diagrams, travel time curves and amplitude-distance curves
in the program ANRAYPL are stored in it. If the file LU2 is specified,
output data for computing ray synthetic seismograms in the program
SYNTAN are stored in it. If the file LU3 is specified, output data for
plotting sections of velocity surfaces are stored in it. If the files
LUBRD, LUGRD, LUIND and LURAY are specified, output data for plotting
model box, interfaces, medium parameter distribution and rays in 3-D
graphics mode are stored in them for the use in the program ANRAYWRM.

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

1) MTEXT

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

2) INULL,ISURF

   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 INULL=4.
   ISURF...    dummy parameter

3) MPRINT,NINT,N1,N2

   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. Default value.
         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 density-
         normalized elastic parameters.

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

   N1,N2... number of points along the first and second horizontal axes
         for triangulation of interfaces for 3-D images, see routine
         FACETS. Default N1=10, N2=10.

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

4a)MX,MY

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

4b)SX(1),...,SX(MX)

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

4c)SY(1),...SY(MY)

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

4d)Z(1),...,Z(MX*MY)

   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

   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.

         ZMIN of the shallowest interface and ZMAX of the deepest
         interface are used as the z-coordinates of the upper and the
         bottom horizontal planes limiting the model box in 3-D
         images.

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

6a) ISQRT,IRHO

   ISQRT... controls the approximation of density-normalized elastic
            parameters.
         ISQRT=0:  Approximation in density-normalized elastic parameters.
         Default value.
         ISQRT=1:  Approximation in square roots of density-normalized
         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 density-normalized elastic
         parameter A11 (which in an isotropic layer corresponds
         to the P-wave velocity). Default value.
         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)

   RHO(I)... density in the I-th layer. It has meaning only if
         IRHO=1. RHO(I)=1 by default.

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

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

   The density-normalized 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

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

   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. Default SIGMA=0.

7b)NX,NY,NZ

   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. Default NX=1.
   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. Default NY=1.
   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. Default NZ=1.

7c)X1(1),X1(2),...,X1(NX)

   X1(I)...  coordinate of the I-th grid line in the x-direction.
             Default X1(1)=0.

7c)Y1(1),Y1(2),...,Y1(NY)

   Y1(I)...  coordinate of the I-th grid line in the y-direction.
             Default Y1(1)=0.

7c)Z1(1),Z1(2),...,Z1(NZ)

   Z1(I)...  coordinate of the I-th grid line in the z-direction.
             Default Z1(1)=0.

7d) Values of density-normalized 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 density-normalized elastic
   parameters are read in successively. The density-normalized 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
   density-normalized elastic parameter is read in in the following way

   DO I=1,NX
   DO J=1,NY
   W(I,J,K),K=1,NZ

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

   Note1: It is possible to specify a water layer by specifying
   S-wave velocity equal zero. This option is applicable only in
   the B-spline approximation (modbs.for).

   Note2: The maximum number of grid lines to specify the distribution
   of density-normalized 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
   (NW1 and NW2 are the first two dimensions of the array W in MODBS.FOR;
   they should be chosen so that NW1.ge.NX and NW2.ge.NY. 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 density-normalized 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)

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

   ANGU(1),...,ANGU(3)... the rotation angles, by which the "elasticity
         tensor" coordinate system is to be rotated at the upper (U)
         interface bounding the layer. The angles are specified in
         degrees. ANGU(1) is a rotation angle around the z-axis of the
         "elasticity tensor" coordinate system. ANGU(1) rotates the
         x- and y-axes of the "elasticity tensor" system into new
         positions x' and y'. ANGU(1) is positive if the rotation is
         anticlockwise. ANGU(2) is a rotation angle around the y'-axis,
         which rotates the "elasticity tensor" coordinates x' and z
         into new positions x'' and z'. ANGU(2) is positive if the
         rotation is anticlockwise. ANGU(3) is a rotation angle around
         the z'-axis. Default ANGU(1)=0., ANGU(2)=0., ANGU(3)=0.

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

         Note 1: The values of rotation angles inside a layer are
         determined in the same manner as the density-normalized
         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)

   A66U... density-normalized elastic parameters 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)

   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)

   A66L... density-normalized elastic parameters 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)

   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 the file LU3 is specified (generation of
   data for plots of sections of velocity surfaces).

   NPAR,LAY

   NPAR... switch indicating which wave surface is to be plotted.
         NPAR=1: slowness surface. Default value.
         NPAR=2: phase velocity surface.
         NPAR=3: group velocity surface.
   LAY... number of the layer, in which the point specified on line 9
         is situated. Default LAY=1.

9) This line is read only if the file LU3 is specified (generation of
   data for plots of velocity surfaces).

   X0,Y0,Z0,DDELTA,AZIM

   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. Deafault X0=1., Y0=1.,
         Z0=1.

   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 DDELTA=0.05.

   AZIM... azimuth (in radians) of the vertical plane intersecting
         the studied surface. For AZIM=0., a section by the (x,z)
         plane is generated. For AZIM=1.570796, a section by the
         (y,z) plane is generated. Default AZIM=0.

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

   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. Default value.
         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.
         Default value.
         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.1: 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. Default value.
         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. Default value.
         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. Default value.
         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.

   MREG...  controls the evaluation of amplitudes at the termination
         points of the rays.
         MREG=0:   Elastic case: displacement vector calculated.
           The conversion coefficients are applied. Default value.
         MREG=1:   Elastic case: displacement vector calculated.
           The conversion coefficients are not applied.
         MREG=2:   Acoustic case: pressure calculated. The conversion
           coefficients are applied.
         MREG=3:   Acoustic case: pressure calculated. The conversion
           coefficients are not applied.
         Note: MREG is set to be equal to 1 automatically
         for VSP computations (when ILOC.EQ.1). In this case
         elastic case is considered and conversion coefficients are
         not used.

   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.
         Default value.
         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. Default value.
         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.
         Default values IPRINT=0, IAMP=0, MTRNS=0, ICOEF=0, IRT=0.

   ILOC... controls the orientation and position of the profile, on
         which receivers are situated.
         ILOC=0:  Receivers along a surface profile. Default value.
         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. Default value.
         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.
         Default value.
         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

   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.
         PROF has no meaning when receivers are distributed along a
         vertical profile.

   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

   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). The azimuth
         and the epicentral distance of the receiver from the source
         are determined automatically and appear in the file LOU.
         If BMIN, BSTEP and BMAX in line 15 or 17 are not specified,
         they are specified automatically from the values of the
         azimuth and the epicentral distance.

   For MEP.GT.1, NDST=MEP.

   PROF,RMIN,RSTEP,XPRF,YPRF

   PROF...  azimuthal angle (in radians) 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.
         PROF has no meaning when receivers are distributed along a
         vertical profile.

   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

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

13) XSOUR,YSOUR,ZSOUR,TSOUR,DT,AC,REPS,PREPS

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

   TSOUR... initial time. Default TSOUR=0.

   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 in the length
         units used for specification of coordinates. 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 in the length
         units used for specification of coordinates. 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

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

   BMIN,BSTEP,BMAX

   AMIN,ASTEP,AMAX,BMIN,BSTEP,BMAX... define the basic system of
         initial angles (in radians) for ray tracing.
         Default values BMIN=PROF(1)-.3, BSTEP=.6, BMAX=PROF(1)+.4
         The initial 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)

   KC... KC.EQ.0:  Ray of direct unconverted wave with no reflections
         and no conversions at interfaces. The ray terminates on
         boundaries of the model.
         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 must be used
         for the 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

   BMIN,BSTEP,BMAX

         The meaning of the above parameters is the same as in lines 14
    and 15.


Example of data LIN for anisotropic model 'an'
Example of data LIN for background isotropic model 'qi' for QI computations
Example of data LIN for anisotropic model 'chs' (isosurface approximation)
Example of data LIN for generating file for plotting phase velocity sections 'vel'
Example of data LIN for generating VRML file 'rfr'

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 for plotting in program ANRAYPL
(if ISPACE=0) or for plotting rays in a graphical mode (ISPACE=1),
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)

      RO... 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... specifies the number of the receiver at which the ray
                terminates.


      Note: Maximum value of N is 2000.

8) (T(J),X(J),Y(J),Z(J),P1(J),P2(J),P3(J),A11(J),A44(J),RHO(J),
   (E(1,K,J),K=1,3),(E(2,L,J),L=1,3),J = 1,N)    FORMAT(16E15.5)

      Values of important quantities at the points along rays,
      starting from the source (J=1) and ending with the termination
      point (J=N). T - travel time, X,Y,Z coordinates of a point
      on a ray, P1,P2,P3 - components of the slowness vector,
      A11,A44 - density-normalized elastic parameters (squares of P
      and S wave velocities in isotropic media). RHO - density,
      E(1),E(2) - basis vectors of the ray-centred coordinate system.

(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. If pressure is
         recorded, it is stored in AMPX.

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(A)

      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. If pressure is
         recorded, it is stored in AMPX.

      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 LU3.NE.0, 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.

      INUM...   number of samples in the plot.

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 density-
         normalized 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 2000 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 quantities related to derivatives of
      the Cartesian coordinates with respect to the ray coordinates,
      the second one quantities related to derivatives of components
      of the slowness vector with respect to the ray coordinates. In
      fact, XDYN and YDYN are obtained from the dynamic ray tracing
      equations, see Appendix A of [5], with the initial conditions
      (A.5), at the same reference.

      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,2000),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: density-normalized 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: the density
              I=17-19: components of the e_1 vector
              I=20-22: components of the e_2 vector
              I=23-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.

[5] Psencik,I., Teles,T.N., 1996. Point source radiation in
    inhomogeneous anisotropic structures, PAGEOPH, 148, 591-623.