C 
C Subroutine file 'hder.for' to calculate the phase-space derivatives
C of the Hamiltonian in anisotropic media.
C
C Version: 7.20
C Date: 2015, May 19
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             Default: DEGREE=-1
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             General anisotropy (TIA1=0, TIA2=0 and TIA3=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             Approximately transversely isotropic medium:
C               KSWAVE=1:  Anisotropic-ray-theory SH-wave rays.
C               KSWAVE=2:  Anisotropic-ray-theory SV-wave rays.
C               Else (default):  Common S-wave rays.
C             Default: KSWAVE=0
C     TIA1=real, TIA2=real, TIA3=real:  Symmetry axis of a transversely
C             isotropic or approximately transversely isotropic (TI)
C             medium.  It can be specified only if DEGREE=2.
C             If specified when tracing the anisotropic-ray-theory
C             S-wave rays instead of common S-wave rays, the rays are
C             separated into SH rays (KSWAVE=1) and SV rays (KSWAVE=2).
C             Default: TIA1=0, TIA2=0, TIA3=0 (no TI symmetry specified)
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: 2015, May 19
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Temporary storage locations:
      LOGICAL LANIRT,LTI
      REAL HSN,HS2,E11,E21,E31,E12,E22,E32,E13,E23,E33,PP1,PP2,PP3
      REAL G11,G22,G33,GSV
      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,G1211,G2211,G1311,G2311,G3311
      REAL G1112,G1212,G2212,G1312,G2312,G3312
      REAL G1122,G1222,G2222,G1322,G2322,G3322
      REAL G1113,G1213,G2213,G1313,G2313,G3313
      REAL G1123,G1223,G2223,G1323,G2323,G3323
      REAL G1133,G1233,G2233,G1333,G2333,G3333
      REAL G1114,G1214,G2214,G1314,G2314,G3314
      REAL G1124,G1224,G2224,G1324,G2324,G3324
      REAL G1134,G1234,G2234,G1334,G2334,G3334
      REAL G1144,G1244,G2244,G1344,G2344,G3344
      REAL G1115,G1215,G2215,G1315,G2315,G3315
      REAL G1125,G1225,G2225,G1325,G2325,G3325
      REAL G1135,G1235,G2235,G1335,G2335,G3335
      REAL G1145,G1245,G2245,G1345,G2345,G3345
      REAL G1155,G1255,G2255,G1355,G2355,G3355
      REAL G1116,G1216,G2216,G1316,G2316,G3316
      REAL G1126,G1226,G2226,G1326,G2326,G3326
      REAL G1136,G1236,G2236,G1336,G2336,G3336
      REAL G1146,G1246,G2246,G1346,G2346,G3346
      REAL G1156,G1256,G2256,G1356,G2356,G3356
      REAL G1166,G1266,G2266,G1366,G2366,G3366
C     Weighting factors for P-wave rays and common S-wave rays
      REAL W11,W22,W33,W1111,W1212,W2222,W1313,W2323,WH
C     Additional weighting factors for SH and SV reference rays
      REAL W12,W13,W23
      REAL W1113,W1213,W2213,W1123,W1223,W2223,W1323,W1333,W2333
      REAL TE1,TE2,TE3,TETE,TE11,TE12,TE22,TE13,TE23,TE33
      REAL AUX
C
C     Degree of the homogeneous Hamiltonian and other input data:
      INTEGER KSWAVE
      REAL DEGREE,DEGRE2,TIA1,TIA2,TIA3,DSWAVE
      SAVE DEGREE,DEGRE2,TIA1,TIA2,TIA3,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('TIA1'  ,TIA1  ,0.)
        CALL RSEP3R('TIA2'  ,TIA2  ,0.)
        CALL RSEP3R('TIA3'  ,TIA3  ,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
      LTI=.FALSE.
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       Anisotropic-ray-theory S-wave ray tracing (general or TI)
        IF(TIA1.EQ.0..AND.TIA2.EQ.0..AND.TIA3.EQ.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
C         SH reference ray in an approximately TI medium
          IF(DEGREE.NE.2) THEN
C           5A4
            CALL ERROR('5A4 in HDER: DEGREE must be 2 in a TI medium')
          END IF
          LANIRT=.FALSE.
          LTI=.TRUE.
          TE1=TIA1*E11+TIA2*E21+TIA3*E31
          TE2=TIA1*E12+TIA2*E22+TIA3*E32
          TE3=TIA1*E13+TIA2*E23+TIA3*E33
          TETE=TE1*TE1+TE2*TE2
          TE11=TE1*TE1/TETE
          TE12=TE1*TE2/TETE
          TE22=TE2*TE2/TETE
          TE13=TE1*TE3/TETE
          TE23=TE2*TE3/TETE
          TE33=TE3*TE3/TETE
          GSV=TE11*G11+TE22*G22
C         Half the weighting factors by Klimes (2015) times multiplicity
          W11=0.5*TE22
          W22=0.5*TE11
          W12=   -TE12
          W13=   -TE13*(GSV-G11)/(G33-G11)
          W23=   -TE23*(GSV-G22)/(G33-G22)
          HS=2.*(W11*G11+W22*G22)
        END IF
      ELSE IF(KSWAVE.EQ.2.AND.ICB.LT.0) THEN
C       Anisotropic-ray-theory S-wave ray tracing (general or TI)
        IF(TIA1.EQ.0..AND.TIA2.EQ.0..AND.TIA3.EQ.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         SV reference ray in an approximately TI medium
          IF(DEGREE.NE.2) THEN
C           5A5
            CALL ERROR('5A5 in HDER: DEGREE must be 2 in a TI medium')
          END IF
          LANIRT=.FALSE.
          LTI=.TRUE.
          TE1=TIA1*E11+TIA2*E21+TIA3*E31
          TE2=TIA1*E12+TIA2*E22+TIA3*E32
          TE3=TIA1*E13+TIA2*E23+TIA3*E33
          TETE=TE1*TE1+TE2*TE2
          TE11=TE1*TE1/TETE
          TE12=TE1*TE2/TETE
          TE22=TE2*TE2/TETE
          TE13=TE1*TE3/TETE
          TE23=TE2*TE3/TETE
          TE33=TE3*TE3/TETE
          GSV=TE11*G11+TE22*G22
C         Half the weighting factors by Klimes (2015) times multiplicity
          W11=0.5*TE11
          W22=0.5*TE22
          W12=    TE12
          W13=    TE13*(GSV-G11)/(G33-G11)
          W23=    TE23*(GSV-G22)/(G33-G22)
          HS=GSV
        END IF
      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
        W33=0.5
        H1=W33*G331
        H2=W33*G332
        H3=W33*G333
        H4=W33*G334
        H5=W33*G335
        H6=W33*G336
C
        IF(IERR.EQ.0) THEN
          W1313=G33/(G33-G11)
          W2323=G33/(G33-G22)
        ELSE
          W1313=0.0
          W2323=0.0
        END IF
        H11=W33*G3311+W1313*G131*G131+W2323*G231*G231
        H12=W33*G3312+W1313*G131*G132+W2323*G231*G232
        H22=W33*G3322+W1313*G132*G132+W2323*G232*G232
        H13=W33*G3313+W1313*G131*G133+W2323*G231*G233
        H23=W33*G3323+W1313*G132*G133+W2323*G232*G233
        H33=W33*G3333+W1313*G133*G133+W2323*G233*G233
        H14=W33*G3314+W1313*G131*G134+W2323*G231*G234
        H24=W33*G3324+W1313*G132*G134+W2323*G232*G234
        H34=W33*G3334+W1313*G133*G134+W2323*G233*G234
        H44=W33*G3344+W1313*G134*G134+W2323*G234*G234
        H15=W33*G3315+W1313*G131*G135+W2323*G231*G235
        H25=W33*G3325+W1313*G132*G135+W2323*G232*G235
        H35=W33*G3335+W1313*G133*G135+W2323*G233*G235
        H45=W33*G3345+W1313*G134*G135+W2323*G234*G235
        H55=W33*G3355+W1313*G135*G135+W2323*G235*G235
        H16=W33*G3316+W1313*G131*G136+W2323*G231*G236
        H26=W33*G3326+W1313*G132*G136+W2323*G232*G236
        H36=W33*G3336+W1313*G133*G136+W2323*G233*G236
        H46=W33*G3346+W1313*G134*G136+W2323*G234*G236
        H56=W33*G3356+W1313*G135*G136+W2323*G235*G236
        H66=W33*G3366+W1313*G136*G136+W2323*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
        IF(LTI) THEN
C         SH or SV reference ray in an approximately TI medium.
C         Weighting factors for SV reference rays (Klimes, 2015, sec.4):
C         Weighting factors also used for common S-wave rays
          W1111=0.
          W2222=0.
          W1212=0.
          WH=0.
          W1313=((2.*GSV-G11-G11)*(TE11-0.25)*TE33-0.5*(G33-GSV)*TE11)
     *          *HS/((G33-G11)*(G33-G11))
          W2323=((2.*GSV-G22-G22)*(TE22-0.25)*TE33-0.5*(G33-GSV)*TE22)
     *          *HS/((G33-G22)*(G33-G22))
C         Additional weighting factors for SH and SV reference rays
          W1323=((2.*GSV-G11-G22)*(TE12     )*TE33-0.5*(G33-GSV)*TE12)
     *          *HS/((G33-G11)*(G33-G22))
          W1113=(TE11*TE13-    TE13*(G33-GSV)/(G33-G11))*HS/(G33-G11)
          W1213=(TE12*TE13-0.5*TE23*(G33-GSV)/(G33-G22))*HS/(G33-G11)
          W2213= TE22*TE13                              *HS/(G33-G11)
          W1123= TE11*TE23                              *HS/(G33-G22)
          W1223=(TE12*TE23-0.5*TE13*(G33-GSV)/(G33-G11))*HS/(G33-G22)
          W2223=(TE22*TE23-    TE23*(G33-GSV)/(G33-G22))*HS/(G33-G22)
          W1333=-W13*HS/(G33-G11)
          W2333=-W23*HS/(G33-G22)
          IF(KSWAVE.EQ.1) THEN
C           Weighting factors for SH reference rays:
C           Weighting factors also used for common S-wave rays
            W1313=-0.5*HS/(G33-G11)-W1313
            W2323=-0.5*HS/(G33-G22)-W2323
C           Additional weighting factors for SH and SV reference rays
            W1323=-W1323
            W1113=-W1113
            W1213=-W1213
            W2213=-W2213
            W1123=-W1123
            W1223=-W1223
            W2223=-W2223
          END IF
C         Half the weighting factors by Klimes (2015) times multiplicity
          W1313=2.*W1313
          W2323=2.*W2323
          W1323=2.*W1323
C         W1113=1.*W1113
          W1213=2.*W1213
C         W2213=1.*W2213
C         W1123=1.*W1123
          W1223=2.*W1223
C         W2223=1.*W2223
C         W1333=1.*W1333
C         W2333=1.*W2333
        ELSE
C         Common S-wave ray
          W11=0.25*G11**(DEGRE2-1.)*HS2/HSN
          W22=0.25*G22**(DEGRE2-1.)*HS2/HSN
          W1313=2.*W11*HS2/(G11-G33)
          W2323=2.*W22*HS2/(G22-G33)
          W1111=(DEGRE2-1.)*W11*HS2/G11
          W2222=(DEGRE2-1.)*W22*HS2/G22
          IF(DEGREE.EQ.-2) THEN
            W1212=-(G11+G22)/(G11*G22)**2
          ELSE IF(DEGREE.EQ.-1) THEN
            W1212=-(G11+SQRT(G11*G22)+G22)/(G11*SQRT(G22)+G22*SQRT(G11))
     *          /(G11*G22)
          ELSE IF(DEGREE.EQ.1) THEN
            W1212=-1./(G11*SQRT(G22)+G22*SQRT(G11))
          ELSE IF(DEGREE.EQ.2) THEN
            W1212=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
          W1212=0.5*W1212*HS2*HS2/HSN
          WH=(2.-DEGREE)
        END IF
C
        H1=W11*G111+W22*G221
        H2=W11*G112+W22*G222
        H3=W11*G113+W22*G223
        H4=W11*G114+W22*G224
        H5=W11*G115+W22*G225
        H6=W11*G116+W22*G226
        IF(LTI) THEN
          H1=W12*G121+W13*G131+W23*G231+H1
          H2=W12*G122+W13*G132+W23*G232+H2
          H3=W12*G123+W13*G133+W23*G233+H3
          H4=W12*G124+W13*G134+W23*G234+H4
          H5=W12*G125+W13*G135+W23*G235+H5
          H6=W12*G126+W13*G136+W23*G236+H6
        END IF
C
        H11=W11*G1111+W1111*G111*G111+W1313*G131*G131+W1212*G121*G121
        H12=W11*G1112+W1111*G111*G112+W1313*G131*G132+W1212*G121*G122
        H22=W11*G1122+W1111*G112*G112+W1313*G132*G132+W1212*G122*G122
        H13=W11*G1113+W1111*G111*G113+W1313*G131*G133+W1212*G121*G123
        H23=W11*G1123+W1111*G112*G113+W1313*G132*G133+W1212*G122*G123
        H33=W11*G1133+W1111*G113*G113+W1313*G133*G133+W1212*G123*G123
        H14=W11*G1114+W1111*G111*G114+W1313*G131*G134+W1212*G121*G124
        H24=W11*G1124+W1111*G112*G114+W1313*G132*G134+W1212*G122*G124
        H34=W11*G1134+W1111*G113*G114+W1313*G133*G134+W1212*G123*G124
        H44=W11*G1144+W1111*G114*G114+W1313*G134*G134+W1212*G124*G124
        H15=W11*G1115+W1111*G111*G115+W1313*G131*G135+W1212*G121*G125
        H25=W11*G1125+W1111*G112*G115+W1313*G132*G135+W1212*G122*G125
        H35=W11*G1135+W1111*G113*G115+W1313*G133*G135+W1212*G123*G125
        H45=W11*G1145+W1111*G114*G115+W1313*G134*G135+W1212*G124*G125
        H55=W11*G1155+W1111*G115*G115+W1313*G135*G135+W1212*G125*G125
        H16=W11*G1116+W1111*G111*G116+W1313*G131*G136+W1212*G121*G126
        H26=W11*G1126+W1111*G112*G116+W1313*G132*G136+W1212*G122*G126
        H36=W11*G1136+W1111*G113*G116+W1313*G133*G136+W1212*G123*G126
        H46=W11*G1146+W1111*G114*G116+W1313*G134*G136+W1212*G124*G126
        H56=W11*G1156+W1111*G115*G116+W1313*G135*G136+W1212*G125*G126
        H66=W11*G1166+W1111*G116*G116+W1313*G136*G136+W1212*G126*G126
        H11=W22*G2211+W2222*G221*G221+W2323*G231*G231+WH*H1*H1+H11
        H12=W22*G2212+W2222*G221*G222+W2323*G231*G232+WH*H1*H2+H12
        H22=W22*G2222+W2222*G222*G222+W2323*G232*G232+WH*H2*H2+H22
        H13=W22*G2213+W2222*G221*G223+W2323*G231*G233+WH*H1*H3+H13
        H23=W22*G2223+W2222*G222*G223+W2323*G232*G233+WH*H2*H3+H23
        H33=W22*G2233+W2222*G223*G223+W2323*G233*G233+WH*H3*H3+H33
        H14=W22*G2214+W2222*G221*G224+W2323*G231*G234+WH*H1*H4+H14
        H24=W22*G2224+W2222*G222*G224+W2323*G232*G234+WH*H2*H4+H24
        H34=W22*G2234+W2222*G223*G224+W2323*G233*G234+WH*H3*H4+H34
        H44=W22*G2244+W2222*G224*G224+W2323*G234*G234+WH*H4*H4+H44
        H15=W22*G2215+W2222*G221*G225+W2323*G231*G235+WH*H1*H5+H15
        H25=W22*G2225+W2222*G222*G225+W2323*G232*G235+WH*H2*H5+H25
        H35=W22*G2235+W2222*G223*G225+W2323*G233*G235+WH*H3*H5+H35
        H45=W22*G2245+W2222*G224*G225+W2323*G234*G235+WH*H4*H5+H45
        H55=W22*G2255+W2222*G225*G225+W2323*G235*G235+WH*H5*H5+H55
        H16=W22*G2216+W2222*G221*G226+W2323*G231*G236+WH*H1*H6+H16
        H26=W22*G2226+W2222*G222*G226+W2323*G232*G236+WH*H2*H6+H26
        H36=W22*G2236+W2222*G223*G226+W2323*G233*G236+WH*H3*H6+H36
        H46=W22*G2246+W2222*G224*G226+W2323*G234*G236+WH*H4*H6+H46
        H56=W22*G2256+W2222*G225*G226+W2323*G235*G236+WH*H5*H6+H56
        H66=W22*G2266+W2222*G226*G226+W2323*G236*G236+WH*H6*H6+H66
        IF(LTI) THEN
          CALL GABX(A,PP1,PP2,PP3,E13,E23,E33,E13,E23,E33,
     *            G331,G332,G333)
          CALL GABP(H,            E13,E23,E33,E13,E23,E33,
     *            G334,G335,G336)
          CALL GABXX(A,PP1,PP2,PP3,E11,E21,E31,E12,E22,E32,
     *            G1211,G1212,G1222,G1213,G1223,G1233)
          CALL GABXP(A,PP1,PP2,PP3,E11,E21,E31,E12,E22,E32,
     *            G1214,G1224,G1234,G1215,G1225,G1235,G1216,G1226,G1236)
          CALL GABPP(              E11,E21,E31,E12,E22,E32,
     *            G1244,G1245,G1255,G1246,G1256,G1266)
          CALL GABXX(A,PP1,PP2,PP3,E11,E21,E31,E13,E23,E33,
     *            G1311,G1312,G1322,G1313,G1323,G1333)
          CALL GABXP(A,PP1,PP2,PP3,E11,E21,E31,E13,E23,E33,
     *            G1314,G1324,G1334,G1315,G1325,G1335,G1316,G1326,G1336)
          CALL GABPP(              E11,E21,E31,E13,E23,E33,
     *            G1344,G1345,G1355,G1346,G1356,G1366)
          CALL GABXX(A,PP1,PP2,PP3,E12,E22,E32,E13,E23,E33,
     *            G2311,G2312,G2322,G2313,G2323,G2333)
          CALL GABXP(A,PP1,PP2,PP3,E12,E22,E32,E13,E23,E33,
     *            G2314,G2324,G2334,G2315,G2325,G2335,G2316,G2326,G2336)
          CALL GABPP(              E12,E22,E32,E13,E23,E33,
     *            G2344,G2345,G2355,G2346,G2356,G2366)
C
          H11=W12*G1211+W13*G1311+W23*G2311+H11
          H12=W12*G1212+W13*G1312+W23*G2312+H12
          H22=W12*G1222+W13*G1322+W23*G2322+H22
          H13=W12*G1213+W13*G1313+W23*G2313+H13
          H23=W12*G1223+W13*G1323+W23*G2323+H23
          H33=W12*G1233+W13*G1333+W23*G2333+H33
          H14=W12*G1214+W13*G1314+W23*G2314+H14
          H24=W12*G1224+W13*G1324+W23*G2324+H24
          H34=W12*G1234+W13*G1334+W23*G2334+H34
          H44=W12*G1244+W13*G1344+W23*G2344+H44
          H15=W12*G1215+W13*G1315+W23*G2315+H15
          H25=W12*G1225+W13*G1325+W23*G2325+H25
          H35=W12*G1235+W13*G1335+W23*G2335+H35
          H45=W12*G1245+W13*G1345+W23*G2345+H45
          H55=W12*G1255+W13*G1355+W23*G2355+H55
          H16=W12*G1216+W13*G1316+W23*G2316+H16
          H26=W12*G1226+W13*G1326+W23*G2326+H26
          H36=W12*G1236+W13*G1336+W23*G2336+H36
          H46=W12*G1246+W13*G1346+W23*G2346+H46
          H56=W12*G1256+W13*G1356+W23*G2356+H56
          H66=W12*G1266+W13*G1366+W23*G2366+H66
C
          H11=(W1113*G111+W1213*G121+W2213*G221+W1333*G331)*G131+H11
          H12=(W1113*G111+W1213*G121+W2213*G221+W1333*G331)*G132+H12
          H22=(W1113*G112+W1213*G122+W2213*G222+W1333*G332)*G132+H22
          H13=(W1113*G111+W1213*G121+W2213*G221+W1333*G331)*G133+H13
          H23=(W1113*G112+W1213*G122+W2213*G222+W1333*G332)*G133+H23
          H33=(W1113*G113+W1213*G123+W2213*G223+W1333*G333)*G133+H33
          H14=(W1113*G111+W1213*G121+W2213*G221+W1333*G331)*G134+H14
          H24=(W1113*G112+W1213*G122+W2213*G222+W1333*G332)*G134+H24
          H34=(W1113*G113+W1213*G123+W2213*G223+W1333*G333)*G134+H34
          H44=(W1113*G114+W1213*G124+W2213*G224+W1333*G334)*G134+H44
          H15=(W1113*G111+W1213*G121+W2213*G221+W1333*G331)*G135+H15
          H25=(W1113*G112+W1213*G122+W2213*G222+W1333*G332)*G135+H25
          H35=(W1113*G113+W1213*G123+W2213*G223+W1333*G333)*G135+H35
          H45=(W1113*G114+W1213*G124+W2213*G224+W1333*G334)*G135+H45
          H55=(W1113*G115+W1213*G125+W2213*G225+W1333*G335)*G135+H55
          H16=(W1113*G111+W1213*G121+W2213*G221+W1333*G331)*G136+H16
          H26=(W1113*G112+W1213*G122+W2213*G222+W1333*G332)*G136+H26
          H36=(W1113*G113+W1213*G123+W2213*G223+W1333*G333)*G136+H36
          H46=(W1113*G114+W1213*G124+W2213*G224+W1333*G334)*G136+H46
          H56=(W1113*G115+W1213*G125+W2213*G225+W1333*G335)*G136+H56
          H66=(W1113*G116+W1213*G126+W2213*G226+W1333*G336)*G136+H66
C
          H11=(W1113*G111+W1213*G121+W2213*G221+W1333*G331)*G131+H11
          H12=(W1113*G112+W1213*G122+W2213*G222+W1333*G332)*G131+H12
          H22=(W1113*G112+W1213*G122+W2213*G222+W1333*G332)*G132+H22
          H13=(W1113*G113+W1213*G123+W2213*G223+W1333*G333)*G131+H13
          H23=(W1113*G113+W1213*G123+W2213*G223+W1333*G333)*G132+H23
          H33=(W1113*G113+W1213*G123+W2213*G223+W1333*G333)*G133+H33
          H14=(W1113*G114+W1213*G124+W2213*G224+W1333*G334)*G131+H14
          H24=(W1113*G114+W1213*G124+W2213*G224+W1333*G334)*G132+H24
          H34=(W1113*G114+W1213*G124+W2213*G224+W1333*G334)*G133+H34
          H44=(W1113*G114+W1213*G124+W2213*G224+W1333*G334)*G134+H44
          H15=(W1113*G115+W1213*G125+W2213*G225+W1333*G335)*G131+H15
          H25=(W1113*G115+W1213*G125+W2213*G225+W1333*G335)*G132+H25
          H35=(W1113*G115+W1213*G125+W2213*G225+W1333*G335)*G133+H35
          H45=(W1113*G115+W1213*G125+W2213*G225+W1333*G335)*G134+H45
          H55=(W1113*G115+W1213*G125+W2213*G225+W1333*G335)*G135+H55
          H16=(W1113*G116+W1213*G126+W2213*G226+W1333*G336)*G131+H16
          H26=(W1113*G116+W1213*G126+W2213*G226+W1333*G336)*G132+H26
          H36=(W1113*G116+W1213*G126+W2213*G226+W1333*G336)*G133+H36
          H46=(W1113*G116+W1213*G126+W2213*G226+W1333*G336)*G134+H46
          H56=(W1113*G116+W1213*G126+W2213*G226+W1333*G336)*G135+H56
          H66=(W1113*G116+W1213*G126+W2213*G226+W1333*G336)*G136+H66
C
          H11=(W1123*G111+W1223*G121+W2223*G221+W2333*G331)*G231+H11
          H12=(W1123*G111+W1223*G121+W2223*G221+W2333*G331)*G232+H12
          H22=(W1123*G112+W1223*G122+W2223*G222+W2333*G332)*G232+H22
          H13=(W1123*G111+W1223*G121+W2223*G221+W2333*G331)*G233+H13
          H23=(W1123*G112+W1223*G122+W2223*G222+W2333*G332)*G233+H23
          H33=(W1123*G113+W1223*G123+W2223*G223+W2333*G333)*G233+H33
          H14=(W1123*G111+W1223*G121+W2223*G221+W2333*G331)*G234+H14
          H24=(W1123*G112+W1223*G122+W2223*G222+W2333*G332)*G234+H24
          H34=(W1123*G113+W1223*G123+W2223*G223+W2333*G333)*G234+H34
          H44=(W1123*G114+W1223*G124+W2223*G224+W2333*G334)*G234+H44
          H15=(W1123*G111+W1223*G121+W2223*G221+W2333*G331)*G235+H15
          H25=(W1123*G112+W1223*G122+W2223*G222+W2333*G332)*G235+H25
          H35=(W1123*G113+W1223*G123+W2223*G223+W2333*G333)*G235+H35
          H45=(W1123*G114+W1223*G124+W2223*G224+W2333*G334)*G235+H45
          H55=(W1123*G115+W1223*G125+W2223*G225+W2333*G335)*G235+H55
          H16=(W1123*G111+W1223*G121+W2223*G221+W2333*G331)*G236+H16
          H26=(W1123*G112+W1223*G122+W2223*G222+W2333*G332)*G236+H26
          H36=(W1123*G113+W1223*G123+W2223*G223+W2333*G333)*G236+H36
          H46=(W1123*G114+W1223*G124+W2223*G224+W2333*G334)*G236+H46
          H56=(W1123*G115+W1223*G125+W2223*G225+W2333*G335)*G236+H56
          H66=(W1123*G116+W1223*G126+W2223*G226+W2333*G336)*G236+H66
C
          H11=(W1123*G111+W1223*G121+W2223*G221+W2333*G331)*G231+H11
          H12=(W1123*G112+W1223*G122+W2223*G222+W2333*G332)*G231+H12
          H22=(W1123*G112+W1223*G122+W2223*G222+W2333*G332)*G232+H22
          H13=(W1123*G113+W1223*G123+W2223*G223+W2333*G333)*G231+H13
          H23=(W1123*G113+W1223*G123+W2223*G223+W2333*G333)*G232+H23
          H33=(W1123*G113+W1223*G123+W2223*G223+W2333*G333)*G233+H33
          H14=(W1123*G114+W1223*G124+W2223*G224+W2333*G334)*G231+H14
          H24=(W1123*G114+W1223*G124+W2223*G224+W2333*G334)*G232+H24
          H34=(W1123*G114+W1223*G124+W2223*G224+W2333*G334)*G233+H34
          H44=(W1123*G114+W1223*G124+W2223*G224+W2333*G334)*G234+H44
          H15=(W1123*G115+W1223*G125+W2223*G225+W2333*G335)*G231+H15
          H25=(W1123*G115+W1223*G125+W2223*G225+W2333*G335)*G232+H25
          H35=(W1123*G115+W1223*G125+W2223*G225+W2333*G335)*G233+H35
          H45=(W1123*G115+W1223*G125+W2223*G225+W2333*G335)*G234+H45
          H55=(W1123*G115+W1223*G125+W2223*G225+W2333*G335)*G235+H55
          H16=(W1123*G116+W1223*G126+W2223*G226+W2333*G336)*G231+H16
          H26=(W1123*G116+W1223*G126+W2223*G226+W2333*G336)*G232+H26
          H36=(W1123*G116+W1223*G126+W2223*G226+W2333*G336)*G233+H36
          H46=(W1123*G116+W1223*G126+W2223*G226+W2333*G336)*G234+H46
          H56=(W1123*G116+W1223*G126+W2223*G226+W2333*G336)*G235+H56
          H66=(W1123*G116+W1223*G126+W2223*G226+W2333*G336)*G236+H66
C
          H11=W1323*(G131*G231+G131*G231)+H11
          H12=W1323*(G131*G232+G132*G231)+H12
          H22=W1323*(G132*G232+G132*G232)+H22
          H13=W1323*(G131*G233+G133*G231)+H13
          H23=W1323*(G132*G233+G133*G232)+H23
          H33=W1323*(G133*G233+G133*G233)+H33
          H14=W1323*(G131*G234+G134*G231)+H14
          H24=W1323*(G132*G234+G134*G232)+H24
          H34=W1323*(G133*G234+G134*G233)+H34
          H44=W1323*(G134*G234+G134*G234)+H44
          H15=W1323*(G131*G235+G135*G231)+H15
          H25=W1323*(G132*G235+G135*G232)+H25
          H35=W1323*(G133*G235+G135*G233)+H35
          H45=W1323*(G134*G235+G135*G234)+H45
          H55=W1323*(G135*G235+G135*G235)+H55
          H16=W1323*(G131*G236+G136*G231)+H16
          H26=W1323*(G132*G236+G136*G232)+H26
          H36=W1323*(G133*G236+G136*G233)+H36
          H46=W1323*(G134*G236+G136*G234)+H46
          H56=W1323*(G135*G236+G136*G235)+H56
          H66=W1323*(G136*G236+G136*G236)+H66
        END IF
      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    *            GAB11,GAB12,GAB22,GAB13,GAB23,GAB33)
C     Subroutine calculating the second spatial derivatives of the
C     diagonal elements of the Christoffel matrix.
C
C     ENTRY GAAXP(A,P1,P2,P3,E1A,E2A,E3A,
C    *            GAB14,GAB24,GAB34,GAB15,GAB25,GAB35,GAB16,GAB26,GAB36)
C     Subroutine calculating the second mixed derivatives of the
C     diagonal elements of the Christoffel matrix.
C
C     ENTRY GAAPP(E1A,E2A,E3A,GAB44,GAB45,GAB55,GAB46,GAB56,GAB66)
C     Subroutine calculating the second slowness derivatives of the
C     diagonal elements of the Christoffel matrix.
C
C     ENTRY GABXX(A,P1,P2,P3,E1A,E2A,E3A,E1B,E2B,E3B,
C    *            GAB11,GAB12,GAB22,GAB13,GAB23,GAB33)
C     Subroutine calculating the second spatial derivatives of the
C     off-diagonal elements of the Christoffel matrix.
C
C     ENTRY GABXP(A,P1,P2,P3,E1A,E2A,E3A,E1B,E2B,E3B,
C    *            GAB14,GAB24,GAB34,GAB15,GAB25,GAB35,GAB16,GAB26,GAB36)
C     Subroutine calculating the second mixed derivatives of the
C     off-diagonal elements of the Christoffel matrix.
C
C     ENTRY GABPP(E1A,E2A,E3A,E1B,E2B,E3B,
c    *            GAB44,GAB45,GAB55,GAB46,GAB56,GAB66)
C     Subroutine calculating the second slowness derivatives of the
C     off-diagonal elements of the 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 GAB11,GAB12,GAB22,GAB13,GAB23,GAB33
      REAL GAB14,GAB24,GAB34,GAB15,GAB25,GAB35,GAB16,GAB26,GAB36
      REAL GAB44,GAB45,GAB55,GAB46,GAB56,GAB66
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     GAB11,GAB12,GAB22,GAB13,GAB23,GAB33,GAB14,GAB24,GAB34,GAB15,GAB25,
C     GAB35,GAB16,GAB26,GAB36,GAB44,GAB45,GAB55,GAB46,GAB56,GAB66...
C             Second partial phase-space derivatives of derivatives of
C             the Christoffel matrix, twice multiplied by the given
C             eigenvector EA, or multiplied by the given eigenvectors EA
C             and EB.
C
C Subroutines and external functions required:
      EXTERNAL EIGEN
C     EIGEN.. File 'eigen.for'.
C
C Date: 2015, April 29
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,E4B,E5B,E6B,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,
     *            GAB11,GAB12,GAB22,GAB13,GAB23,GAB33)
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
      GAB11=G(5)
      GAB12=G(6)
      GAB22=G(7)
      GAB13=G(8)
      GAB23=G(9)
      GAB33=G(10)
      RETURN
C
C-----------------------------------------------------------------------
C
      ENTRY GAAXP(A,P1,P2,P3,E1A,E2A,E3A,
     *            GAB14,GAB24,GAB34,GAB15,GAB25,GAB35,GAB16,GAB26,GAB36)
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     E4A=dW4/dPi, E5A=dW5/dPi, E6A=dW6/dPi
C
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
      GAB14=2.*G(2)
      GAB24=2.*G(3)
      GAB34=2.*G(4)
C
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
      GAB15=2.*G(2)
      GAB25=2.*G(3)
      GAB35=2.*G(4)
C
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
      GAB16=2.*G(2)
      GAB26=2.*G(3)
      GAB36=2.*G(4)
      RETURN
C
C-----------------------------------------------------------------------
C
      ENTRY GAAPP(E1A,E2A,E3A,GAB44,GAB45,GAB55,GAB46,GAB56,GAB66)
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     GAB(i+3,l+3)=[A(i,j,k,l)+A(i,k,j,l)]*EA(j)*EA(k):
      GAB44=A1111*W11+(A1121+A1211)*W12+(A1131+A1311)*W13
     *                      +A1221 *W22+(A1231+A1321)*W23+A1331*W33
      GAB45=A1112*W11+(A1122+A1212)*W12+(A1132+A1312)*W13
     *                      +A1222 *W22+(A1232+A1322)*W23+A1332*W33
      GAB55=A2112*W11+(A2122+A2212)*W12+(A2132+A2312)*W13
     *                      +A2222 *W22+(A2232+A2322)*W23+A2332*W33
      GAB46=A1113*W11+(A1123+A1213)*W12+(A1133+A1313)*W13
     *                      +A1223 *W22+(A1233+A1323)*W23+A1333*W33
      GAB56=A2113*W11+(A2123+A2213)*W12+(A2133+A2313)*W13
     *                      +A2223 *W22+(A2233+A2323)*W23+A2333*W33
      GAB66=A3113*W11+(A3123+A3213)*W12+(A3133+A3313)*W13
     *                      +A3223 *W22+(A3233+A3323)*W23+A3333*W33
      RETURN
C
C-----------------------------------------------------------------------
C
      ENTRY GABXX(A,P1,P2,P3,E1A,E2A,E3A,E1B,E2B,E3B,
     *            GAB11,GAB12,GAB22,GAB13,GAB23,GAB33)
C
C     Subroutine calculating the spatial derivatives of the Christoffel
C     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 70 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)
   70 CONTINUE
      GAB11=G(5)
      GAB12=G(6)
      GAB22=G(7)
      GAB13=G(8)
      GAB23=G(9)
      GAB33=G(10)
      RETURN
C
C-----------------------------------------------------------------------
C
      ENTRY GABXP(A,P1,P2,P3,E1A,E2A,E3A,E1B,E2B,E3B,
     *            GAB14,GAB24,GAB34,GAB15,GAB25,GAB35,GAB16,GAB26,GAB36)
C
C     Subroutine calculating the mixed derivatives of the Christoffel
C     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
C     E4A=dW4A/dPi, E5A=dW5A/dPi, E6A=dW6A/dPi
C     E4B=dW4B/dPi, E5B=dW5B/dPi, E6B=dW6B/dPi
C
C     Derivative with respect to P1 (nonzero E1A,E5A,E6A,E1B,E5B,E6B):
      E5A=E3A
      E6A=E2A
      E5B=E3B
      E6B=E2B
      W11=E1A*W1B+W1A*E1B
      W12=E1A*W2B+W2A*E1B
      W13=E1A*W3B+W3A*E1B
      W14=E1A*W4B+W4A*E1B
      W15=E1A*W5B+W5A*E1B+E5A*W1B+W1A*E5B
      W25=E5A*W2B+W2A*E5B
      W35=E5A*W3B+W3A*E5B
      W45=E5A*W4B+W4A*E5B
      W55=E5A*W5B+W5A*E5B
      W16=E1A*W6B+W6A*E1B+E6A*W1B+W1A*E6B
      W26=E6A*W2B+W2A*E6B
      W36=E6A*W3B+W3A*E6B
      W46=E6A*W4B+W4A*E6B
      W56=E5A*W6B+W6A*E5B+E6A*W5B+W5A*E6B
      W66=E6A*W6B+W6A*E6B
      DO 81 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)
   81 CONTINUE
      GAB14=G(2)
      GAB24=G(3)
      GAB34=G(4)
