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.
C
C Version: 5.80
C Date: 2004, June 11
C
C Coded by: Petr Bulant
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: bulant@seis.karlov.mff.cuni.cz
C
C=======================================================================
C
      SUBROUTINE WAN(LU1,LU2,LU3,GREEN,MGREEN,RPTS,IPTS,MAXPTS,NQ,
     *               ERR,LUQI)
C
C----------------------------------------------------------------------
      INTEGER LU1,LU2,LU3,LUQI
      INTEGER MGREEN,MAXPTS,NQ
      INTEGER    NQPTS
      PARAMETER (NQPTS=43)
      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
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,WACHAN,WAPROJ,
     *WAPERT,WASUM,WASUM4,WASUM5,
     *ERROR,RSEP1,RSEP3T,RSEP3R,RSEP3I,FORM2,LENGTH,EIGEN,
     *HIVD2,BLOCK,PARM3,AP00,AP21,TRANS
      INTEGER LENGTH
      REAL WACHAN
C     WAPTS,WAREAD,WACHRI,WAMAT,WACHAN,WAPROJ,WAPERT,WASUM,WASUM4,WASUM5
C           ... This file.
C     ERROR ... File error.for.
C     RSEP1,RSEP3I,RSEP3T,RSEP3R ...
C               File sep.for.
C     FORM2 ... File forms.for.
C     LENGTH... File length.for.
C     EIGEN ... File eigen.for.
C     HIVD2 ... File means.for
C     BLOCK ... File model.for
C     PARM3 ... File parm.for
C     AP00,AP21 ... File ap.for.
C     TRANS ... File trans.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
C     Auxiliary storage locations:
      INTEGER MPTS
      REAL PIGRA(2,2),PIGIA(2,2),PIGR(2,2),PIGI(2,2)
      REAL AR11,AI11,AR21,AI21,AR12,AI12,AR22,AI22
      REAL G11,G21,G12,G22
      REAL PI
      PARAMETER (PI=3.1415926)
      REAL GREENA(32)
      REAL TTCOR,FREQ,GAMA,AUX0,AUX1,AUX2
      INTEGER NPTS
      INTEGER I,I1,I2,I3,J,IRT
      INTEGER NFFT,NF,KQIPV,KQITT,KQIRAY
      REAL DT,FMIN,FMAX,DF,OF
      INTEGER IPTTMP,NYFTMP,ISRFFT
      REAL YI1TMP,YF1TMP,STEP,TISO,DTLIN,DTAU(3),QICOR(8)
      LOGICAL LPWAVE,LSWAVE
      CHARACTER*71 FORMQI
      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
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,'
        CALL RSEP3I('QIRAY',KQIRAY,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
      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.
        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.
        ENDIF
   30 CONTINUE
C
      IF (LPWAVE) THEN
C       P wave along the whole ray:
        Y(1)=TTCOR+YI(1)
        CALL AP21(KQIPV,RPTS(25,1),RPTS(25,NPTS),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(2)=0.
        IPTTMP=IPT
        YI1TMP=YI(1)
        NYFTMP=NYF
        ISRFFT=ISRFF
        YF1TMP=YF(1)
        YI(1) =RPTS(2,1)
        NYF   =1
        DTLIN=0.
        DO 100, I1=1,NPTS
          IF (IPTS(1,I1).EQ.2) STEP=RPTS(2,I1)-RPTS(2,1)
          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(2)
          ENDIF
          CALL WAPERT(RPTS(1,I1),IPTS(1,I1),NQPTS,STEP,DTAU)
          IF (KQIRAY.NE.0) THEN
            RPTS(21,I1)=RPTS(21,I1)+DTAU(2)-DTAU(1)
          ENDIF
 100    CONTINUE
C
        IF (LUQI.NE.0) THEN
          DTLIN=DTLIN/2.
          QICOR(1)=TISO
          QICOR(2)=TTCOR-DTLIN+DTAU(1)
          QICOR(3)=TTCOR+DTLIN+DTAU(2)
          QICOR(4)=TTCOR-DTLIN-TISO
          QICOR(5)=TTCOR+DTLIN-TISO
          QICOR(6)=DTAU(1)
          QICOR(7)=DTAU(2)
          QICOR(8)=0.25*(DTAU(1)-2.*DTAU(3)+DTAU(2))
          CALL FORM2(8,QICOR,QICOR,FORMQI(8:71))
          WRITE(LUQI,FORMQI) IREC,(' ',QICOR(I),I=1,8),' /'
        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(2)+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:
      DO 90, I2=0,NF-1
        FREQ=OF+I2*DF
        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:
            IRT=IRT+1
            IF (IRT.GT.NRT) THEN
C             WAN-17
              CALL ERROR('WAN-17: 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.
              ELSE
                AUX2=SIN(AUX0)/AUX0
              ENDIF
C             Matrix for this step along the ray:
              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)
            ELSE
C             P wave:
              GOTO 39
            ENDIF
          ENDIF
C         Matrix for all the steps along the ray to current point:
          CALL WAMAT(PIGRA,PIGIA,PIGR,PIGI)
   39     CONTINUE
   40   CONTINUE
C
C       Computing the Green function:
C       Transforming the propagator matrix from eigenvectors
C       to polarization vectors:
        Y(1)=TTCOR+YI(1)
        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,RPTS(25,1),RPTS(25,NPTS),GREENA)
        IF (I2.EQ.0) THEN
          DO 50 I=1,14
            GREEN(I)=GREENA(I)
   50     CONTINUE
          NQ=14
        ENDIF
        J=I2*18
        DO 60 I=15,32
          GREEN(J+I)=GREENA(I)
   60   CONTINUE
        NQ=NQ+18
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. If the value of DELTFI is greater than
C prescribed limit new points are added using subroutine HIVD2 of
C the file 'means.for'.
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     NPTS ... Number of points stored in RPTS (IPTS).
C
C External functions required:
      EXTERNAL WACHAN
      REAL WACHAN
C.......................................................................
C
C     Auxiliary storage locations:
      INTEGER NPTSE
      INTEGER MQANT,MNEWP,NNEWP
      PARAMETER (MNEWP=5)
      PARAMETER (MQANT=26)
      REAL ROLDP(MQANT),RNEWP(MQANT,MNEWP)
      INTEGER KNEWP(MNEWP)
      INTEGER MAXDIV,KQIDT
      PARAMETER (MAXDIV=MNEWP-1)
      REAL AA(10,21),RHO,QQ(21)
      REAL EE(9),DER(9)
      REAL ERRWAN,DEFI,DELTFI,PI
      PARAMETER (PI=3.1415926)
      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
      INTEGER I1,I2
      SAVE NNEWP
      DATA NNEWP/0/,KNEWP/MNEWP*0/
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     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     MAXDIV .. The distance between two points on the rays must not
C               be smaller than original distance between the points
C               computed by CRT divided by MAXDIV  when adding new
C               points on the ray.
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
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.
      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))
      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)
      NNEWP=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: Array RPTS small')
