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 without interfaces.
C
C Version: 5.20
C Date: 1998, November 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)
C
C----------------------------------------------------------------------
      INTEGER LU1,LU2,LU3
      INTEGER MGREEN,MAXPTS,NQ
      INTEGER    NQPTS
      PARAMETER (NQPTS=21)
      REAL GREEN(MGREEN),RPTS(NQPTS,MAXPTS/NQPTS)
      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 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     MGREEN..Dimension of an output array GREEN.
C     MAXPTS..Maximum number of records in arrays RPTS and IPTS.
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
C.......................................................................
C Common block /POINTC/ to store the results of complete ray tracing:
      INCLUDE 'pointc.inc'
C.......................................................................
C
C     Auxiliary storage locations:
      INTEGER MPTS
      REAL PIGRA(2,2),PIGIA(2,2),PIGR(2,2),PIGI(2,2)
      REAL PI
      PARAMETER (PI=3.1415926)
      REAL GREENA(32)
      REAL TTCOR,FREQ,SIGRAY,GAMA,AUX0,AUX1,AUX2,COOR(3)
      INTEGER NPTS
      INTEGER I,I1,I2,J
      INTEGER NFFT,NF
      REAL DT,FMIN,FMAX,DF,OF
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 qS2 (faster) wave.
C        RPTS(15,I)... Eigenvalue of Christoffel matrix in I-th point
C                      corresponding to qS1 (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 qS2 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
C-----------------------------------------------------------------------
C
      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)
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,RPTS,IPTS,MPTS,NPTS)
      IF (NPTS.EQ.0) THEN
C       End of rays:
        NQ=0
        RETURN
      ENDIF
C
C     Computing the values of travel time corrections along the ray:
      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))
        ELSE
C         S wave:
          TTCOR=TTCOR+(1./SQRT((RPTS(14,I2)+RPTS(15,I2))*0.5)+
     *                 1./SQRT((RPTS(14,I1)+RPTS(15,I1))*0.5) )*0.5*
     *          (RPTS(2,I2)-RPTS(2,I1))
        ENDIF
   30 CONTINUE
C
      IF (IPTS(12,I2).GT.0.) THEN
C       P wave: !!! OPRAVIT PO ZARAZENI KONVERTOVANYCH VLN !!!
        Y(1)=TTCOR+YI(1)
        CALL AP21(GREEN)
        NQ=32
        GOTO 91
      ENDIF
C
C     S wave:
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.
        SIGRAY=0.
C       Computing the propagator matrix PiGe along the ray:
C       Loop along points on the ray:
        DO 40, I1=2,NPTS
          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)
C         Matrix for all the steps along the ray to current point:
          CALL WAMAT(PIGRA,PIGIA,PIGR,PIGI)
   40   CONTINUE
C
C       Computing the Green function:
        Y(1)=TTCOR+YI(1)
        G11=RPTS(16,1)
        G12=RPTS(17,1)
        G21=RPTS(18,1)
        G22=RPTS(19,1)
        PIR11=PIGR(1,1)*G11+PIGR(1,2)*G21
        PIR21=PIGR(2,1)*G11+PIGR(2,2)*G21
        PIR12=PIGR(1,1)*G12+PIGR(1,2)*G22
        PIR22=PIGR(2,1)*G12+PIGR(2,2)*G22
        PII11=PIGI(1,1)*G11+PIGI(1,2)*G21
        PII21=PIGI(2,1)*G11+PIGI(2,2)*G21
        PII12=PIGI(1,1)*G12+PIGI(1,2)*G22
        PII22=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*PIR11+G12*PIR21
        Y(29)=G11*PII11+G12*PII21
        Y(30)=G21*PIR11+G22*PIR21
        Y(31)=G21*PII11+G22*PII21
        Y(32)=G11*PIR12+G12*PIR22
        Y(33)=G11*PII12+G12*PII22
        Y(34)=G21*PIR12+G22*PIR22
        Y(35)=G21*PII12+G22*PII22
