C
C Subroutine file 'hder.for' to calculate the phase-space derivatives
C of the Hamiltonian in anisotropic media.
C
C Version: 7.10
C Date: 2014, June 18
C Coded by Ludek Klimes
C
C.......................................................................
C
C This file consists of the following external procedures:
C     HDER... Subroutine designed to calculate the phase-space
C             derivatives of the Hamiltonian in anisotropic media.
C     GAA...  Subroutine calculating the eigenvalues and eigenvectors
C             of the Christoffel matrix, with entries GABX, GABP, GAAXX,
C             GAAXP and GAAPP calculating the first-order and
C             second-order phase-space derivatives of the Christoffel
C             matrix, transformed into the given eigenvectors of the
C             Christoffel matrix.
C
C=======================================================================
C
      SUBROUTINE HDER(ICB,A,P1,P2,P3,HP,HS,H,H1,H2,H3,H4,H5,H6,
     *                H11,H12,H22,H13,H23,H33,H14,H24,H34,H44,
     *                H15,H25,H35,H45,H55,H16,H26,H36,H46,H56,H66,
     *                E,IERR)
      INTEGER ICB,IERR
      REAL A(10,21),P1,P2,P3,HP,HS,H,H1,H2,H3,H4,H5,H6
      REAL H11,H12,H22,H13,H23,H33,H14,H24,H34,H44
      REAL H15,H25,H35,H45,H55,H16,H26,H36,H46,H56,H66,E(9)
C
C Subroutine designed to calculate the phase-space derivatives of the
C Hamiltonian.
C
C Input:
C     ICB...  Wave type:
C             ICB.GT.0: P wave,
C             ICB.LT.0: S wave.
C     A...    Density-reduced elastic moduli and their first and second
C             spatial derivatives.
C             The order of the value, first and second partial
C             derivatives of each parameter Aij (second array index)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     P1,P2,P3... Slowness vector.
C
C Output:
C     HP,HS.. Twice the value of the homogeneous Hamiltonian of the
C             second degree corresponding to P and S waves,
C             respectively.  HP equals the P-wave eigenvalue of the
C             Christoffel matrix, definition of HS depends on the value
C             of the input SEP parameter DEGREE.
C     H...    Twice the value of the homogeneous Hamiltonian of the
C             second degree corresponding to the given wave (for P wave,
C             H equals the P-wave eigenvalue of the Christoffel matrix).
C             H should be equal to 1 if P1,P2,P3 is the slowness vector.
C     H1,H2,H3,H4,H5,H6... First partial phase-space derivatives of the
C             second degree Hamiltonian.  The derivatives are calculated
C             for the normalized slowness vector (divided by the square
C             root of the value of the Hamiltonian).
C     H11,H12,H22,H13,H23,H33,H14,H24,H34,H44,H15,H25,H35,H45,H55,H16,
C     H26,H36,H46,H56,H66... Second partial phase-space derivatives of
C             the second degree Hamiltonian. The derivatives are
C             calculated for the normalized slowness vector (divided by
C             the square root of the value of the Hamiltonian).
C     E...    Components of the eigenvectors of the Christoffel matrix:
C             E(1:3),E(4:6): S-wave eigenvectors.
C             E(7:9): P-wave eigenvector.
C     IERR... Error indicator when tracing the anisotropic-ray-theory
C             S-wave rays.
C             If KSWAVE=1 or KSWAVE=2, and the relative difference
C               between the S-wave eigenvalues of the Christoffel matrix
C               is smaller than input SEP parameter DSWAVE, IERR=1.
C               In this case, H1,H2,...,H6 are considerably inaccurate,
C               and H11,H12,...,H66 are defined but meaningless.
C             Otherwise, IERR=0.
C                                                     
C Input SEP parameters:
C     DEGREE='real'... Degree of the homogeneous Hamiltonian to be
C             arithmetically averaged for common S-wave ray tracing.
C     KSWAVE=integer... Selection of the anisotropic-ray-theory S-wave
C             rays instead of common S-wave rays.
C             Affects only S waves (ICB.LT.0).
C             KSWAVE=1:  Slower anisotropic-ray-theory S-wave rays.
C             KSWAVE=2:  Faster anisotropic-ray-theory S-wave rays.
C             Else (default):  Common S-wave rays.
C             Default: KSWAVE=0
C     DSWAVE='real'... Minimum relative difference between the S-wave
C             eigenvalues G1 and G2 of the Christoffel matrix. If the
C             difference is smaller, output error indicator IERR is
C             set to 1 in order to have a possibility to terminate ray
C             tracing.
C             Used only if ICB.LT.0 and (KSWAVE=1 or KSWAVE=2).
C             The relative difference is defined as ABS(G1-G2)/GS,
C             where GS=(G1+G2)/2.
C             The maximum angular numerical error of the eigenvectors
C             of the Christoffel matrix in radians is roughly equal to
C             (relative rounding error)/DSWAVE.
C             Note that for DSWAVE smaller than its default value
C             the KMAH index and the matrix of geometrical spreading
C             may be considerably erroneous.
C             Default: DSWAVE=0.00001
C
C Subroutines and external functions required:
      EXTERNAL GAA,GABX,GABP,GAAXX,GAAXP,GAAPP,ERROR,RSEP3I,RSEP3R