C         The dimension of the array RPTS is given by the dimension MRAM
C         in the include file ram.inc.
        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)
        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 ((ROLDP(14).EQ.1.).AND.(ROLDP(15).EQ.0.).AND.
     *        (ROLDP(16).EQ.0.).AND.(ROLDP(17).EQ.1.)) THEN
C           New point is the first anisotropic point on the ray,
C           angular change is not to be computed:
            DELTFI=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
        DELTFI=WACHAN(ROLDP(14),RNEWP(11,I1))
        DEFI=1./SQRT(RNEWP(12,I1))-1./SQRT(RNEWP(13,I1))
        DEFI=DEFI-1./SQRT(ROLDP(12))+1./SQRT(ROLDP(13))
        DEFI=DEFI*PI*(OF+FLOAT(NF-1)*DF)*RPTS(2,NPTS)/6.
        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 (KNEWP(NNEWP).GE.MAXDIV-1) 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 than MAXDIV divisions of this interval is
C         needed to keep the change under the prescribed limit. Try to
C         decrease the parameter STORE,
C         or to increase MAXDIV or 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       Material parameters:
        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
          DELTFI=WACHAN(ROLDP(14),RNEWP(11,NNEWP))
          DEFI=1./SQRT(RNEWP(12,NNEWP))-1./SQRT(RNEWP(13,NNEWP))
          DEFI=DEFI-1./SQRT(ROLDP(12))+1./SQRT(ROLDP(13))
          DEFI=DEFI*PI*(OF+FLOAT(NF-1)*DF)*RPTS(2,NPTS)/6.
          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         of in the include file 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)
      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     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
C     Auxiliary storage locations:
      INTEGER ISRFA,ISRF2,KODNEW,ICBNEW,IEN,ISB2A,ISB2B,ICB2A,ICB2B,I1
      INTEGER IYFA(12)
      REAL YLFA(6),YFA(35),YYFA(5),HUI(5)
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     Derivatives of the velocity:
      RPTS(22,NPTS)=YLI(4)
      RPTS(23,NPTS)=YLI(5)
      RPTS(24,NPTS)=YLI(6)
C     S-wave velocity:
      RPTS(34,NPTS)=YLI(2)
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
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 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       Derivatives of the velocity:
        RPTS(22,NPTS)=YLF(4)
        RPTS(23,NPTS)=YLF(5)
        RPTS(24,NPTS)=YLF(6)
C       S-wave velocity:
        RPTS(34,NPTS)=YLF(2)
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
        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
          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 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 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     Derivatives of the velocity:
      RPTS(22,NPTS)=YL(4)
      RPTS(23,NPTS)=YL(5)
      RPTS(24,NPTS)=YL(6)
C     S-wave velocity:
      RPTS(34,NPTS)=YL(2)
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
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE WACHRI(P1,P2,P3,B11,B12,B22,B13,B23,B33,
     *                  B14,B24,B34,B44,B15,B25,B35,B45,B55,
     *                  B16,B26,B36,B46,B56,B66,
     *                  G1,G2,G3,EE,DER)
C
C----------------------------------------------------------------------
C Subroutine to compute the Christoffel matrix, its eigenvalues
C and eigenvectors.
      REAL P1,P2,P3
      REAL B11,B12,B22,B13,B23,B33,B14,B24,B34,B44
      REAL B15,B25,B35,B45,B55,B16,B26,B36,B46,B56,B66
      REAL G1,G2,G3,EE(9),DER(9)
