C
C Subroutine file 'wan.for' to compute quantities along a ray necessary
C for computation of the Green function by means of coupling ray-theory
C in weakly anisotropic models with interfaces, and to calculate second-
C order travel time perturbations in models without interfaces.
C
C Version: 7.20
C Date: 2015, April 8
C
C Coded by: Petr Bulant
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     http://sw3d.cz/staff/bulant.htm
C
C=======================================================================
C
      SUBROUTINE WAN(LU1,LU2,LU3,GREEN,MGREEN,RPTS,IPTS,MAXPTS,NQ,
     *               ERR,LUQI,LUEIG)
C
C----------------------------------------------------------------------
      INTEGER LU1,LU2,LU3,LUQI,LUEIG
      INTEGER MGREEN,MAXPTS,NQ
      INTEGER    NQPTS
      PARAMETER (NQPTS=61)
      REAL GREEN(MGREEN),RPTS(NQPTS,MAXPTS/NQPTS),ERR(2)
      INTEGER            IPTS(NQPTS,MAXPTS/NQPTS)
C     ------------------------------------------------------------------
C Input:
C     LU1 ... Number of the logical unit connected to the CRT output
C             file with quantities along rays.
C     LU2 ... Number of the logical unit connected to the CRT output
C             file with quantities at the storing surface.
C     LU3 ... Number of the logical unit connected to the CRT output
C             file with quantities at the initial points of rays.
C     MGREEN..Dimension of an output array GREEN.
C     MAXPTS..Maximum number of records in arrays RPTS and IPTS.
C     ERR...  Output of previous invocations of WAN, zero for the first
C             invocation.
C     LUQI... Number of the logical unit connected to the output file
C             QILST with the second-order perturbations of anisotropic
C             travel times of S waves.
C             Zero if the file is not to be generated.
C     LUEIG.. Number of the logical unit connected to the output file
C             EIGLST with the eigenvectors of Christoffel matrix.
C             Zero if the file is not to be generated.
C
C Output:
C     GREEN...Array containing the Green function for the given ray
C             and for all the frequencies:
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
C             of the source.
C        GREEN(9:14)... Derivatives of the travel time with respect
C             to the coordinates of the receiver and coordinates of the
C             source.
C        GREEN(15:(14+NF*18)) ...
C             for P wave once, for S wave NF times following 18 numbers,
C             specifying 1 000 000 times enlarged amplitude of the
C             Green function: contravariant components of the complex-
C             valued 3*3 matrix Gij in model coordinates, where the
C             first subscript corresponds to the receiver and the second
C             subscript corresponds to the source.  The components are
C             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     NQ ... Number of records stored in the array GREEN.
C     IPTS,RPTS... Auxiliary arrays, redimensioned in each invocation.
C     ERR...  Maximum of calculated errors of the quasi-isotropic
C             projection of the polarization vectors.
C             ERR(1): at the source,
C             ERR(2): at the receiver.
C
C.......................................................................
C Subroutines and external functions required:
      EXTERNAL WAPTS,WAREAD,WACHRI,WAMAT,WAANGL,WAVECR,WAPROJ,
     *WAPERT,WASUM,WASUM3,WASUM4,WASUM5,
     *ERROR,RSEP1,RSEP3T,RSEP3R,RSEP3I,FORM2,LENGTH,EIGEN,
     *HIVD2,BLOCK,VELOC,PARM3,AP00,AP21,TRANS,HDER
      INTEGER LENGTH
      REAL WAANGL
C     WAPTS,WAREAD,WACHRI,WAMAT,WAANGL,WAVECR,WAPROJ,WAPERT,WASUM,
C       WASUM3,WASUM4,WASUM5 ... This file.
C     ERROR ... File
C     error.for.
C     RSEP1,RSEP3I,RSEP3T,RSEP3R ... File
C     sep.for.
C     FORM2 ... File
C     forms.for.
C     LENGTH... File
C     length.for.
C     EIGEN ... File
C     eigen.for.
C     HIVD2 ... File
C     means.for.
C     BLOCK,VELOC ... File
C     model.for.
C     PARM3 ... File
C     parm.for.
C     AP00,AP21 ... File
C     ap.for.
C     TRANS ... File
C     trans.for.
C     HDER... File hder.for.
C
C Common block /POINTC/ to store the results of complete ray tracing:
      INCLUDE 'pointc.inc'
C ...........................
C Common block /RTMAT/ to store the matrices of reflection-transmittion
C coefficients.
      INTEGER MRT,NRT
      PARAMETER (MRT=100)
      REAL PIRTR(2,2,MRT),PIRTI(2,2,MRT)
      COMMON/RTMAT/NRT,PIRTR,PIRTI
      SAVE /RTMAT/
C     MRT ... Dimension of arrays.
C     NRT ... Number of stored R-T matrices.
C     PIRTR .. Real part of R-T matrices.
C     PIRTI .. Imaginary part of R-T matrices.
C     The matrices are expressed in terms of polarization vectors.
C ...........................
C Common block /WASEP/ to store the input SEP parameters used
C in subroutine WAN and their subroutines.
      REAL CRTANI,DEGREE,SINGLF
      COMMON/WASEP/CRTANI,DEGREE,SINGLF
      SAVE /WASEP/
C     CRTANI ... Switch between isotropic and anisotropic ray tracing.
C     DEGREE ... Degree of the homogeneous Hamiltonian.
C     SINGLF ... Frequency of the prevailing-frequency approximation.
C.......................................................................
C
C     Auxiliary storage locations:
      INTEGER MPTS
      REAL PIGRA(2,2),PIGIA(2,2),PIGR(2,2),PIGI(2,2)
      REAL PIGRB(2,2),PIGIB(2,2)
      REAL DPIGRA(2,2),DPIGIA(2,2),DPIGR(2,2),DPIGI(2,2)
      REAL AR11,AI11,AR21,AI21,AR12,AI12,AR22,AI22
      REAL G11,G21,G12,G22,C,S,DPIG
      REAL PI
      PARAMETER (PI=3.1415926)
      REAL GREENA(32)
      REAL TTCOR,FREQ,GAMA,AUX0,AUX1,AUX2,AUX3,AUX4
      INTEGER NPTS
      INTEGER I,I1,I2,I3,J,IRT,ISINGL,NSINGL
      INTEGER NFFT,NF,KQIPV,KQITT,KQIRAY,NF1
      SAVE    NFFT,NF,KQIPV,KQITT,KQIRAY
      REAL DT,FMIN,FMAX,DF,OF
      SAVE DT,FMIN,FMAX,DF,OF
      INTEGER IPTTMP,NYFTMP,ISRFFT
      REAL YI1TMP,YF1TMP,STEP,TISO,DTLIN,DTAU(7),QICOR(16),EIGOUT(9)
      LOGICAL LPWAVE,LSWAVE
      CHARACTER*135 FORMQI
      CHARACTER*79  FORMEI
      CHARACTER*256 TXTERR
      DATA NFFT/0/
C                                                    
C     IPTS,RPTS... Quantities in the points on the ray:
C        IPTS(1,I)...  Index of the I-th point, zero for points
C                      added to the ray by interpolation.
C        RPTS(2,I)...  Travel time in I-th point.
C        RPTS(3-5,I).. Coordinates of the point.
C        RPTS(6-8,I).. Slowness vector in the point.
C        RPTS(9-11,I)  Polarization vector in the point.
C        IPTS(12,I)... Index of complex block.
C        RPTS(13,I)... Eigenvalue of Christoffel matrix in I-th point
C                      corresponding to qP wave.
C        RPTS(14,I)... Eigenvalue of Christoffel matrix in I-th point
C                      corresponding to qS1 (faster) wave.
C        RPTS(15,I)... Eigenvalue of Christoffel matrix in I-th point
C                      corresponding to qS2 (slower) wave.
C        RPTS(16-19,I) Eigenvectors projected to the plane defined by
C                      polarization vectors.
C        RPTS(20,I)... Angular difference of eigenvectors of Christoffel
C                      matrix corresponding to qS1 wave in I-th and
C                      in (I-1)-st point.
C        RPTS(21,I)... Time integral of the difference between qS1 and
C                      qS2 eigenvalues along the ray from the (I-1)-st
C                      to the I-th point.
C        RPTS(22-24,I) Derivatives of the velocity.
C        RPTS(25-33,I) Eigenvectors of qS1, qS2 and P wave.
C        RPTS(34,I)... S-wave velocity.
C        RPTS(35-38,I) Matrix of geometrical spreading.
C        RPTS(39-42,I) Transformation matrix P.
C        IPTS(43,I)... Index of surface, zero inside a complex block.
C        RPTS(44-52,I) Covariant components of the basis vectors
C                      of the ray-centred coordinate system.
C        RPTS(53-61,I) Contravariant components of the basis vectors
C                      of the ray-centred coordinate system.
C
C-----------------------------------------------------------------------
C
      IF (NFFT.EQ.0) THEN
        CALL RSEP3I('NFFT',NFFT,1)
        CALL RSEP3R('DT  ',DT  ,1.)
        CALL RSEP3R('FMIN',FMIN,0.)
        CALL RSEP3R('FMAX',FMAX,0.5/DT)
        CALL RSEP3R('DF  ',DF  ,1./(FLOAT(NFFT)*DT))
        CALL RSEP3R('OF  ',OF  ,DF*NINT(FMIN/DF))
        CALL RSEP3I('NF  ',NF  ,NINT((FMAX-OF)/DF)+1)
        CALL RSEP3I('QIPV',KQIPV,0)
        CALL RSEP3I('QITT',KQITT,0)
        FORMQI(1:7)='(1I6,A,'
        FORMEI(1:7)='(1I6,A,'
        CALL RSEP3I('QIRAY',KQIRAY,0)
        CALL RSEP3R('CRTANI',CRTANI, 0.)
        CALL RSEP3R('DEGREE',DEGREE,-1.)
        CALL RSEP3R('SINGLF',SINGLF, 0.)
      ENDIF
C
C     Reading in the quantities stored in individual points on
C     a ray, computing all the auxiliary quantities in the points,
C     and, if necessary, adding new points on the ray by interpolation:
      MPTS=MAXPTS/NQPTS
      NPTS=0
      CALL WAPTS(LU1,LU2,LU3,NF,OF,DF,RPTS,IPTS,NQPTS,MPTS,NPTS)
      IF (NPTS.EQ.0) THEN
C       End of rays:
        NQ=0
        RETURN
      ENDIF
C
C     Replacing the reference polarization vectors in common block
C     /POINTC/ by the eigenvectors of the Christoffel matrix:
      DO 9 I=1,9
        PVI(I)=RPTS(24+I,1)
        PV (I)=RPTS(24+I,NPTS)
    9 CONTINUE
C
C     Output to the file EIGLST:
      IF (LUEIG.NE.0) THEN
        WRITE(LUEIG,'(A,1I6,A)') 'IREC =',IREC,' /'
        DO 10, I1=1,NPTS
          EIGOUT(1)=RPTS( 2,I1)
          EIGOUT(2)=RPTS(25,I1)
          EIGOUT(3)=RPTS(26,I1)
          EIGOUT(4)=RPTS(27,I1)
          EIGOUT(5)=RPTS(28,I1)
          EIGOUT(6)=RPTS(29,I1)
          EIGOUT(7)=RPTS(30,I1)
          EIGOUT(8)=RPTS(14,I1)
          EIGOUT(9)=RPTS(15,I1)
          CALL FORM2(9,EIGOUT,EIGOUT,FORMEI(8:79))
          WRITE(LUEIG,FORMEI) IPTS(1,I1),(' ',EIGOUT(I),I=1,9),' /'
   10   CONTINUE
      ENDIF
C
      TISO=Y(1)
C
C     Computing the values of travel time corrections along the ray:
      LPWAVE=.TRUE.
      LSWAVE=.TRUE.
      TTCOR=0.
      DO 30, I1=1,NPTS-1
        I2=I1+1
        IF (IPTS(12,I2).GT.0) THEN
C         P wave:
          TTCOR=TTCOR+(1./SQRT(RPTS(13,I2))+1./SQRT(RPTS(13,I1)))*0.5*
     *          (RPTS(2,I2)-RPTS(2,I1))
          LSWAVE=.FALSE.
C         Updating common block /POINTC/ for subroutine AP21:
          NPVI=3
          NPV=3
        ELSE
C         S wave:
          IF(KQITT.LE.0) THEN
            TTCOR=TTCOR+(0.25/SQRT(RPTS(14,I2))+0.25/SQRT(RPTS(15,I2))
     *                  +0.25/SQRT(RPTS(14,I1))+0.25/SQRT(RPTS(15,I1)))
     *                 *(RPTS(2,I2)-RPTS(2,I1))
          ELSE
            TTCOR=TTCOR+(1.5-0.125*(RPTS(14,I2)+RPTS(15,I2)
     *                             +RPTS(14,I1)+RPTS(15,I1)))
     *                 *(RPTS(2,I2)-RPTS(2,I1))
          END IF
          LPWAVE=.FALSE.
C         Updating common block /POINTC/ for subroutine AP21:
          NPVI=6
          NPV=6
        ENDIF
   30 CONTINUE
C
      IF (LPWAVE) THEN
C       P wave along the whole ray:
        Y(1)=TTCOR+YI(1)
        CALL AP21(KQIPV,GREEN)
        NQ=32
        GOTO 91
      ENDIF
C
      IF (LSWAVE.AND.((LUQI.NE.0).OR.(KQIRAY.NE.0))) THEN
C       S wave along the whole ray,
C       calculating second-order travel-time perturbations:
        DTAU(1)=0.
        DTAU(3)=0.
        IPTTMP=IPT
        YI1TMP=YI(1)
        NYFTMP=NYF
        ISRFFT=ISRFF
        YF1TMP=YF(1)
        YI(1) =RPTS(2,1)
        NYF   =1
        DTLIN=0.
        DO 95, I1=1,NPTS
          IF (IPTS(1,I1).EQ.2) THEN
            STEP=RPTS(2,I1)-RPTS(2,1)
          ENDIF
          IF (IPTS(1,I1).EQ.0) THEN
            STEP=0.
            GOTO 96
          ENDIF
  95    CONTINUE
  96    CONTINUE
        DO 100, I1=1,NPTS
          DTLIN=DTLIN+RPTS(21,I1)
          IPT   =IPTS(1,I1)-1
          ISRFF =IPTS(43,I1)
          YF(1) =RPTS(2,I1)
          IF (KQIRAY.NE.0) THEN
            RPTS(21,I1)=RPTS(21,I1)+DTAU(1)-DTAU(3)
          ENDIF
          CALL WAPERT(RPTS(1,I1),IPTS(1,I1),NQPTS,STEP,DTAU)
          IF (KQIRAY.NE.0) THEN
            RPTS(21,I1)=RPTS(21,I1)+DTAU(3)-DTAU(1)
          ENDIF
 100    CONTINUE