C
        CALL AP21(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
      RETURN
      END
C================================================================================
C
      SUBROUTINE WAPTS(LU1,LU2,LU3,RPTS,IPTS,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 NPTS,NQPTS,MPTS
      PARAMETER (NQPTS=21)
      REAL    RPTS(NQPTS,MPTS)
      INTEGER IPTS(NQPTS,MPTS)
C Input:
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 qS2 (faster) wave.
C        RPTS(15,I)... Eigenvalue of Christoffel matrix in I-th point
C                      corresponding to qS1 (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 qS2 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    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=17)
      REAL ROLDP(MQANT),RNEWP(MQANT,MNEWP)
      INTEGER KNEWP(MNEWP)
      INTEGER MAXDIV
      PARAMETER (MAXDIV=MNEWP-1)
      REAL AA(10,21),RHO,QQ(21)
      REAL EE(9),DER(9)
      REAL DEFI,DELTFI
c     PARAMETER (DEFI=0.01)
      PARAMETER (DEFI=0.03)
c     PARAMETER (DEFI=100.)
      REAL TTOLD,CROLD(3),DCROLD(3)
      REAL TTNEW,CRNEW(3),DCRNEW(3)
      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     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.
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     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     CROLD...Coordinates of the old point stored for interpolation.
C     DCROLD..Derivatives in the old point stored for interpolation.
C     TTNEW...Travel time in the new point stored for interpolation.
C     CRNEW...Coordinates of the new point stored for interpolation.
C     DCRNEW..Derivatives in the new point stored for interpolation.
C     NPTSE...Number of points along the ray, where all the quantities
C             have been checked.
C
C-----------------------------------------------------------------------
C     Reading the quantities computed by CRT for one ray:
      CALL WAREAD(LU1,LU2,LU3,RPTS,IPTS,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)
      IF (ABS(RPTS(15,1)-RPTS(14,1)).LT.0.000001) THEN
C       Isotropic case, no projection of eigenvectors:
        RPTS(16,1)=0.
        RPTS(17,1)=0.
      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)
      CROLD(1)   =RPTS(3,1)
      CROLD(2)   =RPTS(4,1)
      CROLD(3)   =RPTS(5,1)
      IF (IPTS(12,1).GT.0) THEN
C       P waves:
        DCROLD(1)=DER(1)
        DCROLD(2)=DER(2)
        DCROLD(3)=DER(3)
      ELSE
C       S waves:
        DCROLD(1)=(DER(4)+DER(7))*0.5
        DCROLD(2)=(DER(5)+DER(8))*0.5
        DCROLD(3)=(DER(6)+DER(9))*0.5
      ENDIF
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)
      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         of 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       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=RNEWP(1,1)
        CRNEW(1)=RNEWP(2,1)
        CRNEW(2)=RNEWP(3,1)
        CRNEW(3)=RNEWP(4,1)
        IF (IPTS(12,NPTSE+1).GT.0) THEN
C         P waves:
          DCRNEW(1)=DER(1)
          DCRNEW(2)=DER(2)
          DCRNEW(3)=DER(3)
        ELSE
C         S waves:
          DCRNEW(1)=(DER(4)+DER(7))*0.5
          DCRNEW(2)=(DER(5)+DER(8))*0.5
          DCRNEW(3)=(DER(6)+DER(9))*0.5
        ENDIF
C
        IF (IPTS(12,NPTSE).NE.IPTS(12,NPTSE+1)) THEN
C         The ray is crossing an interface:
C         WAN-04
          CALL ERROR('WAN-04: A ray is crossing an interface')
C         This part is not yet coded in this version.
        ENDIF
C
        IF (ABS(RNEWP(12,1)-RNEWP(13,1)).LT.0.000001) THEN
C         New point is isotropic, qS eigenvalues are the same,
C         no change in eigenvectors:
          RNEWP(14,1)=ROLDP(14)
          RNEWP(15,1)=ROLDP(15)
          RNEWP(16,1)=ROLDP(16)
          RNEWP(17,1)=ROLDP(17)
          DELTFI=0.
          GOTO 20
        ENDIF
        IF (ABS(ROLDP(12)-ROLDP(13)).LT.0.000001) THEN
C         Old point is isotropic (new point is anizotropic):
          IF ((ROLDP(14).EQ.0.).AND.(ROLDP(15).EQ.0.)) THEN
C           New point is the first anizotropic point on the ray,
C           angular change is not to be computed:
            DELTFI=0.
            GOTO 20
          ENDIF
        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))
        IF (DELTFI.LE.DEFI) THEN
cc      IF (DELTFI.LE.( 24/(2.*PI*(OF+NF*DF)*RPTS(2,NPTS)*
CC          (RNEWP(12,NNEWP)-RNEWP(13,NNEWP)+ROLDP(12)-ROLDP(13)) ) 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) 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 STEP, or to
C         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,CROLD,DCROLD,TTNEW,CRNEW,DCRNEW,
     *  RNEWP(1,NNEWP),RNEWP(2,NNEWP),DER)