C
C Input:
C     P1,P2,P3... Slowness vector.
C     Bii ... Values of real parts of 21 reduced
C             (divided by the density) elastic parameters.
C Output:
C     G1,G2,G3 ... Eigenvalues of the Christoffel matrix.
C     EE  ... Eigenvectors of the Christoffel matrix.
C     DER ... Derivatives dx/dt=dH/dp=Aijkl Ej Ek pl stored columnwise
C             for qP, qS1 and qS2 waves.
C
C-----------------------------------------------------------------------
C
      REAL A11,A12,A13,A14,A21,A22,A23,A24,A31,A32,A33,A34,A44
      REAL A15,A25,A35,A45,A55,A16,A26,A36,A46,A56,A66
      REAL A1111,A2111,A3111,A1211,A2211,A3211,A1311,A2311,A3311
      REAL A1121,A2121,A3121,A1221,A2221,A3221,A1321,A2321,A3321
      REAL A1131,A2131,A3131,A1231,A2231,A3231,A1331,A2331,A3331
      REAL A1112,A2112,A3112,A1212,A2212,A3212,A1312,A2312,A3312
      REAL A1122,A2122,A3122,A1222,A2222,A3222,A1322,A2322,A3322
      REAL A1132,A2132,A3132,A1232,A2232,A3232,A1332,A2332,A3332
      REAL A1113,A2113,A3113,A1213,A2213,A3213,A1313,A2313,A3313
      REAL A1123,A2123,A3123,A1223,A2223,A3223,A1323,A2323,A3323
      REAL A1133,A2133,A3133,A1233,A2233,A3233,A1333,A2333,A3333
      EQUIVALENCE (A11,A1111)
      EQUIVALENCE (A22,A2222)
      EQUIVALENCE (A33,A3333)
      EQUIVALENCE (A16,A1112,A1121,A1211,A2111)
      EQUIVALENCE (A26,A2221,A2212,A2122,A1222)
      EQUIVALENCE (A15,A1113,A1131,A1311,A3111)
      EQUIVALENCE (A35,A3331,A3313,A3133,A1333)
      EQUIVALENCE (A24,A2223,A2232,A2322,A3222)
      EQUIVALENCE (A34,A3332,A3323,A3233,A2333)
      EQUIVALENCE (A23,A2233,A3322)
      EQUIVALENCE (A13,A1133,A3311)
      EQUIVALENCE (A12,A1122,A2211)
      EQUIVALENCE (A44,            A2323,A3232,A2332,A3223)
      EQUIVALENCE (A55,            A1313,A3131,A1331,A3113)
      EQUIVALENCE (A66,            A1212,A2121,A1221,A2112)
      EQUIVALENCE (A14,A1123,A1132,A2311,A3211)
      EQUIVALENCE (A25,A2213,A2231,A1322,A3122)
      EQUIVALENCE (A36,A3312,A3321,A1233,A2133)
      EQUIVALENCE (A56,A1321,A3112,A2113,A1231,A1213,A2131,A1312,A3121)
      EQUIVALENCE (A46,A2312,A3221,A1223,A2132,A2123,A1232,A2321,A3212)
      EQUIVALENCE (A45,A3213,A2331,A1332,A3123,A3132,A1323,A3231,A2313)
      REAL GAMMA(6),E11,E21,E31,E12,E22,E32,E13,E23,E33
C     EQUIVALENCE (GAMMA(1),G1),(GAMMA(3),G2),(GAMMA(6),G3)
C     EQUIVALENCE (EE(1),E11),(EE(4),E12),(EE(7),E13)
C     EQUIVALENCE (EE(2),E21),(EE(5),E22),(EE(8),E23)
C     EQUIVALENCE (EE(3),E31),(EE(6),E32),(EE(9),E33)
      REAL A111,A112,A121,A122,A113,A123,A131,A132,A133
      REAL A211,A212,A221,A222,A213,A223,A231,A232,A233
      REAL A311,A312,A322,A313,A321,A323,A331,A332,A333
      REAL AUX
      INTEGER KQICHM
      DATA    KQICHM/-1/
C
C     GAMMA,G1,G2,G3...Christoffel matrix, later its eigenvalues.
C        (E11,E12,E13)
C     EE=(E21,E22,E23)... Eigenvectors of the christoffel matrix.
C        (E31,E32,E33)
C     A111,A211,A311,A112,A212,A312,A122,A222,A322,A113,A213,A313,A123,
C     A223,A323,A133,A233,A333... A(I,J,K,L)*P(L) summed over L.
C     A11,A21,A31,A12,A22,A32,A13,A23,A33 ... Aijk*Ek
C
C.......................................................................
C
      A11=B11
      A22=B22
      A33=B33
      A16=B16
      A26=B26
      A15=B15
      A35=B35
      A24=B24
      A34=B34
      A23=B23
      A13=B13
      A12=B12
      A44=B44
      A55=B55
      A66=B66
      A14=B14
      A25=B25
      A36=B36
      A56=B56
      A46=B46
      A45=B45