C
        IF (LUQI.NE.0) THEN
          DTLIN=DTLIN/2.
          QICOR(1)=TISO
          QICOR(2)=TTCOR-DTLIN-TISO
          QICOR(3)=TTCOR+DTLIN-TISO
          QICOR(4)=DTAU(1)
          QICOR(5)=DTAU(3)
          QICOR(6)=TTCOR-DTLIN+DTAU(1)
          QICOR(7)=QICOR(6)
          QICOR(8)=TTCOR+DTLIN+DTAU(3)
          QICOR(9)=QICOR(8)
          IF(CRTANI.EQ.0.) THEN
            QICOR(10)=TTCOR-TISO
            QICOR(11)=0.25*(DTAU(1)+DTAU(2)+DTAU(3))
            QICOR(12)=TISO+QICOR(10)+QICOR(11)
            QICOR(13)=QICOR(12)
            QICOR(14)=0.25*(DTAU(1)-DTAU(2)+DTAU(3))
            QICOR(15)=QICOR(14)
            QICOR(16)=DTAU(1)-DTAU(2)+DTAU(3)
          ELSE
            QICOR(10)=DTAU(7)-TISO
            QICOR(11)=DTAU(6)
            QICOR(12)=DTAU(7)+DTAU(6)
            QICOR(13)=QICOR(12)
            QICOR(14)=DTAU(1)-DTAU(4)+DTAU(6)
            QICOR(15)=DTAU(3)-DTAU(5)+DTAU(6)
            QICOR(16)=DTAU(1)-DTAU(2)+DTAU(3)
          END IF
          CALL FORM2(16,QICOR,QICOR,FORMQI(8:135))
          WRITE(LUQI,FORMQI) IREC,(' ',QICOR(I),I=1,16),' /'
        ENDIF
C
        IPT   =IPTTMP
        YI(1) =YI1TMP
        NYF   =NYFTMP
        ISRFF =ISRFFT
        YF(1) =YF1TMP
C
        IF (KQIRAY.NE.0) THEN
          TTCOR=TTCOR+0.5*(DTAU(3)+DTAU(1))
        ENDIF
      ELSEIF ((.NOT.LSWAVE).AND.((LUQI.NE.0).OR.(KQIRAY.NE.0))) THEN
C       WAN-20
        WRITE(TXTERR,'(2A,1I6,A)')
     *  'WAN-20: Second-order perturbations of anisotropic travel ',
     *  'times may not be calculated for the ray ',IRAY,
     *  ', the ray crossed an interface.'
        CALL WARN(TXTERR(1:LENGTH(TXTERR)))
C       The perturbations may be calculated only along the rays
C       which do not cross any interface.
      ENDIF
C
C     S wave in some part of the ray:
C     Loop over the frequencies:
      NF1=NF-1
      IF (SINGLF.NE.0.) THEN
        NF1=0
      ENDIF
      DO 90, I2=0,NF1
        FREQ=OF+I2*DF
        IF (SINGLF.NE.0.) THEN
          FREQ=SINGLF
          DPIGR(1,1)=0.
          DPIGR(2,1)=0.
          DPIGR(1,2)=0.
          DPIGR(2,2)=0.
          DPIGI(1,1)=0.
          DPIGI(2,1)=0.
          DPIGI(1,2)=0.
          DPIGI(2,2)=0.
        ENDIF
        PIGR(1,1)=1.
        PIGR(2,1)=0.
        PIGR(1,2)=0.
        PIGR(2,2)=1.
        PIGI(1,1)=0.
        PIGI(2,1)=0.
        PIGI(1,2)=0.
        PIGI(2,2)=0.
        IRT=0
C       Computing the propagator matrix PiGe along the ray:
C       Loop along points on the ray:
        DO 40, I1=2,NPTS
          I3=I1-1
          IF (RPTS(2,I1).EQ.RPTS(2,I3)) THEN
C           The ray crossed an interface:
C           WAN-17
            CALL ERROR('WAN-17: Ray crossed an interface.')
C           Not yet coded.
*           IRT=IRT+1
*           IF (IRT.GT.NRT) THEN
C             WAN-XX
*             CALL ERROR('WAN-XX: Disorder in number of R-T.')
C             This error should not appear. Number of reflections-
C             -transmittions should be the same in WAREAD and here.
*           ENDIF
C           Transforming R-T matrix from polarization vectors to
C           eigenvectors:
*           G11=RPTS(16,I3)
*           G21=RPTS(17,I3)
*           G12=RPTS(18,I3)
*           G22=RPTS(19,I3)
*           AR11=PIRTR(1,1,IRT)*G11+PIRTR(1,2,IRT)*G21
*           AI11=PIRTI(1,1,IRT)*G11+PIRTI(1,2,IRT)*G21
*           AR21=PIRTR(2,1,IRT)*G11+PIRTR(2,2,IRT)*G21
*           AI21=PIRTI(2,1,IRT)*G11+PIRTI(2,2,IRT)*G21
*           AR12=PIRTR(1,1,IRT)*G12+PIRTR(1,2,IRT)*G22
*           AI12=PIRTI(1,1,IRT)*G12+PIRTI(1,2,IRT)*G22
*           AR22=PIRTR(2,1,IRT)*G12+PIRTR(2,2,IRT)*G22
*           AI22=PIRTI(2,1,IRT)*G12+PIRTI(2,2,IRT)*G22
*           G11=RPTS(16,I1)
*           G21=RPTS(18,I1)
*           G12=RPTS(17,I1)
*           G22=RPTS(19,I1)
*           PIGRA(1,1)=G11*AR11+G12*AR21
*           PIGIA(1,1)=G11*AI11+G12*AI21
*           PIGRA(2,1)=G21*AR11+G22*AR21
*           PIGIA(2,1)=G21*AI11+G22*AI21
*           PIGRA(1,2)=G11*AR12+G12*AR22
*           PIGIA(1,2)=G11*AI12+G12*AI22
*           PIGRA(2,2)=G21*AR12+G22*AR22
*           PIGIA(2,2)=G21*AI12+G22*AI22
          ELSE
C         Element of the ray inside a complex block:
            IF (IPTS(12,I1).LT.0) THEN
C             S wave:
              GAMA=PI*FREQ*RPTS(21,I1)
              AUX0=SQRT(RPTS(20,I1)**2 + GAMA**2)
              AUX1=COS(AUX0)
              IF (AUX0.EQ.0.) THEN
                AUX2=1.
                AUX3=-1./3.
              ELSE
                AUX2=SIN(AUX0)/AUX0
                AUX3=(COS(AUX0)-AUX2)/AUX0**2
              ENDIF
C             Matrix for this step along the ray:
C             R22, Eq.(49):
              PIGRA(1,1)= AUX1
              PIGRA(2,1)=-RPTS(20,I1)*AUX2
              PIGRA(1,2)=-PIGRA(2,1)
              PIGRA(2,2)= AUX1
              PIGIA(1,1)=-GAMA*AUX2
              PIGIA(2,1)= 0.
              PIGIA(1,2)= 0.
              PIGIA(2,2)=-PIGIA(1,1)
              IF (SINGLF.NE.0.) THEN
C               R22, Eq.(56):
                AUX4=GAMA/(2.*PI*FREQ)
                DPIGRA(1,1)=-AUX2*GAMA*AUX4
                DPIGRA(2,1)=-AUX3*RPTS(20,I1)*GAMA*AUX4
                DPIGRA(1,2)=-DPIGRA(2,1)
                DPIGRA(2,2)= DPIGRA(1,1)
                DPIGIA(1,1)=-(AUX3*GAMA**2+AUX2)*AUX4
                DPIGIA(2,1)= 0.
                DPIGIA(1,2)= 0.
                DPIGIA(2,2)=-DPIGIA(1,1)
              ENDIF
            ELSE
C             P wave:
              GOTO 39
            ENDIF
          ENDIF
C         Matrix for all the steps along the ray to current point:
          IF (SINGLF.NE.0.) THEN
C           R22, Eq.(53):
            CALL WAMAT(PIGRA,PIGIA,DPIGR,DPIGI)
            PIGRB(1,1)=PIGR(1,1)
            PIGRB(2,1)=PIGR(2,1)
            PIGRB(1,2)=PIGR(1,2)
            PIGRB(2,2)=PIGR(2,2)
            PIGIB(1,1)=PIGI(1,1)
            PIGIB(2,1)=PIGI(2,1)
            PIGIB(1,2)=PIGI(1,2)
            PIGIB(2,2)=PIGI(2,2)
            CALL WAMAT(DPIGRA,DPIGIA,PIGRB,PIGIB)
            DPIGR(1,1)=DPIGR(1,1)+PIGRB(1,1)
            DPIGR(2,1)=DPIGR(2,1)+PIGRB(2,1)
            DPIGR(1,2)=DPIGR(1,2)+PIGRB(1,2)
            DPIGR(2,2)=DPIGR(2,2)+PIGRB(2,2)
            DPIGI(1,1)=DPIGI(1,1)+PIGIB(1,1)
            DPIGI(2,1)=DPIGI(2,1)+PIGIB(2,1)
            DPIGI(1,2)=DPIGI(1,2)+PIGIB(1,2)
            DPIGI(2,2)=DPIGI(2,2)+PIGIB(2,2)
          ENDIF
C         R22, Eq.(47):
          CALL WAMAT(PIGRA,PIGIA,PIGR,PIGI)
   39     CONTINUE
   40   CONTINUE
C
C       Prevailing-frequency approximation
        NSINGL=0
        IF (SINGLF.NE.0.) THEN
C         R22, Eq.(21):
          DPIG=DPIGR(1,1)**2+DPIGI(1,1)**2+DPIGR(2,1)**2+DPIGI(2,1)**2
          DPIG=SQRT(DPIG)
          PIGIA(1,1)=PIGI(1,1)
          PIGRA(1,1)=PIGR(1,1)
          PIGIA(2,1)=PIGI(2,1)
          PIGRA(2,1)=PIGR(2,1)
          PIGIA(1,2)=PIGI(1,2)
          PIGRA(1,2)=PIGR(1,2)
          PIGIA(2,2)=PIGI(2,2)
          PIGRA(2,2)=PIGR(2,2)
          NSINGL=32
        ENDIF
C
C       Computing the Green function:
        DO 80 ISINGL=0,NSINGL,32
C         Transforming the propagator matrix from eigenvectors
C         to polarization vectors:
          Y(1)=TTCOR+YI(1)
          IF (SINGLF.NE.0.) THEN
            IF (ISINGL.EQ.0) THEN
C             R22, Eq.(41):
              Y(1)=TTCOR+YI(1)-DPIG
C             R22, Eq.(22):
              IF (DPIG.NE.0) THEN
                PIGIB(1,1)=(PIGIA(1,1)+DPIGR(1,1)/DPIG)/2.
                PIGRB(1,1)=(PIGRA(1,1)-DPIGI(1,1)/DPIG)/2.
                PIGIB(2,1)=(PIGIA(2,1)+DPIGR(2,1)/DPIG)/2.
                PIGRB(2,1)=(PIGRA(2,1)-DPIGI(2,1)/DPIG)/2.
                PIGIB(1,2)=(PIGIA(1,2)+DPIGR(1,2)/DPIG)/2.
                PIGRB(1,2)=(PIGRA(1,2)-DPIGI(1,2)/DPIG)/2.
                PIGIB(2,2)=(PIGIA(2,2)+DPIGR(2,2)/DPIG)/2.
                PIGRB(2,2)=(PIGRA(2,2)-DPIGI(2,2)/DPIG)/2.
              ELSE
                PIGIB(1,1)=PIGIA(1,1)
                PIGRB(1,1)=PIGRA(1,1)
                PIGIB(2,1)=PIGIA(2,1)
                PIGRB(2,1)=PIGRA(2,1)
                PIGIB(1,2)=0.
                PIGRB(1,2)=0.
                PIGIB(2,2)=0.
                PIGRB(2,2)=0.
              ENDIF
              C=COS(2.*PI*SINGLF*DPIG)
              S=SIN(2.*PI*SINGLF*DPIG)
            ELSE
C             R22, Eq.(42):
              Y(1)=TTCOR+YI(1)+DPIG
C             R22, Eq.(23):
              IF (DPIG.NE.0) THEN
                PIGIB(1,1)=(PIGIA(1,1)-DPIGR(1,1)/DPIG)/2.
                PIGRB(1,1)=(PIGRA(1,1)+DPIGI(1,1)/DPIG)/2.
                PIGIB(2,1)=(PIGIA(2,1)-DPIGR(2,1)/DPIG)/2.
                PIGRB(2,1)=(PIGRA(2,1)+DPIGI(2,1)/DPIG)/2.
                PIGIB(1,2)=(PIGIA(1,2)-DPIGR(1,2)/DPIG)/2.
                PIGRB(1,2)=(PIGRA(1,2)+DPIGI(1,2)/DPIG)/2.
                PIGIB(2,2)=(PIGIA(2,2)-DPIGR(2,2)/DPIG)/2.
                PIGRB(2,2)=(PIGRA(2,2)+DPIGI(2,2)/DPIG)/2.
              ELSE
                PIGIB(1,1)=0.
                PIGRB(1,1)=0.
                PIGIB(2,1)=0.
                PIGRB(2,1)=0.
                PIGIB(1,2)=PIGIA(1,2)
                PIGRB(1,2)=PIGRA(1,2)
                PIGIB(2,2)=PIGIA(2,2)
                PIGRB(2,2)=PIGRA(2,2)
              ENDIF
              S=-S
            ENDIF
C           R22, Eq.(43):
            PIGI(1,1)=C*PIGIB(1,1)+S*PIGRB(1,1)
            PIGR(1,1)=C*PIGRB(1,1)-S*PIGIB(1,1)
            PIGI(2,1)=C*PIGIB(2,1)+S*PIGRB(2,1)
            PIGR(2,1)=C*PIGRB(2,1)-S*PIGIB(2,1)
            PIGI(1,2)=C*PIGIB(1,2)+S*PIGRB(1,2)
            PIGR(1,2)=C*PIGRB(1,2)-S*PIGIB(1,2)
            PIGI(2,2)=C*PIGIB(2,2)+S*PIGRB(2,2)
            PIGR(2,2)=C*PIGRB(2,2)-S*PIGIB(2,2)
          ENDIF
C
          IF(KQIPV.EQ.0) THEN
            Y(28)=PIGR(1,1)
            Y(29)=PIGI(1,1)
            Y(30)=PIGR(2,1)
            Y(31)=PIGI(2,1)
            Y(32)=PIGR(1,2)
            Y(33)=PIGI(1,2)
            Y(34)=PIGR(2,2)
            Y(35)=PIGI(2,2)
          ELSE
            G11=RPTS(16,1)
            G12=RPTS(17,1)
            G21=RPTS(18,1)
            G22=RPTS(19,1)
            AR11=PIGR(1,1)*G11+PIGR(1,2)*G21
            AI11=PIGI(1,1)*G11+PIGI(1,2)*G21
            AR21=PIGR(2,1)*G11+PIGR(2,2)*G21
            AI21=PIGI(2,1)*G11+PIGI(2,2)*G21
            AR12=PIGR(1,1)*G12+PIGR(1,2)*G22
            AI12=PIGI(1,1)*G12+PIGI(1,2)*G22
            AR22=PIGR(2,1)*G12+PIGR(2,2)*G22
            AI22=PIGI(2,1)*G12+PIGI(2,2)*G22
            G11=RPTS(16,NPTS)
            G21=RPTS(17,NPTS)
            G12=RPTS(18,NPTS)
            G22=RPTS(19,NPTS)
            Y(28)=G11*AR11+G12*AR21
            Y(29)=G11*AI11+G12*AI21
            Y(30)=G21*AR11+G22*AR21
            Y(31)=G21*AI11+G22*AI21
            Y(32)=G11*AR12+G12*AR22
            Y(33)=G11*AI12+G12*AI22
            Y(34)=G21*AR12+G22*AR22
            Y(35)=G21*AI12+G22*AI22
          ENDIF