C     Subroutine GAA and its entries GABX,GABP,GAAXX,GAAXP,GAAPP...
C             This file.
C     ERROR.. File 'error.for'.
C     RSEP3I,RSEP3R... File 'sep.for'.
C     EIGEN.. File 'eigen.for'.
C
C Date: 2014, June 18
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Temporary storage locations:
      LOGICAL LANIRT
      REAL HSN,HS2,E11,E21,E31,E12,E22,E32,E13,E23,E33,PP1,PP2,PP3
      REAL G11,G22,G33
      REAL G111,G121,G221,G131,G231,G331
      REAL G112,G122,G222,G132,G232,G332
      REAL G113,G123,G223,G133,G233,G333
      REAL G114,G124,G224,G134,G234,G334
      REAL G115,G125,G225,G135,G235,G335
      REAL G116,G126,G226,G136,G236,G336
      REAL G1111,G2211,G3311,G1112,G2212,G3312,G1122,G2222,G3322
      REAL G1113,G2213,G3313,G1123,G2223,G3323,G1133,G2233,G3333
      REAL G1114,G2214,G3314,G1124,G2224,G3324,G1134,G2234,G3334
      REAL G1144,G2244,G3344,G1115,G2215,G3315,G1125,G2225,G3325
      REAL G1135,G2235,G3335,G1145,G2245,G3345,G1155,G2255,G3355
      REAL G1116,G2216,G3316,G1126,G2226,G3326,G1136,G2236,G3336
      REAL G1146,G2246,G3346,G1156,G2256,G3356,G1166,G2266,G3366
      REAL W1,W2,W3,W11,W12,W22,W13,W23,WH,AUX
C
C     Degree of the homogeneous Hamiltonian:
      INTEGER KSWAVE
      REAL DEGREE,DEGRE2,DSWAVE
      SAVE DEGREE,DEGRE2,DSWAVE,KSWAVE
      DATA DEGREE/0./
C
C.......................................................................
C
      IF(DEGREE.EQ.0.) THEN
        CALL RSEP3R('DEGREE',DEGREE,-1.)
        DEGRE2=0.5*DEGREE
        CALL RSEP3I('KSWAVE',KSWAVE,0)
        CALL RSEP3R('DSWAVE',DSWAVE,0.00001)
        IF(KSWAVE.GE.1.AND.KSWAVE.LE.2.AND.DSWAVE.LE.0.0) THEN
C         5A3
          CALL ERROR('5A3 in HDER: SEP parameter DSWAVE not positive')
        END IF
      END IF
      IERR=0
C
C     Eigenvalues and eigenvectors of the Christoffel matrix:
      CALL GAA(A,P1,P2,P3,G11,G22,G33,
     *         E11,E21,E31,E12,E22,E32,E13,E23,E33)
      E(1)=E11
      E(2)=E21
      E(3)=E31
      E(4)=E12
      E(5)=E22
      E(6)=E32
      E(7)=E13
      E(8)=E23
      E(9)=E33
C
C     Twice the homogeneous Hamiltonian of the second degree:
      HP=G33
      IF(KSWAVE.EQ.1.AND.ICB.LT.0) THEN
C       Slower anisotropic-ray-theory S-wave ray
C       Setting anisotropic-ray-theory ray tracing
        LANIRT=.TRUE.
        IF(ABS(G11-G22).LE.DSWAVE*0.5*(G11+G22)) THEN
          IERR=1
        END IF
C       Exchanging P wave and slower S wave
        HS=G11
        G11=G33
        G33=HS
        AUX=E11
        E11=E13
        E13=AUX
        AUX=E21
        E21=E23
        E23=AUX
        AUX=E31
        E31=E33
        E33=AUX
      ELSE IF(KSWAVE.EQ.2.AND.ICB.LT.0) THEN
C       Faster anisotropic-ray-theory S-wave ray
C       Setting anisotropic-ray-theory ray tracing
        LANIRT=.TRUE.
        IF(ABS(G11-G22).LE.DSWAVE*0.5*(G11+G22)) THEN
          IERR=1
        END IF
C       Exchanging P wave and faster S wave
        HS=G22
        G22=G33
        G33=HS
        AUX=E12
        E12=E13
        E13=AUX
        AUX=E22
        E22=E23
        E23=AUX
        AUX=E32
        E32=E33
        E33=AUX
      ELSE
C       Setting S-wave common ray tracing
        LANIRT=.FALSE.
C       Averaging S-wave eigenvalues of the Christoffel matrix
        HSN=0.5*(G11**DEGRE2+G22**DEGRE2)
        HS2=HSN**(1./DEGRE2)
        HS=HS2
      END IF
      IF(ICB.GT.0) THEN
C       P wave:
C       Setting anisotropic-ray-theory ray tracing
        LANIRT=.TRUE.
        H=G33
      ELSE IF(ICB.LT.0) THEN
C       S wave:
        H=HS
      ELSE
C       5A2
        CALL ERROR('5A2 in HDER: Zero wave type')
C       This error should not appear, contact the authors.
      END IF
C
C     Normalized slowness vector:
      PP1=P1/SQRT(H)
      PP2=P2/SQRT(H)
      PP3=P3/SQRT(H)
C
C     Phase-space derivatives of the Christoffel matrix transformed into
C     its own eigenvectors required both for P and S waves:
      CALL GABX(A,PP1,PP2,PP3,E11,E21,E31,E13,E23,E33,G131,G132,G133)
      CALL GABP(H,            E11,E21,E31,E13,E23,E33,G134,G135,G136)
      CALL GABX(A,PP1,PP2,PP3,E12,E22,E32,E13,E23,E33,G231,G232,G233)
      CALL GABP(H,            E12,E22,E32,E13,E23,E33,G234,G235,G236)
C
      IF(LANIRT) THEN
