C 
C Subroutine file 'ap.for': Applications and processing of the results
C of complete ray tracing.
C
C Date: 1999, March 20
C Coded by Ludek Klimes
C
C This file contains subroutines related to the Chapter 7 of the paper
C on C.R.T., designed to read from the files the quantities describing
C the points of rays, and to evaluate many other quantities used in
C seismology and discussed in the Chapter C.R.T.7.  These subroutines
C may be included in user's application programs following the complete
C ray tracing program.  Individual subroutines correspond to the
C individual sections of the Chapter C.R.T.7, and may call other
C subroutines of the files of the packages 'MODEL' and 'CRT' composing
C the complete ray tracing program.
C
C Attention:  The lines dealing with curvilinear coordinates are denoted
C by '*' in the first column.  That is why this code works properly only
C in Cartesian coordinates.  To consider curvilinear coordinates, each
C '*' in the first column has to be replaced by a space.  In such a
C case, subroutine METR1 has to be called first to specify the
C coordinate system.
C
C This is a preliminary version, containing only some of the routines
C corresponding to sections of the Chapter C.R.T.7.  The routines
C denoted by * in the second column will likely be coded in the near
C future, others later on.
C
C This file consists of the following external procedures:
C     POINTB..Block data subroutine defining common block /POINTC/ to
C             store the quantities describing a point on a ray.
C             POINTB
C     AP00... Subroutine designed to read from the files the quantities
C             describing a point on a ray, and to store them into the
C             common block /POINTC/.  The input files are assumed to
C             have the structure of the output files of the complete ray
C             tracing program 'CRT', written by the subroutines of the
C             file 'writ.for'.
C             AP00
C     AP01... Subroutine designed to evaluate the travel time from the
C             initial point of a ray to a point situated on a ray, see
C             C.R.T.7.1.
C             AP01
C     AP02... Subroutine designed to evaluate the components of the
C             slowness vector either at the initial point of a ray or at
C             a point situated on a ray, see C.R.T.7.2.
C             AP02
C     AP03... Subroutine designed to evaluate the covariant components
C             of the basis vectors of the ray-centred coordinate system
C             at the initial point of a ray or at a point situated on a
C             ray, see C.R.T.7.3.
C             AP03
C     AP03A...Auxiliary subroutine to AP03, evaluating the basis of the
C             intrinsic ray-centred coordinate system at the given
C             point.
C             AP04
C*    AP04
C     AP05... Subroutine designed to evaluate the components of the
C             matrix of geometrical spreading at a point situated on a
C             ray, see C.R.T.7.5.
C             AP05
C     AP06... Subroutine designed to evaluate the components of the
C             transformation matrix P at a point situated on a ray, see
C             C.R.T.7.6.
C             AP06
C     AP07... Subroutine designed to evaluate the geometrical spreading
C             at a point situated on a ray, see C.R.T.7.7.
C             AP07
C     AP08... Subroutine designed to evaluate the components of the
C             symmetric 3*3 matrices M and N of second derivatives of
C             the travel-time field at a point situated on a ray, see
C             C.R.T.7.8.  Subroutine AP03 should be called before the
C             invocation of AP08 to define the basis of R.C.C.S.
C             AP08
C*    AP09
C*    AP10(XI,X)
C     AP11... Subroutine designed to evaluate two ray-centred
C             coordinates of a given paraxial ray.
C             AP11
C*    AP12(XI,X)
C*    AP13(XI,XF,X)
C*    AP14
C     AP15... Subroutine designed to evaluate the ray amplitudes at a
C             point situated on a ray, see C.R.T.7.15.  Subroutine AP03
C             should be called before the invocation of AP15 to define
C             the basis of R.C.C.S.
C             AP15
C*    AP16
C     AP21... Subroutine designed to evaluate the ray-theory
C             elastodynamic Green function according to C.R.T.7.21.
C             AP21
C
C Note: There are no application routines corresponding to the sections
C 7.18, 7.23 and 7.27 of the paper on C.R.T.  AP28 and subsequent
C applications, located in subroutine file 'apvar.for' do not correspond
C to any section of the paper on C.R.T.
C
C Storage in the memory:
C     When processing of the results of complete ray tracing, the
C     quantities describing some points on a ray are required to be
C     known.  In the Chapter 7 of the paper on C.R.T., three different
C     points situated on a ray are introduced:
C     O/O (O subscript O)... Initial point of the ray.
C     O/S (O subscript S)... Another point situated on the ray.  In some
C             applications it may be treated as the endpoint of the ray.
C     O/F (O subscript F)... Another point situated on the ray, usually
C             situated between the points O/O and O/S, see C.R.T.7.13:
C             Fresnel volumes.  This point is required just by few
C             applications and usually need not be defined.
C     The quantities describing the points O/O, O/S (and, possibly, O/F)
C     of a ray are stored by the subroutine AP00 into the common block
C     /POINTC/ defined in the following subroutine:
C     ------------------------------------------------------------------
C     
C
      BLOCK DATA POINTB
      INCLUDE 'pointc.inc'
C     pointc.inc
      DATA IWAVE/0/,IRAY/0/
      END
C     ------------------------------------------------------------------
C
C=======================================================================
C
C     
C
      SUBROUTINE AP00(LU1,LU2,LU3)
      INTEGER LU1,LU2,LU3
C
C This subroutine reads from the given files the quantities describing
C some points on a ray, and stores them into the common block /POINTC/.
C Quantities at one point are stored during one invocation of AP00.
C The input files are assumed to have the structure of the output files
C of the complete ray tracing program 'CRT', written by the subroutines
C of the file 'writ.for'.  When all points of the given files are read
C over, IWAVE and IRAY of the common block /POINTC/ are set to zeros.
C The locations IWAVE and IRAY of the common block /POINTC/ should not
C be changed but with the following exception: They must be set to zeros
C before this subroutine is called with altered (reopened) files
C corresponding to the input parameters LU1, LU2 and LU3.
C In the Chapter 7 of the paper on C.R.T., the initial point of the ray
C is denoted O/O (O subscript O) and the points on the ray are denoted
C O/F and O/S.  After the invocation of this subroutine, a user may want
C to call his own subroutine deciding whether the point on a ray is to
C be ignored (e.g. when it is too far from the receivers) or processed
C on.
C
C Input:
C     LU1...  Zero if the point O/F of the ray need not be defined.
C             Otherwise the logical unit number of the external input
C             device containing a file with the quantities along rays
C             (see C.R.T.5.5.1) or a file with the quantities at a
C             specified surface (see C.R.T.5.5.2).  Points O/F are read
C             from file LU1.  If there is a point O/S of the same ray in
C             file LU2, the quantities corresponding to these points are
C             stored in common block /POINTC/ and subroutine AP00 is
C             exited.  A ray with its points stored only in one of files
C             LU1 and LU2 is skipped.
C     LU2...  Logical unit number of the external input device
C             containing a file with the quantities along rays (see
C             C.R.T.5.5.1, just for LU1=0) or a file with the quantities
C             at a specified surface (see C.R.T.5.5.2).  For LU1=0, all
C             points O/S from this file will be successively stored into
C             the common block /POINTC/.
C             If LU1=0 and LU2=0, only initial points of rays will be
C             read from LU3 and stored in common block /POINTC/, one
C             initial point per each invocation of AP00.
C     LU3...  Logical unit number of the external input device
C             containing the file with the quantities at the initial
C             points of rays, corresponding to the above file LU1 or LU2
C             (see C.R.T.6.1).  If LU1=0, file LU3 must contain initial
C             points of all rays contained in file LU2 and may contain
C             initial points of other rays.  Otherwise, file LU3 must
C             contain initial points of all rays common to files LU1 and
C             LU2 and may contain initial points of other rays.
C The input parameters are not altered.
C
C No output.
C
C Common block /POINTC/:
      INCLUDE 'pointc.inc'
C     pointc.inc
C All the storage locations of the common blocks are defined in this
C subroutine.
C
C Subroutines and external functions required:
      EXTERNAL POINTB
C     POINTB..Block data subroutine of this file.
C
C Date: 1997, September 7
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER IWAVEF,IRAYF,IWAVEI,IRAYI,I
C     IWAVEF,IRAYF... Values of IWAVE,IRAY corresponding to the last
C             point O/F of a ray read in during the current invocation.
C     IWAVEI,IRAYI... Values of IWAVE,IRAY corresponding to the last
C             read in initial point of a ray.
C     I...    Auxiliary loop variable.
C
      IF(LU2.NE.0) THEN
C
C       Reading the results of the complete ray tracing:
        IWAVEF=-1
        IRAYF=0
        NYF=0
        IWAVEI=IWAVE
        IRAYI=IRAY
