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: 2011, April 5
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     INISRC..Subroutine reading data for a point source from already
C             opened file SRC.
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     INIE1, INIE2, INIE3
C     INIPAR
C     ADVANC
C are described in file
C 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
C       for each source point.
C     If optional input SEP parameter RPAREW=' ' (default), data in data
C       set RPAR are not
C       repeated, they should be specified for all elementary waves
C       of all source points.
C     If optional input SEP parameter RPAREW is specified and is not
C       blank, data RPAR-(4) for a single source should be located in
C       file RPAREW.
C       If there are several source points, data RPAR-(4) are repeatedly
C       used for all source points.
C     If optional input SEP parameter WRITEW=' ' (default), data in data
C       set WRIT are not
C       repeated, they should be specified for all elementary waves
C       of all source points.
C     If optional input SEP parameter WRITEW is specified and is not
C       blank, data WRIT-(6) for a single source should be located in
C       file WRITEW.
C       If there are several source points, data WRIT-(6) are repeatedly
C       used for 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) or (2.2):
C For option INIDIR.EQ.0:
C (2.1) 'NAME',X1INI,X2INI,X3INI,TTINI,/
C For option INIDIR.NE.0:
C (2.2) 'NAME',X1INI,X2INI,X3INI,TTINI,E1DIR,E2DIR,E3DIR,/
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     E1DIR,E2DIR,E3DIR... Directional vector determining the axis
C             of the coordinate system for ray take-off parameters.
C     Default: X2INI=0, X3INI=0, TTINI=0, E1DIR=0, E2DIR=0, E3DIR=0.
C             In case (2.2), the directional vector has to be specified
C             and non-zero.
C (3) / (a slash) or end of file.
C 
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:
C             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     Description of
C     INIDIM and INIPAR.
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)
C     are stored in the common block /INITC/.  These common blocks
C     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             C.R.T.5.5.4.
C             
C             Description of YL.
C     Y...    Array containing basic quantities computed along the ray
C             at the initial point of the ray, see
C             C.R.T.5.2.1.
C             Description of Y.
C     YY...   Array containing real auxiliary quantities computed along
C             the ray at the initial point of the ray, see
C             C.R.T.5.2.2.
C             
C             Description of YY.
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.
C             
C             Description of IY.
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: 2011, April 5
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,ICB2,I,NY,NFUNC
      PARAMETER (NFUNC=4)
      REAL A(10,21),Q(21)
      REAL HP,HS,H,H1,H2,H3,H4,H5,H6
      REAL H11,H12,H22,H13,H23,H33,H14,H24,H34,H44
      REAL H15,H25,H35,H45,H55,H16,H26,H36,H46,H56,H66
      REAL YI3,YI4,YI5
      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     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     A,Q...  Elastic moduli and their spatial derivatives.
C     HP,HS,H,H1,H2,H3,H4,H5,H6,H11,H12,H22,H13,H23,H33,H14,H24,H34,H44,
C     H15,H25,H35,H45,H55,H16,H26,H36,H46,H56,H66... Hamiltonian and its
C             phase-space derivatives.
C     YI3,YI4,YI5... Temporary copies of the coordinates of the initial
C             point of the ray.
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
C     Anisotropic ray tracing:
      REAL CRTANI
      SAVE CRTANI
      DATA CRTANI/-1./
      IF(CRTANI.EQ.-1.) THEN
        CALL RSEP3R('CRTANI',CRTANI,0.)
        IF(CRTANI.NE.0..AND.KOOR().NE.0) THEN
C         6A1
          CALL ERROR('6A1 in INIT2: Anisotropy in curvilinear coord.')
C         Anisotropic ray tracing is now applicable only to Cartesian
C         coordinates.
        END IF
      END IF
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
C     Values for a non-existing ray
      DO 1 I=1,6
        YLI(I)=0.
    1 CONTINUE
      DO 2 I=1,19
        YI(I)=0.
    2 CONTINUE
C
C     Take-off parameters
      YI(20)=PAR1
      YI(21)=PAR2
C
      CALL INIS2(NFUNC,PAR1,PAR2,E,F,IEND)
      IF(IEND.NE.0) THEN
C       Non-existing ray
        ICB1I=0
        RETURN
      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
      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
C       Initial point of the ray outside the computational volume
        ICB1I=0
        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     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