C       P wave or anisotropic-ray-theory S wave:
C
C       Phase-space derivatives of the Christoffel matrix transformed
C       into its own eigenvectors:
        CALL GAAXX(A,PP1,PP2,PP3,E13,E23,E33,
     *             G3311,G3312,G3322,G3313,G3323,G3333)
        CALL GAAXP(A,PP1,PP2,PP3,E13,E23,E33,
     *            G3314,G3324,G3334,G3315,G3325,G3335,G3316,G3326,G3336)
        CALL GAAPP(              E13,E23,E33,
     *             G3344,G3345,G3355,G3346,G3356,G3366)
        G331=0.5*(G3314*PP1+G3315*PP2+G3316*PP3)
        G332=0.5*(G3324*PP1+G3325*PP2+G3326*PP3)
        G333=0.5*(G3334*PP1+G3335*PP2+G3336*PP3)
        G334=     G3344*PP1+G3345*PP2+G3346*PP3
        G335=     G3345*PP1+G3355*PP2+G3356*PP3
        G336=     G3346*PP1+G3356*PP2+G3366*PP3
C
        W3=0.5
        H1=W3*G331
        H2=W3*G332
        H3=W3*G333
        H4=W3*G334
        H5=W3*G335
        H6=W3*G336
C
        IF(IERR.EQ.0) THEN
          W13=G33/(G33-G11)
          W23=G33/(G33-G22)
        ELSE
          W13=0.0
          W23=0.0
        END IF
        H11=W3*G3311+W13*G131*G131+W23*G231*G231
        H12=W3*G3312+W13*G131*G132+W23*G231*G232
        H22=W3*G3322+W13*G132*G132+W23*G232*G232
        H13=W3*G3313+W13*G131*G133+W23*G231*G233
        H23=W3*G3323+W13*G132*G133+W23*G232*G233
        H33=W3*G3333+W13*G133*G133+W23*G233*G233
        H14=W3*G3314+W13*G131*G134+W23*G231*G234
        H24=W3*G3324+W13*G132*G134+W23*G232*G234
        H34=W3*G3334+W13*G133*G134+W23*G233*G234
        H44=W3*G3344+W13*G134*G134+W23*G234*G234
        H15=W3*G3315+W13*G131*G135+W23*G231*G235
        H25=W3*G3325+W13*G132*G135+W23*G232*G235
        H35=W3*G3335+W13*G133*G135+W23*G233*G235
        H45=W3*G3345+W13*G134*G135+W23*G234*G235
        H55=W3*G3355+W13*G135*G135+W23*G235*G235
        H16=W3*G3316+W13*G131*G136+W23*G231*G236
        H26=W3*G3326+W13*G132*G136+W23*G232*G236
        H36=W3*G3336+W13*G133*G136+W23*G233*G236
        H46=W3*G3346+W13*G134*G136+W23*G234*G236
        H56=W3*G3356+W13*G135*G136+W23*G235*G236
        H66=W3*G3366+W13*G136*G136+W23*G236*G236
      ELSE
C       S wave:
C
C       Phase-space derivatives of the Christoffel matrix transformed
C       into its own eigenvectors:
        CALL GABX(A,PP1,PP2,PP3,E11,E21,E31,E12,E22,E32,G121,G122,G123)
        CALL GABP(H,            E11,E21,E31,E12,E22,E32,G124,G125,G126)
        CALL GAAXX(A,PP1,PP2,PP3,E11,E21,E31,
     *             G1111,G1112,G1122,G1113,G1123,G1133)
        CALL GAAXP(A,PP1,PP2,PP3,E11,E21,E31,
     *            G1114,G1124,G1134,G1115,G1125,G1135,G1116,G1126,G1136)
        CALL GAAPP(           E11,E21,E31,
     *             G1144,G1145,G1155,G1146,G1156,G1166)
        G111=0.5*(G1114*PP1+G1115*PP2+G1116*PP3)
        G112=0.5*(G1124*PP1+G1125*PP2+G1126*PP3)
        G113=0.5*(G1134*PP1+G1135*PP2+G1136*PP3)
        G114=     G1144*PP1+G1145*PP2+G1146*PP3
        G115=     G1145*PP1+G1155*PP2+G1156*PP3
        G116=     G1146*PP1+G1156*PP2+G1166*PP3
        CALL GAAXX(A,PP1,PP2,PP3,E12,E22,E32,
     *            G2211,G2212,G2222,G2213,G2223,G2233)
        CALL GAAXP(A,PP1,PP2,PP3,E12,E22,E32,
     *            G2214,G2224,G2234,G2215,G2225,G2235,G2216,G2226,G2236)
        CALL GAAPP(           E12,E22,E32,
     *             G2244,G2245,G2255,G2246,G2256,G2266)
        G221=0.5*(G2214*PP1+G2215*PP2+G2216*PP3)
        G222=0.5*(G2224*PP1+G2225*PP2+G2226*PP3)
        G223=0.5*(G2234*PP1+G2235*PP2+G2236*PP3)
        G224=     G2244*PP1+G2245*PP2+G2246*PP3
        G225=     G2245*PP1+G2255*PP2+G2256*PP3
        G226=     G2246*PP1+G2256*PP2+G2266*PP3
C
        W1=0.25*G11**(DEGRE2-1.)*HS2/HSN
        W2=0.25*G22**(DEGRE2-1.)*HS2/HSN
        H1=W1*G111+W2*G221
        H2=W1*G112+W2*G222
        H3=W1*G113+W2*G223
        H4=W1*G114+W2*G224
        H5=W1*G115+W2*G225
        H6=W1*G116+W2*G226