C
C       Points O/F and O/S on a ray:
C       For LU1=0 just the point O/S is being read,
C       otherwise, points O/F and O/S must be situated on the same ray.
   10   CONTINUE
        IF(LU1.NE.0.AND.(IWAVEF.LT.IWAVE.OR.
     *                  (IWAVEF.EQ.IWAVE.AND.IRAYF.LT.IRAY))) THEN
C         New point O/F on a ray:
          READ(LU1,END=80)
     *               IWAVEF,IRAYF,NYF,ICB1F,ISRFF,XF,YLF,(YF(I),I=1,NYF)
          IF (IWAVEF.LE.0) THEN
C           701
            CALL ERROR('701 in AP00: Wrong file with computed points')
C           Index of the elementary wave is not positive.  Maybe, the
C           file LU1 has the structure of a file LU3.
          END IF
          GO TO 10
        END IF
        IF(LU1.EQ.0.OR.IWAVEF.GT.IWAVE.OR.
     *                (IWAVEF.EQ.IWAVE.AND.IRAYF.GT.IRAY)) THEN
C         New point O/S on a ray:
          READ(LU2,END=80) IWAVE,IRAY,NY,ICB1,ISRF,X,YL,(Y(I),I=1,NY)
          IF (IWAVE.LE.0) THEN
C           702
            CALL ERROR('702 in AP00: Wrong file with computed points')
C           Index of the elementary wave is not positive.  Maybe, the
C           file LU2 has the structure of a file LU3.
          END IF
          IF(LU1.NE.0) THEN
            GO TO 10
          END IF
        END IF
C
C       Defining IPT:
        IF(IWAVE.NE.IWAVEI) THEN
          IPT=0
        ELSE IF(IRAY.NE.IRAYI) THEN
          IPT=1
        ELSE
          IPT=MAX0(1,IPT)+1
        END IF
C
C       Initial point O/O of a ray:
   20   CONTINUE
        IF(IWAVE.NE.IWAVEI.OR.IRAY.NE.IRAYI) THEN
          IF(IWAVE.LT.IWAVEI) THEN
C           704
            WRITE(*,'('' WAVE:'',I4,'',  RAY:'',I6)') IWAVE,IRAY
            CALL ERROR('704 in AP00: The wave not found')
C           The initial point of a ray from file LU2 is not found in
C           the file LU3.
          ELSE IF(IWAVE.EQ.IWAVEI.AND.IRAY.LT.IRAYI) THEN
C           705
            WRITE(*,'('' WAVE:'',I4,'',  RAY:'',I6)') IWAVE,IRAY
            CALL ERROR
     *               ('705 in AP00: Initial point of the ray not found')
C           The initial point of a ray from file LU2 is not found in
C           the file LU3.
          ELSE
C           Initial point O/O of a ray
            READ(LU3,END=90) IWAVEI,IRAYI,ICB1I,IEND,ISHEET,IREC,YLI,YI
            IF(IWAVEI.LT.0) THEN
              IWAVEI=-IWAVEI
            ELSE
              GO TO 89
            END IF
            GO TO 20
          END IF
        END IF
        RETURN
C
      ELSE
C
C       Reading only initial point of a ray:
        READ(LU3,END=80) IWAVE,IRAY,ICB1I,IEND,ISHEET,IREC,YLI,YI
        IF(IWAVE.LT.0) THEN
          IWAVE=-IWAVE
        ELSE
          GO TO 89
        END IF
        IPT=0
        RETURN
C
      END IF
C
C     End of the file with the computed points:
   80 CONTINUE
      IWAVE=0
      IRAY=0
      IPT=0
      RETURN
C
C     End of the file with the initial points of rays:
   89 CONTINUE
C     703
      CALL ERROR('703 in AP00: Wrong file with initial points')
C     Index of the elementary wave is not supplied with a minus
C     sign. Maybe, the file LU3 has the structure of a file LU2.
C
C     End of the file with the initial points of rays:
   90 CONTINUE
C     706
      CALL ERROR('706 in AP00: End of the file with initial points')
C     The initial point of a ray from file LU2 is not found in file LU3.
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE AP01(TT,TTIM)
      REAL TT,TTIM
C
C This subroutine evaluates the travel time from the initial point of a
C ray to a point situated on a ray, see C.R.T.7.1.
C
C No input.
C
C Output:
C     TT...   Travel time from the initial point O/O of a ray to the
C             point O/S read into the common block /POINTC/ by the last
C             invocation of the subroutine AP00.
C     TTIM... Imaginary travel time from the initial point O/O of a ray
C             to the point O/S (see C.R.T.7.1).
C
C Common block /POINTC/:
      INCLUDE 'pointc.inc'
C     pointc.inc
C None of the storage locations of the common block are altered.
C
C No subroutines and external functions required.
C
C Date: 1994, January 23
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
      TT  =Y(1)-YI(1)
      TTIM=Y(2)-YI(2)
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE AP02(PI,P)
      REAL PI(6),P(6)
C
C This subroutine evaluates the components of the slowness vector either
C at the initial point of a ray or at a point situated on a ray, see
C C.R.T.7.2.
C
C No input.
C
C Output:
C     PI...   Three covariant (PI(1:3)) and three contravariant
C             (PI(4:6)) components of the gradient of the travel time
C             field at the initial point O/O of a ray.
C     P...    Three covariant (P(1:3)) and three contravariant (P(4:6))
C             components of the gradient of the travel time field at the
C             point O/S read into the common block /POINTC/ by the last
C             invocation of the subroutine AP00.
C
C Common block /POINTC/:
      INCLUDE 'pointc.inc'
C     pointc.inc
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
*     INTEGER KOOR
*     EXTERNAL KOOR,METRIC,SMVPRD
C     KOOR,METRIC... File 'metric.for' of the package 'MODEL'.
C     SMVPRD... File 'means.for' of the package 'MODEL'.
C
C Date: 1994, January 23
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
*     REAL PI4,PI5,PI6
*     SAVE PI4,PI5,PI6
*     INTEGER JWAVE,JRAY
*     SAVE    JWAVE,JRAY
*     DATA    JWAVE,JRAY/0,0/
C
C.......................................................................
C
C     Covariant components:
      PI(1)=YI(6)
      PI(2)=YI(7)
      PI(3)=YI(8)
      P (1)=Y (6)
      P (2)=Y (7)
      P (3)=Y (8)
C
C     Contravariant components:
*     IF(KOOR().NE.0) THEN
C       Curvilinear coordinates:
*       IF(IWAVE.NE.JWAVE.OR.IRAY.NE.JRAY) THEN
*         CALL METRIC(YI(3),GSQRD,G,GAMMA)
*         CALL SMVPRD(G(7),YI(6),YI(7),YI(8),PI4,PI5,PI6)
*         JWAVE=IWAVE
*         JRAY=IRAY
*       END IF
*       PI(4)=PI4
*       PI(5)=PI5
*       PI(6)=PI6
*       CALL METRIC(Y (3),GSQRD,G,GAMMA)
*       CALL SMVPRD(G(7),Y(6),Y(7),Y(8),P(4),P(5),P(6))
*     ELSE
C       Cartesian coordinates:
        PI(4)=YI(6)
        PI(5)=YI(7)
        PI(6)=YI(8)
        P (4)=Y (6)
        P (5)=Y (7)
        P (6)=Y (8)
*     END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE AP03(IUSER,HI,H,HUI)
      INTEGER IUSER
      REAL HI(18),H(18),HUI(9)
C
C This subroutine evaluates the covariant components of the basis
C vectors of the ray-centred coordinate system at the initial point of a
C ray or at a point situated on a ray, see C.R.T.7.3.
C
C Input:
C     IUSER...IUSER=0... Intrinsic choice of polarization vectors.
C               Any other input need not be specified.
C             IUSER=1... User's choice of polarization vectors at the
C               initial point O/O of the ray.
C               HI(10:18) and HUI need not be specified.
C               HI(1:9) must be specified if a ray has changed since
C               the last invocation of this subroutine.
C             IUSER=2... User's choice of polarization vectors at the
C               initial point O/O of the ray.
C               HI(1:9) and HUI need not be specified.
C               HI(10:18) must be specified if a ray has changed since
C               the last invocation of this subroutine.
C             IUSER=3... Transformation matrix from the intrinsic
C               ray-centred to the user's ray-centred coordinate system
C               is given.
C               HI(1:18) need not be specified.
C               HUI(1:9) has to be specified if a ray has changed since
C               the last invocation of this subroutine.
C     HI(1:9)... Covariant components of the basis vectors of the user's
C             ray-centred coordinate system at the initial point O/O of
C             a ray.
C     HI(10:18)... Contravariant components of the basis vectors of the
C             user's ray-centred coordinate system at the initial point
C             O/O of a ray.
C     HUI...  Components of the 3*3 transformation matrix from the
C             intrinsic to the user's ray-centred coordinate system.
C
C Output:
C     HI...   Covariant (HI(1:9)) and contravariant (HI(10:18))
C             components of the basis vectors of the user's ray-centred
C             coordinate system at the initial point O/O of a ray.
C     H...    Covariant (H(1:9)) and contravariant (H(10:18)) components
C             of the basis vectors of the user's ray-centred coordinate
C             system at the point O/S read into the common block
C             /POINTC/ by the last invocation of the subroutine AP00.
C     HUI...  Components of the 3*3 transformation matrix from the
C             intrinsic to the user's ray-centred coordinate system.
C
C Common block /POINTC/:
      INCLUDE 'pointc.inc'
