Program WEAKAN
**************

   Program WEAKAN is designed for the calculations in the quasi-
isotropic (QI) approximation. This version calculates QI ray 
amplitudes of the qS waves propagating in continuous, inhomogeneous,
weakly anisotropic media. It solves two coupled ordinary differential
equations providing zero-order QI amplitudes of qS waves. The 
equations are solved along a ray in a background isotropic medium.
   
   The program WEAKAN uses the file LU1 generated by the program
package ANRAY. WEAKAN generates file LU6, which can be used in 
the program FRESAN and following programs for computation of synthetic 
seismograms.

   The program yields regular results in regions, in which the
program ANRAY, based on the standard ray method for anisotropic media
fails, i.e. in weakly anisotropic media and in singular regions like
cross and kiss singularity. It fails in vicinity of caustics.

                                                     
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' 'LU6'/
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 input data file LU1,
          generated by the program ANRAY.
    'LU6' is the name of the output data file LU6.
    / is a slash recomended in batch and script files to enable future
        extensions.
Defaults:
    'LIN'='anray.dat'
    'LOU'='anray.out'
    'LU1'='lu1.dat'
    'LU6'='lu6.dat'
Example of the main input data:
    'weak.qi' /

    Input data consist partially of  the data  generated by the
program ANRAY, stored in the formatted form in the file LU1, and
partially  of the  additional  input  data  read by list-directed 
input (free format), prepared by the user and stored in  the file 
LIN. Output data describing the computations are stored in the file 
LOU. Output data for computing ray synthetic seismograms in the 
program FRESAN are stored in the file LU6.

                                                      
The data stored at LU1
----------------------

    The data are stored in the file LU1 in the formatted form. For
details see the description of the content of the file LU1 in program
ANRAY.

1)  ICONT,NDST,ILOC                               FORMAT(26I3)
2)  ROS                                           FORMAT(8F10.5)
3)  NPN,NPN,NPN                                   FORMAT(26I3)
4)  APN,APN,APN,APN,APN                           FORMAT(5E15.5)
5)  XPRF,YPRF,0.,PROF                             FORMAT(8F10.5)
6)  (DST(I),I=1,NDST)                             FORMAT(8F10.5)
7)  NTOT,II                                       FORMAT(26I3)
8)  (T(J),X(J),Y(J),Z(J),P1(J),P2(J),P3(J),VP(J),VS(J),RHO(J),
    (E(1,K,J),K=1,3),(E(2,L,J),L=1,3),J = 1,N)    FORMAT(16E15.5)
9)  INDEX                                         FORMAT(26I3)
10) (INDI(I),XCOOR(I),ZCOOR(I),TIME(I),AMPX(1,I),AMPY(1,I),
    AMPZ(1,I),I=1,INDEX)                          FORMAT(I5,9E15.5)
11) IF(INDEX.LT.0) ((AMPX(2,I),AMPY(2,I),AMPZ(2,I),I=1,INDEX)
                                                  FORMAT(6E15.5)
12) (P(K),K=1,3)                                  FORMAT(3E15.5)
13) (POL(K),K=1,3)                                FORMAT(3E15.5)
14) IF(INDEX.LT.0) (POL1(K),K=1,3)                FORMAT(3E15.5)

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

                                                      
The additional input data in the file LIN
-----------------------------------------

1) MTEXT

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

2) INULL,INEWB

   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.
   INEWB...    switch deciding if the QI approach is to be used
         with the standard weak anisotropy matrix B_{mn} or with
         its modified version calculated in the routine NEWB
         INEWB=0: standard QI approach.
         INEWB=1: modified weak anisotropy matrix used.
         Default INEWB=0.

3) MPRINT,NINT

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

   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.

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

6a) ISQRT,IRHO

   ISQRT... controls the approximation of elastic parameters.
         ISQRT=0:  Approximation in elastic parameters. Default 
         value.
         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).
         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 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

   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.

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

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

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

   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)

   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.

8) F(1),F(2),F(3)

   F(I)..I-th component of the force vector.

9) FL,FD,NF

   FL... lower limit of the considered frequency range.
   FD... frequency sampling step.
   NF... number of frequency samples (NF may not exceed 1000).


Example of data LIN for isotropic background for QI computations

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

    The computation terminates when NTOT=0 is found on 7th line of LU1.

                                                      
Output to the file LU6
----------------------

    Data necessary for calculation of the frequency response in the
program FRESAN are stored in the file LU6. The data are stored in LU6 
in the following formatted form:

1) FL,FD,NF                                 FORMAT(2F10.5,I5)

   FL... lower limit of the considered frequency range.
   FD... frequency sampling step.
   NF... number of frequency samples.

2) (E(1,K,NTOT),K=1,3),(E(2,L,NTOT),L=1,3)  FORMAT(16E15.5)

   E(I,J,N)... J-th component of the I-th ray-centered basis vector
         at the termination point of the ray (NTOT).

3) FF,W                                     FORMAT(16E15.5)

   FF... frequency  
   W(I)... (I=1,2) complex-valued I-th ray amplitude in the QI 
         approximation for the frequency FF.

The line 3 repeats NF-times for each ray read from LU1. The lines
2 and 3 repeat for each ray.

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

    The file LOU contains a copy of the input data. Additional data
stored in LOU is information about the model. It is controlled by the 
switch MPRINT.

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.