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