C     pointc.inc
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
      EXTERNAL AP03A
C     KOOR,METRIC... File 'metric.for' of the package 'MODEL'.
C     SMVPRD..File 'means.for' of the package 'MODEL'.
C     AP03A...This file.
C
C Date: 1998, October 15
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
      REAL HA(18),HA01,HA02,HA03,HA04,HA05,HA06,HA07,HA08,HA09
      REAL        HA10,HA11,HA12,HA13,HA14,HA15,HA16,HA17,HA18
      EQUIVALENCE (HA(01),HA01),(HA(02),HA02),(HA(03),HA03)
      EQUIVALENCE (HA(04),HA04),(HA(05),HA05),(HA(06),HA06)
      EQUIVALENCE (HA(07),HA07),(HA(08),HA08),(HA(09),HA09)
      EQUIVALENCE (HA(10),HA10),(HA(11),HA11),(HA(12),HA12)
      EQUIVALENCE (HA(13),HA13),(HA(14),HA14),(HA(15),HA15)
      EQUIVALENCE (HA(16),HA16),(HA(17),HA17),(HA(18),HA18)
C
C     HA...   Auxiliary storage location for covariant (HA(1:9)) and
C             contravariant (HA(10:18)) components of the basis vectors
C             of the intrinsic ray-centred coordinate system.
C
      REAL HUI1,HUI2,HUI3,HUI4,HUI5,HUI6,HUI7,HUI8,HUI9
      REAL HI01,HI02,HI03,HI04,HI05,HI06,HI07,HI08,HI09
      REAL HI10,HI11,HI12,HI13,HI14,HI15,HI16,HI17,HI18
      REAL H01,H02,H03,H04,H05,H06,H07,H08,H09
      REAL H10,H11,H12,H13,H14,H15,H16,H17,H18
      SAVE HUI1,HUI2,HUI3,HUI4,HUI5,HUI6,HUI7,HUI8,HUI9
      SAVE HI01,HI02,HI03,HI04,HI05,HI06,HI07,HI08,HI09
      SAVE HI10,HI11,HI12,HI13,HI14,HI15,HI16,HI17,HI18
      SAVE H01,H02,H03,H04,H05,H06,H07,H08,H09
      SAVE H10,H11,H12,H13,H14,H15,H16,H17,H18
C
      INTEGER JWAVE,JRAY,JPT
      SAVE    JWAVE,JRAY,JPT
      DATA    JWAVE,JRAY,JPT/-1,-1,-1/
C
C.......................................................................
C
      IF(IWAVE.NE.JWAVE.OR.IRAY.NE.JRAY) THEN
        IF(IUSER.EQ.0) THEN
          CALL AP03A(YI,HI)
          HUI(1)=1.
          HUI(2)=0.
          HUI(3)=0.
          HUI(4)=0.
          HUI(5)=1.
          HUI(6)=0.
          HUI(7)=0.
          HUI(8)=0.
          HUI(9)=1.
        ELSE
          CALL AP03A(YI,HA)
          IF(IUSER.NE.2) THEN
            IF(IUSER.EQ.1) THEN
              HUI(1)=HI(1)*HA10+HI(2)*HA11+HI(3)*HA12
              HUI(2)=HI(4)*HA10+HI(5)*HA11+HI(6)*HA12
              HUI(3)=HI(7)*HA10+HI(8)*HA11+HI(9)*HA12
              HUI(4)=HI(1)*HA13+HI(2)*HA14+HI(3)*HA15
              HUI(5)=HI(4)*HA13+HI(5)*HA14+HI(6)*HA15
              HUI(6)=HI(7)*HA13+HI(8)*HA14+HI(9)*HA15
              HUI(7)=HI(1)*HA16+HI(2)*HA17+HI(3)*HA18
              HUI(8)=HI(4)*HA16+HI(5)*HA17+HI(6)*HA18
              HUI(9)=HI(7)*HA16+HI(8)*HA17+HI(9)*HA18
            END IF
            HI(10)=HUI(1)*HA10+HUI(4)*HA11+HUI(7)*HA12
            HI(11)=HUI(2)*HA10+HUI(5)*HA11+HUI(8)*HA12
            HI(12)=HUI(3)*HA10+HUI(6)*HA11+HUI(9)*HA12
            HI(13)=HUI(1)*HA13+HUI(4)*HA14+HUI(7)*HA15
            HI(14)=HUI(2)*HA13+HUI(5)*HA14+HUI(8)*HA15
            HI(15)=HUI(3)*HA13+HUI(6)*HA14+HUI(9)*HA15
            HI(16)=HUI(1)*HA16+HUI(4)*HA17+HUI(7)*HA18
            HI(17)=HUI(2)*HA16+HUI(5)*HA17+HUI(8)*HA18
            HI(18)=HUI(3)*HA16+HUI(6)*HA17+HUI(9)*HA18
          END IF
          IF(IUSER.NE.1) THEN
            IF(IUSER.EQ.2) THEN
              HUI(1)=HI(10)*HA01+HI(11)*HA02+HI(12)*HA03
              HUI(2)=HI(13)*HA01+HI(14)*HA02+HI(15)*HA03
              HUI(3)=HI(16)*HA01+HI(17)*HA02+HI(18)*HA03
              HUI(4)=HI(10)*HA04+HI(11)*HA05+HI(12)*HA06
              HUI(5)=HI(13)*HA04+HI(14)*HA05+HI(15)*HA06
              HUI(6)=HI(16)*HA04+HI(17)*HA05+HI(18)*HA06
              HUI(7)=HI(10)*HA07+HI(11)*HA08+HI(12)*HA09
              HUI(8)=HI(13)*HA07+HI(14)*HA08+HI(15)*HA09
              HUI(9)=HI(16)*HA07+HI(17)*HA08+HI(18)*HA09
            END IF
            HI( 1)=HUI(1)*HA01+HUI(4)*HA02+HUI(7)*HA03
            HI( 2)=HUI(2)*HA01+HUI(5)*HA02+HUI(8)*HA03
            HI( 3)=HUI(3)*HA01+HUI(6)*HA02+HUI(9)*HA03
            HI( 4)=HUI(1)*HA04+HUI(4)*HA05+HUI(7)*HA06
            HI( 5)=HUI(2)*HA04+HUI(5)*HA05+HUI(8)*HA06
            HI( 6)=HUI(3)*HA04+HUI(6)*HA05+HUI(9)*HA06
            HI( 7)=HUI(1)*HA07+HUI(4)*HA08+HUI(7)*HA09
            HI( 8)=HUI(2)*HA07+HUI(5)*HA08+HUI(8)*HA09
            HI( 9)=HUI(3)*HA07+HUI(6)*HA08+HUI(9)*HA09
          END IF
        END IF
        HUI1=HUI(1)
        HUI2=HUI(2)
        HUI3=HUI(3)
        HUI4=HUI(4)
        HUI5=HUI(5)
        HUI6=HUI(6)
        HUI7=HUI(7)
        HUI8=HUI(8)
        HUI9=HUI(9)
        HI01=HI( 1)
        HI02=HI( 2)
        HI03=HI( 3)
        HI04=HI( 4)
        HI05=HI( 5)
        HI06=HI( 6)
        HI07=HI( 7)
        HI08=HI( 8)
        HI09=HI( 9)
        HI10=HI(10)
        HI11=HI(11)
        HI12=HI(12)
        HI13=HI(13)
        HI14=HI(14)
        HI15=HI(15)
        HI16=HI(16)
        HI17=HI(17)
        HI18=HI(18)
      ELSE
        HUI(1)=HUI1
        HUI(2)=HUI2
        HUI(3)=HUI3
        HUI(4)=HUI4
        HUI(5)=HUI5
        HUI(6)=HUI6
        HUI(7)=HUI7
        HUI(8)=HUI8
        HUI(9)=HUI9
        HI( 1)=HI01
        HI( 2)=HI02
        HI( 3)=HI03
        HI( 4)=HI04
        HI( 5)=HI05
        HI( 6)=HI06
        HI( 7)=HI07
        HI( 8)=HI08
        HI( 9)=HI09
        HI(10)=HI10
        HI(11)=HI11
        HI(12)=HI12
        HI(13)=HI13
        HI(14)=HI14
        HI(15)=HI15
        HI(16)=HI16
        HI(17)=HI17
        HI(18)=HI18
      END IF