C       Initial point of the ray outside the material complex blocks
        ICB1I=0
        IEND=72
        RETURN
      ENDIF
      IY(2)=0
      IY(6)=0
      IY(8)=ICB2
      CALL CODE(IY,KODIND,ICBNEW,IEND)
      IF(IEND.EQ.22) THEN
C       Initial point of the ray outside the required complex block
        ICB1I=ICB2
        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(CRTANI.EQ.0.) THEN
        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
      ELSE
        IF(F(2,4).NE.0..OR.F(3,4).NE.0.) THEN
C         6A2
          CALL ERROR('6A2 in INIT2: Anisotropy-incompatible source')
C         Anisotropic ray tracing is now applicable to an initial line
C         or initial surface only if travel time is constant along the
C         source.
        END IF
        ICB1I=ICBNEW
        CALL PARM3(IABS(ICBNEW),YI(3),A,YLI(3),Q)
        CALL HDER(ICBNEW,A,Z1,Z2,Z3,HP,HS,H,H1,H2,H3,H4,H5,H6,
     *            H11,H12,H22,H13,H23,H33,H14,H24,H34,H44,
     *            H15,H25,H35,H45,H55,H16,H26,H36,H46,H56,H66)
        V0=SQRT(H)
        YLI(1)=SQRT(HP)
        YLI(2)=SQRT(HS)
        YLI(4)=H1*V0
        YLI(5)=H2*V0
        YLI(6)=H3*V0
      END IF
C
C     Important quantities at the initial point of the ray, see
C     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+F32*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
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, see
C     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 Subroutines and external functions required:
      EXTERNAL INISRC
C
C Date: 2008, February 15
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
      IF(INIDIM.LE.0) THEN
C       Point source
        CALL INISRC(IRAY)
        IF(IRAY.LT.0) THEN
          CLOSE(LUSRC)
        END IF
      ELSE
        IRAY=-1
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE INISRC(IRAY)
      INTEGER IRAY
C
C Subroutine INISRC reads the input data for a point source from already
C opened file SRC.  The logical unit number of the
C file is specified in include file
C initd.inc.
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 Some of the storage locations of the common block are defined in this
C subroutine.
C
C Subroutines and external functions required:
      EXTERNAL UARRAY,ERROR
      REAL UARRAY
C     UARRAY...File 'array.for' of the package 'FORMS'.
C     ERROR... File 'error.for' of the package 'FORMS'.
C
C Date: 2008, February 15
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
      CHARACTER*1 TEXT
      REAL E1DIR,E2DIR,E3DIR,AUX1,AUX2
      REAL UNDEF
      SAVE UNDEF
      UNDEF=UARRAY()
C
      X1INI=UNDEF
      X2INI=0.
      X3INI=0.
      TTINI=0.
      IF(INIDIR.EQ.0) THEN
        READ(LUSRC,*,END=19) TEXT,X1INI,X2INI,X3INI,TTINI
        IF(X1INI.EQ.UNDEF) THEN
          GO TO 19
        END IF
      ELSE
        E1DIR=0.
        E2DIR=0.
        E3DIR=0.
        READ(LUSRC,*,END=19) TEXT,X1INI,X2INI,X3INI,TTINI,
     *                            E1DIR,E2DIR,E3DIR
        IF(X1INI.EQ.UNDEF) THEN
          GO TO 19
        END IF
C
C       Rotation vector:
        AUX1=SQRT(E1DIR*E1DIR+E2DIR*E2DIR+E3DIR*E3DIR)
        IF(AUX1.EQ.0.) THEN
C         605
          CALL ERROR('605 in INIS1: Directional vector not specified')
C         For option
C         INIDIR.NE.0,
C         directional vector E1DIR, E2DIR, E3DIR must be specified for
C         each source point in input file SRC.
        END IF
        AUX1=SQRT(E1DIR*E1DIR+E2DIR*E2DIR)
        IF(AUX1.EQ.0.) THEN
          AUX1=ATAN2(AUX1,E3DIR)
          E1INI=0.
          E2INI=AUX1
        ELSE
          AUX1=ATAN2(AUX1,E3DIR)/AUX1
          E1INI=-E2DIR*AUX1
          E2INI= E1DIR*AUX1
        END IF
        E3INI=0.