C
          CALL AP21(KQIPV,GREENA)
          IF (I2.EQ.0) THEN
            DO 50 I=1,14
              GREEN(ISINGL+I)=GREENA(I)
   50       CONTINUE
            NQ=ISINGL+14
          ENDIF
          J=I2*18
          DO 60 I=15,32
            GREEN(ISINGL+J+I)=GREENA(I)
   60     CONTINUE
          NQ=NQ+18
   80   CONTINUE
C
   90 CONTINUE
   91 CONTINUE
C
      IF (KQIPV.EQ.1) THEN
C       Calculating the error of the QI projection of polariz. vectors
        G11=RPTS(16,1)
        G21=RPTS(17,1)
        G12=RPTS(18,1)
        G22=RPTS(19,1)
        ERR(1)=AMAX1(SQRT(2.*ABS(1.-ABS(G11*G22-G12*G21))),ERR(1))
        G11=RPTS(16,NPTS)
        G21=RPTS(17,NPTS)
        G12=RPTS(18,NPTS)
        G22=RPTS(19,NPTS)
        ERR(2)=AMAX1(SQRT(2.*ABS(1.-ABS(G11*G22-G12*G21))),ERR(2))
      ENDIF
C
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE WAPTS(LU1,LU2,LU3,NF,OF,DF,RPTS,IPTS,NQPTS,MPTS,NPTS)
C
C-----------------------------------------------------------------------
C Subroutine to read in the quantities stored in individual points
C on the ray, to compute all the auxiliary quantities in the points,
C and, if necessary, to add new points on the rays by interpolation.
C Reading the files with points on the rays is done by a simple
C invocation of subroutine WAREAD.
C Computation of elastic parameters is completed by invocation
C of subroutine PARM3 of file 'parm.for'.
C Then Christoffel matrix is evaluated and its eigenvalues and
C eigenvectors are computed by subroutine EIGEN of file 'eigen.for'.
C In the next step the angular difference DELTFI is computed for
C each subinterval along the ray.
C If the value of DELTFI is greater than prescribed limit, new points
C are added on the ray.  Coordinates, slowness vector and polarization
C vector are interpolated by subroutine HIVD2 of the file 'means.for'.
C Material parameters, Christoffel matrix and their eigenvalues and
C eigenvectors are calculated in the new point. Matrix of geometrical
C spreading, transformation matrix P and covariant and contravariant
C components of the basis vectors of the RCCS are interpolated linearly.
C
      INTEGER LU1,LU2,LU3,NF
      REAL OF,DF
      INTEGER NPTS,NQPTS,MPTS
      REAL    RPTS(NQPTS,MPTS)
      INTEGER IPTS(NQPTS,MPTS)
C Input:
C     LU1 ... Number of the logical unit connected to the CRT output
C             file with quantities along rays.
C     LU2 ... Number of the logical unit connected to the CRT output
C             file with quantities at the initial points of rays.
C     LU3 ... Number of the logical unit connected to the CRT output
C             file with quantities at the storing surface.
C     NF  ... Number of frequencies.
C     OF  ... First frequency.
C     DF  ... Frequency increment.
C     NQPTS.. Dimension of arrays RPTS and IPTS.
C     MPTS... Dimension of arrays RPTS and IPTS.
C Output:
C     RPTS,IPTS ... The arrays are filled with all the quantities for
C                  single two-point ray during one invocation of this
C                  subroutine:
C        IPTS(1,I) ... Index of the I-th point, zero for points
C                      added to the ray by interpolation.
C        RPTS(2,I) ... Travel time in I-th point.
C        RPTS(3-5,I) . Coordinates of the point.
C        RPTS(6-8,I) . Slowness vector in the point.
C        RPTS(9-11,I)  Polarization vector in the point.
C        IPTS(12,I) .. Index of complex block.
C        RPTS(13,I)... Eigenvalue of Christoffel matrix in I-th point
C                      corresponding to qP wave.
C        RPTS(14,I)... Eigenvalue of Christoffel matrix in I-th point
C                      corresponding to qS1 (faster) wave.
C        RPTS(15,I)... Eigenvalue of Christoffel matrix in I-th point
C                      corresponding to qS2 (slower) wave.
C        RPTS(16-19,I) Eigenvectors projected to the plane defined by
C                      polarization vectors.
C        RPTS(20,I) .. Angular difference of eigenvectors of Christoffel
C                      matrix corresponding to qS1 wave in I-th and
C                      in (I-1)-st point.
C        RPTS(21,I) .. Time integral of the difference between qS1 and
C                      qS2 eigenvalues along the ray from the (I-1)-st
C                      to the I-th point.
C        RPTS(22-24,I) Derivatives of the velocity.
C        RPTS(25-33,I) Eigenvectors of qS1, qS2 and P wave.
C        RPTS(34,I)... S-wave velocity.
C        RPTS(35-38,I) Matrix of geometrical spreading.
C        RPTS(39-42,I) Transformation matrix P.
C        IPTS(43,I)... Index of surface, zero inside a complex block.
C        RPTS(44-52,I) Covariant components of the basis vectors
C                      of the ray-centred coordinate system.
C        RPTS(53-61,I) Contravariant components of the basis vectors
C                      of the ray-centred coordinate system.
C     NPTS ... Number of points stored in RPTS (IPTS).
C
C ...........................
C Common block /WASEP/ to store the input SEP parameters used
C in subroutine WAN and their subroutines.
      REAL CRTANI,DEGREE,SINGLF
      COMMON/WASEP/CRTANI,DEGREE,SINGLF
      SAVE /WASEP/
C     CRTANI ... Switch between isotropic and anisotropic ray tracing.
C     DEGREE ... Degree of the homogeneous Hamiltonian.
C     SINGLF ... Frequency of the prevailing-frequency approximation.
C.......................................................................
C External functions required:
      EXTERNAL WAVECR,WAANGL
      REAL WAANGL
C.......................................................................
C
C     Auxiliary storage locations:
      INTEGER NPTSE,LPTISO
      INTEGER KQIDT
      SAVE    KQIDT
      REAL UP(10),US(10),RHO,QP,QS,VP,VS,VD(10),QL,AA(10,21),QQ(21)
      REAL EE(9),DER(9)
      REAL ERRWAN,DEFI,DELTFI,PI
      SAVE ERRWAN
      PARAMETER (PI=3.1415926)
      INTEGER I1,I2
      INTEGER MQANT,MNEWP,NNEWP
      PARAMETER (MNEWP=50)
      PARAMETER (MQANT=56)
      REAL ROLDP(MQANT),RNEWP(MQANT,MNEWP)
      INTEGER KNEWP(MNEWP)
      REAL TTOLD,XOLD(3),DXOLD(3),POLD(3),DPOLD(3),EOLD(3),DEOLD(3)
      REAL TTNEW,XNEW(3),DXNEW(3),PNEW(3),DPNEW(3),ENEW(3),DENEW(3),AUX
C
C     ROLDP(I),RNEWP(I,J)
C           I=1    ...   Travel time.
C             2-4  ...   Coordinates.
C             5-7  ...   Slowness vector.
C             8-10 ...   Polarization vector E1.
C             11   ...   qP eigenvalue of Christoffel matrix.
C             12   ...   qS1 eigenvalue of Christoffel matrix.
C             13   ...   qS2 eigenvalue of Christoffel matrix.
C             14   ...   First component of the qS1 eigenvector
C                  projected to the plane perpendicular to the ray,
C                  defined by two basis vectors of ray-centered
C                  coordinate system.
C             15   ...   Second component of the qS1 eigenvector.
C             16   ...   First component of the qS2 eigenvector.
C             17   ...   Second component of the qS2 eigenvector.
C             18-26...   Eigenvectors of qS1, qS2 and P wave.
C             27   ...   S-wave velocity.
C             28-30...   Derivatives of S-wave velocity.
C             31-34...   Matrix of geometrical spreading.
C             34-38...   Transformation matrix P.
C             39-56...   Covariant and contravariant components of
C                        the basis vectors of the RCCS.
C     KNEWP...Division index of points in RNEWP.
C     AA..    Values, first and second partial derivatives of real
C             parts of 21 reduced (divided by the density) elastic
C             parameters.  The order of the value, first and second
C             partial derivatives of each parameter Aij is:
C             Aij,Aij1,Aij2,Aij3,Aij11,Aij12,Aij22,Aij13,Aij23,Aij33.
C             The order of parameters (second array index) is:
C             A11,A12,A22,A13,A23,A33,A14,A24,A34,A44,A15,A25,A35,A45,
C             A55,A16,A26,A36,A46,A56,A66.
C     RHO...  Density at the given point.
C     QQ ...  Imaginary parts of 21 reduced elastic parameters at the
C             given point, ordered as
C             Q11,Q12,Q22,Q13,Q23,Q33,Q14,Q24,Q34,Q44,Q15,Q25,Q35,Q45,
C             Q55,Q16,Q26,Q36,Q46,Q56,Q66.
C     EE ...  Eigenvectors of the Christoffel matrix.
C     DER ... Derivatives dx/dt (anisotropic).
C     ERRWAN..Maximum error in propagator matrix along whole ray.
C     DEFI .. Maximum allowed angular change for eigenvectors of
C             Christoffel matrix between neighboring points on the ray.
C     TTOLD...Travel time in the old point stored for interpolation.
C     XOLD ...Coordinates of the old point stored for interpolation.
C     DXOLD ..Derivatives in the old point stored for interpolation.
C     POLD ...Slowness vector in the old point stored for interpolation.
C     DPOLD ..Derivatives of the slowness vector in the old point.
C     EOLD ...Polarization vector in the old point.
C     DEOLD ..Derivatives of the polarization vector in the old point.
C     TTNEW,XNEW,DXNEW,PNEW,DPNEW,ENEW,DENEW ... The same quantities
C             in the new point. Here old point and new point mean two
C             consecutive points in which the output of CRT was
C             recorded.
C     NPTSE...Number of points along the ray, where all the quantities
C             have been checked.
C     LPTISO..LPTISO=1 in case that all already checked points were
C             isotropic.
C
C-----------------------------------------------------------------------
C
      CALL RSEP3R('ERRWAN',ERRWAN,0.001)
      CALL RSEP3I('QIDT',KQIDT,0)
C
C     Reading the quantities computed by CRT for one ray:
      CALL WAREAD(LU1,LU2,LU3,RPTS,IPTS,NQPTS,MPTS,NPTS)
      IF (NPTS.EQ.0)
C       End of rays:
     *  RETURN
      IF (NPTS.LT.2) THEN
C       WAN-01
        CALL ERROR('WAN-01: A ray formed by single point')
C       This error should not appear.
C       Each ray should be represented by at least two points.
      ENDIF
C
C
C     Reading the material parameters in the initial point:
      CALL PARM3(IPTS(12,1),RPTS(3,1),AA,RHO,QQ)
C     Computing eigenvectors and eigenvalues of
C     the Christoffel matrix in the initial point:
      CALL WACHRI(RPTS(6,1),RPTS(7,1),RPTS(8,1),
     * AA(1,1),AA(1,2),AA(1,3),AA(1,4),AA(1,5),AA(1,6),AA(1,7),
     * AA(1,8),AA(1,9),AA(1,10),AA(1,11),AA(1,12),AA(1,13),AA(1,14),
     * AA(1,15),AA(1,16),AA(1,17),AA(1,18),AA(1,19),AA(1,20),AA(1,21),
     * RPTS(13,1),RPTS(14,1),RPTS(15,1),EE,DER)
C     Storing the eigenvectors of the Christoffel matrix
      RPTS(25,1)=EE(4)
      RPTS(26,1)=EE(5)
      RPTS(27,1)=EE(6)
      RPTS(28,1)=EE(7)
      RPTS(29,1)=EE(8)
      RPTS(30,1)=EE(9)
      RPTS(31,1)=EE(1)
      RPTS(32,1)=EE(2)
      RPTS(33,1)=EE(3)
      IF (ABS(RPTS(15,1)-RPTS(14,1)).LT.0.00001) THEN
C       Isotropic case, no projection of eigenvectors:
        RPTS(16,1)=1.
        RPTS(17,1)=0.
        RPTS(18,1)=0.
        RPTS(19,1)=1.
        LPTISO=1
      ELSE
C       Projecting the qS eigenvectors to the plane
C       perpendicular to the ray:
        CALL WAPROJ(RPTS(6,1),RPTS(7,1),RPTS(8,1),
     *              RPTS(9,1),RPTS(10,1),RPTS(11,1),
     *              EE(4),EE(5),EE(6),EE(7),EE(8),EE(9),
     *              RPTS(16,1),RPTS(17,1),RPTS(18,1),RPTS(19,1))
        LPTISO=0
      ENDIF
      RPTS(20,1)=0.
      RPTS(21,1)=0.
C
C     Quantities for future possible interpolation:
      TTOLD=RPTS(2,1)
C     Coordinates:
      XOLD(1)=RPTS(3,1)
      XOLD(2)=RPTS(4,1)
      XOLD(3)=RPTS(5,1)
C     dx/dt=p*v**2=p/sum(pi pi):
      AUX=(RPTS(6,1)**2+RPTS(7,1)**2+RPTS(8,1)**2)
      DXOLD(1)=RPTS(6,1)/AUX
      DXOLD(2)=RPTS(7,1)/AUX
      DXOLD(3)=RPTS(8,1)/AUX
C     Slowness:
      POLD(1)=RPTS(6,1)
      POLD(2)=RPTS(7,1)
      POLD(3)=RPTS(8,1)
C     dpi/dt=-1/v*dv/dxi:
      AUX=SQRT(AUX)
      DPOLD(1)=-RPTS(22,1)/AUX
      DPOLD(2)=-RPTS(23,1)/AUX
      DPOLD(3)=-RPTS(24,1)/AUX
C     Polarization:
      EOLD(1)=RPTS( 9,1)
      EOLD(2)=RPTS(10,1)
      EOLD(3)=RPTS(11,1)
C     dei/dt=v*sum(dv/dxk ek)*pi
      AUX=AUX*
     *(RPTS(22,1)*RPTS(9,1)+RPTS(23,1)*RPTS(10,1)+RPTS(24,1)*RPTS(11,1))
      DEOLD(1)=AUX*RPTS(6,1)
      DEOLD(2)=AUX*RPTS(7,1)
      DEOLD(3)=AUX*RPTS(8,1)