C
      IF(IWAVE.NE.JWAVE.OR.IRAY.NE.JRAY.OR.IPT.NE.JPT) THEN
        IF(IUSER.EQ.0) THEN
          CALL AP03A(Y,H)
        ELSE
          CALL AP03A(Y,HA)
          H( 1)=HUI(1)*HA01+HUI(4)*HA02+HUI(7)*HA03
          H( 2)=HUI(2)*HA01+HUI(5)*HA02+HUI(8)*HA03
          H( 3)=HUI(3)*HA01+HUI(6)*HA02+HUI(9)*HA03
          H( 4)=HUI(1)*HA04+HUI(4)*HA05+HUI(7)*HA06
          H( 5)=HUI(2)*HA04+HUI(5)*HA05+HUI(8)*HA06
          H( 6)=HUI(3)*HA04+HUI(6)*HA05+HUI(9)*HA06
          H( 7)=HUI(1)*HA07+HUI(4)*HA08+HUI(7)*HA09
          H( 8)=HUI(2)*HA07+HUI(5)*HA08+HUI(8)*HA09
          H( 9)=HUI(3)*HA07+HUI(6)*HA08+HUI(9)*HA09
          H(10)=HUI(1)*HA10+HUI(4)*HA11+HUI(7)*HA12
          H(11)=HUI(2)*HA10+HUI(5)*HA11+HUI(8)*HA12
          H(12)=HUI(3)*HA10+HUI(6)*HA11+HUI(9)*HA12
          H(13)=HUI(1)*HA13+HUI(4)*HA14+HUI(7)*HA15
          H(14)=HUI(2)*HA13+HUI(5)*HA14+HUI(8)*HA15
          H(15)=HUI(3)*HA13+HUI(6)*HA14+HUI(9)*HA15
          H(16)=HUI(1)*HA16+HUI(4)*HA17+HUI(7)*HA18
          H(17)=HUI(2)*HA16+HUI(5)*HA17+HUI(8)*HA18
          H(18)=HUI(3)*HA16+HUI(6)*HA17+HUI(9)*HA18
        END IF
        H01=H( 1)
        H02=H( 2)
        H03=H( 3)
        H04=H( 4)
        H05=H( 5)
        H06=H( 6)
        H07=H( 7)
        H08=H( 8)
        H09=H( 9)
        H10=H(10)
        H11=H(11)
        H12=H(12)
        H13=H(13)
        H14=H(14)
        H15=H(15)
        H16=H(16)
        H17=H(17)
        H18=H(18)
      ELSE
        H( 1)=H01
        H( 2)=H02
        H( 3)=H03
        H( 4)=H04
        H( 5)=H05
        H( 6)=H06
        H( 7)=H07
        H( 8)=H08
        H( 9)=H09
        H(10)=H10
        H(11)=H11
        H(12)=H12
        H(13)=H13
        H(14)=H14
        H(15)=H15
        H(16)=H16
        H(17)=H17
        H(18)=H18
      END IF
C
      JWAVE=IWAVE
      JRAY=IRAY
      JPT=IPT
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE AP03A(YA,HA)
      REAL YA(11),HA(18)
C
C Auxiliary subroutine to AP03, evaluates the basis of the intrinsic
C ray-centred coordinate system at the given point.
C
C Input:
C     YA...   Quantities describing a point of a ray
C
C Output:
C     HA...   Covariant (HA(1:9)) and contravariant (HA(10:18))
C             components of the basis vectors of the intrinsic
C             ray-centred coordinate system at the given point.
C
C Subroutines and external functions required:
*     INTEGER KOOR
*     EXTERNAL KOOR,METRIC,SMVPRD
C     KOOR,METRIC... File 'metric.for' of the package 'MODEL'.
C     SMVPRD...File 'means.for' of the package 'MODEL'.
C
C Date: 1994, January 15
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
      REAL AUX
C
C.......................................................................
C
      HA(01)=YA( 9)
      HA(02)=YA(10)
      HA(03)=YA(11)
*     IF(KOOR().NE.0) THEN
C       Curvilinear coordinates:
*       CALL METRIC(YA(3),GSQRD,G,GAMMA)
*       CALL SMVPRD(G(7),HA(01),HA(02),HA(03),HA(10),HA(11),HA(12))
*       CALL SMVPRD(G(7),YA(6),YA(7),YA(8),HA(16),HA(17),HA(18))
*       AUX=SQRT(YA(6)*HA(16)+YA(7)*HA(17)+YA(8)*HA(18))
*       HA(07)=YA(6)/AUX
*       HA(08)=YA(7)/AUX
*       HA(09)=YA(8)/AUX
*       HA(16)=HA(16)/AUX
*       HA(17)=HA(17)/AUX
*       HA(18)=HA(18)/AUX
*       HA(04)=(HA(17)*HA(12)-HA(18)*HA(11))*GSQRD
*       HA(05)=(HA(18)*HA(10)-HA(16)*HA(12))*GSQRD
*       HA(06)=(HA(16)*HA(11)-HA(17)*HA(10))*GSQRD
*       HA(13)=(HA(08)*HA(03)-HA(09)*HA(02))/GSQRD
*       HA(14)=(HA(09)*HA(01)-HA(07)*HA(03))/GSQRD
*       HA(15)=(HA(07)*HA(02)-HA(06)*HA(01))/GSQRD
*     ELSE
C       Cartesian coordinates:
        AUX=SQRT(YA(6)*YA(6)+YA(7)*YA(7)+YA(8)*YA(8))
        HA(07)=YA(6)/AUX
        HA(08)=YA(7)/AUX
        HA(09)=YA(8)/AUX
        HA(10)=HA(01)
        HA(11)=HA(02)
        HA(12)=HA(03)
        HA(16)=HA(07)
        HA(17)=HA(08)
        HA(18)=HA(09)
        HA(04)=HA(17)*HA(12)-HA(18)*HA(11)
        HA(05)=HA(18)*HA(10)-HA(16)*HA(12)
        HA(06)=HA(16)*HA(11)-HA(17)*HA(10)
        HA(13)=HA(04)
        HA(14)=HA(05)
        HA(15)=HA(06)
*     END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE AP05(IUSER,HUI,Q11,Q21,Q12,Q22)
      INTEGER IUSER
      REAL HUI(5),Q11,Q21,Q12,Q22
C
C This subroutine evaluates the components of the matrix of geometrical
C spreading at a point situated on a ray, see C.R.T.7.5.
C
C Input:
C     IUSER...IUSER=0... Intrinsic choice of polarization vectors.
C               Any other input need not be specified.
C             Otherwise, user's choice of polarization vectors at the
C               initial point O/O of the ray.
C               HUI(1:9) has to be specified.
C     HUI(1),HUI(2),HUI(4),HUI(5)...  Components HUI11,HUI21,HUI12,HUI22
C             of the 2*2 transformation matrix from the intrinsic to the
C             user's ray-centred coordinate system.
C
C Output:
C     Q11,Q21,Q12,Q22... Components of the matrix of geometrical
C             spreading in the user's ray-centred coordinate system at
C             the point O/S read into the common block /POINTC/ by the
C             last invocation of the subroutine AP00.
C
C Common block /POINTC/:
      INCLUDE 'pointc.inc'
C     pointc.inc
C None of the storage locations of the common block are altered.
C
C No subroutines and external functions required.
C
C Date: 1994, January 23
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
      REAL AUX
C
      Q11=Y(12)*YI(12)+Y(16)*YI(13)+Y(20)*YI(14)+Y(24)*YI(15)
      Q21=Y(13)*YI(12)+Y(17)*YI(13)+Y(21)*YI(14)+Y(25)*YI(15)
      Q12=Y(12)*YI(16)+Y(16)*YI(17)+Y(20)*YI(18)+Y(24)*YI(19)
      Q22=Y(13)*YI(16)+Y(17)*YI(17)+Y(21)*YI(18)+Y(25)*YI(19)
      IF(IUSER.NE.0) THEN
        AUX=Q11
        Q11=HUI(1)*AUX+HUI(4)*Q21
        Q21=HUI(2)*AUX+HUI(5)*Q21
        AUX=Q12
        Q12=HUI(1)*AUX+HUI(4)*Q22
        Q22=HUI(2)*AUX+HUI(5)*Q22
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE AP06(IUSER,HUI,P11,P21,P12,P22)
      INTEGER IUSER
      REAL HUI(5),P11,P21,P12,P22