C
        W13=2.*W1*HS2/(G11-G33)
        W23=2.*W2*HS2/(G22-G33)
        W11=(DEGRE2-1.)*W1*HS2/G11
        W22=(DEGRE2-1.)*W2*HS2/G22
        IF(DEGREE.EQ.-2) THEN
          W12=-(G11+G22)/(G11*G22)**2
        ELSE IF(DEGREE.EQ.-1) THEN
          W12=-(G11+SQRT(G11*G22)+G22)/(G11*SQRT(G22)+G22*SQRT(G11))
     *        /(G11*G22)
        ELSE IF(DEGREE.EQ.1) THEN
          W12=-1./(G11*SQRT(G22)+G22*SQRT(G11))
        ELSE IF(DEGREE.EQ.2) THEN
          W12=0.
        ELSE
C         5A1
          CALL ERROR('5A1 in HDER: Wrong degree of Hamiltonian')
C         Only homogeneous Hamiltonians of degrees -2, -1, 1 and 2 are
C         now allowed.
        END IF
        W12=0.5*W12*HS2*HS2/HSN
        WH=(2.-DEGREE)
        H11=W1*G1111+W11*G111*G111+W13*G131*G131+W12*G121*G121
        H12=W1*G1112+W11*G111*G112+W13*G131*G132+W12*G121*G122
        H22=W1*G1122+W11*G112*G112+W13*G132*G132+W12*G122*G122
        H13=W1*G1113+W11*G111*G113+W13*G131*G133+W12*G121*G123
        H23=W1*G1123+W11*G112*G113+W13*G132*G133+W12*G122*G123
        H33=W1*G1133+W11*G113*G113+W13*G133*G133+W12*G123*G123
        H14=W1*G1114+W11*G111*G114+W13*G131*G134+W12*G121*G124
        H24=W1*G1124+W11*G112*G114+W13*G132*G134+W12*G122*G124
        H34=W1*G1134+W11*G113*G114+W13*G133*G134+W12*G123*G124
        H44=W1*G1144+W11*G114*G114+W13*G134*G134+W12*G124*G124
        H15=W1*G1115+W11*G111*G115+W13*G131*G135+W12*G121*G125
        H25=W1*G1125+W11*G112*G115+W13*G132*G135+W12*G122*G125
        H35=W1*G1135+W11*G113*G115+W13*G133*G135+W12*G123*G125
        H45=W1*G1145+W11*G114*G115+W13*G134*G135+W12*G124*G125
        H55=W1*G1155+W11*G115*G115+W13*G135*G135+W12*G125*G125
        H16=W1*G1116+W11*G111*G116+W13*G131*G136+W12*G121*G126
        H26=W1*G1126+W11*G112*G116+W13*G132*G136+W12*G122*G126
        H36=W1*G1136+W11*G113*G116+W13*G133*G136+W12*G123*G126
        H46=W1*G1146+W11*G114*G116+W13*G134*G136+W12*G124*G126
        H56=W1*G1156+W11*G115*G116+W13*G135*G136+W12*G125*G126
        H66=W1*G1166+W11*G116*G116+W13*G136*G136+W12*G126*G126
        H11=W2*G2211+W22*G221*G221+W23*G231*G231+WH*H1*H1+H11
        H12=W2*G2212+W22*G221*G222+W23*G231*G232+WH*H1*H2+H12
        H22=W2*G2222+W22*G222*G222+W23*G232*G232+WH*H2*H2+H22
        H13=W2*G2213+W22*G221*G223+W23*G231*G233+WH*H1*H3+H13
        H23=W2*G2223+W22*G222*G223+W23*G232*G233+WH*H2*H3+H23
        H33=W2*G2233+W22*G223*G223+W23*G233*G233+WH*H3*H3+H33
        H14=W2*G2214+W22*G221*G224+W23*G231*G234+WH*H1*H4+H14
        H24=W2*G2224+W22*G222*G224+W23*G232*G234+WH*H2*H4+H24
        H34=W2*G2234+W22*G223*G224+W23*G233*G234+WH*H3*H4+H34
        H44=W2*G2244+W22*G224*G224+W23*G234*G234+WH*H4*H4+H44
        H15=W2*G2215+W22*G221*G225+W23*G231*G235+WH*H1*H5+H15
        H25=W2*G2225+W22*G222*G225+W23*G232*G235+WH*H2*H5+H25
        H35=W2*G2235+W22*G223*G225+W23*G233*G235+WH*H3*H5+H35
        H45=W2*G2245+W22*G224*G225+W23*G234*G235+WH*H4*H5+H45
        H55=W2*G2255+W22*G225*G225+W23*G235*G235+WH*H5*H5+H55
        H16=W2*G2216+W22*G221*G226+W23*G231*G236+WH*H1*H6+H16
        H26=W2*G2226+W22*G222*G226+W23*G232*G236+WH*H2*H6+H26
        H36=W2*G2236+W22*G223*G226+W23*G233*G236+WH*H3*H6+H36
        H46=W2*G2246+W22*G224*G226+W23*G234*G236+WH*H4*H6+H46
        H56=W2*G2256+W22*G225*G226+W23*G235*G236+WH*H5*H6+H56
        H66=W2*G2266+W22*G226*G226+W23*G236*G236+WH*H6*H6+H66
      END IF
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE GAA(A,P1,P2,P3,
     *               G11,G22,G33,E11,E21,E31,E12,E22,E32,E13,E23,E33)
