C 
C Subroutine file 'hder.for' to calculate the phase-space derivatives
C of the Hamiltonian in anisotropic media.
C
C Date: 2003, May 15
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=======================================================================
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)
      INTEGER ICB
      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
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...    Elastic moduli and their first and second spatial
C             derivatives.
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
C Input SEP parameter:
C     DEGREE='real'... Degree of the homogeneous Hamiltonian to be
C             arithmetically averaged for common S-wave ray tracing.
C
C Subroutines and external functions required:
      EXTERNAL GAA,GABX,GABP,GAAXX,GAAXP,GAAPP,EIGEN,RSEP3R
C     Subroutine GAA and its entries GABX,GABP,GAAXX,GAAXP,GAAPP...
C             This file.
C     EIGEN.. File 'eigen.for'.
C     RSEP3R. File 'sep.for'.
C
C Date: 2003, May 6
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Temporary storage locations:
      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
C
C     Degree of the homogeneous Hamiltonian:
      REAL DEGREE,DEGRE2
      SAVE DEGREE,DEGRE2
      DATA DEGREE/0./
      IF(DEGREE.EQ.0.) THEN
        CALL RSEP3R('DEGREE',DEGREE,-1.)
        DEGRE2=0.5*DEGREE
      END IF
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)
C
C     Twice the homogeneous Hamiltonian of the second degree:
      HP=G33
      HSN=0.5*(G11**DEGRE2+G22**DEGRE2)
      HS2=HSN**(1./DEGRE2)
      HS=HS2
      IF(ICB.GT.0) THEN
C       P wave:
        H=G33
      ELSE
C       S wave:
        H=HS
      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(ICB.GT.0) THEN
C       P 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
        W13=G33/(G33-G11)
        W23=G33/(G33-G22)
        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 IF(ICB.LT.0) THEN
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
      ELSE
C       5A2
        CALL ERROR('5A2 in HDER: Zero wave type')
C       This error should not appear, contact the authors.
      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 spatial derivatives of the Christoffel
C     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 mixed derivatives of the Christoffel
C     matrix.
C
C     ENTRY GAAPP(E1A,E2A,E3A,GAA44,GAA45,GAA55,GAA46,GAA56,GAA66)
C     Subroutine calculating the slowness derivatives of the Christoffel
C     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