C
C This subroutine evaluates the components of the transformation matrix
C P at a point situated on a ray, see C.R.T.7.6.
C
C Input:
C     IUSER...IUSER=0... Intrinsic choice of polarization vectors.
C               Any other input need not be specified.
C             Otherwise, user's choice of polarization vectors at the
C               initial point O/O of the ray.
C               HUI(1:9) has to be specified.
C     HUI(1),HUI(2),HUI(4),HUI(5)...  Components HUI11,HUI21,HUI12,HUI22
C             of the 2*2 transformation matrix from the intrinsic to the
C             user's ray-centred coordinate system.
C
C Output:
C     P11,P21,P12,P22... Components of the transformation matrix P from
C             ray coordinates to the user's ray-centred components of
C             the slowness vector at the point O/S read into the common
C             block /POINTC/ by the last invocation of the subroutine
C             AP00.
C
C Common block /POINTC/:
      INCLUDE 'pointc.inc'
C     pointc.inc
C None of the storage locations of the common block are altered.
C
C No subroutines and external functions required.
C
C Date: 1994, January 23
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
      REAL AUX
C
      P11=Y(14)*YI(12)+Y(18)*YI(13)+Y(22)*YI(14)+Y(26)*YI(15)
      P21=Y(15)*YI(12)+Y(19)*YI(13)+Y(23)*YI(14)+Y(27)*YI(15)
      P12=Y(14)*YI(16)+Y(18)*YI(17)+Y(22)*YI(18)+Y(26)*YI(19)
      P22=Y(15)*YI(16)+Y(19)*YI(17)+Y(23)*YI(18)+Y(27)*YI(19)
      IF(IUSER.NE.0) THEN
        AUX=P11
        P11=HUI(1)*AUX+HUI(4)*P21
        P21=HUI(2)*AUX+HUI(5)*P21
        AUX=P12
        P12=HUI(1)*AUX+HUI(4)*P22
        P22=HUI(2)*AUX+HUI(5)*P22
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE AP07(QDETI,QDET,VI,V,RHOI,RHO,INIDIM)
      REAL QDETI,QDET,VI,V,RHOI,RHO
      INTEGER INIDIM
C
C This subroutine evaluates the geometrical spreading at a point
C situated on a ray, see C.R.T.7.7.
C
C No input.
C
C Output:
C     QDETI.. For a regular surface source, geometrical spreading at the
C               initial point O/O of the ray.
C             For a line source, geometrical spreading at the distance
C               epsilon from the initial point O/O, measured along the
C               ray, divided by the square root from distance epsilon
C               (limit for epsilon approaching zero from the right).
C             For a point source, geometrical spreading at the distance
C               epsilon from the initial point O/O, divided by distance
C               epsilon (limit for epsilon approaching zero from the
C               right).  Refer to equation (7.47) divided by equation
C               (7.49) in C.R.T.7.15.
C     QDET... Geometrical spreading at the point O/S read into common
C             block /POINTC/ by the last invocation of subroutine AP00.
C     VI...   Velocity at the initial point O/O of the ray.
C     V...    Velocity at the point O/S read into common block /POINTC/
C             by the last invocation of subroutine AP00.
C     RHOI... Density at the initial point O/O of the ray.
C     RHO...  Density at the point O/S read into common block /POINTC/
C             by the last invocation of subroutine AP00.
C     INIDIM..0: Point source
C             1: Line source
C             2: Regular surface source
C
C Common block /POINTC/:
      INCLUDE 'pointc.inc'
C     pointc.inc
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
      EXTERNAL AP05
C     AP05... This file.
C
C Date: 1995, August 11
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
      REAL DUMMY(5),Q11,Q21,Q12,Q22
C
C     Velocities and densities:
      VI=1./SQRT(YI(6)**2+YI(7)**2+YI(8)**2)
      V =1./SQRT(Y (6)**2+Y (7)**2+Y (8)**2)
      RHOI=YLI(3)
      RHO =YL (3)
C
C     Geometrical spreading at the initial point:
      IF(YI(12).EQ.0..AND.YI(13).EQ.0..AND.
     *   YI(16).EQ.0..AND.YI(17).EQ.0.) THEN
C       Point source
        INIDIM=0
        QDETI=ABS((YI(14)*YI(19)-YI(15)*YI(18)))*VI*VI
      ELSE
C       Surface source
        INIDIM=2
        QDETI=ABS(YI(12)*YI(17)-YI(13)*YI(16))
        IF(QDETI.LT.0.000001*ABS(YI(12)*YI(17))) THEN
C         Line source
          INIDIM=1
          QDETI=ABS((YI(12)*YI(19)-YI(13)*YI(18)
     *              +YI(14)*YI(17)-YI(15)*YI(16)))*VI
        END IF
      END IF
C
C     Geometrical spreading at the ray point:
      CALL AP05(0,DUMMY,Q11,Q21,Q12,Q22)
      QDET=ABS(Q11*Q22-Q12*Q21)
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE AP08(IUSER,H,HUI,RM,RN)
      INTEGER IUSER
      REAL H(9),HUI(9),RM(6),RN(6)
C
C This subroutine evaluates the components of the symmetric 3*3 matrices
C M and N of second derivatives of the travel-time field at a point
C situated on a ray, see C.R.T.7.8.  Subroutine AP03 should be called
C before the invocation of AP08 to define input arguments H and HUI, see
C below.
C
C Input:
C     IUSER...IUSER=0... Intrinsic choice of polarization vectors.
C               HUI(1:9) need not be specified.
C             Otherwise, user's choice of polarization vectors at the
C               initial point O/O of the ray.
C               HUI(1:9) has to be specified.
C     H...    Covariant components of the basis vectors of the user's
C             ray-centred coordinate system at the point O/S read into
C             the common block /POINTC/ by the last invocation of the
C             subroutine AP00.
C     HUI...  Components of the 3*3 transformation matrix from the
C             intrinsic to the user's ray-centred coordinate system.
C
C Output:
C     RM...   Components M11,M12,M22,M13,M23,M33 of the second covariant
C             derivatives of the travel-time field in the user's
C             ray-centred coordinate system, at the point O/S read into
C             the common block /POINTC/ by the last invocation of the
C             subroutine AP00.
C     RN...   Components N11,N12,N22,N13,N23,N33 of the second partial
C             derivatives of the travel-time field in the general model
C             coordinates, at the point O/S read into the common block
C             /POINTC/ by the last invocation of the subroutine AP00.
C
C Common block /POINTC/:
      INCLUDE 'pointc.inc'
C     pointc.inc
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
*     INTEGER KOOR
*     EXTERNAL KOOR,METRIC
C     KOOR,METRIC... File 'metric.for' of the package 'MODEL'.
C
C Date: 1997, November 13
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
      REAL QDET,Q11,Q21,Q12,Q22,P11,P21,P12,P22
      REAL RM1,RM2,RM3,AUX1,AUX2,AUX3
C
C     Matrix M in the intrinsic R.C.C.S.
      CALL AP05(0,HUI,Q11,Q21,Q12,Q22)
      CALL AP06(0,HUI,P11,P21,P12,P22)
      QDET=Q11*Q22-Q12*Q21
      IF(QDET.EQ.0.) THEN
        QDET=1.0E-12*(Q11+Q22)**2
      END IF
      IF(QDET.EQ.0.) THEN
        RM(1)=0.
        RM(2)=0.
        RM(3)=0.
      ELSE
        RM(1)=( P11*Q22-P12*Q21)/QDET
        RM(2)=( P21*Q22-P22*Q21)/QDET
        RM(3)=(-P21*Q12+P22*Q11)/QDET
      END IF
      IF(ICB1.GE.0) THEN
        AUX2=YL(1)**2
      ELSE
        AUX2=YL(2)**2
      END IF
C     Slowness gradient in R.C.C.S.
      RM(4)=-(H(1)*YL(4)+H(2)*YL(5)+H(3)*YL(6))/AUX2
      RM(5)=-(H(4)*YL(4)+H(5)*YL(5)+H(6)*YL(6))/AUX2
      RM(6)=-(H(7)*YL(4)+H(8)*YL(5)+H(9)*YL(6))/AUX2