C     Subroutine calculating the eigenvalues and eigenvectors of the
C     Christoffel matrix.
C
C     ENTRY GABX(A,P1,P2,P3,E1A,E2A,E3A,E1B,E2B,E3B,GAB1,GAB2,GAB3)
C     Subroutine calculating the first spatial derivatives of the
C     Christoffel matrix.
C
C     ENTRY GABP(H,E1A,E2A,E3A,E1B,E2B,E3B,GAB4,GAB5,GAB6)
C     Subroutine calculating the first partial derivatives of the
C     Christoffel matrix with respect to the slowness vector.
C
C     ENTRY GAAXX(A,P1,P2,P3,E1A,E2A,E3A,
C    *            GAA11,GAA12,GAA22,GAA13,GAA23,GAA33)
C     Subroutine calculating the second spatial derivatives of the
C     Christoffel matrix.
C
C     ENTRY GAAXP(A,P1,P2,P3,E1A,E2A,E3A,
C    *            GAA14,GAA24,GAA34,GAA15,GAA25,GAA35,GAA16,GAA26,GAA36)
C     Subroutine calculating the second mixed derivatives of the
C     Christoffel matrix.
C
C     ENTRY GAAPP(E1A,E2A,E3A,GAA44,GAA45,GAA55,GAA46,GAA56,GAA66)
C     Subroutine calculating the second slowness derivatives of the
C     Christoffel matrix.
C
      REAL A(10,21),H,P1,P2,P3,E1A,E2A,E3A,E1B,E2B,E3B
      REAL G11,G22,G33,E11,E21,E31,E12,E22,E32,E13,E23,E33
      REAL GAB1,GAB2,GAB3,GAB4,GAB5,GAB6
      REAL GAA11,GAA12,GAA22,GAA13,GAA23,GAA33
      REAL GAA14,GAA24,GAA34,GAA15,GAA25,GAA35,GAA16,GAA26,GAA36
      REAL GAA44,GAA45,GAA55,GAA46,GAA56,GAA66
C
C Input:
C     A...    Elastic moduli and their first and second spatial
C             derivatives.
C     H...    Value of the homogeneous Hamiltonian of the second degree
C             corresponding to the slowness vector before normalization
C             (division by the square root of the value of the
C             Hamiltonian).  H is used to rescale the quantities
C             calculated previously by subroutine GAA and used by entry
C             GAAP.  If there is no rescaling of the slowness vector
C             between invocations of subroutines GAA and GAAP, set H=1.
C     P1,P2,P3... Slowness vector.
C     E1A,E2A,E3A... One of the eigenvectors of the Christoffel matrix,
C             previously calculated by subroutine GAA.
C     E1B,E2B,E3B... Another eigenvector of the Christoffel matrix.
C
C Output:
C     G11,G22,G33... Eigenvalues of the Christoffel matrix.
C             Eigenvalue G33 corresponds to the P wave.
C     E11,E21,E31,E12,E22,E32,E13,E23,E33... Eigenvectors of the
C             Christoffel matrix.
C     GAB1,GAB2,GAB3,GAB4,GAB5,GAB6... First partial phase-space
C             derivatives of the Christoffel matrix, multiplied by
C             the given eigenvectors EA and EB.
C     GAA11,GAA12,GAA22,GAA13,GAA23,GAA33,GAA14,GAA24,GAA34,GAA15,GAA25,
C     GAA35,GAA16,GAA26,GAA36,GAA44,GAA45,GAA55,GAA46,GAA56,GAA66...
C             Second partial phase-space derivatives of derivatives of
C             the Christoffel matrix, twice multiplied by the given
C             eigenvector EA.
C
C Subroutines and external functions required:
      EXTERNAL EIGEN
C     EIGEN.. File 'eigen.for'.
C
C Date: 2003, May 5
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Quantities to be saved:
C
C     Elastic moduli A(i,j,k,l):
      REAL AIJKL(21)
      SAVE AIJKL
      REAL A11,A12,A22,A13,A23,A33,A14,A24,A34,A44
      REAL A15,A25,A35,A45,A55,A16,A26,A36,A46,A56,A66
      REAL A1111
      REAL A1122,A2211
      REAL A2222
      REAL A1133,A3311
      REAL A2233,A3322
      REAL A3333
      REAL A1123,A1132,A2311,A3211
      REAL A2223,A2232,A2322,A3222
      REAL A3323,A3332,A2333,A3233
      REAL A2323,A3223,A2332,A3232
      REAL A1113,A1131,A1311,A3111
      REAL A2213,A2231,A1322,A3122
      REAL A3313,A3331,A1333,A3133
      REAL A2313,A3213,A2331,A3231,A1323,A3123,A1332,A3132
      REAL A1313,A3113,A1331,A3131
      REAL A1112,A1121,A1211,A2111
      REAL A2212,A2221,A1222,A2122
      REAL A3312,A3321,A1233,A2133
      REAL A2312,A3212,A2321,A3221,A1223,A2123,A1232,A2132
      REAL A1312,A3112,A1321,A3121,A1213,A2113,A1231,A2131
      REAL A1212,A2112,A1221,A2121
      EQUIVALENCE (AIJKL( 1),A11,A1111)
      EQUIVALENCE (AIJKL( 2),A12,A1122,A2211)
      EQUIVALENCE (AIJKL( 3),A22,A2222)
      EQUIVALENCE (AIJKL( 4),A13,A1133,A3311)
      EQUIVALENCE (AIJKL( 5),A23,A2233,A3322)
      EQUIVALENCE (AIJKL( 6),A33,A3333)
      EQUIVALENCE (AIJKL( 7),A14,A1123,A1132,A2311,A3211)
      EQUIVALENCE (AIJKL( 8),A24,A2223,A2232,A2322,A3222)
      EQUIVALENCE (AIJKL( 9),A34,A3323,A3332,A2333,A3233)
      EQUIVALENCE (AIJKL(10),A44,A2323,A3223,A2332,A3232)
      EQUIVALENCE (AIJKL(11),A15,A1113,A1131,A1311,A3111)
      EQUIVALENCE (AIJKL(12),A25,A2213,A2231,A1322,A3122)
      EQUIVALENCE (AIJKL(13),A35,A3313,A3331,A1333,A3133)
      EQUIVALENCE (AIJKL(14),A45,A2313,A3213,A2331,A3231,
     *                           A1323,A1332,A3123,A3132)
      EQUIVALENCE (AIJKL(15),A55,A1313,A3113,A1331,A3131)
      EQUIVALENCE (AIJKL(16),A16,A1112,A1121,A1211,A2111)
      EQUIVALENCE (AIJKL(17),A26,A2212,A2221,A1222,A2122)
      EQUIVALENCE (AIJKL(18),A36,A3312,A3321,A1233,A2133)
      EQUIVALENCE (AIJKL(19),A46,A2312,A3212,A2321,A3221,
     *                           A1223,A1232,A2123,A2132)
      EQUIVALENCE (AIJKL(20),A56,A1312,A3112,A1321,A3121,
     *                           A1213,A1231,A2113,A2131)
      EQUIVALENCE (AIJKL(21),A66,A1212,A2112,A1221,A2121)