C
C       Rotational matrix:
        EINI=SQRT(E1INI**2+E2INI**2+E3INI**2)
        IF(EINI.LE.0.03) THEN
          AUX1=0.5
        ELSE
          AUX1=(1.-COS(EINI))/(EINI*EINI)
        END IF
        IF(EINI.EQ.0.) THEN
          AUX2=1.0
        ELSE
          AUX2=SIN(EINI)/EINI
        END IF
        E11INI=AUX1*E1INI*E1INI+COS(EINI)
        E21INI=AUX1*E2INI*E1INI+AUX2*E3INI
        E31INI=AUX1*E3INI*E1INI-AUX2*E2INI
        E12INI=AUX1*E1INI*E2INI-AUX2*E3INI
        E22INI=AUX1*E2INI*E2INI+COS(EINI)
        E32INI=AUX1*E3INI*E2INI+AUX2*E1INI
        E13INI=AUX1*E1INI*E3INI+AUX2*E2INI
        E23INI=AUX1*E2INI*E3INI-AUX2*E1INI
        E33INI=AUX1*E3INI*E3INI+COS(EINI)
      END IF
      IRAY=0
      RETURN
C
   19 CONTINUE
      IRAY=-1
      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,INISRC
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: 2008, February 15
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER LFUNC,MFUNC,I
      REAL AUX1,AUX2
      CHARACTER*80 TEXTI
      CHARACTER*3 TFUNCT
      DATA TFUNCT/'   '/
C
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 RSEP3I('INIDIR',INIDIR,0)
      IF(INIPAR.GE.0) THEN
        CALL RSEP3R('INIE1' ,E1INI, 0.)
      ELSE
        CALL RSEP3R('INIE1' ,E1INI, 3.141592654)
      END IF
      CALL RSEP3R('INIE2' ,E2INI, 0.)
      CALL RSEP3R('INIE3' ,E3INI, 0.)
      CALL RSEP3R('ADVANC',ADVANC,0.)
C
      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)
        CALL INISRC(I)
        IF(I.LT.0) THEN
C         604
          CALL ERROR('604 in INIS1: No source point')
C         For INIDIM.LE.0, there must be specified at least one source
C         point in file SRC.
        END IF
      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 
C         main input data 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
C
C     Optional rotational matrix:
      EINI=SQRT(E1INI**2+E2INI**2+E3INI**2)
      IF(EINI.LE.0.03) THEN
        AUX1=0.5
      ELSE
        AUX1=(1.-COS(EINI))/(EINI*EINI)
      END IF
      IF(EINI.EQ.0.) THEN
        AUX2=1.0
      ELSE
        AUX2=SIN(EINI)/EINI
      END IF
      E11INI=AUX1*E1INI*E1INI+COS(EINI)
      E21INI=AUX1*E2INI*E1INI+AUX2*E3INI
      E31INI=AUX1*E3INI*E1INI-AUX2*E2INI
      E12INI=AUX1*E1INI*E2INI-AUX2*E3INI
      E22INI=AUX1*E2INI*E2INI+COS(EINI)
      E32INI=AUX1*E3INI*E2INI+AUX2*E1INI
      E13INI=AUX1*E1INI*E3INI+AUX2*E2INI
      E23INI=AUX1*E2INI*E3INI-AUX2*E1INI
      E33INI=AUX1*E3INI*E3INI+COS(EINI)
C
      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(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(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(PAR1,PAR2,FF,IEND)
      INTEGER 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     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               INIPAR.EQ.1: polar-like spherical coordinates
C                 (colatitude, longitude).
C                 If SIN(colatitude).LT.0, the ray is reported out of
C                 the 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: Azimuthal equidistant projection of a unit
C                 sphere.
C                 If SQRT(PAR1*PAR1+PAR2*PAR2).GT.PI, the ray is
C                 reported out of the ray-parameter domain: IEND=71.
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 No subroutines and external functions required.
C
C Date: 2003, May 5
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I
      REAL C1,C2,S1,S2,R,R1,R2,R11,R12,R22,R111,R112,R122,R222,F1,F2,F3
C
C     I...    Loop variable.
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(EINI.NE.0.) THEN
        DO 11 I=1,6
          F1=FF(I,1)
          F2=FF(I,2)
          F3=FF(I,3)
          FF(I,1)=E11INI*F1+E12INI*F2+E13INI*F3
          FF(I,2)=E21INI*F1+E22INI*F2+E23INI*F3
          FF(I,3)=E31INI*F1+E32INI*F2+E33INI*F3
   11   CONTINUE
      END IF
C
      RETURN
      END
C
C=======================================================================
C