C     Christoffel matrix:
      A111=A1111*P1+A1112*P2+A1113*P3
      A112=A1121*P1+A1122*P2+A1123*P3
      A121=A1211*P1+A1212*P2+A1213*P3
      A122=A1221*P1+A1222*P2+A1223*P3
      A113=A1131*P1+A1132*P2+A1133*P3
      A123=A1231*P1+A1232*P2+A1233*P3
      A131=A1311*P1+A1312*P2+A1313*P3
      A132=A1321*P1+A1322*P2+A1323*P3
      A133=A1331*P1+A1332*P2+A1333*P3
      A211=A2111*P1+A2112*P2+A2113*P3
      A212=A2121*P1+A2122*P2+A2123*P3
      A221=A2211*P1+A2212*P2+A2213*P3
      A222=A2221*P1+A2222*P2+A2223*P3
      A213=A2131*P1+A2132*P2+A2133*P3
      A223=A2231*P1+A2232*P2+A2233*P3
      A231=A2311*P1+A2312*P2+A2313*P3
      A232=A2321*P1+A2322*P2+A2323*P3
      A233=A2331*P1+A2332*P2+A2333*P3
      A311=A3111*P1+A3112*P2+A3113*P3
      A312=A3121*P1+A3122*P2+A3123*P3
      A322=A3221*P1+A3222*P2+A3223*P3
      A313=A3131*P1+A3132*P2+A3133*P3
      A321=A3211*P1+A3212*P2+A3213*P3
      A323=A3231*P1+A3232*P2+A3233*P3
      A331=A3311*P1+A3312*P2+A3313*P3
      A332=A3321*P1+A3322*P2+A3323*P3
      A333=A3331*P1+A3332*P2+A3333*P3
      GAMMA(1)=P1*A111+P2*A211+P3*A311
      GAMMA(2)=P1*A112+P2*A212+P3*A312
      GAMMA(3)=P1*A122+P2*A222+P3*A322
      GAMMA(4)=P1*A113+P2*A213+P3*A313
      GAMMA(5)=P1*A123+P2*A223+P3*A323
      GAMMA(6)=P1*A133+P2*A233+P3*A333
C
C     Quasi-isotropic modification of the Christoffel matrix:
      IF(KQICHM.LT.0) THEN
        CALL RSEP3I('QICHM',KQICHM,0)
      END IF
      IF(KQICHM.GT.0) THEN
        AUX=SQRT(P1*P1+P2*P2+P3*P3)
        E11=P1/AUX
        E21=P2/AUX
        E31=P3/AUX
        E12=GAMMA(1)*E11+GAMMA(2)*E21+GAMMA(4)*E31
        E22=GAMMA(2)*E11+GAMMA(3)*E21+GAMMA(5)*E31
        E32=GAMMA(4)*E11+GAMMA(5)*E21+GAMMA(6)*E31
        AUX=2.*(E11*E12+E21*E22+E31*E32)
        GAMMA(1)=GAMMA(1)-E11*E12-E12*E11+AUX*E11*E11
        GAMMA(2)=GAMMA(2)-E11*E22-E12*E21+AUX*E11*E21
        GAMMA(3)=GAMMA(3)-E21*E22-E22*E21+AUX*E21*E21
        GAMMA(4)=GAMMA(4)-E11*E32-E12*E31+AUX*E11*E31
        GAMMA(5)=GAMMA(5)-E21*E32-E22*E31+AUX*E21*E31
        GAMMA(6)=GAMMA(6)-E31*E32-E32*E31+AUX*E31*E31
      END IF
C
C     Selecting eigenvalue and eigenvector of the Christoffel matrix:
      CALL EIGEN(GAMMA,EE,3,0)
      G1=GAMMA(1)
      G2=GAMMA(3)
      G3=GAMMA(6)
      E11=EE(1)
      E21=EE(2)
      E31=EE(3)
      E12=EE(4)
      E22=EE(5)
      E32=EE(6)
      E13=EE(7)
      E23=EE(8)
      E33=EE(9)
      IF (G3.LT.0.) THEN
C       WAN-11
        CALL ERROR('WAN-11: Negative eigenvalue of Christoffel matrix')
C       This error should not appear.
      END IF
      AUX=E11*E11+E21*E21+E31*E31
      IF (ABS(AUX-1.).GT.0.00001) THEN
C       WAN-12
        CALL ERROR('WAN-12: Eigenvector is not normalized')
C       This error should not appear.
      ENDIF
      AUX=E12*E12+E22*E22+E32*E32
      IF (ABS(AUX-1.).GT.0.00001) THEN
C       WAN-13
        CALL ERROR('WAN-13: Eigenvector is not normalized')
C       This error should not appear.
      ENDIF
      AUX=E13*E13+E23*E23+E33*E33
      IF (ABS(AUX-1.).GT.0.00001) THEN
C       WAN-14
        CALL ERROR('WAN-14: Eigenvector is not normalized')
C       This error should not appear.
      ENDIF