C
C     A(i,j,k)=A(i,j,k,l)*P(l):
      REAL A111,A211,A311,A121,A221,A321,A131,A231,A331
      REAL A112,A212,A312,A122,A222,A322,A132,A232,A332
      REAL A113,A213,A313,A123,A223,A323,A133,A233,A333
      SAVE A111,A121,A221,A131,A231,A331
      SAVE A112,A122,A222,A132,A232,A332
      SAVE A113,A123,A223,A133,A233,A333
      EQUIVALENCE (A121,A211),(A131,A311),(A231,A321)
      EQUIVALENCE (A122,A212),(A132,A312),(A232,A322)
      EQUIVALENCE (A123,A213),(A133,A313),(A233,A323)
C
C-----------------------------------------------------------------------
C
C     Temporary storage locations:
C
      INTEGER I
      REAL G(10),E(9),E4A,E5A,E6A,W1,W2,W3,W4,W5,W6
      REAL W1A,W2A,W3A,W4A,W5A,W6A,W1B,W2B,W3B,W4B,W5B,W6B
      REAL W11,W12,W22,W13,W23,W33,W14,W24,W34,W44
      REAL W15,W25,W35,W45,W55,W16,W26,W36,W46,W56,W66
C
C     I...    Loop variable.
C     G...    Christoffel matrix, its eigenvalues, or their derivatives.
C     E...    Eigenvectors of the Christoffel matrix.
C     EA4,EA5,EA6... Copy of a component of the given eigenvector to
C             clearly arrange the Voigt notation in some equations.
C     W**...  Weighting factors.  Indices correspond to Voigt notation.
C
C-----------------------------------------------------------------------
C
C     SUBROUTINE GAA(A,P1,P2,P3,G11,G22,G33,
C    *               E11,E21,E31,E12,E22,E32,E13,E23,E33)
C
C     Storing the elastic moduli:
      DO 10 I=1,21
        AIJKL(I)=A(1,I)
   10 CONTINUE
C
C     Christoffel matrix:
      A111=A1111*P1+A1112*P2+A1113*P3
      A121=A1211*P1+A1212*P2+A1213*P3
      A221=A2211*P1+A2212*P2+A2213*P3
      A131=A1311*P1+A1312*P2+A1313*P3
      A231=A2311*P1+A2312*P2+A2313*P3
      A331=A3311*P1+A3312*P2+A3313*P3
      A112=A1121*P1+A1122*P2+A1123*P3
      A122=A1221*P1+A1222*P2+A1223*P3
      A222=A2221*P1+A2222*P2+A2223*P3
      A132=A1321*P1+A1322*P2+A1323*P3
      A232=A2321*P1+A2322*P2+A2323*P3
      A332=A3321*P1+A3322*P2+A3323*P3
      A113=A1131*P1+A1132*P2+A1133*P3
      A123=A1231*P1+A1232*P2+A1233*P3
      A223=A2231*P1+A2232*P2+A2233*P3
      A133=A1331*P1+A1332*P2+A1333*P3
      A233=A2331*P1+A2332*P2+A2333*P3
      A333=A3331*P1+A3332*P2+A3333*P3
      G(1)=P1*A111+P2*A211+P3*A311
      G(2)=P1*A112+P2*A212+P3*A312
      G(3)=P1*A122+P2*A222+P3*A322
      G(4)=P1*A113+P2*A213+P3*A313
      G(5)=P1*A123+P2*A223+P3*A323
      G(6)=P1*A133+P2*A233+P3*A333
C
C     Eigenvalues and eigenvectors of the Christoffel matrix:
      CALL EIGEN(G,E,3,0)
      G33=G(1)
      G22=G(3)
      G11=G(6)
      E13=E(1)
      E23=E(2)
      E33=E(3)
      E12=E(4)
      E22=E(5)
      E32=E(6)
      E11=E(7)
      E21=E(8)
      E31=E(9)
      RETURN
C
C-----------------------------------------------------------------------
C
      ENTRY GABX(A,P1,P2,P3,E1A,E2A,E3A,E1B,E2B,E3B,GAB1,GAB2,GAB3)