C
      ROLDP( 1)=RPTS( 2,1)
      ROLDP( 2)=RPTS( 3,1)
      ROLDP( 3)=RPTS( 4,1)
      ROLDP( 4)=RPTS( 5,1)
      ROLDP( 5)=RPTS( 6,1)
      ROLDP( 6)=RPTS( 7,1)
      ROLDP( 7)=RPTS( 8,1)
      ROLDP( 8)=RPTS( 9,1)
      ROLDP( 9)=RPTS(10,1)
      ROLDP(10)=RPTS(11,1)
      ROLDP(11)=RPTS(13,1)
      ROLDP(12)=RPTS(14,1)
      ROLDP(13)=RPTS(15,1)
      ROLDP(14)=RPTS(16,1)
      ROLDP(15)=RPTS(17,1)
      ROLDP(16)=RPTS(18,1)
      ROLDP(17)=RPTS(19,1)
      ROLDP(18)=RPTS(25,1)
      ROLDP(19)=RPTS(26,1)
      ROLDP(20)=RPTS(27,1)
      ROLDP(21)=RPTS(28,1)
      ROLDP(22)=RPTS(29,1)
      ROLDP(23)=RPTS(30,1)
      ROLDP(24)=RPTS(31,1)
      ROLDP(25)=RPTS(32,1)
      ROLDP(26)=RPTS(33,1)
      ROLDP(27)=RPTS(34,1)
      ROLDP(28)=RPTS(22,1)
      ROLDP(29)=RPTS(23,1)
      ROLDP(30)=RPTS(24,1)
      ROLDP(31)=RPTS(35,1)
      ROLDP(32)=RPTS(36,1)
      ROLDP(33)=RPTS(37,1)
      ROLDP(34)=RPTS(38,1)
      ROLDP(35)=RPTS(39,1)
      ROLDP(36)=RPTS(40,1)
      ROLDP(37)=RPTS(41,1)
      ROLDP(38)=RPTS(42,1)
      ROLDP(39)=RPTS(44,1)
      ROLDP(40)=RPTS(45,1)
      ROLDP(41)=RPTS(46,1)
      ROLDP(42)=RPTS(47,1)
      ROLDP(43)=RPTS(48,1)
      ROLDP(44)=RPTS(49,1)
      ROLDP(45)=RPTS(50,1)
      ROLDP(46)=RPTS(51,1)
      ROLDP(47)=RPTS(52,1)
      ROLDP(48)=RPTS(53,1)
      ROLDP(49)=RPTS(54,1)
      ROLDP(50)=RPTS(55,1)
      ROLDP(51)=RPTS(56,1)
      ROLDP(52)=RPTS(57,1)
      ROLDP(53)=RPTS(58,1)
      ROLDP(54)=RPTS(59,1)
      ROLDP(55)=RPTS(60,1)
      ROLDP(56)=RPTS(61,1)
      NNEWP=0
      KNEWP(1)=0
      NPTSE=1
C
C
C     Loop along the ray:
  5   CONTINUE
      IF (NNEWP.LE.0) THEN
C       Reading the material parameters in a new point on the ray:
        IF (NNEWP.LT.0) THEN
C         WAN-02
          CALL ERROR('WAN-02: Negative number of points')
C         This error should not appear.
C         The number of new points should be zero or positive integer.
        ENDIF
        IF (NPTSE+1.GT.NPTS) THEN
C         WAN-03
          CALL ERROR('WAN-03: Point along ray missing')
C         Here we assume that we are checking the interval of points
C         NPTSE to NPTSE+1, and thus NPTS should be at least NPTSE+1
C         or more.  This error should not appear.
        ENDIF
        I1=NPTSE+1
        RNEWP( 1,1)=RPTS( 2,I1)
        RNEWP( 2,1)=RPTS( 3,I1)
        RNEWP( 3,1)=RPTS( 4,I1)
        RNEWP( 4,1)=RPTS( 5,I1)
        RNEWP( 5,1)=RPTS( 6,I1)
        RNEWP( 6,1)=RPTS( 7,I1)
        RNEWP( 7,1)=RPTS( 8,I1)
        RNEWP( 8,1)=RPTS( 9,I1)
        RNEWP( 9,1)=RPTS(10,I1)
        RNEWP(10,1)=RPTS(11,I1)
        RNEWP(27,1)=RPTS(34,I1)
        RNEWP(28,1)=RPTS(22,I1)
        RNEWP(29,1)=RPTS(23,I1)
        RNEWP(30,1)=RPTS(24,I1)
        RNEWP(31,1)=RPTS(35,I1)
        RNEWP(32,1)=RPTS(36,I1)
        RNEWP(33,1)=RPTS(37,I1)
        RNEWP(34,1)=RPTS(38,I1)
        RNEWP(35,1)=RPTS(39,I1)
        RNEWP(36,1)=RPTS(40,I1)
        RNEWP(37,1)=RPTS(41,I1)
        RNEWP(38,1)=RPTS(42,I1)
        RNEWP(39,1)=RPTS(44,I1)
        RNEWP(40,1)=RPTS(45,I1)
        RNEWP(41,1)=RPTS(46,I1)
        RNEWP(42,1)=RPTS(47,I1)
        RNEWP(43,1)=RPTS(48,I1)
        RNEWP(44,1)=RPTS(49,I1)
        RNEWP(45,1)=RPTS(50,I1)
        RNEWP(46,1)=RPTS(51,I1)
        RNEWP(47,1)=RPTS(52,I1)
        RNEWP(48,1)=RPTS(53,I1)
        RNEWP(49,1)=RPTS(54,I1)
        RNEWP(50,1)=RPTS(55,I1)
        RNEWP(51,1)=RPTS(56,I1)
        RNEWP(52,1)=RPTS(57,I1)
        RNEWP(53,1)=RPTS(58,I1)
        RNEWP(54,1)=RPTS(59,I1)
        RNEWP(55,1)=RPTS(60,I1)
        RNEWP(56,1)=RPTS(61,I1)
        NNEWP=1
        CALL PARM3(IPTS(12,I1),RNEWP(2,1),AA,RHO,QQ)
C       Computing the Christoffel matrix in the new point:
        CALL WACHRI(RNEWP(5,1),RNEWP(6,1),RNEWP(7,1),
     *  AA(1,1),AA(1,2),AA(1,3),AA(1,4),AA(1,5),AA(1,6),AA(1,7),
     *  AA(1,8),AA(1,9),AA(1,10),AA(1,11),AA(1,12),AA(1,13),AA(1,14),
     *  AA(1,15),AA(1,16),AA(1,17),AA(1,18),AA(1,19),AA(1,20),
     *  AA(1,21),RNEWP(11,1),RNEWP(12,1),RNEWP(13,1),EE,DER)
C       Storing the eigenvectors of the Christoffel matrix
        RNEWP(18,1)=EE(4)
        RNEWP(19,1)=EE(5)
        RNEWP(20,1)=EE(6)
        RNEWP(21,1)=EE(7)
        RNEWP(22,1)=EE(8)
        RNEWP(23,1)=EE(9)
        RNEWP(24,1)=EE(1)
        RNEWP(25,1)=EE(2)
        RNEWP(26,1)=EE(3)
C       Projecting the qS eigenvectors to the plane
C       perpendicular to the ray:
        CALL WAPROJ(RNEWP(5,1),RNEWP(6,1),RNEWP(7,1),RNEWP(8,1),
     *  RNEWP(9,1),RNEWP(10,1),EE(4),EE(5),EE(6),EE(7),EE(8),EE(9),
     *  RNEWP(14,1),RNEWP(15,1),RNEWP(16,1),RNEWP(17,1))
C
C       Quantities for future possible interpolation:
        TTNEW=RPTS(2,I1)
C       Coordinates:
        XNEW(1)=RPTS(3,I1)
        XNEW(2)=RPTS(4,I1)
        XNEW(3)=RPTS(5,I1)
C       dx/dt=p*v**2=p/sum(pi pi):
        AUX=(RPTS(6,I1)**2+RPTS(7,I1)**2+RPTS(8,I1)**2)
        DXNEW(1)=RPTS(6,I1)/AUX
        DXNEW(2)=RPTS(7,I1)/AUX
        DXNEW(3)=RPTS(8,I1)/AUX
C       Slowness:
        PNEW(1)=RPTS(6,I1)
        PNEW(2)=RPTS(7,I1)
        PNEW(3)=RPTS(8,I1)
C       dpi/dt=-1/v*dv/dxi:
        AUX=SQRT(AUX)
        DPNEW(1)=-RPTS(22,I1)/AUX
        DPNEW(2)=-RPTS(23,I1)/AUX
        DPNEW(3)=-RPTS(24,I1)/AUX
C       Polarization:
        ENEW(1)=RPTS( 9,I1)
        ENEW(2)=RPTS(10,I1)
        ENEW(3)=RPTS(11,I1)
C       dei/dt=v*sum(dv/dxk ek)*pi
        AUX=AUX*(RPTS(22,I1)*RPTS(9,I1)+RPTS(23,I1)*RPTS(10,I1)+
     *  RPTS(24,I1)*RPTS(11,I1))
        DENEW(1)=AUX*RPTS(6,I1)
        DENEW(2)=AUX*RPTS(7,I1)
        DENEW(3)=AUX*RPTS(8,I1)
C
        IF (ABS(RNEWP(12,1)-RNEWP(13,1)).LT.0.00001) THEN
C         New point is isotropic, qS eigenvalues are the same.
C         No change in eigenvalues and eigenvectors:
          RNEWP(14,1)=ROLDP(14)
          RNEWP(15,1)=ROLDP(15)
          RNEWP(16,1)=ROLDP(16)
          RNEWP(17,1)=ROLDP(17)
          RNEWP(18,1)=ROLDP(18)
          RNEWP(19,1)=ROLDP(19)
          RNEWP(20,1)=ROLDP(20)
          RNEWP(21,1)=ROLDP(21)
          RNEWP(22,1)=ROLDP(22)
          RNEWP(23,1)=ROLDP(23)
          RNEWP(24,1)=ROLDP(24)
          RNEWP(25,1)=ROLDP(25)
          RNEWP(26,1)=ROLDP(26)
          DELTFI=0.
          GOTO 20
        ENDIF
        IF (ABS(ROLDP(12)-ROLDP(13)).LT.0.00001) THEN
C         Old point is isotropic, new point is anisotropic:
          IF (LPTISO.EQ.1) THEN
C           New point is the first anisotropic point on the ray,
C           angular change is not to be computed:
            DELTFI=0.
            LPTISO=0
            GOTO 20
          ENDIF
        ENDIF
        IF (RPTS(2,NPTSE).EQ.RPTS(2,NPTSE+1)) THEN
C         The ray is crossing an interface.
          DELTFI=0.
          GOTO 20
        ENDIF
      ENDIF
C
C     Computing the angular change in eigenvectors, adding new points
C     on the ray if necessary:
      DO 10, I1=1,NNEWP
        CALL WAVECR(ROLDP(18),ROLDP(21),RNEWP(18,I1),RNEWP(21,I1),
     *              RNEWP(12,I1),RNEWP(13,I1))
        DELTFI=WAANGL(ROLDP(18),ROLDP(21),RNEWP(18,I1),RNEWP(21,I1))
        DEFI=1./SQRT(RNEWP(12,I1))-1./SQRT(RNEWP(13,I1))
        DEFI=DEFI-1./SQRT(ROLDP(12))+1./SQRT(ROLDP(13))
        IF (SINGLF.NE.0.) THEN
          DEFI=DEFI*PI*SINGLF*RPTS(2,NPTS)/6.
        ELSE
          DEFI=DEFI*PI*(OF+FLOAT(NF-1)*DF)*RPTS(2,NPTS)/6.
        ENDIF
        DEFI=ABS(DEFI)
        IF (DEFI*PI/8..GT.ERRWAN) THEN
          DEFI=ERRWAN/DEFI
        ELSE
          DEFI=PI/8.
        END IF
        IF (ABS(DELTFI).LE.DEFI) THEN
C         Angular change is less than prescribed limit, the I1-th point
C         of array RNEWP will be used as the next point on the ray:
          NNEWP=I1
          GOTO 20
        ENDIF
   10 CONTINUE
C     Angular change is greater than prescribed limit for all the points
C     of array RNEWP.
   15 CONTINUE
C     Loop for adding new points on the ray until the angular change is
C     small enough:
        IF (NNEWP.GE.MNEWP-2) THEN
C         WAN-05
          CALL ERROR('WAN-05: Maximum number of divisions exceeded')
C         The angular change of eigenvectors in two consecutive points
C         is too big. More divisions than fits into the memory is
C         needed to keep the change under the prescribed limit. Try to
C         decrease the parameter
C         STORE,
C         or to increase the dimension MNEWP or the angular change DEFI.
        ENDIF
C       Adding new point to the ray:
        KNEWP(NNEWP)=KNEWP(NNEWP)+1
        NNEWP=NNEWP+1
        IF (NNEWP.GT.MNEWP) THEN
C         WAN-06
          CALL ERROR('WAN-06: Array KNEWP too small')
C         This error should not appear, error05 should appear instead.
        ENDIF
C       Travel time:
        RNEWP(1,NNEWP)=(ROLDP(1)+RNEWP(1,NNEWP-1))*0.5
C       Coordinates:
        CALL HIVD2(3,TTOLD,XOLD,DXOLD,TTNEW,XNEW,DXNEW,
     *  RNEWP(1,NNEWP),RNEWP(2,NNEWP),DER)
C       Slowness vector:
        CALL HIVD2(3,TTOLD,POLD,DPOLD,TTNEW,PNEW,DPNEW,
     *  RNEWP(1,NNEWP),RNEWP(5,NNEWP),DER)
C       Polarization vector:
        CALL HIVD2(3,TTOLD,EOLD,DEOLD,TTNEW,ENEW,DENEW,
     *  RNEWP(1,NNEWP),RNEWP(8,NNEWP),DER)
C       Matrices Q and P, basis vectors of RCCS:
        DO 17, I1=31,56
          RNEWP(I1,NNEWP)=(ROLDP(I1)+RNEWP(I1,NNEWP-1))*0.5
  17    CONTINUE
C       Material parameters:
        CALL PARM2(IPTS(12,NPTSE+1),RNEWP(2,NNEWP),UP,US,RHO,QP,QS)
        CALL VELOC(IPTS(12,NPTSE+1),UP,US,QP,QS,VP,VS,VD,QL)
        RNEWP(27,NNEWP)=VD(1)
        RNEWP(28,NNEWP)=VD(2)
        RNEWP(29,NNEWP)=VD(3)
        RNEWP(30,NNEWP)=VD(4)
        CALL PARM3(IPTS(12,NPTSE+1),RNEWP(2,NNEWP),AA,RHO,QQ)
C       Christoffel matrix and eigenvalues:
        CALL WACHRI(RNEWP(5,NNEWP),RNEWP(6,NNEWP),RNEWP(7,NNEWP),
     *  AA(1,1),AA(1,2),AA(1,3),AA(1,4),AA(1,5),AA(1,6),AA(1,7),
     *  AA(1,8),AA(1,9),AA(1,10),AA(1,11),AA(1,12),AA(1,13),AA(1,14),
     *  AA(1,15),AA(1,16),AA(1,17),AA(1,18),AA(1,19),AA(1,20),AA(1,21),
     *  RNEWP(11,NNEWP),RNEWP(12,NNEWP),RNEWP(13,NNEWP),EE,DER)
C       Storing the eigenvectors of the Christoffel matrix
        RNEWP(18,NNEWP)=EE(4)
        RNEWP(19,NNEWP)=EE(5)
        RNEWP(20,NNEWP)=EE(6)
        RNEWP(21,NNEWP)=EE(7)
        RNEWP(22,NNEWP)=EE(8)
        RNEWP(23,NNEWP)=EE(9)
        RNEWP(24,NNEWP)=EE(1)
        RNEWP(25,NNEWP)=EE(2)
        RNEWP(26,NNEWP)=EE(3)
C       Projection of the qS eigenvectors to the plane
C       perpendicular to the ray:
        CALL WAPROJ(RNEWP(5,NNEWP),RNEWP(6,NNEWP),RNEWP(7,NNEWP),
     *  RNEWP(8,NNEWP),RNEWP(9,NNEWP),RNEWP(10,NNEWP),
     *  EE(4),EE(5),EE(6),EE(7),EE(8),EE(9),
     *  RNEWP(14,NNEWP),RNEWP(15,NNEWP),RNEWP(16,NNEWP),RNEWP(17,NNEWP))