C
C     Computation of derivatives dx/dt:
      A11=  A111*E11+A112*E21+A113*E31
      A21=  A211*E11+A212*E21+A213*E31
      A31=  A311*E11+A312*E21+A313*E31
      A12=  A121*E11+A122*E21+A123*E31
      A22=  A221*E11+A222*E21+A223*E31
      A32=  A321*E11+A322*E21+A323*E31
      A13=  A131*E11+A132*E21+A133*E31
      A23=  A231*E11+A232*E21+A233*E31
      A33=  A331*E11+A332*E21+A333*E31
      DER(1)=A11*E11+ A12*E21+ A13*E31
      DER(2)=A21*E11+ A22*E21+ A23*E31
      DER(3)=A31*E11+ A32*E21+ A33*E31
      A11=  A111*E12+A112*E22+A113*E32
      A21=  A211*E12+A212*E22+A213*E32
      A31=  A311*E12+A312*E22+A313*E32
      A12=  A121*E12+A122*E22+A123*E32
      A22=  A221*E12+A222*E22+A223*E32
      A32=  A321*E12+A322*E22+A323*E32
      A13=  A131*E12+A132*E22+A133*E32
      A23=  A231*E12+A232*E22+A233*E32
      A33=  A331*E12+A332*E22+A333*E32
      DER(4)=A11*E12+ A12*E22+ A13*E32
      DER(5)=A21*E12+ A22*E22+ A23*E32
      DER(6)=A31*E12+ A32*E22+ A33*E32
      A11=  A111*E13+A112*E23+A113*E33
      A21=  A211*E13+A212*E23+A213*E33
      A31=  A311*E13+A312*E23+A313*E33
      A12=  A121*E13+A122*E23+A123*E33
      A22=  A221*E13+A222*E23+A223*E33
      A32=  A321*E13+A322*E23+A323*E33
      A13=  A131*E13+A132*E23+A133*E33
      A23=  A231*E13+A232*E23+A233*E33
      A33=  A331*E13+A332*E23+A333*E33
      DER(7)=A11*E13+ A12*E23+ A13*E33
      DER(8)=A21*E13+ A22*E23+ A23*E33
      DER(9)=A31*E13+ A32*E23+ A33*E33
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE WAMAT(A,B,C,D)
C
C----------------------------------------------------------------------
C Subroutine to compute the product of two 2x2 complex matrices.
C The second matrix (C+iD) is destroyed in the computation.
      REAL A(2,2),B(2,2),C(2,2),D(2,2)
C Input:
C     A,B,C,D ... Real and imaginary parts of the two matrices.
C Output:
C     C,D ... Real and imaginary parts of resulting matrix.
C
C.......................................................................
C     Auxiliary storage locations:
      REAL CR11,CR21,CR12,CR22,CI11,CI21,CI12,CI22
C.......................................................................
      CR11=A(1,1)*C(1,1)-B(1,1)*D(1,1)+A(1,2)*C(2,1)-B(1,2)*D(2,1)
      CR21=A(2,1)*C(1,1)-B(2,1)*D(1,1)+A(2,2)*C(2,1)-B(2,2)*D(2,1)
      CR12=A(1,1)*C(1,2)-B(1,1)*D(1,2)+A(1,2)*C(2,2)-B(1,2)*D(2,2)
      CR22=A(2,1)*C(1,2)-B(2,1)*D(1,2)+A(2,2)*C(2,2)-B(2,2)*D(2,2)
C
      CI11=A(1,1)*D(1,1)+B(1,1)*C(1,1)+A(1,2)*D(2,1)+B(1,2)*C(2,1)
      CI21=A(2,1)*D(1,1)+B(2,1)*C(1,1)+A(2,2)*D(2,1)+B(2,2)*C(2,1)
      CI12=A(1,1)*D(1,2)+B(1,1)*C(1,2)+A(1,2)*D(2,2)+B(1,2)*C(2,2)
      CI22=A(2,1)*D(1,2)+B(2,1)*C(1,2)+A(2,2)*D(2,2)+B(2,2)*C(2,2)
C
      C(1,1)=CR11
      C(2,1)=CR21
      C(1,2)=CR12
      C(2,2)=CR22
      D(1,1)=CI11
      D(2,1)=CI21
      D(1,2)=CI12
      D(2,2)=CI22
C
      RETURN
      END
C
C=======================================================================
C
      REAL FUNCTION WACHAN(G,H)
C
C----------------------------------------------------------------------
C Subroutine to compute the smallest angle between the two-dimensional
C vector G and vectors U,V,-U,-V.
C The subroutine reorganizes the vectors U and V in such way,
C that the pair U,V is equal to pair G,H rotated with angle WACHAN.
C The real numbers R1 (associated with U) and R2 (associated with V),
C and the vectors PU (associated with U) and PV (associated with V),
C are reorganized in the same way.
C
      REAL G(4),H(13)
C     G(1:4)=G1,G2,H1,H2
C     H(1:7)=R1,R2,R3,U1,U2,V1,V2,PU1,PU2,PU3,PV1,PV2,PV3
C Input:
C     G1,G2,H1,H2 ... A pair of two-dimensional orthonormal vectors.
C     U1,U2,V1,V2 ... A pair of two-dimensional orthonormal vectors.
C     PU1,PU2,PU3,PV1,PV2,PV3 ... A pair of three-dimensional vectors
C                     associated to U and V.
C     R1,R2  ...  Real numbers associated to U and V.
C Output:
C     WACHAN  ... The smallest one from the angles between vector G
C                 and vectors U,V,-U,-V.
C     U1,U2,V1,V2 ... Selection from U,V,-U,-V in such way, that
C                     the pair U,V is equal to pair G,H rotated
C                     with angle WACHAN.
C     PU1,PU2,PU3,PV1,PV2,PV3 ... A pair of three-dimensional vectors
C                     associated to U and V.
C     R1,R2  ...  Real numbers associeted to U and V.
C
C.......................................................................
C
C     Auxiliary storage locations:
      REAL SP1,SP2,SP3,SP4,A1,A2,A3,AUX
      REAL G1,G2,H1,H2,U1,U2,V1,V2,PU1,PU2,PU3,PV1,PV2,PV3,R1,R2