C
      IF(IUSER.NE.0) THEN
        AUX1=RM(1)*HUI(1)+RM(2)*HUI(4)+RM(4)*HUI(7)
        AUX2=RM(2)*HUI(1)+RM(3)*HUI(4)+RM(5)*HUI(7)
        AUX3=RM(4)*HUI(1)+RM(5)*HUI(4)+RM(6)*HUI(7)
        RM1  =HUI(1)*AUX1+HUI(4)*AUX2+HUI(7)*AUX3
        AUX1=RM(1)*HUI(2)+RM(2)*HUI(5)+RM(4)*HUI(8)
        AUX2=RM(2)*HUI(2)+RM(3)*HUI(5)+RM(5)*HUI(8)
        AUX3=RM(4)*HUI(2)+RM(5)*HUI(5)+RM(6)*HUI(8)
        RM2  =HUI(1)*AUX1+HUI(4)*AUX2+HUI(7)*AUX3
        RM3  =HUI(2)*AUX1+HUI(5)*AUX2+HUI(8)*AUX3
        AUX1=RM(1)*HUI(3)+RM(2)*HUI(6)+RM(4)*HUI(9)
        AUX2=RM(2)*HUI(3)+RM(3)*HUI(6)+RM(5)*HUI(9)
        AUX3=RM(4)*HUI(3)+RM(5)*HUI(6)+RM(6)*HUI(9)
        RM(4)=HUI(1)*AUX1+HUI(4)*AUX2+HUI(7)*AUX3
        RM(5)=HUI(2)*AUX1+HUI(5)*AUX2+HUI(8)*AUX3
        RM(6)=HUI(3)*AUX1+HUI(6)*AUX2+HUI(9)*AUX3
        RM(1)=RM1
        RM(2)=RM2
        RM(3)=RM3
      END IF
      AUX1=RM(1)*H(1)+RM(2)*H(4)+RM(4)*H(7)
      AUX2=RM(2)*H(1)+RM(3)*H(4)+RM(5)*H(7)
      AUX3=RM(4)*H(1)+RM(5)*H(4)+RM(6)*H(7)
      RN(1)=H(1)*AUX1+H(4)*AUX2+H(7)*AUX3
      AUX1=RM(1)*H(2)+RM(2)*H(5)+RM(4)*H(8)
      AUX2=RM(2)*H(2)+RM(3)*H(5)+RM(5)*H(8)
      AUX3=RM(4)*H(2)+RM(5)*H(5)+RM(6)*H(8)
      RN(2)=H(1)*AUX1+H(4)*AUX2+H(7)*AUX3
      RN(3)=H(2)*AUX1+H(5)*AUX2+H(8)*AUX3
      AUX1=RM(1)*H(3)+RM(2)*H(6)+RM(4)*H(9)
      AUX2=RM(2)*H(3)+RM(3)*H(6)+RM(5)*H(9)
      AUX3=RM(4)*H(3)+RM(5)*H(6)+RM(6)*H(9)
      RN(4)=H(1)*AUX1+H(4)*AUX2+H(7)*AUX3
      RN(5)=H(2)*AUX1+H(5)*AUX2+H(8)*AUX3
      RN(6)=H(3)*AUX1+H(6)*AUX2+H(9)*AUX3
*     IF(KOOR().NE.0) THEN
C       curvilinear coordinates:
*       CALL METRIC(Y(3),GSQRD,G,GAMMA)
*       RN(1)=RN(1)+GAMMA(1)*Y(6)+GAMMA( 7)*Y(6)+GAMMA(13)*Y(8)
*       RN(2)=RN(2)+GAMMA(2)*Y(6)+GAMMA( 8)*Y(6)+GAMMA(14)*Y(8)
*       RN(3)=RN(3)+GAMMA(3)*Y(6)+GAMMA( 9)*Y(6)+GAMMA(15)*Y(8)
*       RN(4)=RN(4)+GAMMA(4)*Y(6)+GAMMA(10)*Y(6)+GAMMA(16)*Y(8)
*       RN(5)=RN(5)+GAMMA(5)*Y(6)+GAMMA(11)*Y(6)+GAMMA(17)*Y(8)
*       RN(6)=RN(6)+GAMMA(6)*Y(6)+GAMMA(12)*Y(6)+GAMMA(18)*Y(8)
*     END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE AP11(IUSER,HUI,DPAR1,DPAR2,DQ1,DQ2,DP1,DP2)
      INTEGER IUSER
      REAL HUI(5),DPAR1,DPAR2,DQ1,DQ2,DP1,DP2
C
C Subroutine designed to evaluate two ray-centred coordinates of a given
C paraxial ray, see C.R.T.7.11.
C
C Input:
C     IUSER...IUSER=0... Intrinsic choice of polarization vectors.
C               Any other input need not be specified.
C             Otherwise, user's choice of polarization vectors at the
C               initial point O/O of the ray.
C               HUI(1:9) has to be specified.
C     HUI(1),HUI(2),HUI(4),HUI(5)...  Components HUI11,HUI21,HUI12,HUI22
C             of the 2*2 transformation matrix from the intrinsic to the
C             user's ray-centred coordinate system.
C     DPAR1,DPAR2... Increment in take-off ray parameters of a paraxial
C             ray.
C
C Output:
C     DQ1,DQ2... Ray-centred coordinates of a point of the paraxial ray.
C     DP1,DP2... Ray-centred components of the slowness vector the
C             paraxial ray.
C
C Subroutines and external functions required:
      EXTERNAL AP05,AP06
C     AP05... This file.
C     AP06... This file.
C
C Date: 1995, August 11
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
      REAL Q11,Q21,Q12,Q22,P11,P21,P12,P22
C
      CALL AP05(IUSER,HUI,Q11,Q21,Q12,Q22)
      CALL AP06(IUSER,HUI,P11,P21,P12,P22)
      DQ1=Q11*DPAR1+Q12*DPAR2
      DQ2=Q21*DPAR1+Q22*DPAR2
      DP1=P11*DPAR1+P12*DPAR2
      DP2=P21*DPAR1+P22*DPAR2
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE AP15(IUSER,HI,H,HUI,ZI,Z,AMPLI,AMPL)
      INTEGER IUSER
      REAL HI(18),H(18),HUI(9),ZI(9),Z(9),AMPLI(6),AMPL(6)
C
C This subroutine evaluates the ray amplitudes at a point situated on a
C ray, see C.R.T.7.15.  Subroutine AP03 should be called before the
C invocation of AP15 to define the input arguments HI and H, see below.
C
C Input:
C     IUSER...IUSER=0... Intrinsic choice of polarization vectors.
C               HUI(1:9) need not be specified.
C             Otherwise, user's choice of polarization vectors at the
C               initial point O/O of the ray.
C               HUI(1:9) has to be specified.
C     HI...   Covariant (HI(1:9)) and contravariant (HI(10:18))
C             components of the basis vectors of the user's ray-centred
C             coordinate system at the initial point O/O of a ray.
C     H...    Covariant (H(1:9)) and contravariant (H(10:18)) components
C             of the basis vectors of the user's ray-centred coordinate
C             system at the point O/S read into the common block
C             /POINTC/ by the last invocation of the subroutine AP00.
C     HUI...  Components of the 3*3 transformation matrix from the
C             intrinsic to the user's ray-centred coordinate system.
C     ZI...   Contravariant components of the basis vectors of the local
C             coordinate system at the initial point O/O of a ray.
C             See C.R.T.7.15, eq.(7.44).
C     Z...    Contravariant components of the basis vectors of the local
C             recording coordinate system at the point O/S read into the
C             common block /POINTC/ by the last invocation of the
C             subroutine AP00.  See C.R.T.7.15, eq.(7.43).
C     AMPLI...Components Re(A1),Im(A1),Re(A2),Im(A2),Re(A3),Im(A3) of
C             the ray amplitude in the local coordinate system
C             multiplied by the square root of the geometrical
C             spreading, at the initial point O/O of a ray.
C             See C.R.T.7.15, eq.(7.47).
C
C Output:
C     AMPL... Components Re(A1),Im(A1),Re(A2),Im(A2),Re(A3),Im(A3) of
C             the ray amplitude in the local recording coordinate system
C             at the point O/S read into the common block /POINTC/ by
C             the last invocation of the subroutine AP00.
C             See C.R.T.7.15, eq.(7.45).
C
C Common block /POINTC/:
      INCLUDE 'pointc.inc'
C     pointc.inc
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
      EXTERNAL AP07
C     AP05... This file.
C     AP07... This file.
C
C Date: 1996, September 30
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
      INTEGER INIDIM,I
      REAL QDETI,QDET,VI,V,RHOI,RHO
      REAL AUX01,AUX02,AUX03,AUX04,AUX05,AUX06,AUX07,AUX08,AUX09,AUX10
      REAL AUX11,AUX12,AUX13,AUX14,AUX15,AUX16,AUX17,AUX18,AUX19,SCALAR
C
      IF(IUSER.NE.0) THEN
C       715
        CALL ERROR
     *           ('715 in AP15: User''s choice of R.C.C.S. not allowed')