C       Index of the division:
        KNEWP(NNEWP)=KNEWP(NNEWP-1)
C
        IF (ABS(RNEWP(12,NNEWP)-RNEWP(13,NNEWP)).LT.0.00001) THEN
C         Isotropic case, qS eigenvalues are the same,
C         no change in eigenvectors:
          RNEWP(14,NNEWP)=ROLDP(14)
          RNEWP(15,NNEWP)=ROLDP(15)
          RNEWP(16,NNEWP)=ROLDP(16)
          RNEWP(17,NNEWP)=ROLDP(17)
          RNEWP(18,NNEWP)=ROLDP(18)
          RNEWP(19,NNEWP)=ROLDP(19)
          RNEWP(20,NNEWP)=ROLDP(20)
          RNEWP(21,NNEWP)=ROLDP(21)
          RNEWP(22,NNEWP)=ROLDP(22)
          RNEWP(23,NNEWP)=ROLDP(23)
          RNEWP(24,NNEWP)=ROLDP(24)
          RNEWP(25,NNEWP)=ROLDP(25)
          RNEWP(26,NNEWP)=ROLDP(26)
          DELTFI=0.
        ELSE
          CALL WAVECR(ROLDP(18),ROLDP(21),
     *                RNEWP(18,NNEWP),RNEWP(21,NNEWP),
     *                RNEWP(12,NNEWP),RNEWP(13,NNEWP))
          DELTFI=WAANGL(ROLDP(18),ROLDP(21),
     *                  RNEWP(18,NNEWP),RNEWP(21,NNEWP))
          DEFI=1./SQRT(RNEWP(12,NNEWP))-1./SQRT(RNEWP(13,NNEWP))
          DEFI=DEFI-1./SQRT(ROLDP(12))+1./SQRT(ROLDP(13))
          IF (SINGLF.NE.0.) THEN
            DEFI=DEFI*PI*SINGLF*RPTS(2,NPTS)/6.
          ELSE
            DEFI=DEFI*PI*(OF+FLOAT(NF-1)*DF)*RPTS(2,NPTS)/6.
          ENDIF
          DEFI=ABS(DEFI)
          IF (DEFI*PI/8..GT.ERRWAN) THEN
            DEFI=ERRWAN/DEFI
          ELSE
            DEFI=PI/8.
          END IF
        ENDIF
        IF (ABS(DELTFI).LE.DEFI) THEN
C         Angular change is less than prescribed limit, this point
C         of array RNEWP will be used as the next point on the ray:
          GOTO 20
        ELSE
C         Angular change is greater than prescribed limit, adding
C         a new point to the ray:
          GOTO 15
        ENDIF
C     End of the loop.
C
   20 CONTINUE
C     Angular change DELTFI for points ROLDP, RNEWP(i,NNEWP) is less
C     than prescribed limit. Recording the computed quantities.
      NPTSE=NPTSE+1
      IF (NNEWP.NE.1) THEN
C       The new point was computed by interpolation.
C       Shifting the array RPTS:
        NPTS=NPTS+1
        IF (NPTS.GT.MPTS) THEN
C         WAN-07
          CALL ERROR('WAN-07: Array RPTS small')
C         The dimension of the array RPTS is given by the dimension MRAM
C         in the include file
C         ram.inc.
        ENDIF
        DO 31, I1=NPTS,NPTSE+1,-1
          DO 30, I2=1,NQPTS
            IF ((I2.EQ.1).OR.(I2.EQ.12).OR.(I2.EQ.43)) THEN
              IPTS(I2,I1)=IPTS(I2,I1-1)
            ELSE
              RPTS(I2,I1)=RPTS(I2,I1-1)
            ENDIF
  30      CONTINUE
  31    CONTINUE
C       Recording interpolated quantities:
        IPTS( 1,NPTSE)=0
        RPTS( 2,NPTSE)=RNEWP( 1,NNEWP)
        RPTS( 3,NPTSE)=RNEWP( 2,NNEWP)
        RPTS( 4,NPTSE)=RNEWP( 3,NNEWP)
        RPTS( 5,NPTSE)=RNEWP( 4,NNEWP)
        RPTS( 6,NPTSE)=RNEWP( 5,NNEWP)
        RPTS( 7,NPTSE)=RNEWP( 6,NNEWP)
        RPTS( 8,NPTSE)=RNEWP( 7,NNEWP)
        RPTS( 9,NPTSE)=RNEWP( 8,NNEWP)
        RPTS(10,NPTSE)=RNEWP( 9,NNEWP)
        RPTS(11,NPTSE)=RNEWP(10,NNEWP)
        IPTS(12,NPTSE)=IPTS(12,NPTSE-1)
        RPTS(34,NPTSE)=RNEWP(27,NNEWP)
        RPTS(22,NPTSE)=RNEWP(28,NNEWP)
        RPTS(23,NPTSE)=RNEWP(29,NNEWP)
        RPTS(24,NPTSE)=RNEWP(30,NNEWP)
        RPTS(35,NPTSE)=RNEWP(31,NNEWP)
        RPTS(36,NPTSE)=RNEWP(32,NNEWP)
        RPTS(37,NPTSE)=RNEWP(33,NNEWP)
        RPTS(38,NPTSE)=RNEWP(34,NNEWP)
        RPTS(39,NPTSE)=RNEWP(35,NNEWP)
        RPTS(40,NPTSE)=RNEWP(36,NNEWP)
        RPTS(41,NPTSE)=RNEWP(37,NNEWP)
        RPTS(42,NPTSE)=RNEWP(38,NNEWP)
        IPTS(43,NPTSE)=0
        RPTS(44,NPTSE)=RNEWP(39,NNEWP)
        RPTS(45,NPTSE)=RNEWP(40,NNEWP)
        RPTS(46,NPTSE)=RNEWP(41,NNEWP)
        RPTS(47,NPTSE)=RNEWP(42,NNEWP)
        RPTS(48,NPTSE)=RNEWP(43,NNEWP)
        RPTS(49,NPTSE)=RNEWP(44,NNEWP)
        RPTS(50,NPTSE)=RNEWP(45,NNEWP)
        RPTS(51,NPTSE)=RNEWP(46,NNEWP)
        RPTS(52,NPTSE)=RNEWP(47,NNEWP)
        RPTS(53,NPTSE)=RNEWP(48,NNEWP)
        RPTS(54,NPTSE)=RNEWP(49,NNEWP)
        RPTS(55,NPTSE)=RNEWP(50,NNEWP)
        RPTS(56,NPTSE)=RNEWP(51,NNEWP)
        RPTS(57,NPTSE)=RNEWP(52,NNEWP)
        RPTS(58,NPTSE)=RNEWP(53,NNEWP)
        RPTS(59,NPTSE)=RNEWP(54,NNEWP)
        RPTS(60,NPTSE)=RNEWP(55,NNEWP)
        RPTS(61,NPTSE)=RNEWP(56,NNEWP)
      ENDIF
C     Recording quantities for computation of anisotropic corrections:
      RPTS(13,NPTSE)=RNEWP(11,NNEWP)
      RPTS(14,NPTSE)=RNEWP(12,NNEWP)
      RPTS(15,NPTSE)=RNEWP(13,NNEWP)
      RPTS(16,NPTSE)=RNEWP(14,NNEWP)
      RPTS(17,NPTSE)=RNEWP(15,NNEWP)
      RPTS(18,NPTSE)=RNEWP(16,NNEWP)
      RPTS(19,NPTSE)=RNEWP(17,NNEWP)
      RPTS(20,NPTSE)=DELTFI
      IF (KQIDT.LE.0) THEN
        RPTS(21,NPTSE)=
     *     (0.5/SQRT(RNEWP(13,NNEWP))-0.5/SQRT(RNEWP(12,NNEWP))
     *     +0.5/SQRT(ROLDP(13))      -0.5/SQRT(ROLDP(12)))
     *    *(RNEWP(1,NNEWP)-ROLDP(1))
      ELSE
        RPTS(21,NPTSE)=0.25*(RNEWP(12,NNEWP)-RNEWP(13,NNEWP)
     *                      +ROLDP(12)      -ROLDP(13))
     *                     *(RNEWP(1,NNEWP)-ROLDP(1))
      END IF
      RPTS(25,NPTSE)=RNEWP(18,NNEWP)
      RPTS(26,NPTSE)=RNEWP(19,NNEWP)
      RPTS(27,NPTSE)=RNEWP(20,NNEWP)
      RPTS(28,NPTSE)=RNEWP(21,NNEWP)
      RPTS(29,NPTSE)=RNEWP(22,NNEWP)
      RPTS(30,NPTSE)=RNEWP(23,NNEWP)
      RPTS(31,NPTSE)=RNEWP(24,NNEWP)
      RPTS(32,NPTSE)=RNEWP(25,NNEWP)
      RPTS(33,NPTSE)=RNEWP(26,NNEWP)
C
      IF (NPTSE.LT.NPTS) THEN
C       Continuing with the next point on the ray:
        DO 100, I1=1,MQANT
          ROLDP(I1)=RNEWP(I1,NNEWP)
  100   CONTINUE
        KNEWP(NNEWP)=0
        NNEWP=NNEWP-1
        IF (NNEWP.LE.0) THEN
          TTOLD=TTNEW
          DO 101, I1=1,3
            XOLD(I1)=XNEW(I1)
            DXOLD(I1)=DXNEW(I1)
            POLD(I1)=PNEW(I1)
            DPOLD(I1)=DPNEW(I1)
            EOLD(I1)=ENEW(I1)
            DEOLD(I1)=DENEW(I1)
  101     CONTINUE
        ENDIF
        GOTO 5
      ENDIF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE WAREAD(LU1,LU2,LU3,RPTS,IPTS,NQPTS,MPTS,NPTS)
C
C----------------------------------------------------------------------
C Subroutine reads the unformatted output of program CRT.
C Reading the output files is completed by a simple invocation of
C subroutine AP00 of file 'ap.for'. Subroutine reads the quantities
C stored in individual points on the ray into array RPTS.
C
      INTEGER LU1,LU2,LU3
      INTEGER NPTS,NQPTS,MPTS
      REAL    RPTS(NQPTS,MPTS)
      INTEGER IPTS(NQPTS,MPTS)
C Input:
C     LU1 ... Number of the logical unit connected to the CRT output
C             file with quantities along rays.
C     LU2 ... Number of the logical unit connected to the CRT output
C             file with quantities at the initial points of rays.
C     LU3 ... Number of the logical unit connected to the CRT output
C             file with quantities at the storing surface.
C     NQPTS.. Dimension of arrays RPTS and IPTS.
C     MPTS... Dimension of arrays RPTS and IPTS.
C Output:
C     RPTS,IPTS ... The arrays are filled with the quantities stored in
C                  CRT output files for single two-point ray during one
C                  invocation of this subroutine:
C        IPTS(1,I) ... Index of the I-th point, zero for points
C                      added to the ray by interpolation.
C        RPTS(2,I) ... Travel time in I-th point.
C        RPTS(3-5,I) . Coordinates of the point.
C        RPTS(6-8,I) . Slowness vector in the point.
C        RPTS(9-11,I)  Polarization vector in the point.
C        IPTS(12,I) .. Index of complex block.
C        RPTS(22-24,I) Derivatives of the velocity.
C        RPTS(34,I)... S-wave velocity.
C        RPTS(35-38,I) Matrix of geometrical spreading.
C        RPTS(39-42,I) Transformation matrix P.
C        IPTS(43,I)... Index of surface, zero inside a complex block.
C        RPTS(44-52,I) Covariant components of the basis vectors
C                      of the ray-centred coordinate system.
C        RPTS(53-61,I) Contravariant components of the basis vectors
C                      of the ray-centred coordinate system.
C     NPTS ... Number of points stored in RPTS (IPTS).
C
C ...........................
C Common block /POINTC/ to store the results of complete ray tracing:
      INCLUDE 'pointc.inc'
C ...........................
C Common block /RTMAT/ to store the matrices of reflection-transmittion
C coefficients.
      INTEGER MRT,NRT
      PARAMETER (MRT=100)
      REAL PIRTR(2,2,MRT),PIRTI(2,2,MRT)
      COMMON/RTMAT/NRT,PIRTR,PIRTI
      SAVE /RTMAT/
C     MRT ... Dimension of arrays.
C     NRT ... Number of stored R-T matrices.
C     PIRTR .. Real part of R-T matrices.
C     PIRTI .. Imaginary part of R-T matrices.
C     The matrices are defined here in terms of polarization vectors.
C ...........................
C Common block /WASEP/ to store the input SEP parameters used
C in subroutine WAN and their subroutines.
      REAL CRTANI,DEGREE,SINGLF
      COMMON/WASEP/CRTANI,DEGREE,SINGLF
      SAVE /WASEP/
C     CRTANI ... Switch between isotropic and anisotropic ray tracing.
C     DEGREE ... Degree of the homogeneous Hamiltonian.
C     SINGLF ... Frequency of the prevailing-frequency approximation.
C.......................................................................
C
C     Auxiliary storage locations:
      INTEGER ISRFA,ISRF2,KODNEW,ICBNEW,IEN,ISB2A,ISB2B,ICB2A,ICB2B,I1
      INTEGER IYFA(13)
      REAL UP(10),US(10),RHO,QP,QS,VP,VS,VD(10),QL
      REAL YLFA(6),YFA(35),YYFA(5)
      REAL HI(18),H(18),HUI(9)
C
C-----------------------------------------------------------------------
C
      ISRFA=0
      NRT=0
  5   CONTINUE
      CALL AP00(LU1,LU2,LU3)
      IF (IWAVE.LT.1) THEN
C       End of rays:
        NPTS=0
        RETURN
      ENDIF
      IF (IREC.LE.0)
C       Only two-point rays are to be written into RPTS (IPTS).
     *  GOTO 5
      IF ((IPT.NE.0).AND.(IPT.NE.1))
C       End of previous ray behind the reference surface:
     *  GOTO 5
C
C     Initial point of a new two-point ray:
      NPTS=1
C     Index of a point:
      IPTS(1,NPTS)=NPTS
C     Travel time:
      RPTS(2,NPTS)=YI(1)
C     Coordinates:
      RPTS(3,NPTS)=YI(3)
      RPTS(4,NPTS)=YI(4)
      RPTS(5,NPTS)=YI(5)
C     Slowness vector:
      RPTS(6,NPTS)=YI(6)
      RPTS(7,NPTS)=YI(7)
      RPTS(8,NPTS)=YI(8)
C     Polarization vector:
      RPTS(9,NPTS)= YI(9)
      RPTS(10,NPTS)=YI(10)
      RPTS(11,NPTS)=YI(11)
C     Index of complex block:
      IPTS(12,NPTS)=ICB1I
C     S-wave velocity and their derivatives:
      IF (CRTANI.EQ.0) THEN
        RPTS(34,NPTS)=YLI(2)
        RPTS(22,NPTS)=YLI(4)
        RPTS(23,NPTS)=YLI(5)
        RPTS(24,NPTS)=YLI(6)
      ELSE
        CALL PARM2(IPTS(12,NPTS),RPTS(3,NPTS),UP,US,RHO,QP,QS)
        CALL VELOC(IPTS(12,NPTS),UP,US,QP,QS,VP,VS,VD,QL)
        RPTS(34,NPTS)=VD(1)
        RPTS(22,NPTS)=VD(2)
        RPTS(23,NPTS)=VD(3)
        RPTS(24,NPTS)=VD(4)
      ENDIF