C
C     Subroutine calculating the first spatial derivatives of the
C     Christoffel matrix.
C
      W1A=E1A*P1
      W2A=E2A*P2
      W3A=E3A*P3
      W4A=E2A*P3+E3A*P2
      W5A=E1A*P3+E3A*P1
      W6A=E1A*P2+E2A*P1
      W1B=E1B*P1
      W2B=E2B*P2
      W3B=E3B*P3
      W4B=E2B*P3+E3B*P2
      W5B=E1B*P3+E3B*P1
      W6B=E1B*P2+E2B*P1
      W11=W1A*W1B
      W12=W1A*W2B+W1B*W2A
      W22=W2A*W2B
      W13=W1A*W3B+W1B*W3A
      W23=W2A*W3B+W2B*W3A
      W33=W3A*W3B
      W14=W1A*W4B+W1B*W4A
      W24=W2A*W4B+W2B*W4A
      W34=W3A*W4B+W3B*W4A
      W44=W4A*W4B
      W15=W1A*W5B+W1B*W5A
      W25=W2A*W5B+W2B*W5A
      W35=W3A*W5B+W3B*W5A
      W45=W4A*W5B+W4B*W5A
      W55=W5A*W5B
      W16=W1A*W6B+W1B*W6A
      W26=W2A*W6B+W2B*W6A
      W36=W3A*W6B+W3B*W6A
      W46=W4A*W6B+W4B*W6A
      W56=W5A*W6B+W5B*W6A
      W66=W6A*W6B
      DO 20 I=2,4
        G(I)=W11*A(I, 1)+W12*A(I, 2)+W22*A(I, 3)+W13*A(I, 4)+W23*A(I, 5)
     *      +W33*A(I, 6)+W14*A(I, 7)+W24*A(I, 8)+W34*A(I, 9)+W44*A(I,10)
     *      +W15*A(I,11)+W25*A(I,12)+W35*A(I,13)+W45*A(I,14)+W55*A(I,15)
     *      +W16*A(I,16)+W26*A(I,17)+W36*A(I,18)+W46*A(I,19)+W56*A(I,20)
     *      +W66*A(I,21)
   20 CONTINUE
      GAB1=G(2)
      GAB2=G(3)
      GAB3=G(4)
      RETURN
C
C-----------------------------------------------------------------------
C
      ENTRY GABP(H,E1A,E2A,E3A,E1B,E2B,E3B,GAB4,GAB5,GAB6)
C
C     Subroutine calculating the first partial derivatives of the
C     Christoffel matrix with respect to the slowness vector.
C
      W1=SQRT(H)
      W11=E1A*E1B*2.
      W12=E1A*E2B+E2A*E1B
      W22=E2A*E2B*2.
      W13=E1A*E3B+E3A*E1B
      W23=E2A*E3B+E3A*E2B
      W33=E3A*E3B*2.
C     GAB(i+3)=[A(i,j,k,l)+A(i,k,j,l)]*EA(j)*EB(k)*P(l):
      GAB4=A111*W11+(A112+A121)*W12+(A113+A131)*W13
     *                   +A122 *W22+(A123+A132)*W23+A133*W33
      GAB5=A211*W11+(A212+A221)*W12+(A213+A231)*W13
     *                   +A222 *W22+(A223+A232)*W23+A233*W33
      GAB6=A311*W11+(A312+A321)*W12+(A313+A331)*W13
     *                   +A322 *W22+(A323+A332)*W23+A333*W33
      GAB4=GAB4/W1
      GAB5=GAB5/W1
      GAB6=GAB6/W1
      RETURN
C
C-----------------------------------------------------------------------
C
      ENTRY GAAXX(A,P1,P2,P3,E1A,E2A,E3A,
     *            GAA11,GAA12,GAA22,GAA13,GAA23,GAA33)
C
C     Subroutine calculating the spatial derivatives of the Christoffel
C     matrix.
C
      W1=E1A*P1
      W2=E2A*P2
      W3=E3A*P3
      W4=E2A*P3+E3A*P2
      W5=E1A*P3+E3A*P1
      W6=E1A*P2+E2A*P1
      W11=W1*W1
      W12=W1*W2*2.
      W22=W2*W2
      W13=W1*W3*2.
      W23=W2*W3*2.
      W33=W3*W3
      W14=W1*W4*2.
      W24=W2*W4*2.
      W34=W3*W4*2.
      W44=W4*W4
      W15=W1*W5*2.
      W25=W2*W5*2.
      W35=W3*W5*2.
      W45=W4*W5*2.
      W55=W5*W5
      W16=W1*W6*2.
      W26=W2*W6*2.
      W36=W3*W6*2.
      W46=W4*W6*2.
      W56=W5*W6*2.
      W66=W6*W6
      DO 40 I=5,10
        G(I)=W11*A(I, 1)+W12*A(I, 2)+W22*A(I, 3)+W13*A(I, 4)+W23*A(I, 5)
     *      +W33*A(I, 6)+W14*A(I, 7)+W24*A(I, 8)+W34*A(I, 9)+W44*A(I,10)
     *      +W15*A(I,11)+W25*A(I,12)+W35*A(I,13)+W45*A(I,14)+W55*A(I,15)
     *      +W16*A(I,16)+W26*A(I,17)+W36*A(I,18)+W46*A(I,19)+W56*A(I,20)
     *      +W66*A(I,21)
   40 CONTINUE
      GAA11=G(5)
      GAA12=G(6)
      GAA13=G(7)
      GAA22=G(8)
      GAA23=G(9)
      GAA33=G(10)
      RETURN
C
C-----------------------------------------------------------------------
C
      ENTRY GAAXP(A,P1,P2,P3,E1A,E2A,E3A,
     *            GAA14,GAA24,GAA34,GAA15,GAA25,GAA35,GAA16,GAA26,GAA36)