C
C     Derivative with respect to P2 (nonzero EA2,EA4,EA6,EB2,EB4,EB6):
      E4A=E3A
      E6A=E1A
      E4B=E3B
      E6B=E1B
      W12=E2A*W1B+W1A*E2B
      W22=E2A*W2B+W2A*E2B
      W23=E2A*W3B+W3A*E2B
      W14=E4A*W1B+W1A*E4B
      W24=E2A*W4B+W4A*E2B+E4A*W2B+W2A*E4B
      W34=E4A*W3B+W3A*E4B
      W44=E4A*W4B+W4A*E4B
      W25=E2A*W5B+W5A*E2B
      W45=E4A*W5B+W5A*E4B
      W16=E6A*W1B+W1A*E6B
      W26=E2A*W6B+W6A*E2B+E6A*W2B+W2A*E6B
      W36=E6A*W3B+W3A*E6B
      W46=E4A*W6B+W6A*E4B+E6A*W4B+W4A*E6B
      W56=E6A*W5B+W5A*E6B
      W66=E6A*W6B+W6A*E6B
      DO 82 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)
   82 CONTINUE
      GAB15=G(2)
      GAB25=G(3)
      GAB35=G(4)
C
C     Derivative with respect to P3 (nonzero EA3,EA4,EA5,EB3,EB4,EB5):
      E4A=E2A
      E5A=E1A
      E4B=E2B
      E5B=E1B
      W13=E3A*W1B+W1A*E3B
      W23=E3A*W2B+W2A*E3B
      W33=E3A*W3B+W3A*E3B
      W14=E4A*W1B+W1A*E4B
      W24=E4A*W2B+W2A*E4B
      W34=E3A*W4B+W4A*E3B+E4A*W3B+W3A*E4B
      W44=E4A*W4B+W4A*E4B
      W15=E5A*W1B+W1A*E5B
      W25=E5A*W2B+W2A*E5B
      W35=E3A*W5B+W5A*E3B+E5A*W3B+W3A*E5B
      W45=E4A*W5B+W5A*E4B+E5A*W4B+W4A*E5B
      W55=E5A*W5B+W5A*E5B
      W36=E3A*W6B+W6A*E3B
      W46=E4A*W6B+W6A*E4B
      W56=E5A*W6B+W6A*E5B
      DO 83 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)
   83 CONTINUE
      GAB16=G(2)
      GAB26=G(3)
      GAB36=G(4)
      RETURN