C
C-----------------------------------------------------------------------
      G1=G(1)
      G2=G(2)
      H1=G(3)
      H2=G(4)
      R1=H(2)
      R2=H(3)
      U1=H(4)
      U2=H(5)
      V1=H(6)
      V2=H(7)
      PU1=H(8)
      PU2=H(9)
      PU3=H(10)
      PV1=H(11)
      PV2=H(12)
      PV3=H(13)
C
      SP1=ABS(1 -( G1*U1+G2*U2))
      SP2=ABS(1 -(-G1*U1-G2*U2))
      SP3=ABS(1 -( G1*V1+G2*V2))
      SP4=ABS(1 -(-G1*V1-G2*V2))
      AUX=AMIN1(SP1,SP2,SP3,SP4)
      IF (AUX.EQ.SP1) THEN
C       No action.
      ELSEIF (AUX.EQ.SP2) THEN
        U1=-U1
        U2=-U2
        PU1=-PU1
        PU2=-PU2
        PU3=-PU3
      ELSEIF (AUX.EQ.SP3) THEN
        A1=U1
        A2=U2
        U1=V1
        U2=V2
        V1=A1
        V2=A2
        A1=PU1
        A2=PU2
        A3=PU3
        PU1=PV1
        PU2=PV2
        PU3=PV3
        PV1=A1
        PV2=A2
        PV3=A3
        AUX=R1
        R1=R2
        R2=AUX
      ELSEIF (AUX.EQ.SP4) THEN
        A1=U1
        A2=U2
        U1=-V1
        U2=-V2
        V1=A1
        V2=A2
        A1=PU1
        A2=PU2
        A3=PU3
        PU1=-PV1
        PU2=-PV2
        PU3=-PV3
        PV1=A1
        PV2=A2
        PV3=A3
        AUX=R1
        R1=R2
        R2=AUX
      ENDIF
      SP1=ABS(1 - ( H1*V1+H2*V2))
      SP2=ABS(1 - (-H1*V1-H2*V2))
      AUX=AMIN1(SP1,SP2)
      IF (AUX.EQ.SP1) THEN
C     No action.
      ELSEIF (AUX.EQ.SP2) THEN
        V1=-V1
        V2=-V2
        PV1=-PV1
        PV2=-PV2
        PV3=-PV3
      ENDIF
      WACHAN=ASIN(0.5*((G1+U1)*(H1-V1)+(G2+U2)*(H2-V2)))
      H(2)=R1
      H(3)=R2
      H(4)=U1
      H(5)=U2
      H(6)=V1
      H(7)=V2
      H(8)=PU1
      H(9)=PU2
      H(10)=PU3
      H(11)=PV1
      H(12)=PV2
      H(13)=PV3
      RETURN
      END
C
C
C=======================================================================
C
      SUBROUTINE WAPROJ(P1,P2,P3,E1,E2,E3,G1,G2,G3,H1,H2,H3,
     *                  GE1,GE2,HE1,HE2)
C
C----------------------------------------------------------------------
C Subroutine to project vectors G and H to the plane defined by vector E
C and vector PxE.
C
      REAL P1,P2,P3,E1,E2,E3,G1,G2,G3,H1,H2,H3,GE1,GE2,HE1,HE2
C Input:
C     P1,P2,P3,E1,E2,E3 ... Vectors defining the plane.
C     G1,G2,G3,H1,H2,H3 ... Vectors to be projected.
C Output:
C     GE1,GE2,HE1,HE2 ... Projected vectors.
C
C.......................................................................
C
C     Auxiliary storage locations:
      REAL F1,F2,F3,AUX
C
C-----------------------------------------------------------------------
C     Second vector defining the plane:
      F1=P2*E3-E2*P3
      F2=P3*E1-E3*P1
      F3=P1*E2-E1*P2
      AUX=SQRT(F1*F1+F2*F2+F3*F3)
      F1=F1/AUX
      F2=F2/AUX
      F3=F3/AUX
C     Projecting vectors G and H to the plane defined by E and F:
      GE1=E1*G1+E2*G2+E3*G3
      GE2=F1*G1+F2*G2+F3*G3
      HE1=E1*H1+E2*H2+E3*H3
      HE2=F1*H1+F2*H2+F3*H3
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(3)
      INTEGER IPTS(NQPTS)
C
C Input:
C     IPTS,RPTS,NPTS... 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 STORE.
C Output:
C     DTAU... Second-order travel-time corrections.
C
C.......................................................................
C Subroutines and external functions required:
      EXTERNAL WASUM,WASUM4,WASUM5
      REAL WASUM,WASUM4,WASUM5
