C 
C Subroutine file 'init.for' to read the input data for the initial
C surface, and to define initial values for complete ray tracing.
C
C Date: 2000, May 19
C Coded by Ludek Klimes
C
C.......................................................................
C
C This file consists of the following subroutines:
C     INIT1...Interface subroutine to INIS1 reading the input data for
C             the four functions specifying the initial conditions for
C             the computed rays.
C             INIT1
C     INIT2...Subroutine evaluating, for given ray take-off parameters,
C             the values of the computed quantities at the initial point
C             of the ray, and storing the important quantities at the
C             initial point of the ray in the common block /INITC/.
C             INIT2
C     INIT5...Subroutine determining whether ray tracing should be
C             repeated for another source point.
C     INIS1...Sample subroutine designed to read the input data for the
C             initial points of rays.  A two-parametric system of rays
C             of each elementary wave is assumed.  A ray of the
C             elementary wave is specified by its two take-off
C             parameters.  The computed rays may start from a single
C             initial point common to all rays, from a curve along which
C             an initial travel time is defined, from an initial surface
C             along which an initial travel time is defined, etc.
C             INIS1
C     INIS2...Sample subroutine returning the functional values and
C             their first and second derivatives, of the functions
C             describing the initial surface.
C             INIS2
C     SQRT3...Subroutine evaluating the square root of the given
C             real-valued positive-semidefinite symmetric 3*3 matrix.
C             SQRT3
C     SPHERE..Subroutine transforming spherical coordinates PAR1, PAR2
C             into the Cartesian coordinates of the corresponding point
C             on the unit sphere.  It also evaluates the first and
C             second derivatives of the Cartesian coordinates with
C             respect to PAR1 and PAR2.
C             SPHERE
C Subroutines INIS1 and INIS2, defining the common initial point,
C initial curve or initial surface, call subroutines VAL1 and VAL2 which
C must be appended.  In addition, subroutines CURVN1 (or its alternative
C CURVB1), CURV2D (or its alternative CURVBD), SURFB1, SURFBD, VAL3B1,
C VAL3BD, VGEN, TERMS, SNHCSH, TRIDEC, TRISOL, DSPLNZ, INTRVL from the
C subroutine package 'FITPACK' by Alan Kaylor Cline, Department of
C Computer Sciences, University of Texas at Austin, are used.  In the
C complete ray tracing, subroutines INIS1 and INIS2 may be replaced by
C any user-defined package containing subroutines INIS1 and INIS2 with
C the same number, type and meaning of their parameters as in this
C file.
C
C.......................................................................
C
C Description of data files:
C                                                     
C Input parameters taken from the input SEP parameter file for the
C Complete Ray Tracing program
C     INIDIM
C     INIPAR
C     ADVANC
C are described in file crtin.for.
C
C                                                     
C Input file SRC containing the coordinates of the initial point:
C     The data are read only if INIDIM.LE.0.
C     If there are several source points, ray tracing is repeated for
C     all source points, repeating all elementary waves specified in
C     data set CODE for each source point.
C     Data in data sets RPAR and
C     WRIT are not repeated, they should be
C     specified for all elementary waves of all source points.
C (1) Several strings terminated by / (a slash).
C             The simplest way is to submit just the /.
C (2) For each source point data (2.1):
C (2.1) 'NAME',X1INI,X2INI,X3INI,TTINI,/
C     'NAME'... String reserved for the name of the initial point.
C             No meaning here, anything in apostrophes may be submitted.
C     X1INI,X2INI,X3INI... Coordinates of a single initial point.
C     TTINI... Initial value of the travel time.
C     Default: X2INI=0,  X3INI=0,  TTINI=0.
C (3) / (a slash) or end of file.
C Example of data set SRC
C
C                                                    
C Input data INIT for the initial points of rays:
C     The data are read only if INIDIM.GT.0.
C     The data are read in by the list directed input (free format).  In
C     the list of input data below, each numbered paragraph indicates
C     the beginning of a new input operation (new READ statement).  If
C     the first letter of the symbolic name of the input variable is
C     I-N, the corresponding value in input data must be of the type
C     INTEGER.  Otherwise (except TEXTI), the input parameter is of the
C     type real.
C (1) TEXTI
C     String describing the data.  Only the first 80 characters of the
C     string are significant.
C (2) Data describing the initial point, curve or surface.
C     For INIDIM=1:
C             (2B) Input data for NFUNC=4 functions X1(.,.),X2(.,.),
C               X3(.,.), TT(.,.) of two variables.  The INIPAR-th one of
C               the 2 variables being simultaneously the ray parameter).
C     For INIDIM=2, INIPAR.LE.0:
C             (2B) Input data for NFUNC=4 functions X1(.,.),X2(.,.),
C               X3(.,.), TT(.,.) of two variables (two initial surface
C               parameters, being simultaneously the ray parameters).
C     For INIDIM=2, INIPAR.GE.1: the following two data sets (2A), (2B):
C             (2A) X1INI,X2INI,X3INI... Coordinates of a given point,
C               see the description of the input data (1).
C             (2B) Input data for NFUNC=2 functions R(.,.),TT(.,.) of
C               two variables (two initial surface parameters, being
C               simultaneously the ray parameters).
C               R(.,.) describes the radius, see input data (1),
C               TT(.,.) is the initial travel time.
C The structure of the input data (2B) is given by the subroutine VAL1
C and is described below.
C Examples of data set INIT:
C 
C Plane wave incident at the model bottom
C 
C Zero-offset rays (exploding reflector)
C
C Above mentioned input data (2B) for the initial curve or for the
C initial surface are read in by the subroutine VAL1 and have the
C following structure:
C     These input data define at least NFUNC individual functions
C     describing the initial conditions.  They are read in by subroutine
C     VAL1 called by INIS1.  The number MFUNC of all functions specified
C     in the input data may be greater or equal to NFUNC.  The data are
C     read in by the list directed input (free format).
C (1) MFUNC
C     The number of all input functions.  It must be greater or equal to
C     the number NFUNC of the functions required to describe the
C     coordinates and travel time along the initial curve or surface.
C     The functions indexed 1 to NFUNC must be the functions describing
C     the coordinates and travel time along the initial curve or
C     surface.
C (2) NFUNC-times (i.e. once for each function) input data (2A)+(2B):
C (2A) TEXTF,IFUNC
C     Identification of the function.
C     TEXTF...Any string. Its first 3 characters must differ from 'END'.
C     IFUNC...Index of the function:
C             1 to 3... coordinates and 4... travel time, or
C             1... radius and 2... travel time.  Amplitude and/or other
C             quantities may follow.
C (2B) 'Input data for one function', see below.
C     Input data for one function
C (3) TEXTE,AUX
C     End of data.
C     TEXTE...String, the first 3 characters of which must be upper-case
C             'END'.
C     AUX...  Any number or a slash.
C Remark:
C     If the initial surface coincides with a structural interface
C     (e.g., exploding-reflector initial conditions), it is reasonable
C     to slightly shift the initial surface from the structural
C     interface into the complex block into which the wave propagates.
C     Otherwise, because of rounding errors, there is a danger that some
C     parts of the initial surface are situated at the opposite side
C     of the interface.
C
C                                                  
C Input data for one function:
C     The data are read in by the list directed input (free format).  In
C     the list of input data below, each numbered paragraph indicates
C     the beginning of a new input operation (new READ statement).  If
C     the first letter of the symbolic name of the input variable is
C     I-N, the corresponding value in input data must be of the type
C     INTEGER.  Otherwise, the input parameter is of the type REAL.
C (1) IVAR1,IVAR2,0,SIGMA
C     The form of the function.
C     IVAR1,IVAR2... Denote the form of the function. The function must
C             be of the form
C               F(G1,G2) = W(A1,A2)-B1-B2  .
C             G1, G2 are the ray parameters. Each of A1, A2, B1, B2 must
C             be either:  (a) one of ray parameters G1, G2,  or (b) must
C             be left out.  At most 2 of parameters A1-B2 may be of kind
C             (a).  Note that IVAR1 controls the type of A1 and B1,
C             IVAR2 controls the type of A2 and B2.
C             For IVAR1.EQ.0: A1, B1 are empty (left out),
C             for IVAR1.EQ.1: A1=G1, B1 is empty,
C             for IVAR1.EQ.2: A1=G2, B1 is empty,
C             for IVAR1.EQ.-1: B1=G1, A1 is empty,
C             for IVAR1.EQ.-2: B1=G2, A1 is empty,
C             the meaning of the parameter IVAR2 is similar.
C             Examples:
C             IVAR1: IVAR2: the form of the function:
C               1      2     F(G1,G2)=W(G1,G2)
C               2      1     F(G1,G2)=W(G2,G1)
C               1      0     F(G1,G2)=W(G1)
C               1     -2     F(G1,G2)=W(G1)-G2
C             Function W is interpolated by means of splines under
C             tension.
C     SIGMA...Is the tension factor (its sign is ignored).  This value
C             indicates the curviness desired.  If ABS(SIGMA) is nearly
C             zero (e.g. 0.001), the resulting surface is approximately
C             the tensor product of cubic splines.  If ABS(SIGMA) is
C             large (e.g. 50.), the resulting surface is approximately
C             tri-linear.  If SIGMA equals zero, tensor products of
C             cubic splines result.  A recommended value for SIGMA is
C             approximately 1. In absolute value.
C (2) NX(1),...,NX(NVAR)
C     The numbers of grid coordinates for the interpolation.
C     This input is performed if at least one of IVAR1, IVAR2 is
C     positive.
C     Each of NX(1),...,NX(NVAR) corresponds to one positive value of
C     IVAR1, IVAR2 and specifies the number of grid coordinates
C     corresponding to that independent variable of function W, see (1).
C     The sign of NX(1),...,NX(NVAR) is ignored. NVAR (.LE.2) is the
C     number of positive values of the above quantities IVAR1, IVAR2,
C     i.e. the number of independent variables of function W, see (1).
C (3) X1(1),...,X1(NX(1))
C     The grid coordinates corresponding to the first independent
C     variable of function W, see (1).
C     This input is performed if NX(1) is specified, see (2), and is not
C     zero. The grid coordinates may be specified in any order.
C (4) X2(1),...,X2(NX(2))
C     The grid coordinates corresponding to the second independent
C     variable of function W, see (1).
C     This input is performed if NX(2) is specified, see (2), and is not
C     zero. The grid coordinates may be specified in any order.
C (5) (((W(I,J),I=1,MAX(NX(1),1)),J=1,MAX(NX(2),1))
C     The values of function W at grid points. Function value W(I,J)
C     corresponds to point (X1(I),X2(J)).
C
C.......................................................................
C
C Storage in the memory:
C     Input data INIT (2A) are stored in the common block
C     /INISC/.  The important quantities at the initial point of the ray
C     (see C.R.T.6.1) are stored in the common block /INITC/.  These
C     common blocks are defined in the include file 'initd.inc'.
C     initd.inc
C
C=======================================================================
C
C     
C
      SUBROUTINE INIT1(LUN)
      INTEGER LUN
C
C Subroutine INIT1 calls the subroutine INIS1 to read the input data for
C the initial points of rays.
C
C Input:
C     LUN...  Logical unit number of the external input device
C             containing the input data.
C The input parameter is not altered.
C
C No output.
C
C Common block /INITC/:
      INCLUDE 'initc.inc'
C     initc.inc
C None of the storage locations of the common block, except ICB1I, are
C altered.  ICB1I is set to zero.
C
C Subroutines and external functions required:
      EXTERNAL INIS1
C     INIS1... This file.
C
C Date: 1993, December 18
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER NFUNC
      PARAMETER (NFUNC=4)
C
      CALL INIS1(LUN,NFUNC)
      ICB1I=0
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE INIT2(PAR1,PAR2,YL,Y,YY,IY,IEND,IWAVE0,IKODE)
      REAL PAR1,PAR2,YL(6),Y(35),YY(5)
      INTEGER IY(12),IEND,IWAVE0,IKODE
C
C Subroutine INIT2 evaluates, for given ray take-off parameters, the
C values of the computed quantities at the initial point of the ray, and
C stores the important quantities at the initial point of the ray in the
C common block /INITC/.
C
C Input:
C     PAR1,PAR2... Ray take-off parameters.
C     YL,Y,YY,IY... Arrays of dimensions at least 6, 35, 5, 12,
C             respectively.
C     IWAVE0..Index of the already computed elementary wave having the
C             most numerous common elements with the current elementary
C             wave.  Need not be defined if IKODE=0.
C     IKODE...The length of the common part of the codes of the IWAVE-th
C             and IWAVE0-th elementary waves.
C The input parameters PAR1, PAR2 are not altered.
C
C Output.
C     YL...   Array containing local quantities at the initial point of
C             the ray (see C.R.T.5.5.4).  The quantities are listed in
C             the subroutine file 'ray.for'.
C     Y...    Array containing basic quantities computed along the ray
C             at the initial point of the ray (see C.R.T.5.2.1).  The
C             quantities are listed in the subroutine file 'ray.for'.
C     YY...   Array containing real auxiliary quantities computed along
C             the ray at the initial point of the ray (see C.R.T.5.2.2).
C             The quantities are listed in the subroutine file
C             'ray.for'.
C     IY...   Array containing integer auxiliary quantities computed
C             along the ray at the initial point of the ray (see
C             C.R.T.5.2.2).  The quantities are listed in the subroutine
C             file 'ray.for'.
C     IEND... Information on the initial point of the ray:
C             IEND
C             0...  Computation of the ray may follow.
C             71... There is no ray corresponding to the given ray
C               parameters.  E.g., the given parameters do not belong to
C               the domain of the initial surface.
C             72... Initial point of the ray is not situated in the
C               required complex block.
C             73... Initial point of the ray is not situated in the
C               computational volume.
C             74... Ray of the generated wave is not real-valued.
C
C Common block /DCRT/ (see subroutine file 'ray.for'):
      INCLUDE 'dcrt.inc'
C     dcrt.inc
C None of the storage locations of the common block are altered.
C
C Common block /INISC/:
      INCLUDE 'initd.inc'
C     initd.inc
C None of the storage locations of the common block are altered.
C
C Common block /INITC/:
      INCLUDE 'initc.inc'
C     initc.inc
C All the storage locations of the common block are defined in this
C subroutine.  ICB1I must be zero before the first invocation of this
C subroutine, other storage locations may be undefined.
C
C Subroutines and external functions required:
      EXTERNAL NSRFC,BLOCK,VELOC,KOOR,METRIC,PARM2,SMVPRD,CODE,INIS2
      INTEGER NSRFC,KOOR
C     NSRFC,BLOCK,VELOC... File 'model.for' of the package 'MODEL'.
C     KOOR,METRIC... File 'metric.for' of the package 'MODEL'.
C     PARM2... File 'parm.for' of the package 'MODEL'.
C     SMVPRD... File 'means.for' of the package 'MODEL'.
C     CODE... File 'code.for'.
C     INIS2... This file.
C
C Date: 2000, May 19
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations for local model parameters:  FAUX(10),
C       G(12),GAMMA(18),GSQRD,  UP(10),US(10),RO,QP,QS, VP,VS,VD(10),QL:
      INCLUDE 'auxmod.inc'
C     auxmod.inc
C
C.......................................................................
C
C     Auxiliary storage locations:
      INTEGER KODIND,ICBNEW,ISB2,ICB2,I,NY,NFUNC
      PARAMETER (NFUNC=4)
      REAL E11,F21,F31,C11,B11,FF11,FFE11,CB11,ECB11,U1,T1,BT1
      REAL E12,F22,F32,C12,B12,FF12,FFE12,CB12,ECB12,U2,T2,BT2
      REAL E22,F23,F33,C22,B22,FF21,FFE21,CB21,ECB21
      REAL E(3),F(6,NFUNC),    FF22,FFE22,CB22,ECB22
      EQUIVALENCE (E(1),E11),(E(2),E12),(E(3),E22)
      REAL V0,DETC,TRECE,Z1,Z2,Z3,AUX1,AUX2,AUX3
      SAVE V0
C
C     KODIND..Position in the code of elementary wave.
C     ICBNEW..The index of the complex block in which the initial point
C             of the ray should be situated, supplemented by the sign
C             '+' for P wave or '-' for S wave.  Output from subroutine
C             CODE.
C     ISB2... Index of the simple block in which the initial point is
C             situated.
C     ICB2... Index of the complex block in which the initial point is
C             situated.
C     I...    Auxiliary and loop variable.
C     NY...   Number of defined locations of array Y.
C     INIDIM..Distinguishes between initial point, line, or surface.
C     E=(E11,E12,E22)... Projection matrix, see the subroutine INIS2.
C     F...Functions describing the initial surface and their
C             derivatives, see the subroutine INIS2.
C     F21,F22,F23, F31,F32,F33... Covariant components of the vectors
C             F(2,1),F(2,2),F(2,3), F(3,1),F(3,2),F(3,3) tangent to the
C             initial surface.
C     C11,C12,C22... Components of matrix C of space or space-time
C             scalar products, eventually of its square root.
C     B11,B12,B22... Components of inverse square root of matrix C.
C     FF11,FF12,FF21,FF22... Components of the auxiliary matrix FF.
C     FFE11,FFE21,FFE12,FFE22... Components of the matrix FF*E.
C     CB11,CB21,CB12,CB22... Components of the matrix 1-C*B.
C     ECB11,ECB21,ECB12,ECB22... Components of the matrix E*CB/Tr(E*C).
C     U1,U2... Slowness derivatives with respect to PAR1,PAR2.
C     T1,T2... Velocity*derivatives of the travel time along the initial
C              surface with respect to PAR1,PAR2.
C     BT1,BT2... Components of the vector (B+ECB)*T.
C     V0...   Propagation velocity at the initial point.
C     DETC... Determinant of matrix C, eventually of its square root.
C     TRECE... Trace of the matrix E*C*E.
C     Z1,Z2,Z3... Unit normal to the initial surface - covariant comp.
C     AUX1,AUX2,AUX3... Auxiliary storage locations.
C
C.......................................................................
C
      IY(11)=0
      IF(MODCRT.GE.3) THEN
C       613
        CALL ERROR('613 in INIT2: This program branch is not coded')
C       Interpolation mode of the complete ray tracing program has
C       not been enabled yet.  See 'ray.for', description of input
C       data, (2), MODCRT.
      END IF
C
      CALL INIS2(NFUNC,PAR1,PAR2,E,F,IEND)
      IF(IEND.NE.0) THEN
        RETURN
      END IF
      IF(F(1,1).LT.BOUNDR(1).OR.BOUNDR(2).LT.F(1,1).OR.
     *   F(1,2).LT.BOUNDR(3).OR.BOUNDR(4).LT.F(1,2).OR.
     *   F(1,3).LT.BOUNDR(5).OR.BOUNDR(6).LT.F(1,3).OR.
     *                          BOUNDR(7).LT.F(1,4)) THEN
        IEND=73
        RETURN
      END IF
      IF(E11.EQ.0..AND.E12.EQ.0..AND.E22.EQ.0.) THEN
C       Initial point:
        IF(INIDIM.GT.0) THEN
C         617
          CALL ERROR('617 in INIT2: Strange common initial point')
C         Initial conditions determined by subroutine INIS2 do not
C         resemble common initial point.  Contact the authors.
        END IF
      ELSE IF(E11*E22-E12*E12.LE.0.001) THEN
C       Initial line:
        IF(INIDIM.NE.1) THEN
C         618
          CALL ERROR('618 in INIT2: Strange initial line')
C         Initial conditions determined by subroutine INIS2 do not
C         resemble initial line.  Contact the authors.
        END IF
      ELSE
C       Initial surface:
        IF(INIDIM.NE.2) THEN
C         619
          CALL ERROR('619 in INIT2: Strange initial surface')
C         Initial conditions determined by subroutine INIS2 do not
C         resemble initial surface.  Contact the authors.
        END IF
      END IF
C
C     Initial travel time
      YI(1)=F(1,4)
C     Initial imaginary travel time
      YI(2)=0.
C
C     Coordinates of the initial point of the ray
      YI3=YI(3)
      YI4=YI(4)
      YI5=YI(5)
      YI(3)=F(1,1)
      YI(4)=F(1,2)
      YI(5)=F(1,3)
C
C     Local coordinate system:
C     Covariant components of vectors tangent to the initial surface
      IF(KOOR().NE.0) THEN
        CALL METRIC(YI(3),GSQRD,G,GAMMA)
        CALL SMVPRD(G,F(2,1),F(2,2),F(2,3),F21,F22,F23)
        CALL SMVPRD(G,F(3,1),F(3,2),F(3,3),F31,F32,F33)
      ELSE
        GSQRD=1.
        F21=F(2,1)
        F22=F(2,2)
        F23=F(2,3)
        F31=F(3,1)
        F32=F(3,2)
        F33=F(3,3)
      END IF
C     Scalar products of vectors tangent to the initial surface
      C11=F(2,1)*F21+F(2,2)*F22+F(2,3)*F23
      C12=F(2,1)*F31+F(2,2)*F32+F(2,3)*F33
      C22=F(3,1)*F31+F(3,2)*F32+F(3,3)*F33
      DETC=C11*C22-C12*C12
C     Unit normal to the initial surface - covariant components
      AUX1=GSQRD/SQRT(DETC)
      Z1=(F(2,2)*F(3,3)-F(2,3)*F(3,2))*AUX1
      Z2=(F(2,3)*F(3,1)-F(2,1)*F(3,3))*AUX1
      Z3=(F(2,1)*F(3,2)-F(2,2)*F(3,1))*AUX1
C
C     Modification of the coordinates of the initial point of the ray:
      IF(ADVANC.NE.0.) THEN
C       Shifting the initial point in the direction of the
C       unit normal to the initial surface - contravariant components
        IF(KOOR().NE.0) THEN
          YI(3)=F(1,1)+ADVANC*(G(1)*Z1+G(2)*Z2+G(4)*Z3)
          YI(4)=F(1,2)+ADVANC*(G(2)*Z1+G(3)*Z2+G(5)*Z3)
          YI(5)=F(1,3)+ADVANC*(G(4)*Z1+G(5)*Z2+G(6)*Z3)
        ELSE
          YI(3)=F(1,1)+ADVANC*Z1
          YI(4)=F(1,2)+ADVANC*Z2
          YI(5)=F(1,3)+ADVANC*Z3
        END IF
      END IF
C
C     Material parameters (defining ISB1I, ICB1I, YL(1) to YL(6)):
      IF(ICB1I.NE.0) THEN
        IF(YI(3).NE.YI3.OR.YI(4).NE.YI4.OR.YI(5).NE.YI5) THEN
          ICB1I=0
        END IF
      END IF
      IF(ICB1I.EQ.0) THEN
        CALL BLOCK(YI(3),0,0,I,ISB1I,ICB2)
      ELSE
        ICB2=IABS(ICB1I)
      ENDIF
C     Defining locations of the array FSRFCA of the common /INITC/:
      DO 11 I=1,NSRFCA
        CALL SRFC2(NSRFC()+I,YI(3),FAUX)
        FSRFCA(I)=FAUX(1)
   11 CONTINUE
      IF(ICB2.EQ.0) THEN
        IEND=72
        RETURN
      ENDIF
      IY(2)=0
      IY(6)=0
      IY(8)=ICB2
      CALL CODE(IY,KODIND,ICBNEW,IEND)
      IF(IEND.EQ.22) THEN
        IEND=72
        RETURN
      ELSE IF(IEND.NE.0) THEN
C       611
        WRITE(*,'(A,I2)') ' Subroutine CODE reports IEND=',IEND
        CALL ERROR('611 in INIT2: Wrong function of subroutine CODE')
C       Other reason of the termination of the ray computation
C       than 22 should not be reported by the subroutine CODE when
C       referenced by the subroutine INIT2.  Contact the authors.
      END IF
      IF(ICB2.NE.IABS(ICBNEW)) THEN
C       612
        CALL ERROR('612 in INIT2: Wrong function of subroutine CODE')
C       Subroutine CODE requires the first ray element to be
C       situated in another complex block than the initial point.
C       This error should not appear.  Contact the authors.
      END IF
      IF(ICB1I.EQ.0.OR.ICB1I.NE.ICBNEW) THEN
        ICB1I=ICBNEW
        CALL PARM2(IABS(ICBNEW),YI(3),UP,US,YLI(3),QP,QS)
        CALL VELOC(ICBNEW,UP,US,QP,QS,YLI(1),YLI(2),VD,QL)
        V0=VD(1)
        YLI(4)=VD(2)
        YLI(5)=VD(3)
        YLI(6)=VD(4)
      ENDIF
C
C     Important quantities at the initial point of the ray (C.R.T.6.1):
C     Slowness derivatives with respect to ray parameters
      AUX1=-(YLI(4)*F(2,1)+YLI(5)*F(2,2)+YLI(6)*F(2,3))/(V0*V0)
      AUX2=-(YLI(4)*F(3,1)+YLI(5)*F(3,2)+YLI(6)*F(3,3))/(V0*V0)
      U1=E11*AUX1+E12*AUX2
      U2=E12*AUX1+E22*AUX2
C     Slowness vector
      AUX1=( C22*F(2,4)-C12*F(3,4))/DETC
      AUX2=(-C12*F(2,4)+C11*F(3,4))/DETC
      AUX3=V0**(-2)-AUX1*F(2,4)-AUX2*F(3,4)
      IF(AUX3.LE.0.) THEN
C       Evanescent wave
        IEND=74
        RETURN
      END IF
      AUX3=SQRT(AUX3)
      YI(6)=F21*AUX1+F31*AUX2+Z1*AUX3
      YI(7)=F22*AUX1+F33*AUX2+Z2*AUX3
      YI(8)=F23*AUX1+F33*AUX2+Z3*AUX3
C     Space-time scalar products of vectors tangent to the surface
      T1=F(2,4)*V0
      T2=F(3,4)*V0
      C11=C11-T1*T1
      C12=C12-T1*T2
      C22=C22-T2*T2
      DETC=SQRT(C11*C22-C12*C12)
      IF(INIDIM.NE.1) THEN
C       Initial surface or initial point:
C       Square root of the matrix C
        AUX1=SQRT(C11+C22+DETC+DETC)
        C11=(C11+DETC)/AUX1
        C12=C12/AUX1
        C22=(C22+DETC)/AUX1
C       Inverse square root of the matrix C
        B11= C22/DETC
        B12=-C12/DETC
        B22= C11/DETC
C       First basis vector of ray-centred coordinate system
        AUX3=V0*(B11*T1+B12*T2)
        YI(9) =F21*B11+F31*B12-YI(6)*AUX3
        YI(10)=F22*B11+F32*B12-YI(7)*AUX3
        YI(11)=F23*B11+F33*B12-YI(8)*AUX3
C       Geometrical spreading matrix Q
        YI(12)=C11*E11+C12*E12
        YI(13)=C12*E11+C22*E12
        YI(16)=C11*E12+C12*E22
        YI(17)=C12*E12+C22*E22
C       Matrix P
        FF11=F(4,4)-YI(6)*F(4,1)-YI(7)*F(4,2)-YI(8)*F(4,3)-T1*U1
        FF12=F(5,4)-YI(6)*F(5,1)-YI(7)*F(5,2)-YI(8)*F(5,3)
        FF21=FF12-T2*U1
        FF12=FF12-T1*U2
        FF22=F(6,4)-YI(6)*F(6,1)-YI(7)*F(6,2)-YI(8)*F(6,3)-T2*U2
        YI(14)=B11*FF11+B12*FF21
        YI(15)=B12*FF11+B22*FF21
        YI(18)=B11*FF12+B12*FF22
        YI(19)=B12*FF12+B22*FF22
      ELSE
C       Initial line:
C       Infinite part of the inverse square root of the matrix C
        B11=(1.-E11)/DETC
        B12=   -E12 /DETC
        B22=(1.-E22)/DETC
C       Matrix CB=1-C*B
        CB11=1.-C11*B11+C12*B12
        CB21=  -C12*B11+C22*B12
        CB12=  -C11*B12+C12*B22
        CB22=1.-C12*B12+C22*B22
C       E-projection of the finite part of the inverse square root of C
        TRECE=SQRT(E11*C11+2.*E12*C12+E22*C22)
        ECB11=(E11*CB11+E12*CB21)/TRECE
        ECB21=(E12*CB11+E22*CB21)/TRECE
        ECB12=(E11*CB12+E12*CB22)/TRECE
        ECB22=(E12*CB12+E22*CB22)/TRECE
C       First basis vector of ray-centred coordinate system
        AUX1=B11+ECB11
        AUX2=B12+ECB12
        BT1=AUX1*T1+AUX2*T2
        BT2=(B12+ECB21)*T1+(B22+ECB22)*T2
        AUX3=V0*BT1
        YI(9) =F21*AUX1+F31*AUX2-YI(6)*AUX3
        YI(10)=F22*AUX1+F32*AUX2-YI(7)*AUX3
        YI(11)=F23*AUX1+F33*AUX2-YI(8)*AUX3
C       Geometrical spreading matrix Q
        YI(12)=E11*TRECE
        YI(13)=E12*TRECE
        YI(16)=E12*TRECE
        YI(17)=E22*TRECE
C       Matrix P
        FF11=F(4,4)-YI(6)*F(4,1)-YI(7)*F(4,2)-YI(8)*F(4,3)
        FF12=F(5,4)-YI(6)*F(5,1)-YI(7)*F(5,2)-YI(8)*F(5,3)
        FF22=F(6,4)-YI(6)*F(6,1)-YI(7)*F(6,2)-YI(8)*F(6,3)
        FFE11=FF11*E11+FF12*E12
        FFE21=FF12*E11+FF22*E12
        FFE12=FF11*E12+FF12*E22
        FFE22=FF12*E12+FF22*E22
        YI(14)=B11*FF11+B12*FF12+ECB11*FFE11+ECB12*FFE21-BT1*U1
        YI(15)=B12*FF11+B22*FF12+ECB12*FFE11+ECB22*FFE21-BT2*U1
        YI(18)=B11*FF12+B12*FF22+ECB11*FFE12+ECB12*FFE22-BT1*U2
        YI(19)=B12*FF12+B22*FF22+ECB12*FFE12+ECB22*FFE22-BT2*U2
      END IF
C     Take-off parameters
      YI(20)=PAR1
      YI(21)=PAR2
C
C     Modification of the initial travel time:
      IF(ADVANC.NE.0.) THEN
        YI(1)=YI(1)+(YI(3)-F(1,1))*YI(6)
     *             +(YI(4)-F(1,2))*YI(7)
     *             +(YI(5)-F(1,3))*YI(8)
      END IF
C
C
C     Initial values for the complete ray tracing (C.R.T.6.2):
      DO 21 I=1,6
        YL(I)=YLI(I)
   21 CONTINUE
      DO 22 I=1,11
        Y(I)=YI(I)
   22 CONTINUE
      IF(ICB1I.GE.0) THEN
        NY=27+2
      ELSE
        NY=27+8
      ENDIF
      DO 23 I=12,NY
        Y(I)=0.0
   23 CONTINUE
      Y(12)=1.0
      Y(17)=1.0
      Y(22)=1.0
      Y(27)=1.0
      Y(28)=1.0
      IF(NY.GE.34) Y(34)=1.0
      YY(1)=0.0
      YY(2)=UEB
      YY(3)=0.0
      YY(4)=0.0
      YY(5)=0.0
      IY(1)=NY
      IY(2)=KODIND
      IY(3)=0
      IY(4)=ISB1I
      IY(5)=ICB1I
      IY(6)=0
C     Note: IY(7),IY(8) may be undefined
      IY(7)=0
      IY(8)=0
      IY(9)=0
      IY(10)=0
C     IY(11) is initiallized in the beginning of this subroutine.
      IY(12)=0
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE INIT5(IRAY)
      INTEGER IRAY
C
C Subroutine INIT5 determines whether ray tracing should be repeated for
C another source point (common initial point of rays).
C
C No input.
C
C Output:
C     IRAY... IRAY=0: There is another source point.  The source point
C               has replaced the previous source point in locations
C               X1INI, X2INI, X3INI, TTINI of common block /INISC/
C             IRAY=-1: There is not another source point.
C
C Common block /INISC/:
      INCLUDE 'initd.inc'
C     initd.inc
C
C No subroutines and external functions required.
C
C Date: 1999, September 23
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I
C
      IF(NSRC.GT.0) THEN
C       There is another source point
        X1INI=X1SRC(1)
        X2INI=X2SRC(1)
        X3INI=X3SRC(1)
        TTINI=TTSRC(1)
        NSRC=NSRC-1
        DO 10 I=1,NSRC
          X1SRC(I)=X1SRC(I+1)
          X2SRC(I)=X2SRC(I+1)
          X3SRC(I)=X3SRC(I+1)
          TTSRC(I)=TTSRC(I+1)
   10   CONTINUE
        IRAY=0
      ELSE
        IRAY=-1
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE INIS1(LUN,NFUNC)
      INTEGER LUN,NFUNC
C
C Subroutine INIS1 reads the input data for the initial points of rays
C and stores them in common block /INISC/, and if required, calls the
C subroutine VAL1 to read the input data for the interpolated functions
C of two variables (ray parameters), to determine the coefficients
C necessary to compute an interpolatory function on a two dimensional
C rectangular grid, and to store them in the memory.  The functions
C determined can be represented as a tensor product of splines under
C tension.  For actual mapping of points it is necessary to call the
C subroutine INIS2, which also returns the first and second partial
C derivatives.
C
C Input:
C     LUN...  Logical unit number of the external input device
C             containing the input data.
C             May be zero for a point source.
C     NFUNC...Number of the functions required to be defined during the
C             current invocation of INIS1.  Since the functions
C             specified in the input data do not coincide with the
C             required functions but are transformed to them, NFUNC need
C             not equal the number of functions specified in the input
C             data.
C None of the input parameters are altered.
C
C No output.
C
C Common block /INISC/:
      INCLUDE 'initd.inc'
C     initd.inc
C All the storage locations of the common block are defined in this
C subroutine.
C
C Subroutines and external functions required:
      EXTERNAL RSEP3I,RSEP3R,RSEP3T,VAL1
C     VAL1, SORTV, READV... File 'val.for' of the package 'MODEL'.
C     CURVN1 or CURVB1 (alternatives), SURFB1, VAL3B1, SNHCSH, VGEN,
C              TERMS, TRIDEC, TRISOL... Subroutine package 'FITPACK'
C              (file 'fit.for' of the package 'MODEL').
C
C Date: 1999, September 23
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER LUSRC
      PARAMETER (LUSRC=9)
      REAL UNDEF
      PARAMETER (UNDEF=-999999.)
      INTEGER LFUNC,MFUNC,I
      CHARACTER*80 TEXTI
      CHARACTER*3 TFUNCT
      DATA TFUNCT/'   '/
C
C     LUSRC...Logical unit number of the source-point file.  The file is
C             opened and closed here.
C     LFUNC...Difference between the number NFUNC of the required
C             functions and the number of input functions specifying
C             them.
C     MFUNC...Number of functions specified in the input data.
C     TEXTI...String of 80 characters for various purposes.
C
C     Quantities defining the kind of initial conditions
      CALL RSEP3I('INIDIM',INIDIM,0)
      CALL RSEP3I('INIPAR',INIPAR,2)
      CALL RSEP3R('ADVANC',ADVANC,0.)
C
      NSRC=0
      IF(INIDIM.LE.0) THEN
C       Point source
        CALL RSEP3T('SRC',TEXTI,'src.dat')
        OPEN(LUSRC,FILE=TEXTI,STATUS='OLD')
        READ(LUSRC,*) (TEXTI,I=1,20)
        X2INI=0.
        X3INI=0.
        TTINI=0.
        READ(LUSRC,*) TEXTI,X1INI,X2INI,X3INI,TTINI
        DO 18 I=1,MSRC
          X1SRC(I)=UNDEF
          X2SRC(I)=0.
          X3SRC(I)=0.
          TTSRC(I)=0.
          READ(LUSRC,*,END=19) TEXTI,X1SRC(I),X2SRC(I),X3SRC(I),TTSRC(I)
          IF(X1SRC(I).EQ.UNDEF) THEN
            GO TO 19
          END IF
   18   CONTINUE
C         604
          CALL ERROR('604 in INIS1: Too many source points')
C         Dimension MSRC of arrays X1SRC, X2SRC, X3SRC, TTSRC in file
C         initd.inc is smaller than the number
C         of source points.
   19   CONTINUE
        NSRC=I-1
        CLOSE(LUSRC)
      ELSE
        IF(LUN.LE.0) THEN
C         603
          CALL ERROR('603 in INIS1: No input device')
C         For INIDIM.GT.0, there must be specified input file INIT,
C         see the main input data for CRT.
        END IF
C
C       (1) Reading the name of the input data
        READ(LUN,*) TEXTI
C
C       (3) Data describing the initial point, curve or surface.
        IF(INIDIM.EQ.2.AND.INIPAR.GT.0) THEN
          READ(LUN,*) X1INI,X2INI,X3INI
          LFUNC=2
        ELSE
          LFUNC=0
          IF(INIDIM.EQ.1.AND.(INIPAR.LT.1.OR.INIPAR.GT.2)) THEN
C           602
            CALL ERROR('602 in INIS1: Wrong value of INIPAR')
C           For INIDIM=1, there must be INIPAR=1 or INIPAR=2.
          END IF
        END IF
        READ(LUN,*) MFUNC
        IF(NFUNC-LFUNC.LE.MFUNC) THEN
          CALL VAL1(LUN,3,MFUNC,1,TFUNCT)
        ELSE
C         601
          CALL ERROR('601 in INIS1: Small number of input functions')
C         The number of input functions is less than the number of
C         functions necessary to describe coordinates and travel
C         time along the initial surface.
        END IF
C
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE INIS2(NFUNC,PAR1,PAR2,E,F,IEND)
      INTEGER NFUNC,IEND
      REAL PAR1,PAR2,E(3),F(6,NFUNC)
C
C Subroutine INIS2 evaluates the functional values and the derivatives
C of the functions describing the initial surface.  The first three
C functions of given ray parameters are coordinates of the point
C corresponding to the given ray parameters, the fourth function is the
C initial value of the travel time.  The single initial point common to
C all rays or the initial line are treated as singular limiting cases of
C the initial surface.  The input data specifying the functions must
C have been read by the subroutine INIS1.
C
C Input:
C     NFUNC...Number of functions required.  It is assumed to be 4 in
C             this version (three coordinates and the travel time).
C     PAR1,PAR2... Ray parameters.
C     E...    Array of the dimension at least 3.
C     F...    Array of the dimension at least 6*NFUNC.
C None of the input parameters except E and F are altered.
C
C Output:
C     E...    Array containing the components E11, E12, E22 of the 2*2
C             symmetric projection matrix onto the tangent space to the
C             ray parameter's manifold.  For a non-degenerate initial
C             surface, E is the identity matrix.  For a single initial
C             point, E is the zero matrix.  For the initial line, E is
C             the projection matrix of the rank 1.  Note that a
C             projection matrix E satisfies the relation E*E=E.
C     F(1:6,I)... For a non-degenerate initial surface, the value and
C             the first and second partial derivatives F, F1, F2, F11,
C             F12, F22 of the I-th function F(PAR1,PAR2).  Note that
C             F1 = E11,E12 * F1
C             F2   E12,E22   F2 ,
C             and
C             F11,F12 = E11,E12 * F11,F12 * E11,E12
C             F12,F22   E12,E22   F12,F22   E12,E22 .
C             Thus, in a degenerate case (i.e. if E is not the identity
C             matrix), the first derivatives are modified in the
C             following way,
C             F1 = F1 + F31 - E11,E12 * F31
C             F2   F2   F32   E12,E22   F32 ,
C             and second derivatives are modified as follows:
C             F11,F12 = F11,F12 + F311,F312 - E11,E12*F311,F312*E11,E12
C             F12,F22   F12,F22   F312,F322   E12,E22 F312,F322 E12,E22.
C             Here F31, F32, F311, F312 and F322 are the derivatives of
C             F1, F2, F11, F12 and F22 with respect to the small
C             parameter (e.g. a radius) which shrinks to zero upon an
C             initial line or at a single initial point.
C     IEND... Information on the initial point of the ray:
C             0...  Computation of the ray may follow.
C             71... There is no ray corresponding to the given ray
C               parameters.  E.g., the given parameters do not belong to
C               the domain of the initial surface.
C
C Common block /INISC/:
      INCLUDE 'initd.inc'
C     initd.inc
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
      EXTERNAL METRIC,VAL2,SPHERE,SQRT3
C     METRIC..File 'metric.for' of the package 'MODEL'.
C     VAL2... File 'val.for' of the package 'MODEL'.
C     CURV2D OR CURVBD (alternatives), SURFBD, VAL3BD, SNHCSH, DSPLNZ,
C             INTRVL... Subroutine package 'FITPACK' (file 'fit.for' of
C             the package 'MODEL').
C     SPHERE,SQRT3... This file.
C
C Date: 1996, May 9
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I,I1,I2,I11,I22,LFUNC
      REAL AUX1,AUX2,DUMMY,R(3),G(12),FF(6,3)
C
C     I...    Auxiliary loop variable.
C     I1,I11..Array subscripts of the first and second derivatives with
C             respect to INIPAR-th ray parameter.
C     I2,I22..Array subscripts of the first and second derivatives with
C             respect to the other than INIPAR-th ray parameter.
C     LFUNC...Difference between the numbers of required (NFUNC) and
C             interpolated functions
C     AUX1,AUX2... Auxiliary storage locations.
C     DUMMY...Dummy storage location.
C     R...    Array used for general coordinates or ray parameters.
C     G...    G(1)-G(6)... Covariant components of the symmetric 3*3
C               metric tensor, or contravariant components of the
C               symmetric 3*3 matrix of the basis vectors of the local
C               Cartesian coordinate system (i.e. the square root of the
C               contravariant metric tensor).
C             G(7)-G(12)... Contravariant components of the symmetric
C               3*3 metric tensor.
C     FF...   General-purpose auxiliary array.  Used to store local
C             Cartesian coordinates and their derivatives with respect
C             to the ray parameters.  Temporarily used also as dummy
C             storage location for Christoffel symbols, for the last
C             interpolated function, e.t.c.
C
C.......................................................................
C
      IEND=0
      IF(INIDIM.LE.1.OR.INIPAR.GE.1) THEN
C       Matrix E of local Cartesian coordinate system basis vectors
        R(1)=X1INI
        R(2)=X2INI
        R(3)=X3INI
        CALL METRIC(R,DUMMY,G,FF)
      END IF
      IF(INIDIM.LE.0) THEN
C       Single initial point:
C       Projection matrix
        E(1)=0.
        E(2)=0.
        E(3)=0.
C       Contravariant components of the symmetric 3*3 matrix of the
C       basis vectors of the local Cartesian coordinate system
        CALL SQRT3(G(7),G)
C       Mapping of the ray parameters onto a unit sphere
        CALL SPHERE(INIPAR,PAR1,PAR2,FF,IEND)
        IF(IEND.NE.0) THEN
          RETURN
        END IF
C       Required functions (3 general coordinates and a travel time)
        F(1,1)=X1INI
        F(1,2)=X2INI
        F(1,3)=X3INI
        F(1,4)=TTINI
        DO 14 I=2,6
          F(I,1)=G(1)*FF(I,1)+G(2)*FF(I,2)+G(4)*FF(I,3)
          F(I,2)=G(2)*FF(I,1)+G(3)*FF(I,2)+G(5)*FF(I,3)
          F(I,3)=G(4)*FF(I,1)+G(5)*FF(I,2)+G(6)*FF(I,3)
          F(I,4)=0.
   14   CONTINUE
      ELSE
C       Initial line or initial surface:
C       Interpolated functions
        R(1)=PAR1
        R(2)=PAR2
        R(3)=0.
        IF(INIDIM.EQ.1) THEN
          R(3-INIPAR)=0.
        END IF
        IF(INIDIM.EQ.2.AND.INIPAR.GT.0) THEN
          LFUNC=2
        ELSE
          LFUNC=0
        END IF
        DO 21 I=LFUNC+1,NFUNC-1
          CALL VAL2(3,I-LFUNC,1,R,F(1,I),DUMMY)
          F(4,I)=F(5,I)
          F(5,I)=F(6,I)
          F(6,I)=F(1,I+1)
   21   CONTINUE
        CALL VAL2(3,NFUNC-LFUNC,1,R,FF,DUMMY)
        F(1,NFUNC)=FF(1,1)
        F(2,NFUNC)=FF(2,1)
        F(3,NFUNC)=FF(3,1)
        F(4,NFUNC)=FF(5,1)
        F(5,NFUNC)=FF(6,1)
        F(6,NFUNC)=FF(1,2)
        IF(INIDIM.EQ.1) THEN
C         Initial line:
C         Covariant components of the vector tangent to the initial line
          I1=1+INIPAR
          FF(6,1)=G(1)*F(I1,1)+G(2)*F(I1,2)+G(4)*F(I1,3)
          FF(6,2)=G(2)*F(I1,1)+G(3)*F(I1,2)+G(5)*F(I1,3)
          FF(6,3)=G(4)*F(I1,1)+G(5)*F(I1,2)+G(6)*F(I1,3)
C         Contravariant unit vector tangent to the initial line
          AUX2=F(I1,1)*FF(6,1)+F(I1,2)*FF(6,2)+F(I1,3)*FF(6,3)
          AUX1=SQRT(AUX2)
          FF(1,1)=F(I1,1)/AUX1
          FF(1,2)=F(I1,2)/AUX1
          FF(1,3)=F(I1,3)/AUX1
C         Derivative of the unit vector tangent to the initial line
          I11=2*I1
          AUX2=(F(I11,1)*FF(6,1)+F(I11,2)*FF(6,2)+F(I11,3)*FF(6,3))/AUX2
          FF(2,1)=(F(I11,1)-FF(1,1)*AUX2)/AUX1
          FF(2,2)=(F(I11,2)-FF(1,2)*AUX2)/AUX1
          FF(2,3)=(F(I11,3)-FF(1,3)*AUX2)/AUX1
C         Covariant vector normal to the interpolated surface
          I2=5-I1
          FF(3,1)=FF(1,2)*F(I2,3)-FF(1,3)*F(I2,2)
          FF(3,2)=FF(1,3)*F(I2,1)-FF(1,1)*F(I2,3)
          FF(3,3)=FF(1,1)*F(I2,2)-FF(1,2)*F(I2,1)
C         Derivative of the vector normal to the interpolated surface
          FF(4,1)=FF(2,2)*F(I2,3)-FF(2,3)*F(I2,2)+
     *            FF(1,2)*F(5,3) -FF(1,3)*F(5,2)
          FF(4,2)=FF(2,3)*F(I2,1)-FF(2,1)*F(I2,3)+
     *            FF(1,3)*F(5,1) -FF(1,1)*F(5,3)
          FF(4,3)=FF(2,1)*F(I2,2)-FF(2,2)*F(I2,1)+
     *            FF(1,1)*F(5,2) -FF(1,2)*F(5,1)
C         Contravariant components
          FF(5,1)=G(7) *FF(3,1)+G(8) *FF(3,2)+G(10)*FF(3,3)
          FF(5,2)=G(8) *FF(3,1)+G(9) *FF(3,2)+G(11)*FF(3,3)
          FF(5,3)=G(10)*FF(3,1)+G(11)*FF(3,2)+G(12)*FF(3,3)
          FF(6,1)=G(7) *FF(4,1)+G(8) *FF(4,2)+G(10)*FF(4,3)
          FF(6,2)=G(8) *FF(4,1)+G(9) *FF(4,2)+G(11)*FF(4,3)
          FF(6,3)=G(10)*FF(4,1)+G(11)*FF(4,2)+G(12)*FF(4,3)
C         Required functions (3 general coordinates and a travel time)
          E(2)=0.
          IF(INIPAR.LE.1) THEN
            E(1)=1.
            E(3)=0.
            AUX1=COS(PAR2)
            AUX2=SIN(PAR2)
          ELSE
            E(1)=0.
            E(3)=1.
            AUX1=COS(PAR1)
            AUX2=SIN(PAR1)
          END IF
          I22=10-I11
          DO 24 I=1,3
            F(I22,I)=-AUX2*F(I2,I)+AUX1*FF(5,I)
            F(I2,I) = AUX1*F(I2,I)+AUX2*FF(5,I)
            F(5,I)  = AUX1*F(5,I) +AUX2*FF(6,I)
   24     CONTINUE
          F(I22,4)=0.
          F(I2,4)=0.
          F(5,4)=0.
        ELSE
C         Initial surface:
C         Projection matrix
          E(1)=1.
          E(2)=0.
          E(3)=1.
C         Required functions (3 general coordinates and a travel time)
          IF(INIPAR.GE.1) THEN
C           Contravariant components of the symmetric 3*3 matrix of the
C           basis vectors of the local Cartesian coordinate system
            CALL SQRT3(G(7),G)
C           Mapping of the ray parameters onto a unit sphere
            CALL SPHERE(INIPAR,PAR1,PAR2,FF,IEND)
            IF(IEND.NE.0) THEN
              RETURN
            END IF
C           Local Cartesian coordinates
            DO 33 I=1,3
              FF(6,I)=F(6,3)*FF(1,I)+2.*F(3,3)*FF(3,I)+F(1,3)*FF(6,I)
              FF(5,I)=F(5,3)*FF(1,I)+F(3,3)*FF(2,I)+F(3,3)*FF(1,I)
     *                                                +F(1,3)*FF(5,I)
              FF(4,I)=F(4,3)*FF(1,I)+2.*F(2,3)*FF(2,I)+F(1,3)*FF(4,I)
              FF(3,I)=F(3,3)*FF(1,I)+F(1,1)*FF(3,1)
              FF(2,I)=F(2,3)*FF(1,I)+F(1,1)*FF(2,1)
              FF(1,I)=F(1,3)*FF(1,I)
   33       CONTINUE
C           General coordinates
            DO 34 I=1,6
              F(I,1)=G(1)*FF(I,1)+G(2)*FF(I,2)+G(4)*FF(I,3)
              F(I,2)=G(2)*FF(I,1)+G(3)*FF(I,2)+G(5)*FF(I,3)
              F(I,3)=G(4)*FF(I,1)+G(5)*FF(I,2)+G(6)*FF(I,3)
   34       CONTINUE
            F(1,1)=F(1,1)+X1INI
            F(1,2)=F(1,2)+X2INI
            F(1,3)=F(1,3)+X3INI
          END IF
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE SQRT3(B,A)
      REAL B(6),A(6)
C
C Subroutine SQRT3 evaluates the square root A of the given real-valued
C positive-semidefinite symmetric 3*3 matrix B.  The square root is the
C real-valued positive-semidefinite symmetric 3*3 matrix A satisfying
C the equation A*A=B.
C
C Input:
C     B...    Array of dimension at least 6, containing the components
C             B11, B12, B22, B13, B23, B33 of the given symmetric 3*3
C             matrix B.
C     A...    Array of dimension at least 6.
C The input parameter B is not altered.
C
C Output.
C     A...    Array containing the components A11, A12, A22, A13, A23,
C             A33 of the symmetric 3*3 matrix a which is the square root
C             of the given matrix B.
C
C No subroutines and external functions required.
C
C Date: 1990, January 22
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     No auxiliary storage locations
C
      IF(B(2).EQ.0..AND.B(4).EQ.0..AND.B(5).EQ.0.) THEN
C       Diagonal matrix
        IF(B(1).LT.0..OR.B(3).LT.0..OR.B(6).LT.0.) THEN
C         614
          CALL ERROR
     *             ('614 in SQRT3: Matrix is not positive-semidefinite')
C         Input matrix B is not positive-semidefinite.
        ELSE
          A(1)=SQRT(B(1))
          A(2)=0.
          A(3)=SQRT(B(3))
          A(4)=0.
          A(5)=0.
          A(6)=SQRT(B(6))
        END IF
      ELSE
C       General symmetric matrix
C       615
        CALL ERROR('615 in SQRT3: This program branch is not coded')
C       The square root of general symmetric matrix B has not been
C       coded.  At present, only the square root of diagonal matrix can
C       be evaluated.
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE SPHERE(INIPAR,PAR1,PAR2,FF,IEND)
      INTEGER INIPAR,IEND
      REAL PAR1,PAR2,FF(6,3)
C
C Subroutine SPHERE transforms spherical coordinates PAR1, PAR2 into the
C Cartesian coordinates of the corresponding point on the unit sphere.
C It also evaluates the first and second derivatives of the Cartesian
C coordinates with respect to PAR1 and PAR2.
C
C Input:
C     INIPAR..Determines the type of spherical coordinates:
C             INIPAR.LT.0: The same as for IABS(INIPAR), but with the
C               unit vector (T1,T2,T3) tangent to the ray changed to
C               (T1,-T2,-T3).
C             INIPAR.EQ.1: Polar-like spherical coordinates (colatitude,
C               longitude).
C               If SIN(colatitude).LT.0, the ray is reported out of the
C               ray-parameter domain: IEND=71.
C             INIPAR.GE.2: Geographic-like spherical coordinates
C               (longitude, latitude).
C               If COS(latitude).LT.0, the ray is reported out of the
C               ray-parameter domain: IEND=71.
C             INIPAR.EQ.3:  The unit vector tangent to the ray,
C               expressed in the local Cartesian coordinate system
C               which basis vectors are given by the square root of
C               the metric tensor at the initial point, is given by
C                 T1=PAR1*SIN(R)/R
C                 T2=PAR2*SIN(R)/R
C                 T3=     COS(R)
C               with  R=SQRT(PAR1*PAR1+PAR2*PAR2).
C               If R.GT.PI, the ray is reported out of the
C               ray-parameter domain: IEND=71.
C     PAR1,PAR2... Ray parameters.
C     FF...   Array of the dimension at least 6*3.
C None of the input parameters except FF are altered.
C
C Output:
C     FF(1:6,I)...I-th Cartesian coordinate of the point on the unit
C             sphere given by PAR1, PAR2, and its first and second
C             partial derivatives with respect to PAR1 and PAR2 in the
C             order FF, FF1, FF2, FF11, FF12, FF22.
C     IEND... Information on the initial point of the ray:
C             0...  Computation of the ray may follow.
C             71... There is no ray corresponding to the given ray
C               parameters.  I.e., the given parameters are outside the
C               domain allowed.
C
C No subroutines and external functions required.
C
C Date: 1999, April 29
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      REAL C1,C2,S1,S2,R,R1,R2,R11,R12,R22,R111,R112,R122,R222
C
C     C1,C2...Cosines of the take-off angles at a single initial point.
C     S1,S2...Sines of the take-off angles at a single initial point.
C     R,R1,R2,R11,R12,R22,R111,R112,R122,R222...
C             R=SQRT(PAR1*PAR1+PAR2*PAR2) and its partial derivatives.
C
      IEND=0
      IF(IABS(INIPAR).LE.2) THEN
        C1=COS(PAR1)
        S1=SIN(PAR1)
        C2=COS(PAR2)
        S2=SIN(PAR2)
        IF(IABS(INIPAR).LE.1) THEN
C         Polar-like spherical coordinates
          IF(S1.LT.0.) THEN
            IEND=71
            RETURN
          END IF
          FF(1,1)= S1*C2
          FF(1,2)= S1*S2
          FF(1,3)= C1
          IF(S1.EQ.0.) THEN
C           Avoiding null vectors
            S1=0.000001
          END IF
          FF(2,1)= C1*C2
          FF(2,2)= C1*S2
          FF(2,3)=-S1
          FF(3,1)=-S1*S2
          FF(3,2)= S1*C2
          FF(3,3)= 0.
          FF(4,1)=-S1*C2
          FF(4,2)=-S1*S2
          FF(4,3)=-C1
          FF(5,1)=-C1*S2
          FF(5,2)= C1*C2
          FF(5,3)= 0.
          FF(6,1)=-S1*C2
          FF(6,2)=-S1*S2
          FF(6,3)= 0.
        ELSE
C         Geographic-like spherical coordinates
          IF(C2.LT.0.) THEN
            IEND=71
            RETURN
          END IF
          FF(1,1)= C2*C1
          FF(1,2)= C2*S1
          FF(1,3)= S2
          IF(C2.EQ.0.) THEN
C           Avoiding null vectors
            C2=0.000001
          END IF
          FF(2,1)=-C2*S1
          FF(2,2)= C2*C1
          FF(2,3)= 0.
          FF(3,1)=-S2*C1
          FF(3,2)=-S2*S1
          FF(3,3)= C2
          FF(4,1)=-C2*C1
          FF(4,2)=-C2*S1
          FF(4,3)= 0.
          FF(5,1)= S2*S1
          FF(5,2)=-S2*C1
          FF(5,3)= 0.
          FF(6,1)=-C2*C1
          FF(6,2)=-C2*S1
          FF(6,3)=-S2
        END IF
      ELSE
C       Azimuthal equidistant projection
        R=SQRT(PAR1*PAR1+PAR2*PAR2)
        IF(R.GT.3.2) THEN
          IEND=71
          RETURN
        ELSE IF(R.GT.0.) THEN
          S1=SIN(R)
          IF(S1.LT.0.) THEN
            IEND=71
            RETURN
          END IF
          C1=COS(R)
          R1=PAR1/R
          R2=PAR2/R
          FF(1,1)= S1*R1
          FF(1,2)= S1*R2
          FF(1,3)= C1
          IF(S1.LT.0.004.AND.R.GT.3.1) THEN
C           Correction for nearly parallel F(2,*) and F(3,*) near
C           the singularity at R=pi
            S1=0.004
          END IF
          R11= R2*R2/R
          R12=-R1*R2/R
          R22= R1*R1/R
          R111=      -3.*R11*R1 /R
          R112=(-R2/R-3.*R12*R1)/R
          R122=(-R1/R-3.*R12*R2)/R
          R222=      -3.*R22*R2 /R
          FF(2,1)= C1*R1*R1+S1*R11
          FF(2,2)= C1*R1*R2+S1*R12
          FF(2,3)=-FF(1,1)
          FF(3,1)= FF(2,2)
          FF(3,2)= C1*R2*R2+S1*R22
          FF(3,3)=-FF(1,2)
          FF(4,1)=-S1*R1*R1*R1          +3.*C1*R1*R11+S1*R111
          FF(4,2)=-S1*R1*R1*R2+C1*R11*R2+2.*C1*R1*R12+S1*R112
          FF(4,3)=-FF(2,1)
          FF(5,1)= FF(4,2)
          FF(5,2)=-S1*R1*R2*R2+C1*R1*R22+2.*C1*R2*R12+S1*R122
          FF(5,3)=-FF(2,2)
          FF(6,1)= FF(5,2)
          FF(6,2)=-S1*R2*R2*R2          +3.*C1*R2*R22+S1*R222
          FF(6,3)=-FF(3,2)
        ELSE
          FF(1,1)= 0.
          FF(1,2)= 0.
          FF(1,3)= 1.
          FF(2,1)= 1.
          FF(2,2)= 0.
          FF(2,3)= 0.
          FF(3,1)= 0.
          FF(3,2)= 1.
          FF(3,3)= 0.
          FF(4,1)= 0.
          FF(4,2)= 0.
          FF(4,3)=-1.
          FF(5,1)= 0.
          FF(5,2)= 0.
          FF(5,3)= 0.
          FF(6,1)= 0.
          FF(6,2)= 0.
          FF(6,3)=-1.
        END IF
      END IF
C
      IF(INIPAR.LT.0) THEN
      DO 12 J=2,3
        DO 11 I=1,6
          FF(I,J)=-FF(I,J)
   11   CONTINUE
   12 CONTINUE
      END IF
C
      RETURN
      END
C
C=======================================================================
C