C     Matrix of geometrical spreading:
      RPTS(35,NPTS)=YI(12)
      RPTS(36,NPTS)=YI(13)
      RPTS(37,NPTS)=YI(16)
      RPTS(38,NPTS)=YI(17)
C     Transformation matrix P:
      RPTS(39,NPTS)=YI(14)
      RPTS(40,NPTS)=YI(15)
      RPTS(41,NPTS)=YI(18)
      RPTS(42,NPTS)=YI(19)
C     Index of the surface:
      IPTS(43,NPTS)=0
C     Covariant and contravariant components of the basis vectors
C     of the RCCS:
      CALL AP03(0,RPTS(44,NPTS),H,HUI)
      IF (CRTANI.NE.0) THEN
        CALL AN03(ICB1I,YI,RPTS(44,NPTS))
      ENDIF
C
C
  20  CONTINUE
C     New point:
      IF (NYF.EQ.0) THEN
C       WAN-04
        CALL ERROR('WAN-04: Undefined quantities at point O/F.')
C       Quantities at the point O/F are not defined. This may happen
C       if the CRT output quantities are stored at the specified surface
C       only, and not along rays.
      ENDIF
      IF (YF(1).LT.Y(1)) THEN
C       The point along the ray is before the reference surface,
C       recording the point along the ray:
        NPTS=NPTS+1
        IF (NPTS.GT.MPTS) THEN
C         WAN-08
          CALL ERROR('WAN-08: Array RPTS small')
C         The dimension of the array RPTS is given by the dimension MRAM
C         in the include file
C         ram.inc.
        ENDIF
C       Index of a point:
        IPTS(1,NPTS)=NPTS
C       Travel time:
        RPTS(2,NPTS)=YF(1)
C       Coordinates:
        RPTS(3,NPTS)=YF(3)
        RPTS(4,NPTS)=YF(4)
        RPTS(5,NPTS)=YF(5)
C       Slowness vector:
        RPTS(6,NPTS)=YF(6)
        RPTS(7,NPTS)=YF(7)
        RPTS(8,NPTS)=YF(8)
C       Polarization vector:
        RPTS(9,NPTS)= YF(9)
        RPTS(10,NPTS)=YF(10)
        RPTS(11,NPTS)=YF(11)
C       Index of complex block:
        IPTS(12,NPTS)=ICB1F
C       S-wave velocity and their derivatives:
        IF (CRTANI.EQ.0) THEN
          RPTS(34,NPTS)=YLF(2)
          RPTS(22,NPTS)=YLF(4)
          RPTS(23,NPTS)=YLF(5)
          RPTS(24,NPTS)=YLF(6)
        ELSE
          CALL PARM2(IPTS(12,NPTS),RPTS(3,NPTS),UP,US,RHO,QP,QS)
          CALL VELOC(IPTS(12,NPTS),UP,US,QP,QS,VP,VS,VD,QL)
          RPTS(34,NPTS)=VD(1)
          RPTS(22,NPTS)=VD(2)
          RPTS(23,NPTS)=VD(3)
          RPTS(24,NPTS)=VD(4)
        ENDIF
C       Matrix of geometrical spreading:
        RPTS(35,NPTS)=YF(12)*YI(12)+YF(16)*YI(13)+
     *                YF(20)*YI(14)+YF(24)*YI(15)
        RPTS(36,NPTS)=YF(13)*YI(12)+YF(17)*YI(13)+
     *                YF(21)*YI(14)+YF(25)*YI(15)
        RPTS(37,NPTS)=YF(12)*YI(16)+YF(16)*YI(17)+
     *                YF(20)*YI(18)+YF(24)*YI(19)
        RPTS(38,NPTS)=YF(13)*YI(16)+YF(17)*YI(17)+
     *                YF(21)*YI(18)+YF(25)*YI(19)
C       Transformation matrix P:
        RPTS(39,NPTS)=YF(14)*YI(12)+YF(18)*YI(13)+
     *                YF(22)*YI(14)+YF(26)*YI(15)
        RPTS(40,NPTS)=YF(15)*YI(12)+YF(19)*YI(13)+
     *                YF(23)*YI(14)+YF(27)*YI(15)
        RPTS(41,NPTS)=YF(14)*YI(16)+YF(18)*YI(17)+
     *                YF(22)*YI(18)+YF(26)*YI(19)
        RPTS(42,NPTS)=YF(15)*YI(16)+YF(19)*YI(17)+
     *                YF(23)*YI(18)+YF(27)*YI(19)
C       Index of the surface:
        IPTS(43,NPTS)=ISRFF
C       Covariant and contravariant components of the basis vectors
C       of the RCCS:
        CALL APYYF
        CALL AP03(0,HI,RPTS(44,NPTS),HUI)
        CALL APYYF
        IF (CRTANI.NE.0) THEN
          CALL AN03(ICB1F,YF,RPTS(44,NPTS))
        ENDIF
C
        IF ((ISRFF.NE.0).AND.(ISRFA.EQ.0)) THEN
C         Crossing an interface - first point at an interface:
          CALL BLOCK(YF(3),0,0,ISRF2,ISB2A,ICB2A)
          CALL BLOCK(YF(3),ISRFF,ISB2A,ISRF2,ISB2B,ICB2B)
          IF (ICB2A.NE.IABS(ICB1F)) THEN
            I1=ICB2B
            ICB2B=ICB2A
            ICB2A=I1
            I1=ISB2B
            ISB2B=ISB2A
            ISB2A=I1
          ENDIF
          IF (ICB2A.NE.IABS(ICB1F)) THEN
C           WAN-15
            CALL ERROR('WAN-15: Wrong indices of complex blocks.')
C           This error should not appear, at least one of the indices
C           ICB2A, ICB2B should equal IABS(ICB1F).
          ENDIF
          DO 28, I1=1,5
            YYFA(I1)=0.
            YLFA(I1)=YLF(I1)
  28      CONTINUE
          YLFA(6)=YLF(6)
          DO 30, I1=1,27
            YFA(I1)=YF(I1)
  30      CONTINUE
          YFA(28)=1.
          YFA(29)=0.
          YFA(30)=0.
          YFA(31)=0.
          YFA(32)=0.
          YFA(33)=0.
          YFA(34)=1.
          YFA(35)=0.
          IF (ICB1F.GE.0) THEN
            IYFA(1)=31
          ELSE
            IYFA(1)=35
          ENDIF
          IYFA(2)=0
          IYFA(3)=0
          IYFA(4)=ISB2A
          IYFA(5)=ICB1F
          IYFA(6)=ISRFF
          IYFA(7)=ISB2B
          IYFA(8)=ICB2B
          IYFA(9)=0
          IYFA(10)=0
          IYFA(11)=0
          IYFA(12)=0
          IYFA(13)=0
          KODNEW=0
          ISRFA=ISRFF
        ELSEIF ((ISRFF.NE.0).AND.(ISRFA.NE.0)) THEN
C         Crossing an interface - second point at an interface:
          ICBNEW=ICB1F
          CALL TRANS(YLFA,YFA,YYFA,IYFA,KODNEW,ICBNEW,IEN)
          NRT=NRT+1
          IF (NRT.GT.MRT) THEN
C           WAN-16
            CALL ERROR('WAN-16: Small dimension of PIRTR and PIRTI.')
C           The arrays store R-T coefficients along the ray, the
C           dimension MRT should exceed the number of interfaces
C           passed by the ray.
          ENDIF
          PIRTR(1,1,NRT)=YFA(28)
          PIRTI(1,1,NRT)=YFA(29)
          IF (IPTS(12,NPTS-1).GE.0) THEN
C           Incident P wave:
            PIRTR(2,1,NRT)=YFA(30)
            PIRTI(2,1,NRT)=YFA(31)
            PIRTR(1,2,NRT)=0.
            PIRTI(1,2,NRT)=0.
            PIRTR(2,2,NRT)=0.
            PIRTI(2,2,NRT)=0.
          ELSE
C           Incident S wave:
            PIRTR(1,2,NRT)=YFA(30)
            PIRTI(1,2,NRT)=YFA(31)
            PIRTR(2,1,NRT)=YFA(32)
            PIRTI(2,1,NRT)=YFA(33)
            PIRTR(2,2,NRT)=YFA(34)
            PIRTI(2,2,NRT)=YFA(35)
          ENDIF
          RPTS(2,NPTS)=RPTS(2,NPTS-1)
          ISRFA=0
        ENDIF
C
C       Reading the results of the complete ray tracing:
        CALL AP00(LU1,LU2,LU3)
        IF ((IWAVE.LT.1).OR.(IPT.LE.1)) THEN
C         This should not happen, the ray must reach the
C         reference surface.
C         WAN-09
          CALL ERROR('WAN-09: The ray missed the reference surface')
C         As only the two-point rays are considered by the subroutine
C         "WAN", each of the rays should pass the reference surface.
C         Check, whether you have specified right names
C         of the input files with points along rays,
C         points at their initial points and points at the reference
C         surface in file
C         CRTOUT,
C         or whether you have correctly specified its name
C         CRTOUT.
        ENDIF
C
        GOTO 20
      ENDIF
C
C     The point along the ray is at or above the reference surface,
C     recording the point at the reference surface:
      NPTS=NPTS+1
      IF (NPTS.GT.MPTS) THEN
C       WAN-10
        CALL ERROR('WAN-10: Array RPTS small')
C       The dimension of the array RPTS is given by the dimension MRAM
C       of in the include file
C       ram.inc.
      ENDIF
C     Index of a point:
      IPTS(1,NPTS)=NPTS
C     Travel time:
      RPTS(2,NPTS)=Y(1)
C     Coordinates:
      RPTS(3,NPTS)=Y(3)
      RPTS(4,NPTS)=Y(4)
      RPTS(5,NPTS)=Y(5)
C     Slowness vector:
      RPTS(6,NPTS)=Y(6)
      RPTS(7,NPTS)=Y(7)
      RPTS(8,NPTS)=Y(8)
C     Polarization vector:
      RPTS(9,NPTS)= Y(9)
      RPTS(10,NPTS)=Y(10)
      RPTS(11,NPTS)=Y(11)
C     Index of complex block:
      IPTS(12,NPTS)=ICB1
C     S-wave velocity and their derivatives:
      IF (CRTANI.EQ.0) THEN
        RPTS(34,NPTS)=YL(2)
        RPTS(22,NPTS)=YL(4)
        RPTS(23,NPTS)=YL(5)
        RPTS(24,NPTS)=YL(6)
      ELSE
        CALL PARM2(IPTS(12,NPTS),RPTS(3,NPTS),UP,US,RHO,QP,QS)
        CALL VELOC(IPTS(12,NPTS),UP,US,QP,QS,VP,VS,VD,QL)
        RPTS(34,NPTS)=VD(1)
        RPTS(22,NPTS)=VD(2)
        RPTS(23,NPTS)=VD(3)
        RPTS(24,NPTS)=VD(4)
      ENDIF
C     Matrix of geometrical spreading:
      CALL AP05(0,HUI,
     *          RPTS(35,NPTS),RPTS(36,NPTS),RPTS(37,NPTS),RPTS(38,NPTS))
C     Transformation matrix P:
      CALL AP06(0,HUI,
     *          RPTS(39,NPTS),RPTS(40,NPTS),RPTS(41,NPTS),RPTS(42,NPTS))
C     Index of the surface:
      IPTS(43,NPTS)=ISRF
C     Covariant and contravariant components of the basis vectors
C     of the RCCS:
      CALL AP03(0,HI,RPTS(44,NPTS),HUI)
      IF (CRTANI.NE.0) THEN
        CALL AN03(ICB1,Y,RPTS(44,NPTS))
      ENDIF
C
      RETURN
      END
C
C=======================================================================
C
       SUBROUTINE WAPERT(RPTS,IPTS,NQPTS,STEP,DTAU)
C
C----------------------------------------------------------------------
C Subroutine to compute second-order travel-time corrections.
      INTEGER NQPTS
      REAL    RPTS(NQPTS),STEP,DTAU(7)
      INTEGER IPTS(NQPTS)
C
C Input:
C     IPTS,RPTS,NQPTS... Quantities in the point on the ray.
C             See the description above.
C     STEP... Step of the independent variable along the ray, see
C             the CRT input parameter
C             STORE.
C Output:
C     DTAU... Second-order travel-time corrections.
C
C.......................................................................
C
C Subroutines and external functions required:
      EXTERNAL HDER,WASUM,WASUM3,WASUM4,WASUM5
      REAL WASUM,WASUM3,WASUM4,WASUM5
C     HDER... File hder.for.
C     WASUM,WASUM3,WASUM4,WASUM5 ... This file.
C ...........................
C Common block /WASEP/ to store the input SEP parameters used
C in subroutine WAN and their subroutines.
      REAL CRTANI,DEGREE,SINGLF
      COMMON/WASEP/CRTANI,DEGREE,SINGLF
      SAVE /WASEP/
C     CRTANI ... Switch between isotropic and anisotropic ray tracing.
C     DEGREE ... Degree of the homogeneous Hamiltonian.
C     SINGLF ... Frequency of the prevailing-frequency approximation.
C
C.......................................................................
C
      INTEGER IERR
      REAL AA(10,21),RHO,QQ(21),
     *     AAA(3,3,3,3),AAA1(3,3,3,3),AAA2(3,3,3,3),AAA3(3,3,3,3)
      REAL V0,VV0(3),P(3),E1(3),E2(3),H1(3),H2(3),H3(3),
     *     HC1(3),HC2(3),HC3(3),
     *     HHC0(3),HHC0SP(3,3),HHC0S(3,3),HHTC0S(2,2),
     *     QQT(4),PPT(4),QIT(4),GIJT(2,2),
     *     HH0(3),HHT0(3),HH1,HH2,HH3,
     *     HH1L(3),HH2L(3),HH3L(3),HH1U(3),HH2U(3),HH3U(3),
     *     HHT1L(2),HHT2L(2),HHT3L(2),HHT1U(2),HHT2U(2),HHT3U(2),
     *     QTKT(7),TT(7),TAUT(4),TOUT(2,3),VK(6),T(6)
      REAL HDEP,HDES,HDE,HDE1,HDE2,HDE3,HDE4,HDE5,HDE6,
     *     HDE11,HDE12,HDE22,HDE13,HDE23,HDE33,HDE14,HDE24,HDE34,
     *     HDE44,HDE15,HDE25,HDE35,HDE45,HDE55,HDE16,HDE26,HDE36,
     *     HDE46,HDE56,HDE66
      REAL EE(9)
      INTEGER ISRF1T,IFUN1T(7),IFUN2T(7),ISRF1K,IFUN1K(6),IFUN2K(6),
     *     NSUM1,NSUM2,IX,NDER,NFUN11,NFUN12,NFUN21,NFUN22
      REAL AUX,PIPI,X1T,FUN1T(7),X1K,FUN1K(6)
      SAVE NSUM1,IX,NDER,X1T,ISRF1T,NFUN11,NFUN21,IFUN1T,FUN1T,IFUN2T,
     *     NSUM2,X1K,ISRF1K,NFUN12,NFUN22,IFUN1K,FUN1K,IFUN2K