C     WASUM,WASUM4,WASUM5 ... This file.
C-----------------------------------------------------------------------
C
      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),
     *     QQT(4),PPT(4),QIT(4),
     *     HH0(3),HHT0(3),HH1,HH2,
     *     HH1L(3),HH2L(3),HH1U(3),HH2U(3),
     *     HHT1L(3),HHT2L(3),HHT1U(3),HHT2U(3),
     *     QTKT(4),TT(4),TAUT(4),VK(3),T(3)
      INTEGER ISRF1T,IFUN1T(4),IFUN2T(4),ISRF1K,IFUN1K(3),IFUN2K(3),
     *     NSUM1,NSUM2,IX,NDER,NFUN11,NFUN12,NFUN21,NFUN22
      REAL AUX,X1T,FUN1T(4),X1K,FUN1K(3)
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.......................................................................
C
      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)
      H1(1)=RPTS(9)
      H1(2)=RPTS(10)
      H1(3)=RPTS(11)
      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)
      AUX=SQRT(P(1)**2+P(2)**2+P(3)**2)
      IF (AUX.EQ.0) THEN
C       WAN-18
        CALL ERROR('WAN-18: Polarization vector equal zero.')
C       This error should not appear.
      ENDIF
      H3(1)=P(1)/AUX
      H3(2)=P(2)/AUX
      H3(3)=P(3)/AUX
      H2(1)=H3(2)*H1(3)-H3(3)*H1(2)
      H2(2)=H3(3)*H1(1)-H3(1)*H1(3)
      H2(3)=H3(1)*H1(2)-H3(2)*H1(1)
      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)=-1.*(HH1**3)*WASUM4(AAA,E1,E1,P,1)
      HH1U(2)=-1.*(HH1**3)*WASUM4(AAA,E1,E1,P,2)
      HH1U(3)=-1.*(HH1**3)*WASUM4(AAA,E1,E1,P,3)
      HH2U(1)=-1.*(HH2**3)*WASUM4(AAA,E2,E2,P,1)
      HH2U(2)=-1.*(HH2**3)*WASUM4(AAA,E2,E2,P,2)
      HH2U(3)=-1.*(HH2**3)*WASUM4(AAA,E2,E2,P,3)
      HHT1L(1)=WASUM(HH1L,H1)
      HHT1L(2)=WASUM(HH1L,H2)
      HHT1L(3)=WASUM(HH1L,H3)
      HHT2L(1)=WASUM(HH2L,H1)
      HHT2L(2)=WASUM(HH2L,H2)
      HHT2L(3)=WASUM(HH2L,H3)
      HHT1U(1)=WASUM(HH1U,H1)
      HHT1U(2)=WASUM(HH1U,H2)
      HHT1U(3)=WASUM(HH1U,H3)
      HHT2U(1)=WASUM(HH2U,H1)
      HHT2U(2)=WASUM(HH2U,H2)
      HHT2U(3)=WASUM(HH2U,H3)
      QTKT(1)=-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)=-1.*((HHT2L(1)+HH2*HHT0(1))*QQT(1)+
     *             (HHT2L(2)+HH2*HHT0(2))*QQT(2)+
     *             HHT2U(1)*PPT(1)+HHT2U(2)*PPT(2))
      QTKT(3)=-1.*((HHT1L(1)+HH1*HHT0(1))*QQT(3)+
     *             (HHT1L(2)+HH1*HHT0(2))*QQT(4)+
     *             HHT1U(1)*PPT(3)+HHT1U(2)*PPT(4))
      QTKT(4)=-1.*((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)=-1.*HHT1U(1)/AUX
        TAUT(2)=-1.*HHT1U(2)/AUX
        TAUT(3)=-1.*HHT2U(1)/AUX
        TAUT(4)=-1.*HHT2U(2)/AUX
        FUN1K(1)=-1.*(2.*(HHT1U(1)*TAUT(1)+HHT1U(2)*TAUT(2)) +
     *                (TAUT(1)*TAUT(1)+TAUT(2)*TAUT(2)))
        FUN1K(2)=-1.*(2.*(HHT2U(1)*TAUT(3)+HHT2U(2)*TAUT(4)) +
     *                (TAUT(3)*TAUT(3)+TAUT(4)*TAUT(4)))
        FUN1K(3)=-1.*(HHT1U(1)*TAUT(3)+HHT1U(2)*TAUT(4) +
     *                HHT2U(1)*TAUT(1)+HHT2U(2)*TAUT(2) +
     *                (TAUT(1)*TAUT(3)+TAUT(2)*TAUT(4)))
        IFUN2K(1)=1
        IFUN2K(2)=2
        IFUN2K(3)=3
C       Output values:
        DTAU(1)=0.
        DTAU(2)=0.
        DTAU(3)=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-19
        CALL ERROR('WAN-19: 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)=-1.*(2.*(HHT1U(1)*TAUT(1)+HHT1U(2)*TAUT(2)) +
     *           (TAUT(1)*TAUT(1)+TAUT(2)*TAUT(2))*V0**2)
      VK(2)=-1.*(2.*(HHT2U(1)*TAUT(3)+HHT2U(2)*TAUT(4)) +
     *           (TAUT(3)*TAUT(3)+TAUT(4)*TAUT(4))*V0**2)
      VK(3)=-1.*(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)
      DTAU(1)=T(1)/2.
      DTAU(2)=T(2)/2.
      DTAU(3)=T(3)/2.
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE WAIJKL(A,B,I)
C
C----------------------------------------------------------------------
      REAL A(10,21),B(3,3,3,3)
      INTEGER I
      B(1,1,1,1)=A(I,1)
      B(1,1,2,2)=A(I,2)
      B(2,2,1,1)=A(I,2)
      B(2,2,2,2)=A(I,3)
      B(1,1,3,3)=A(I,4)
      B(3,3,1,1)=A(I,4)
      B(2,2,3,3)=A(I,5)
      B(3,3,2,2)=A(I,5)
      B(3,3,3,3)=A(I,6)
      B(1,1,2,3)=A(I,7)
      B(1,1,3,2)=A(I,7)
      B(2,3,1,1)=A(I,7)
      B(3,2,1,1)=A(I,7)
      B(2,2,2,3)=A(I,8)
      B(2,2,3,2)=A(I,8)
      B(2,3,2,2)=A(I,8)
      B(3,2,2,2)=A(I,8)
      B(3,3,3,2)=A(I,9)
      B(3,3,2,3)=A(I,9)
      B(3,2,3,3)=A(I,9)
      B(2,3,3,3)=A(I,9)
      B(2,3,2,3)=A(I,10)
      B(3,2,3,2)=A(I,10)
      B(2,3,3,2)=A(I,10)
      B(3,2,2,3)=A(I,10)
      B(1,1,1,3)=A(I,11)
      B(1,1,3,1)=A(I,11)
      B(1,3,1,1)=A(I,11)
      B(3,1,1,1)=A(I,11)
      B(2,2,1,3)=A(I,12)
      B(2,2,3,1)=A(I,12)
      B(1,3,2,2)=A(I,12)
      B(3,1,2,2)=A(I,12)
      B(3,3,3,1)=A(I,13)
      B(3,3,1,3)=A(I,13)
      B(3,1,3,3)=A(I,13)
      B(1,3,3,3)=A(I,13)
      B(3,2,1,3)=A(I,14)
      B(2,3,3,1)=A(I,14)
      B(1,3,3,2)=A(I,14)
      B(3,1,2,3)=A(I,14)
      B(3,1,3,2)=A(I,14)
      B(1,3,2,3)=A(I,14)
      B(3,2,3,1)=A(I,14)
      B(2,3,1,3)=A(I,14)
      B(1,3,1,3)=A(I,15)
      B(3,1,3,1)=A(I,15)
      B(1,3,3,1)=A(I,15)
      B(3,1,1,3)=A(I,15)
      B(1,1,1,2)=A(I,16)
      B(1,1,2,1)=A(I,16)
      B(1,2,1,1)=A(I,16)
      B(2,1,1,1)=A(I,16)
      B(2,2,2,1)=A(I,17)
      B(2,2,1,2)=A(I,17)
      B(2,1,2,2)=A(I,17)
      B(1,2,2,2)=A(I,17)
      B(3,3,1,2)=A(I,18)
      B(3,3,2,1)=A(I,18)
      B(1,2,3,3)=A(I,18)
      B(2,1,3,3)=A(I,18)
      B(2,3,1,2)=A(I,19)
      B(3,2,2,1)=A(I,19)
      B(1,2,2,3)=A(I,19)
      B(2,1,3,2)=A(I,19)
      B(2,1,2,3)=A(I,19)
      B(1,2,3,2)=A(I,19)
      B(2,3,2,1)=A(I,19)
      B(3,2,1,2)=A(I,19)
      B(1,3,2,1)=A(I,20)
      B(3,1,1,2)=A(I,20)
      B(2,1,1,3)=A(I,20)
      B(1,2,3,1)=A(I,20)
      B(1,2,1,3)=A(I,20)
      B(2,1,3,1)=A(I,20)
      B(1,3,1,2)=A(I,20)
      B(3,1,2,1)=A(I,20)
      B(1,2,1,2)=A(I,21)
      B(2,1,2,1)=A(I,21)
      B(1,2,2,1)=A(I,21)
      B(2,1,1,2)=A(I,21)
      RETURN
      END
C
C=======================================================================
C
      REAL FUNCTION WASUM(A,B)
C
C----------------------------------------------------------------------
      REAL A(3),B(3)
      INTEGER I
      WASUM=0.
      DO 10, I=1,3
        WASUM=WASUM+A(I)*B(I)
  10  CONTINUE
      RETURN
      END
C
C=======================================================================
C
      REAL FUNCTION WASUM5(A,B,C,D,E)
C
C----------------------------------------------------------------------
      REAL A(3,3,3,3),B(3),C(3),D(3),E(3)
      INTEGER I,J,K,L
      WASUM5=0.
      DO 13, I=1,3
        DO 12, J=1,3
          DO 11, K=1,3
            DO 10, L=1,3
              WASUM5=WASUM5+A(I,J,K,L)*B(I)*C(J)*D(K)*E(L)
  10        CONTINUE
  11      CONTINUE
  12    CONTINUE
  13  CONTINUE
      RETURN
      END
C
C=======================================================================
C
      REAL FUNCTION WASUM4(A,B,C,D,J)
C
C----------------------------------------------------------------------
      REAL A(3,3,3,3),B(3),C(3),D(3)
      INTEGER I,J,K,L
      WASUM4=0.
      DO 13, I=1,3
        DO 11, K=1,3
          DO 10, L=1,3
            WASUM4=WASUM4+A(I,J,K,L)*B(I)*C(K)*D(L)
  10      CONTINUE
  11    CONTINUE
  13  CONTINUE
      RETURN
      END
C
C=======================================================================
C