C       Nonzero input parameter IUSER of the subroutine AP15 indicates
C       user's choice of polarization vectors.  This option has not been
C       included in this subroutine yet.  Subroutine AP03 should be
C       called before the invocation of AP15 with IUSER=0.
      END IF
C
C     Scalar multiplication factor in eq.(7.45)
      CALL AP07(QDETI,QDET,VI,V,RHOI,RHO,INIDIM)
      SCALAR=YLI(3)/(YL(3)*QDET)
      IF(ICB1.GT.0) THEN
C       P-wave
        SCALAR=SCALAR/YL(1)
      ELSE
C       S-wave
        SCALAR=SCALAR/YL(2)
      END IF
C     Note: The square root and velocity at the initial point of the ray
C     will be applied later on.
C
C     Transposed inverse matrix to Z multiplied by its determinant
      AUX01=Z(5)*Z(9)-Z(6)*Z(8)
      AUX02=Z(6)*Z(7)-Z(4)*Z(9)
      AUX03=Z(4)*Z(8)-Z(5)*Z(7)
      AUX04=Z(3)*Z(8)-Z(2)*Z(9)
      AUX05=Z(1)*Z(9)-Z(3)*Z(7)
      AUX06=Z(2)*Z(7)-Z(1)*Z(8)
      AUX07=Z(2)*Z(6)-Z(3)*Z(5)
      AUX08=Z(3)*Z(4)-Z(1)*Z(6)
      AUX09=Z(1)*Z(5)-Z(2)*Z(4)
      AUX10=Z(1)*AUX01+Z(2)*AUX02+Z(3)*AUX03
C
C     Transformation matrix from R.C.C.S. to the local recording
C     coordinate system Z
      IF(ICB1.GT.0.OR.NY.EQ.33.OR.NY.EQ.39) THEN
C       P-wave or interface conversion coefficients applied
        AUX17=AUX01*H(16)+AUX04*H(17)+AUX07*H(18)
        AUX18=AUX02*H(16)+AUX05*H(17)+AUX08*H(18)
        AUX19=AUX03*H(16)+AUX06*H(17)+AUX09*H(18)
      END IF
      IF(ICB1.LT.0.OR.NY.EQ.33.OR.NY.EQ.39) THEN
C       S-wave or interface conversion coefficients applied
        AUX11=AUX01*H(10)+AUX04*H(11)+AUX07*H(12)
        AUX12=AUX02*H(10)+AUX05*H(11)+AUX08*H(12)
        AUX13=AUX03*H(10)+AUX06*H(11)+AUX09*H(12)
        AUX14=AUX01*H(13)+AUX04*H(14)+AUX07*H(15)
        AUX15=AUX02*H(13)+AUX05*H(14)+AUX08*H(15)
        AUX16=AUX03*H(13)+AUX06*H(14)+AUX09*H(15)
      END IF
C
      IF(ICB1I.GT.0) THEN
C       P-wave at the initial point:
        SCALAR=SQRT(SCALAR*YLI(1))/AUX10
C       Matrix ZI in R.C.C.S.
        AUX01=HI(16)*ZI(1)+HI(17)*ZI(2)+HI(18)*ZI(3)
        AUX02=HI(16)*ZI(4)+HI(17)*ZI(5)+HI(18)*ZI(6)
        AUX03=HI(16)*ZI(7)+HI(17)*ZI(8)+HI(18)*ZI(9)
C       AMPLI in R.C.C.S. multiplied by scalar
        AUX07=SCALAR*(AUX01*AMPLI(1)+AUX02*AMPLI(3)+AUX03*AMPLI(5))
        AUX08=SCALAR*(AUX01*AMPLI(2)+AUX02*AMPLI(4)+AUX03*AMPLI(6))
C       AMPL in R.C.C.S.
        AUX01=Y(28)*AUX07-Y(29)*AUX08
        AUX02=Y(29)*AUX07+Y(28)*AUX08
        IF(NY.GT.29) THEN
C         P-wave to S-wave or interface conversion coefficients applied
          AUX03=Y(30)*AUX07-Y(31)*AUX08
          AUX04=Y(31)*AUX07+Y(30)*AUX08
          IF(NY.GT.31) THEN
C           Interface conversion coefficients applied
            AUX05=Y(32)*AUX07-Y(33)*AUX08
            AUX06=Y(33)*AUX07+Y(32)*AUX08
          END IF
        END IF
      ELSE
C       S-wave at the initial point:
        SCALAR=SQRT(SCALAR*YLI(2))/AUX10
C       Matrix ZI in R.C.C.S.
        AUX01=HI(10)*ZI(1)+HI(11)*ZI(2)+HI(12)*ZI(3)
        AUX02=HI(10)*ZI(4)+HI(11)*ZI(5)+HI(12)*ZI(6)
        AUX03=HI(10)*ZI(7)+HI(11)*ZI(8)+HI(12)*ZI(9)
        AUX04=HI(13)*ZI(1)+HI(14)*ZI(2)+HI(15)*ZI(3)
        AUX05=HI(13)*ZI(4)+HI(14)*ZI(5)+HI(15)*ZI(6)
        AUX06=HI(13)*ZI(7)+HI(14)*ZI(8)+HI(15)*ZI(9)
C       AMPLI in R.C.C.S. Multiplied by scalar
        AUX07=SCALAR*(AUX01*AMPLI(1)+AUX02*AMPLI(3)+AUX03*AMPLI(5))
        AUX08=SCALAR*(AUX01*AMPLI(2)+AUX02*AMPLI(4)+AUX03*AMPLI(6))
        AUX09=SCALAR*(AUX04*AMPLI(1)+AUX05*AMPLI(3)+AUX06*AMPLI(5))
        AUX10=SCALAR*(AUX04*AMPLI(2)+AUX05*AMPLI(4)+AUX06*AMPLI(6))
C       AMPL in R.C.C.S.
        I=(NY+27)/2
        AUX01=Y(28)*AUX07-Y(29)*AUX08+Y(I+1)*AUX09-Y(I+2)*AUX10
        AUX02=Y(29)*AUX07+Y(28)*AUX08+Y(I+2)*AUX09+Y(I+1)*AUX10
        IF(NY.GT.31) THEN
C         S-wave to S-wave or interface conversion coefficients applied
          AUX03=Y(30)*AUX07-Y(31)*AUX08+Y(I+3)*AUX09-Y(I+4)*AUX10
          AUX04=Y(31)*AUX07+Y(30)*AUX08+Y(I+4)*AUX09-Y(I+3)*AUX10
          IF(NY.GT.35) THEN
C           Interface conversion coefficients applied
            AUX05=Y(32)*AUX07-Y(33)*AUX08+Y(38)*AUX09-Y(39)*AUX10
            AUX06=Y(33)*AUX07+Y(32)*AUX08+Y(39)*AUX09-Y(38)*AUX10
          END IF
        END IF
      END IF
C
C     AMPL in the local recording coordinate system Z
      IF(NY.EQ.33.OR.NY.EQ.39) THEN
C       Interface conversion coefficients applied
        AMPL(1)=AUX11*AUX01+AUX14*AUX03+AUX17*AUX05
        AMPL(2)=AUX11*AUX02+AUX14*AUX04+AUX17*AUX06
        AMPL(3)=AUX12*AUX01+AUX15*AUX03+AUX18*AUX05
        AMPL(4)=AUX12*AUX02+AUX15*AUX04+AUX18*AUX06
        AMPL(5)=AUX13*AUX01+AUX16*AUX03+AUX19*AUX05
        AMPL(6)=AUX13*AUX02+AUX16*AUX04+AUX19*AUX06
      ELSE IF(ICB1.GT.0) THEN
C       P-wave
        AMPL(1)=AUX17*AUX01
        AMPL(2)=AUX17*AUX02
        AMPL(3)=AUX18*AUX01
        AMPL(4)=AUX18*AUX02
        AMPL(5)=AUX19*AUX01
        AMPL(6)=AUX19*AUX02
      ELSE
C       S-wave
        AMPL(1)=AUX11*AUX01+AUX14*AUX03
        AMPL(2)=AUX11*AUX02+AUX14*AUX04
        AMPL(3)=AUX12*AUX01+AUX15*AUX03
        AMPL(4)=AUX12*AUX02+AUX15*AUX04
        AMPL(5)=AUX13*AUX01+AUX16*AUX03
        AMPL(6)=AUX13*AUX02+AUX16*AUX04
      END IF
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE AP21(GREEN)
      REAL GREEN(32)