C     AA ...  Values, first and second partial derivatives of real
C             parts of 21 reduced (divided by the density) elastic
C             parameters.  The order of the value, first and second
C             partial derivatives of each parameter Aij is:
C             Aij,Aij1,Aij2,Aij3,Aij11,Aij12,Aij22,Aij13,Aij23,Aij33.
C             The order of parameters (second array index) is:
C             A11,A12,A22,A13,A23,A33,A14,A24,A34,A44,A15,A25,A35,A45,
C             A55,A16,A26,A36,A46,A56,A66.
C     AAA...  3*3*3*3 matrix of reduced elastic parameters.
C     AAA1... 3*3*3*3 matrix of first partial derivatives of reduced
C             elastic parameters according to X1.
C     AAA2... 3*3*3*3 matrix of first partial derivatives of reduced
C             elastic parameters according to X2.
C     AAA3... 3*3*3*3 matrix of first partial derivatives of reduced
C             elastic parameters according to X3.
C     EE...   Eigenvectors of the Christoffel matrix.  Not used here.
C
C-----------------------------------------------------------------------
C
      IF (CRTANI.EQ.0.) THEN
        V0=RPTS(34)
        VV0(1)=RPTS(22)
        VV0(2)=RPTS(23)
        VV0(3)=RPTS(24)
        P(1)=RPTS(6)
        P(2)=RPTS(7)
        P(3)=RPTS(8)
        E1(1)=RPTS(25)
        E1(2)=RPTS(26)
        E1(3)=RPTS(27)
        E2(1)=RPTS(28)
        E2(2)=RPTS(29)
        E2(3)=RPTS(30)
        QQT(1)=RPTS(35)
        QQT(2)=RPTS(36)
        QQT(3)=RPTS(37)
        QQT(4)=RPTS(38)
        PPT(1)=RPTS(39)
        PPT(2)=RPTS(40)
        PPT(3)=RPTS(41)
        PPT(4)=RPTS(42)
        CALL PARM3(IPTS(12),RPTS(3),AA,RHO,QQ)
        CALL WAIJKL(AA,AAA,1)
        CALL WAIJKL(AA,AAA1,2)
        CALL WAIJKL(AA,AAA2,3)
        CALL WAIJKL(AA,AAA3,4)
        H1(1)=RPTS(44)
        H1(2)=RPTS(45)
        H1(3)=RPTS(46)
        H2(1)=RPTS(47)
        H2(2)=RPTS(48)
        H2(3)=RPTS(49)
        H3(1)=RPTS(50)
        H3(2)=RPTS(51)
        H3(3)=RPTS(52)
        HH0(1)=VV0(1)/V0
        HH0(2)=VV0(2)/V0
        HH0(3)=VV0(3)/V0
        HHT0(1)=WASUM(HH0,H1)
        HHT0(2)=WASUM(HH0,H2)
        HHT0(3)=WASUM(HH0,H3)
        HH1=-1./SQRT(WASUM5(AAA,E1,P,E1,P))
        HH2=-1./SQRT(WASUM5(AAA,E2,P,E2,P))
        HH1L(1)=-0.5*(HH1**3)*WASUM5(AAA1,E1,P,E1,P)
        HH1L(2)=-0.5*(HH1**3)*WASUM5(AAA2,E1,P,E1,P)
        HH1L(3)=-0.5*(HH1**3)*WASUM5(AAA3,E1,P,E1,P)
        HH2L(1)=-0.5*(HH2**3)*WASUM5(AAA1,E2,P,E2,P)
        HH2L(2)=-0.5*(HH2**3)*WASUM5(AAA2,E2,P,E2,P)
        HH2L(3)=-0.5*(HH2**3)*WASUM5(AAA3,E2,P,E2,P)
        HH1U(1)=-(HH1**3)*WASUM4(AAA,E1,E1,P,1)
        HH1U(2)=-(HH1**3)*WASUM4(AAA,E1,E1,P,2)
        HH1U(3)=-(HH1**3)*WASUM4(AAA,E1,E1,P,3)
        HH2U(1)=-(HH2**3)*WASUM4(AAA,E2,E2,P,1)
        HH2U(2)=-(HH2**3)*WASUM4(AAA,E2,E2,P,2)
        HH2U(3)=-(HH2**3)*WASUM4(AAA,E2,E2,P,3)
        HHT1L(1)=WASUM(HH1L,H1)
        HHT1L(2)=WASUM(HH1L,H2)
        HHT2L(1)=WASUM(HH2L,H1)
        HHT2L(2)=WASUM(HH2L,H2)
        HHT1U(1)=WASUM(HH1U,H1)
        HHT1U(2)=WASUM(HH1U,H2)
        HHT2U(1)=WASUM(HH2U,H1)
        HHT2U(2)=WASUM(HH2U,H2)
        QTKT(1)=-((HHT1L(1)+HH1*HHT0(1))*QQT(1)+
     *            (HHT1L(2)+HH1*HHT0(2))*QQT(2)+
     *            HHT1U(1)*PPT(1)+HHT1U(2)*PPT(2))
        QTKT(2)=-((HHT2L(1)+HH2*HHT0(1))*QQT(1)+
     *            (HHT2L(2)+HH2*HHT0(2))*QQT(2)+
     *            HHT2U(1)*PPT(1)+HHT2U(2)*PPT(2))
        QTKT(3)=-((HHT1L(1)+HH1*HHT0(1))*QQT(3)+
     *            (HHT1L(2)+HH1*HHT0(2))*QQT(4)+
     *            HHT1U(1)*PPT(3)+HHT1U(2)*PPT(4))
        QTKT(4)=-((HHT2L(1)+HH2*HHT0(1))*QQT(3)+
     *            (HHT2L(2)+HH2*HHT0(2))*QQT(4)+
     *            HHT2U(1)*PPT(3)+HHT2U(2)*PPT(4))
        IF (IPTS(1).EQ.1) THEN
C         Initial point of the ray:
C         Initiating first integration:
          NSUM1=4
          IX=1
          NDER=1
          NFUN11=4
          NFUN21=4
          IFUN1T(1)=1
          IFUN1T(2)=3
          IFUN1T(3)=2
          IFUN1T(4)=4
          FUN1T(1)=QTKT(1)
          FUN1T(2)=QTKT(2)
          FUN1T(3)=QTKT(3)
          FUN1T(4)=QTKT(4)
          IFUN2T(1)=1
          IFUN2T(2)=3
          IFUN2T(3)=2
          IFUN2T(4)=4
C         Initiating second integration:
          NSUM2=3
          NFUN12=3
          NFUN22=3
          IFUN1K(1)=1
          IFUN1K(2)=2
          IFUN1K(3)=3
          AUX=V0**2
          TAUT(1)=-HHT1U(1)/AUX
          TAUT(2)=-HHT1U(2)/AUX
          TAUT(3)=-HHT2U(1)/AUX
          TAUT(4)=-HHT2U(2)/AUX
          FUN1K(1)=-(2.*(HHT1U(1)*TAUT(1)+HHT1U(2)*TAUT(2)) +
     *               (TAUT(1)*TAUT(1)+TAUT(2)*TAUT(2))*AUX)
          FUN1K(2)=-(2.*(HHT2U(1)*TAUT(3)+HHT2U(2)*TAUT(4)) +
     *               (TAUT(3)*TAUT(3)+TAUT(4)*TAUT(4))*AUX)
          FUN1K(3)=-(HHT1U(1)*TAUT(3)+HHT1U(2)*TAUT(4) +
     *               HHT2U(1)*TAUT(1)+HHT2U(2)*TAUT(2) +
     *               (TAUT(1)*TAUT(3)+TAUT(2)*TAUT(4))*AUX)
          IFUN2K(1)=1
          IFUN2K(2)=2
          IFUN2K(3)=3
C         Output values:
          DTAU(1)=0.
          DTAU(2)=0.
          DTAU(3)=0.
          DTAU(4)=0.
          DTAU(5)=0.
          DTAU(6)=0.
          DTAU(7)=0.
          RETURN
        ENDIF
        CALL AP28(NSUM1,TT,IX,NDER,STEP,X1T,ISRF1T,NFUN11,IFUN1T,FUN1T,
     *            NFUN21,IFUN2T,QTKT)
        AUX=QQT(1)*QQT(4)-QQT(2)*QQT(3)
        IF (AUX.EQ.0) THEN
C         WAN-18
          CALL ERROR('WAN-18: Zero determinant of geom. spread.')
C         This error should not appear.
        ENDIF
        QIT(1)= QQT(4)/AUX
        QIT(4)= QQT(1)/AUX
        QIT(2)=-QQT(2)/AUX
        QIT(3)=-QQT(3)/AUX
        TAUT(1)=QIT(1)*TT(1)+QIT(2)*TT(2)
        TAUT(2)=QIT(3)*TT(1)+QIT(4)*TT(2)
        TAUT(3)=QIT(1)*TT(3)+QIT(2)*TT(4)
        TAUT(4)=QIT(3)*TT(3)+QIT(4)*TT(4)
        VK(1)=-(2.*(HHT1U(1)*TAUT(1)+HHT1U(2)*TAUT(2)) +
     *          (TAUT(1)*TAUT(1)+TAUT(2)*TAUT(2))*V0**2)
        VK(2)=-(2.*(HHT2U(1)*TAUT(3)+HHT2U(2)*TAUT(4)) +
     *          (TAUT(3)*TAUT(3)+TAUT(4)*TAUT(4))*V0**2)
        VK(3)=-(HHT1U(1)*TAUT(3)+HHT1U(2)*TAUT(4) +
     *          HHT2U(1)*TAUT(1)+HHT2U(2)*TAUT(2) +
     *          (TAUT(1)*TAUT(3)+TAUT(2)*TAUT(4))*V0**2)
        CALL AP28(NSUM2,T,IX,NDER,STEP,X1K,ISRF1K,NFUN12,IFUN1K,FUN1K,
     *            NFUN22,IFUN2K,VK)
C       Second order perturbations towards anisotropic rays:
        DTAU(1)=T(1)/2.
        DTAU(3)=T(2)/2.
C       Estimation of the two equal second order perturbations from
C       anisotropic common ray towards anisotropic rays:
        DTAU(2)=T(3)
        DTAU(4)=0.
        DTAU(5)=0.
        DTAU(6)=0.
        DTAU(7)=0.
        RETURN
      ENDIF
C
      IF (CRTANI.NE.0.) THEN
        V0=RPTS(34)
        VV0(1)=RPTS(22)
        VV0(2)=RPTS(23)
        VV0(3)=RPTS(24)
        P(1)=RPTS(6)
        P(2)=RPTS(7)
        P(3)=RPTS(8)
        E1(1)=RPTS(25)
        E1(2)=RPTS(26)
        E1(3)=RPTS(27)
        E2(1)=RPTS(28)
        E2(2)=RPTS(29)
        E2(3)=RPTS(30)
        QQT(1)=RPTS(35)
        QQT(2)=RPTS(36)
        QQT(3)=RPTS(37)
        QQT(4)=RPTS(38)
        PPT(1)=RPTS(39)
        PPT(2)=RPTS(40)
        PPT(3)=RPTS(41)
        PPT(4)=RPTS(42)
        CALL PARM3(IPTS(12),RPTS(3),AA,RHO,QQ)
        CALL WAIJKL(AA,AAA,1)
        CALL WAIJKL(AA,AAA1,2)
        CALL WAIJKL(AA,AAA2,3)
        CALL WAIJKL(AA,AAA3,4)
        H1(1)=RPTS(44)
        H1(2)=RPTS(45)
        H1(3)=RPTS(46)
        H2(1)=RPTS(47)
        H2(2)=RPTS(48)
        H2(3)=RPTS(49)
        H3(1)=RPTS(50)
        H3(2)=RPTS(51)
        H3(3)=RPTS(52)
        HC1(1)=RPTS(53)
        HC1(2)=RPTS(54)
        HC1(3)=RPTS(55)
        HC2(1)=RPTS(56)
        HC2(2)=RPTS(57)
        HC2(3)=RPTS(58)
        HC3(1)=RPTS(59)
        HC3(2)=RPTS(60)
        HC3(3)=RPTS(61)
        CALL HDER(-1,AA,P(1),P(2),P(3),HDEP,HDES,HDE,
     *            HDE1,HDE2,HDE3,HDE4,HDE5,HDE6,
     *            HDE11,HDE12,HDE22,HDE13,HDE23,HDE33,HDE14,HDE24,HDE34,
     *            HDE44,HDE15,HDE25,HDE35,HDE45,HDE55,HDE16,HDE26,HDE36,
     *            HDE46,HDE56,HDE66,EE,IERR)
        HH0(1)=HDE1
        HH0(2)=HDE2
        HH0(3)=HDE3
        HHC0(1)=HDE4
        HHC0(2)=HDE5
        HHC0(3)=HDE6
        HHC0SP(1,1)=HDE44
        HHC0SP(2,1)=HDE45
        HHC0SP(3,1)=HDE46
        HHC0SP(1,2)=HDE45
        HHC0SP(2,2)=HDE55
        HHC0SP(3,2)=HDE56
        HHC0SP(1,3)=HDE46
        HHC0SP(2,3)=HDE56
        HHC0SP(3,3)=HDE66
C       Eq. 10:
        HHC0S(1,1)=HHC0SP(1,1)-3.*HHC0(1)*HHC0(1)
        HHC0S(2,1)=HHC0SP(2,1)-3.*HHC0(2)*HHC0(1)
        HHC0S(3,1)=HHC0SP(3,1)-3.*HHC0(3)*HHC0(1)
        HHC0S(1,2)=HHC0SP(1,2)-3.*HHC0(1)*HHC0(2)
        HHC0S(2,2)=HHC0SP(2,2)-3.*HHC0(2)*HHC0(2)
        HHC0S(3,2)=HHC0SP(3,2)-3.*HHC0(3)*HHC0(2)
        HHC0S(1,3)=HHC0SP(1,3)-3.*HHC0(1)*HHC0(3)
        HHC0S(2,3)=HHC0SP(2,3)-3.*HHC0(2)*HHC0(3)
        HHC0S(3,3)=HHC0SP(3,3)-3.*HHC0(3)*HHC0(3)
C       Eq. 39:
        HHT0(1)=WASUM(HH0,HC1)
        HHT0(2)=WASUM(HH0,HC2)
        HHT0(3)=WASUM(HH0,HC3)
C       Eq. 41:
        HHTC0S(1,1)=WASUM3(HHC0S,H1,H1)
        HHTC0S(1,2)=WASUM3(HHC0S,H1,H2)
        HHTC0S(2,1)=WASUM3(HHC0S,H2,H1)
        HHTC0S(2,2)=WASUM3(HHC0S,H2,H2)
C       Eq. 11:
        HH1=-1./SQRT(WASUM5(AAA,E1,P,E1,P))
        HH2=-1./SQRT(WASUM5(AAA,E2,P,E2,P))
C       Eq. 15:
        PIPI=(P(1)**2+P(2)**2+P(3)**2)
        HH3=-1./V0/SQRT(PIPI)
C       Eq. 12:
        HH1L(1)=-0.5*(HH1**3)*WASUM5(AAA1,E1,P,E1,P)
        HH1L(2)=-0.5*(HH1**3)*WASUM5(AAA2,E1,P,E1,P)
        HH1L(3)=-0.5*(HH1**3)*WASUM5(AAA3,E1,P,E1,P)
        HH2L(1)=-0.5*(HH2**3)*WASUM5(AAA1,E2,P,E2,P)
        HH2L(2)=-0.5*(HH2**3)*WASUM5(AAA2,E2,P,E2,P)
        HH2L(3)=-0.5*(HH2**3)*WASUM5(AAA3,E2,P,E2,P)