C
C     Subroutine calculating the mixed derivatives of the Christoffel
C     matrix.
C
      W1=E1A*P1
      W2=E2A*P2
      W3=E3A*P3
      W4=E2A*P3+E3A*P2
      W5=E1A*P3+E3A*P1
      W6=E1A*P2+E2A*P1
C     Derivative with respect to P1 (nonzero E1A,E5A,E6A):
      E5A=E3A
      E6A=E2A
      W11=E1A*W1
      W12=E1A*W2
      W13=E1A*W3
      W14=E1A*W4
      W15=E1A*W5+E5A*W1
      W25=E5A*W2
      W35=E5A*W3
      W45=E5A*W4
      W55=E5A*W5
      W16=E1A*W6+E6A*W1
      W26=E6A*W2
      W36=E6A*W3
      W46=E6A*W4
      W56=E5A*W6+E6A*W5
      W66=E6A*W6
      DO 51 I=2,4
        G(I)=W11*A(I, 1)+W12*A(I, 2)+W13*A(I, 4)+W14*A(I, 7)+W15*A(I,11)
     *      +W25*A(I,12)+W35*A(I,13)+W45*A(I,14)+W55*A(I,15)+W16*A(I,16)
     *      +W26*A(I,17)+W36*A(I,18)+W46*A(I,19)+W56*A(I,20)+W66*A(I,21)
   51 CONTINUE
      GAA14=2.*G(2)
      GAA24=2.*G(3)
      GAA34=2.*G(4)
C     Derivative with respect to P2 (nonzero EA2,EA4,EA6):
      E4A=E3A
      E6A=E1A
      W12=E2A*W1
      W22=E2A*W2
      W23=E2A*W3
      W14=E4A*W1
      W24=E2A*W4+E4A*W2
      W34=E4A*W3
      W44=E4A*W4
      W25=E2A*W5
      W45=E4A*W5
      W16=E6A*W1
      W26=E2A*W6+E6A*W2
      W36=E6A*W3
      W46=E4A*W6+E6A*W4
      W56=E6A*W5
      W66=E6A*W6
      DO 52 I=2,4
        G(I)=W12*A(I, 2)+W22*A(I, 3)+W23*A(I, 5)+W14*A(I, 7)+W24*A(I, 8)
     *      +W34*A(I, 9)+W44*A(I,10)+W25*A(I,12)+W45*A(I,14)+W16*A(I,16)
     *      +W26*A(I,17)+W36*A(I,18)+W46*A(I,19)+W56*A(I,20)+W66*A(I,21)
   52 CONTINUE
      GAA15=2.*G(2)
      GAA25=2.*G(3)
      GAA35=2.*G(4)
C     Derivative with respect to P3 (nonzero EA3,EA4,EA5):
      E4A=E2A
      E5A=E1A
      W13=E3A*W1
      W23=E3A*W2
      W33=E3A*W3
      W14=E4A*W1
      W24=E4A*W2
      W34=E3A*W4+E4A*W3
      W44=E4A*W4
      W15=E5A*W1
      W25=E5A*W2
      W35=E3A*W5+E5A*W3
      W45=E4A*W5+E5A*W4
      W55=E5A*W5
      W36=E3A*W6
      W46=E4A*W6
      W56=E5A*W6
      DO 53 I=2,4
        G(I)=W13*A(I, 4)+W23*A(I, 5)+W33*A(I, 6)+W14*A(I, 7)+W24*A(I, 8)
     *      +W34*A(I, 9)+W44*A(I,10)+W15*A(I,11)+W25*A(I,12)+W35*A(I,13)
     *      +W45*A(I,14)+W55*A(I,15)+W36*A(I,18)+W46*A(I,19)+W56*A(I,20)
   53 CONTINUE
      GAA16=2.*G(2)
      GAA26=2.*G(3)
      GAA36=2.*G(4)
      RETURN
C
C-----------------------------------------------------------------------
C
      ENTRY GAAPP(E1A,E2A,E3A,GAA44,GAA45,GAA55,GAA46,GAA56,GAA66)
C
C     Subroutine calculating the slowness derivatives of the Christoffel
C     matrix.
C
      W11=E1A*E1A*2.
      W12=E1A*E2A*2.
      W22=E2A*E2A*2.
      W13=E1A*E3A*2.
      W23=E2A*E3A*2.
      W33=E3A*E3A*2.
C     GAA(i+3,l+3)=A(i,j,k,l)*EA(j)*EA(k):
      GAA44=A1111*W11+(A1121+A1211)*W12+(A1131+A1311)*W13
     *                      +A1221 *W22+(A1231+A1321)*W23+A1331*W33
      GAA45=A1112*W11+(A1122+A1212)*W12+(A1132+A1312)*W13
     *                      +A1222 *W22+(A1232+A1322)*W23+A1332*W33
      GAA55=A2112*W11+(A2122+A2212)*W12+(A2132+A2312)*W13
     *                      +A2222 *W22+(A2232+A2322)*W23+A2332*W33
      GAA46=A1113*W11+(A1123+A1213)*W12+(A1133+A1313)*W13
     *                      +A1223 *W22+(A1233+A1323)*W23+A1333*W33
      GAA56=A2113*W11+(A2123+A2213)*W12+(A2133+A2313)*W13
     *                      +A2223 *W22+(A2233+A2323)*W23+A2333*W33
      GAA66=A3113*W11+(A3123+A3213)*W12+(A3133+A3313)*W13
     *                      +A3223 *W22+(A3233+A3323)*W23+A3333*W33
      RETURN
      END
C
C=======================================================================
C