C
C This subroutine evaluates the ray-theory elastodynamic Green function
c according to C.R.T.7.21.
C
C Attention:
C     This subroutine should be applied only to the rays starting from
C     common initial point (point source).  Otherwise, the phase shift
C     due to caustics would be incorrect.
C
C No input.
C
C Output:
C     GREEN(1)... Travel time between receiver and source.
C     GREEN(2)... Imaginary part of the complex-valued travel time
C             between receiver and source due to attenuation.
C     GREEN(3:8)... Coordinates of the receiver and coordinates of the
C             source.
C     GREEN(9:14)... Derivatives of the travel time with respect to the
C             coordinates of the receiver and coordinates of the source.
C     GREEN(15:32)... Amplitude of the Green function: contravariant
C             components of the complex-valued 3*3 matrix Gij in model
C             coordinates, where the first subscript corresponds to the
C             receiver and the second subscript corresponds to the
C             source.  The components are ordered as
C             Re(G11),Im(G11),Re(G21),Im(G21),Re(G31),Im(G31),
C             Re(G12),Im(G12),Re(G22),Im(G22),Re(G32),Im(G32),
C             Re(G13),Im(G13),Re(G23),Im(G23),Re(G33),Im(G33).
C
C Common block /POINTC/:
      INCLUDE 'pointc.inc'
C     pointc.inc
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
      EXTERNAL AP03
C     AP03,AP03A...This file.
C     KOOR,METRIC... File 'metric.for' of the package 'MODEL'.
C     SMVPRD..File 'means.for' of the package 'MODEL'.
C
C Date: 1998, October 6
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
      REAL VI,V,RHOI,RHO,DETQ2,A,HI(18),H(18),HUI(9),AUX1,AUX2,AUX3
      REAL AR11,AR21,AR31,AR12,AR22,AR32,AR13,AR23,AR33
      REAL AI11,AI21,AI31,AI12,AI22,AI32,AI13,AI23,AI33
C
      GREEN( 1)=Y(1)-YI(1)
      GREEN( 2)=Y(2)-YI(2)
      GREEN( 3)= Y (3)
      GREEN( 4)= Y (4)
      GREEN( 5)= Y (5)
      GREEN( 6)= YI(3)
      GREEN( 7)= YI(4)
      GREEN( 8)= YI(5)
      GREEN( 9)= Y (6)
      GREEN(10)= Y (7)
      GREEN(11)= Y (8)
      GREEN(12)=-YI(6)
      GREEN(13)=-YI(7)
      GREEN(14)=-YI(8)
C
C     Velocities, densities, reciprocal geometrical spreading:
      VI=1./SQRT(YI(6)**2+YI(7)**2+YI(8)**2)
      V =1./SQRT(Y (6)**2+Y (7)**2+Y (8)**2)
      RHOI=YLI(3)
      RHO =YL (3)
      DETQ2=ABS(Y(20)*Y(25)-Y(21)*Y(24))
C
C     Scalar multiplication factor (7.71):
      A=1./12.56637/SQRT(VI*RHOI*V*RHO*DETQ2)
C
C     Amplitude of the Green function in R.C.C.S., see 5.2.1 and 5.4.4:
      AR11=0.
      AI11=0.
      AR21=0.
      AI21=0.
      AR31=0.
      AI31=0.
      AR12=0.
      AI12=0.
      AR22=0.
      AI22=0.
      AR32=0.
      AI32=0.
      AR13=0.
      AI13=0.
      AR23=0.
      AI23=0.
      AR33=0.
      AI33=0.
      IF(NY.EQ.29.AND.ICB1I.GT.0.AND.ICB1.GT.0) THEN
        AR33=A*Y(28)
        AI33=A*Y(29)
      ELSE IF(NY.EQ.31) THEN
        IF(ICB1I.GT.0.AND.ICB1.LT.0) THEN
          AR13=A*Y(28)
          AI13=A*Y(29)
          AR23=A*Y(30)
          AI23=A*Y(31)
        ELSE IF(ICB1I.LT.0.AND.ICB1.GT.0) THEN
          AR31=A*Y(28)
          AI31=A*Y(29)
          AR32=A*Y(30)
          AI32=A*Y(31)
        END IF
      ELSE IF(NY.EQ.35.AND.ICB1I.LT.0.AND.ICB1.LT.0) THEN
        AR11=A*Y(28)
        AI11=A*Y(29)
        AR21=A*Y(30)
        AI21=A*Y(31)
        AR12=A*Y(32)
        AI12=A*Y(33)
        AR22=A*Y(34)
        AI22=A*Y(35)
      ELSE IF(NY.EQ.33.AND.ICB1I.GT.0) THEN
        AR13=A*Y(28)
        AI13=A*Y(29)
        AR23=A*Y(30)
        AI23=A*Y(31)
        AR33=A*Y(32)
        AI33=A*Y(33)
      ELSE IF(NY.EQ.39.AND.ICB1I.LT.0) THEN
        AR11=A*Y(28)
        AI11=A*Y(29)
        AR21=A*Y(30)
        AI21=A*Y(31)
        AR31=A*Y(32)
        AI31=A*Y(33)
        AR12=A*Y(34)
        AI12=A*Y(35)
        AR22=A*Y(36)
        AI22=A*Y(37)
        AR32=A*Y(38)
        AI32=A*Y(39)
      END IF
C
C     Basis of R.C.C.S.:
      CALL AP03(0,HI,H,HUI)
C     Contravariant components of the basis vectors: HI(10:18),H(10:18).
C
C     Contravariant components of the amplitude of the Green function:
      AUX1=AR11*HI(10)+AR12*HI(13)+AR13*HI(16)
      AUX2=AR21*HI(10)+AR22*HI(13)+AR23*HI(16)
      AUX3=AR31*HI(10)+AR32*HI(13)+AR33*HI(16)
      GREEN(15)=H(10)*AUX1+H(13)*AUX2+H(16)*AUX3
      GREEN(17)=H(11)*AUX1+H(14)*AUX2+H(17)*AUX3
      GREEN(19)=H(12)*AUX1+H(15)*AUX2+H(18)*AUX3
      AUX1=AI11*HI(10)+AI12*HI(13)+AI13*HI(16)
      AUX2=AI21*HI(10)+AI22*HI(13)+AI23*HI(16)
      AUX3=AI31*HI(10)+AI32*HI(13)+AI33*HI(16)
      GREEN(16)=H(10)*AUX1+H(13)*AUX2+H(16)*AUX3
      GREEN(18)=H(11)*AUX1+H(14)*AUX2+H(17)*AUX3
      GREEN(20)=H(12)*AUX1+H(15)*AUX2+H(18)*AUX3
      AUX1=AR11*HI(11)+AR12*HI(14)+AR13*HI(17)
      AUX2=AR21*HI(11)+AR22*HI(14)+AR23*HI(17)
      AUX3=AR31*HI(11)+AR32*HI(14)+AR33*HI(17)
      GREEN(21)=H(10)*AUX1+H(13)*AUX2+H(16)*AUX3
      GREEN(23)=H(11)*AUX1+H(14)*AUX2+H(17)*AUX3
      GREEN(25)=H(12)*AUX1+H(15)*AUX2+H(18)*AUX3
      AUX1=AI11*HI(11)+AI12*HI(14)+AI13*HI(17)
      AUX2=AI21*HI(11)+AI22*HI(14)+AI23*HI(17)
      AUX3=AI31*HI(11)+AI32*HI(14)+AI33*HI(17)
      GREEN(22)=H(10)*AUX1+H(13)*AUX2+H(16)*AUX3
      GREEN(24)=H(11)*AUX1+H(14)*AUX2+H(17)*AUX3
      GREEN(26)=H(12)*AUX1+H(15)*AUX2+H(18)*AUX3
      AUX1=AR11*HI(12)+AR12*HI(15)+AR13*HI(18)
      AUX2=AR21*HI(12)+AR22*HI(15)+AR23*HI(18)
      AUX3=AR31*HI(12)+AR32*HI(15)+AR33*HI(18)
      GREEN(27)=H(10)*AUX1+H(13)*AUX2+H(16)*AUX3
      GREEN(29)=H(11)*AUX1+H(14)*AUX2+H(17)*AUX3
      GREEN(31)=H(12)*AUX1+H(15)*AUX2+H(18)*AUX3
      AUX1=AI11*HI(12)+AI12*HI(15)+AI13*HI(18)
      AUX2=AI21*HI(12)+AI22*HI(15)+AI23*HI(18)
      AUX3=AI31*HI(12)+AI32*HI(15)+AI33*HI(18)
      GREEN(28)=H(10)*AUX1+H(13)*AUX2+H(16)*AUX3
      GREEN(30)=H(11)*AUX1+H(14)*AUX2+H(17)*AUX3
      GREEN(32)=H(12)*AUX1+H(15)*AUX2+H(18)*AUX3
C
      RETURN
      END
C
C=======================================================================
C