C       Eq. 16:
        HH3L(1)=VV0(1)/V0**2/SQRT(PIPI)
        HH3L(2)=VV0(2)/V0**2/SQRT(PIPI)
        HH3L(3)=VV0(3)/V0**2/SQRT(PIPI)
C       Eq. 13:
        HH1U(1)=-(HH1**3)*WASUM4(AAA,E1,E1,P,1)
        HH1U(2)=-(HH1**3)*WASUM4(AAA,E1,E1,P,2)
        HH1U(3)=-(HH1**3)*WASUM4(AAA,E1,E1,P,3)
        HH2U(1)=-(HH2**3)*WASUM4(AAA,E2,E2,P,1)
        HH2U(2)=-(HH2**3)*WASUM4(AAA,E2,E2,P,2)
        HH2U(3)=-(HH2**3)*WASUM4(AAA,E2,E2,P,3)
C       Eq. 17:
        HH3U(1)=1./V0/SQRT(PIPI**3)*P(1)
        HH3U(2)=1./V0/SQRT(PIPI**3)*P(2)
        HH3U(3)=1./V0/SQRT(PIPI**3)*P(3)
C       Eq. 42:
        HHT1L(1)=WASUM(HH1L,HC1)
        HHT1L(2)=WASUM(HH1L,HC2)
        HHT2L(1)=WASUM(HH2L,HC1)
        HHT2L(2)=WASUM(HH2L,HC2)
        HHT3L(1)=WASUM(HH3L,HC1)
        HHT3L(2)=WASUM(HH3L,HC2)
C       Eq. 43:
        HHT1U(1)=WASUM(HH1U,H1)
        HHT1U(2)=WASUM(HH1U,H2)
        HHT2U(1)=WASUM(HH2U,H1)
        HHT2U(2)=WASUM(HH2U,H2)
        HHT3U(1)=WASUM(HH3U,H1)
        HHT3U(2)=WASUM(HH3U,H2)
C       Eq. 55:
        QTKT(1)=-((HHT1L(1)+HH1*HHT0(1))*QQT(1)+
     *            (HHT1L(2)+HH1*HHT0(2))*QQT(2)+
     *            HHT1U(1)*PPT(1)+HHT1U(2)*PPT(2))
        QTKT(2)=-((HHT2L(1)+HH2*HHT0(1))*QQT(1)+
     *            (HHT2L(2)+HH2*HHT0(2))*QQT(2)+
     *            HHT2U(1)*PPT(1)+HHT2U(2)*PPT(2))
        QTKT(3)=-((HHT3L(1)+HH3*HHT0(1))*QQT(1)+
     *            (HHT3L(2)+HH3*HHT0(2))*QQT(2)+
     *            HHT3U(1)*PPT(1)+HHT3U(2)*PPT(2))
        QTKT(4)=-((HHT1L(1)+HH1*HHT0(1))*QQT(3)+
     *            (HHT1L(2)+HH1*HHT0(2))*QQT(4)+
     *            HHT1U(1)*PPT(3)+HHT1U(2)*PPT(4))
        QTKT(5)=-((HHT2L(1)+HH2*HHT0(1))*QQT(3)+
     *            (HHT2L(2)+HH2*HHT0(2))*QQT(4)+
     *            HHT2U(1)*PPT(3)+HHT2U(2)*PPT(4))
        QTKT(6)=-((HHT3L(1)+HH3*HHT0(1))*QQT(3)+
     *            (HHT3L(2)+HH3*HHT0(2))*QQT(4)+
     *            HHT3U(1)*PPT(3)+HHT3U(2)*PPT(4))
        QTKT(7)=-HH3
        IF (IPTS(1).EQ.1) THEN
C         Initial point of the ray:
C         Initiating first integration:
          NSUM1=7
          IX=1
          NDER=1
          NFUN11=7
          NFUN21=7
          IFUN1T(1)=1
          IFUN1T(2)=3
          IFUN1T(3)=5
          IFUN1T(4)=2
          IFUN1T(5)=4
          IFUN1T(6)=6
          IFUN1T(7)=7
          FUN1T(1)=QTKT(1)
          FUN1T(2)=QTKT(2)
          FUN1T(3)=QTKT(3)
          FUN1T(4)=QTKT(4)
          FUN1T(5)=QTKT(5)
          FUN1T(6)=QTKT(6)
          FUN1T(7)=QTKT(7)
          IFUN2T(1)=1
          IFUN2T(2)=3
          IFUN2T(3)=5
          IFUN2T(4)=2
          IFUN2T(5)=4
          IFUN2T(6)=6
          IFUN2T(7)=7
C         Initiating second integration:
          NSUM2=6
          NFUN12=6
          NFUN22=6
          IFUN1K(1)=1
          IFUN1K(2)=2
          IFUN1K(3)=3
          IFUN1K(4)=4
          IFUN1K(5)=5
          IFUN1K(6)=6
          AUX=HHTC0S(1,1)*HHTC0S(2,2)-HHTC0S(2,1)*HHTC0S(1,2)
          IF (AUX.EQ.0) THEN
C           WAN-19
            CALL ERROR
     *           ('WAN-19: Zero determinant of second der. of Ham.')
C           This error should not appear.
          ENDIF
          GIJT(1,1)= HHTC0S(2,2)/AUX
          GIJT(2,2)= HHTC0S(1,1)/AUX
          GIJT(2,1)=-HHTC0S(2,1)/AUX
          GIJT(1,2)=-HHTC0S(1,2)/AUX
          TOUT(1,1)=-GIJT(1,1)*HHT1U(1)-GIJT(1,2)*HHT1U(2)
          TOUT(2,1)=-GIJT(2,1)*HHT1U(1)-GIJT(2,2)*HHT1U(2)
          TOUT(1,2)=-GIJT(1,1)*HHT2U(1)-GIJT(1,2)*HHT2U(2)
          TOUT(2,2)=-GIJT(2,1)*HHT2U(1)-GIJT(2,2)*HHT2U(2)
          TOUT(1,3)=-GIJT(1,1)*HHT3U(1)-GIJT(1,2)*HHT3U(2)
          TOUT(2,3)=-GIJT(2,1)*HHT3U(1)-GIJT(2,2)*HHT3U(2)
          FUN1K(1)=-(HHT1U(1)*TOUT(1,1)+HHT1U(2)*TOUT(2,1) +
     *               HHT1U(1)*TOUT(1,1)+HHT1U(2)*TOUT(2,1) +
     *    HHTC0S(1,1)*TOUT(1,1)*TOUT(1,1)+
     *    HHTC0S(2,1)*TOUT(2,1)*TOUT(1,1)+
     *    HHTC0S(1,2)*TOUT(1,1)*TOUT(2,1)+
     *    HHTC0S(2,2)*TOUT(2,1)*TOUT(2,1))
          FUN1K(2)=-(HHT1U(1)*TOUT(1,2)+HHT1U(2)*TOUT(2,2) +
     *               HHT2U(1)*TOUT(1,1)+HHT2U(2)*TOUT(2,1) +
     *    HHTC0S(1,1)*TOUT(1,1)*TOUT(1,2)+
     *    HHTC0S(2,1)*TOUT(2,1)*TOUT(1,2)+
     *    HHTC0S(1,2)*TOUT(1,1)*TOUT(2,2)+
     *    HHTC0S(2,2)*TOUT(2,1)*TOUT(2,2))
          FUN1K(3)=-(HHT2U(1)*TOUT(1,2)+HHT2U(2)*TOUT(2,2) +
     *               HHT2U(1)*TOUT(1,2)+HHT2U(2)*TOUT(2,2) +
     *    HHTC0S(1,1)*TOUT(1,2)*TOUT(1,2)+
     *    HHTC0S(2,1)*TOUT(2,2)*TOUT(1,2)+
     *    HHTC0S(1,2)*TOUT(1,2)*TOUT(2,2)+
     *    HHTC0S(2,2)*TOUT(2,2)*TOUT(2,2))
          FUN1K(4)=-(HHT1U(1)*TOUT(1,3)+HHT1U(2)*TOUT(2,3) +
     *               HHT3U(1)*TOUT(1,1)+HHT3U(2)*TOUT(2,1) +
     *    HHTC0S(1,1)*TOUT(1,1)*TOUT(1,3)+
     *    HHTC0S(2,1)*TOUT(2,1)*TOUT(1,3)+
     *    HHTC0S(1,2)*TOUT(1,1)*TOUT(2,3)+
     *    HHTC0S(2,2)*TOUT(2,1)*TOUT(2,3))
          FUN1K(5)=-(HHT2U(1)*TOUT(1,3)+HHT2U(2)*TOUT(2,3) +
     *               HHT3U(1)*TOUT(1,2)+HHT3U(2)*TOUT(2,2) +
     *    HHTC0S(1,1)*TOUT(1,2)*TOUT(1,3)+
     *    HHTC0S(2,1)*TOUT(2,2)*TOUT(1,3)+
     *    HHTC0S(1,2)*TOUT(1,2)*TOUT(2,3)+
     *    HHTC0S(2,2)*TOUT(2,2)*TOUT(2,3))
          FUN1K(6)=-(HHT3U(1)*TOUT(1,3)+HHT3U(2)*TOUT(2,3) +
     *               HHT3U(1)*TOUT(1,3)+HHT3U(2)*TOUT(2,3) +
     *    HHTC0S(1,1)*TOUT(1,3)*TOUT(1,3)+
     *    HHTC0S(2,1)*TOUT(2,3)*TOUT(1,3)+
     *    HHTC0S(1,2)*TOUT(1,3)*TOUT(2,3)+
     *    HHTC0S(2,2)*TOUT(2,3)*TOUT(2,3))
          IFUN2K(1)=1
          IFUN2K(2)=2
          IFUN2K(3)=3
          IFUN2K(4)=4
          IFUN2K(5)=5
          IFUN2K(6)=6
C         Output values:
          DTAU(1)=0.
          DTAU(2)=0.
          DTAU(3)=0.
          DTAU(4)=0.
          DTAU(5)=0.
          DTAU(6)=0.
          DTAU(7)=0.
          RETURN
        ENDIF
C       Eq. 51:
        CALL AP28(NSUM1,TT,IX,NDER,STEP,X1T,ISRF1T,NFUN11,IFUN1T,FUN1T,
     *            NFUN21,IFUN2T,QTKT)
C       Eq. 57:
        AUX=QQT(1)*QQT(4)-QQT(2)*QQT(3)
        IF (AUX.EQ.0) THEN
C         WAN-21
          CALL ERROR('WAN-21: Zero determinant of geom. spread.')
C         This error should not appear.
        ENDIF
        QIT(1)= QQT(4)/AUX
        QIT(4)= QQT(1)/AUX
        QIT(2)=-QQT(2)/AUX
        QIT(3)=-QQT(3)/AUX
        TOUT(1,1)=QIT(1)*TT(1)+QIT(2)*TT(2)
        TOUT(2,1)=QIT(3)*TT(1)+QIT(4)*TT(2)
        TOUT(1,2)=QIT(1)*TT(3)+QIT(2)*TT(4)
        TOUT(2,2)=QIT(3)*TT(3)+QIT(4)*TT(4)
        TOUT(1,3)=QIT(1)*TT(5)+QIT(2)*TT(6)
        TOUT(2,3)=QIT(3)*TT(5)+QIT(4)*TT(6)
C       Eq. 62:
        VK(1)=-(HHT1U(1)*TOUT(1,1)+HHT1U(2)*TOUT(2,1) +
     *          HHT1U(1)*TOUT(1,1)+HHT1U(2)*TOUT(2,1) +
     *  HHTC0S(1,1)*TOUT(1,1)*TOUT(1,1)+HHTC0S(2,1)*TOUT(2,1)*TOUT(1,1)+
     *  HHTC0S(1,2)*TOUT(1,1)*TOUT(2,1)+HHTC0S(2,2)*TOUT(2,1)*TOUT(2,1))
        VK(2)=-(HHT1U(1)*TOUT(1,2)+HHT1U(2)*TOUT(2,2) +
     *          HHT2U(1)*TOUT(1,1)+HHT2U(2)*TOUT(2,1) +
     *  HHTC0S(1,1)*TOUT(1,1)*TOUT(1,2)+HHTC0S(2,1)*TOUT(2,1)*TOUT(1,2)+
     *  HHTC0S(1,2)*TOUT(1,1)*TOUT(2,2)+HHTC0S(2,2)*TOUT(2,1)*TOUT(2,2))
        VK(3)=-(HHT2U(1)*TOUT(1,2)+HHT2U(2)*TOUT(2,2) +
     *          HHT2U(1)*TOUT(1,2)+HHT2U(2)*TOUT(2,2) +
     *  HHTC0S(1,1)*TOUT(1,2)*TOUT(1,2)+HHTC0S(2,1)*TOUT(2,2)*TOUT(1,2)+
     *  HHTC0S(1,2)*TOUT(1,2)*TOUT(2,2)+HHTC0S(2,2)*TOUT(2,2)*TOUT(2,2))
        VK(4)=-(HHT1U(1)*TOUT(1,3)+HHT1U(2)*TOUT(2,3) +
     *          HHT3U(1)*TOUT(1,1)+HHT3U(2)*TOUT(2,1) +
     *  HHTC0S(1,1)*TOUT(1,1)*TOUT(1,3)+HHTC0S(2,1)*TOUT(2,1)*TOUT(1,3)+
     *  HHTC0S(1,2)*TOUT(1,1)*TOUT(2,3)+HHTC0S(2,2)*TOUT(2,1)*TOUT(2,3))
        VK(5)=-(HHT2U(1)*TOUT(1,3)+HHT2U(2)*TOUT(2,3) +
     *          HHT3U(1)*TOUT(1,2)+HHT3U(2)*TOUT(2,2) +
     *  HHTC0S(1,1)*TOUT(1,2)*TOUT(1,3)+HHTC0S(2,1)*TOUT(2,2)*TOUT(1,3)+
     *  HHTC0S(1,2)*TOUT(1,2)*TOUT(2,3)+HHTC0S(2,2)*TOUT(2,2)*TOUT(2,3))
        VK(6)=-(HHT3U(1)*TOUT(1,3)+HHT3U(2)*TOUT(2,3) +
     *          HHT3U(1)*TOUT(1,3)+HHT3U(2)*TOUT(2,3) +
     *  HHTC0S(1,1)*TOUT(1,3)*TOUT(1,3)+HHTC0S(2,1)*TOUT(2,3)*TOUT(1,3)+
     *  HHTC0S(1,2)*TOUT(1,3)*TOUT(2,3)+HHTC0S(2,2)*TOUT(2,3)*TOUT(2,3))
C       Eq. 59:
        CALL AP28(NSUM2,T,IX,NDER,STEP,X1K,ISRF1K,NFUN12,IFUN1K,FUN1K,
     *            NFUN22,IFUN2K,VK)
C       Second order perturbations towards anisotropic rays:
        DTAU(1)=T(1)/2.
        DTAU(2)=T(2)
        DTAU(3)=T(3)/2.
C       Second order perturbation towards isotropic ray:
        DTAU(4)=T(4)
        DTAU(5)=T(5)
        DTAU(6)=T(6)/2.
C       Estimation of the isotropic ray travel time:
        DTAU(7)=TT(7)
        RETURN
      ENDIF
      END
C
C=======================================================================
C