C       Slowness vector:
        RNEWP(5,NNEWP)=(ROLDP(5)+RNEWP(5,NNEWP-1))*0.5
        RNEWP(6,NNEWP)=(ROLDP(6)+RNEWP(6,NNEWP-1))*0.5
        RNEWP(7,NNEWP)=(ROLDP(7)+RNEWP(7,NNEWP-1))*0.5
C       Polarization vector: (could be interpolated by HIVD2 with
C       derivatives obtained from FCT from file 'raycb.for')
        RNEWP(8,NNEWP) =(ROLDP(8) +RNEWP(8,NNEWP-1)) *0.5
        RNEWP(9,NNEWP) =(ROLDP(9) +RNEWP(9,NNEWP-1)) *0.5
        RNEWP(10,NNEWP)=(ROLDP(10)+RNEWP(10,NNEWP-1))*0.5
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       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.000001) 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)
          DELTFI=0.
        ELSE
          DELTFI=WACHAN(ROLDP(14),RNEWP(11,NNEWP))
        ENDIF
        IF (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
            RPTS(I2,I1)=RPTS(I2,I1-1)
  30      CONTINUE
  31    CONTINUE
C       Recording interpolated quantities:
        RPTS( 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
C-LK  RPTS(21,NPTSE)=
C-LK *  (RNEWP(12,NNEWP)-RNEWP(13,NNEWP)+ROLDP(12)-ROLDP(13))*0.25
C-LK *  *(RNEWP(1,NNEWP)-ROLDP(1))
      RPTS(21,NPTSE)=
     *  (1./SQRT(RNEWP(13,NNEWP))-1./SQRT(RNEWP(12,NNEWP))
     *  +1./SQRT(ROLDP(13))-1./SQRT(ROLDP(12)))*0.5
     *  *(RNEWP(1,NNEWP)-ROLDP(1))
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
        GOTO 5
      ENDIF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE WAREAD(LU1,LU2,LU3,RPTS,IPTS,MPTS,NPTS)
      INTEGER LU1,LU2,LU3
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 of common
C block PTS.
C The subroutine expects, that files LU and LUI are already opened
C and that at least one point was already read by calling AP00.
C
      INTEGER NPTS,NQPTS,MPTS
      PARAMETER (NQPTS=21)
      REAL    RPTS(NQPTS,MPTS)
      INTEGER IPTS(NQPTS,MPTS)
C Input:
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    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 None of the storage locations of the common block are altered.
C.......................................................................
C
C     Auxiliary storage locations:
C
C-----------------------------------------------------------------------
C
  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
C
  20  CONTINUE
C     New point:
      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         of 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
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
C         reference surface in file CRT,
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
      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
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     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.000001) 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.000001) 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.000001) 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 are reorganized in the same way.
C
      REAL G(4),H(7)
C     G(1:4)=G1,G2,H1,H2
C     H(1:7)=R1,R2,R3,U1,U2,V1,V2
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     R1,R2  ...  Real numbers associeted 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     R1,R2  ...  Real numbers associeted to U and V.
C
C.......................................................................
C
C     Auxiliary storage locations:
      REAL SP1,SP2,SP3,SP4,A1,A2,AUX
      REAL G1,G2,H1,H2,U1,U2,V1,V2,R1,R2,R3
C
C-----------------------------------------------------------------------
      G1=G(1)
      G2=G(2)
      H1=G(3)
      H2=G(4)
      R3=H(1)
      R2=H(2)
      R1=H(3)
      U1=H(4)
      U2=H(5)
      V1=H(6)
      V2=H(7)
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
      ELSEIF (AUX.EQ.SP3) THEN
        A1=U1
        A2=U2
        U1=V1
        U2=V2
        V1=A1
        V2=A2
        AUX=R1
        R1=R2
        R2=AUX
      ELSEIF (AUX.EQ.SP4) THEN
        A1=U1
        A2=U2
        U1=-V1
        U2=-V2
        V1=A1
        V2=A2
        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
      ENDIF
      WACHAN=ASIN(0.5*((G1+U1)*(H1-V1)+(G2+U2)*(H2-V2)))
      H(1)=R3
      H(2)=R2
      H(3)=R1
      H(4)=U1
      H(5)=U2
      H(6)=V1
      H(7)=V2
      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
      INCLUDE 'eigen.for'
C     eigen.for
      INCLUDE 'means.for'
C     means.for
C
C=======================================================================
C