C
C-----------------------------------------------------------------------
C
      ENTRY GABPP(E1A,E2A,E3A,E1B,E2B,E3B,
     *            GAB44,GAB45,GAB55,GAB46,GAB56,GAB66)
C
C     Subroutine calculating the slowness derivatives of the Christoffel
C     matrix.
C
      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,l+3)=[A(i,j,k,l)+A(i,k,j,l)]*EA(j)*EB(k):
      GAB44=A1111*W11+(A1121+A1211)*W12+(A1131+A1311)*W13
     *                      +A1221 *W22+(A1231+A1321)*W23+A1331*W33
      GAB45=A1112*W11+(A1122+A1212)*W12+(A1132+A1312)*W13
     *                      +A1222 *W22+(A1232+A1322)*W23+A1332*W33
      GAB55=A2112*W11+(A2122+A2212)*W12+(A2132+A2312)*W13
     *                      +A2222 *W22+(A2232+A2322)*W23+A2332*W33
      GAB46=A1113*W11+(A1123+A1213)*W12+(A1133+A1313)*W13
     *                      +A1223 *W22+(A1233+A1323)*W23+A1333*W33
      GAB56=A2113*W11+(A2123+A2213)*W12+(A2133+A2313)*W13
     *                      +A2223 *W22+(A2233+A2323)*W23+A2333*W33
      GAB66=A3113*W11+(A3123+A3213)*W12+(A3133+A3313)*W13
     *                      +A3223 *W22+(A3233+A3323)*W23+A3333*W33
      RETURN
      END
C
C=======================================================================
C