C
C Subroutine file 'hder.for' to calculate the phase-space derivatives C of the Hamiltonian function in anisotropic elastic or viscoelastic C media. C C Version: 8.10 C Date: 2022, December 5 C Coded by Ludek Klimes C C....................................................................... C C This file consists of the following external procedures: C HDER... Interface subroutine to HDERC or HDERR designed to C calculate the phase-space derivatives of the real-valued C reference Hamiltonian function in anisotropic viscoelastic C or elastic media, depending on input SEP parameter KQ. C HDERC.. Subroutine designed to calculate the phase-space C derivatives of the complex-valued Hamiltonian function C in anisotropic viscoelastic media. C CHR0... Auxiliary subroutine to HDERC designed to calculate C the complex-valued Christoffel matrix, with entries CHR1, C CHR11 and CHR2 which calculate the linear combinations of C its first-order and second-order phase-space derivatives. C SQRTC.. Auxiliary subroutine to HDERC designed to calculate the C square root of a complex number. C HDERR.. Subroutine designed to calculate the phase-space C derivatives of the real-valued Hamiltonian function C in anisotropic elastic media. C Identical to subroutine HDER version 7.20. C GAA... Auxiliary subroutine to HDERR calculating the eigenvalues C and eigenvectors of the real-valued Christoffel matrix, C with entries GABX, GABP, GAAXX, GAAXP and GAAPP which C calculate the first-order and second-order phase-space C derivatives of the Christoffel matrix, transformed into C the given eigenvectors of the Christoffel matrix. C C======================================================================= C SUBROUTINE HDER(ICB,A,Q,P1,P2,P3,HP,HS,H,HI,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),Q(10,21),P1,P2,P3,HP,HS,H,HI,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 real-valued reference homogeneous Hamiltonian function of the second C degree in anisotropic viscoelastic or elastic media. C C In viscoelastic media, subroutine HDERC calculates the complex-valued C homogeneous Hamiltonian function of degree DEGREE. Its real part C represents the real-valued reference Hamiltonian function, and is then C converted into the real-valued reference homogeneous Hamiltonian C function of the second degree. C In elastic media, subroutine HDERR calculates directly the real-valued C homogeneous Hamiltonian function of the second degree. C C Input: C ICB... Wave type: C ICB.GT.0: P wave, C ICB.LT.0: S wave. C A... Real parts of the density-reduced viscoelastic moduli C and their first and second spatial derivatives. C The order of the value, first and second partial C derivatives of each parameter Aij (first 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 Q... Imaginary parts of the density-reduced viscoelastic moduli C and their first and second spatial derivatives, C ordered analogoustly to the real parts. C Not used in elastic media. C P1,P2,P3... Slowness vector. C C Output: C HP,HS.. Twice the value of the real-valued reference Hamiltonian C function corresponding to P and S waves, respectively. C The output Hamiltonian function is homogeneous of the C second degree with respect to the slowness vector. C Definition of HP and HS depends on the value of the input C SEP parameter DEGREE. C H... Twice the value of the reference homogeneous Hamiltonian C function of the second degree corresponding to the given C wave, i.e. HP for a P wave or HS for an S wave. C H should be equal to 1 if P1,P2,P3 is the reference C slowness vector. C HI... Imaginary part of the complex-valued homogeneous C Hamiltonian function of degree DEGREE, calculated for the C normalized slowness vector, i.e. divided by H**(DEGREE/2). C Note that the first-order perturbation of the imaginary C part of travel time is represented by the integral of -HI C with respect to the reference travel time. C HI=0 in elastic media. C H1,H2,H3,H4,H5,H6... First partial phase-space derivatives of the C second-degree Hamiltonian function. The derivatives are C calculated for the normalized slowness vector (H1,H2,H3 C divided by the value of the reference Hamiltonian function C H, H4,H5,H6 divided by the square root of the value of the C the reference Hamiltonian function H). 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 function. The derivatives C are calculated for the normalized slowness vector C (H11,H12,H22,H13,H23,H33 divided by the value of the C reference Hamiltonian function H, H14,H24,H34,H15,H25,H35, C H16,H26,H36 divided by the square root of the value of the C the reference Hamiltonian function H). 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 in an elastic medium. 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 KQ=integer... Switch between the viscoelastic and elastic media. C KQ=0: Imaginary part of the stiffness tensor is ignored. C The Hamiltonian function is thus real-valued. C Otherwise: Reference Hamiltonian function equals the real C part of the complex-valued Hamiltonian function. C DEGREE=real... Degree of the homogeneous complex-valued C Hamiltonian function. C Simultaneously the degree of the homogeneous Hamiltonian C function to be arithmetically averaged for common S-wave C ray tracing. C DEGREE must equal 2, 1, -1 or -2. C Default: DEGREE=-1 C KSWAVE=integer... Selection of the anisotropic-ray-theory S-wave C rays instead of common S-wave rays in elastic media. C Affects only S waves (ICB.LT.0). C There must be KSWAVE=0 (common S-wave rays) if KQ.NE.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 elastic 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 Considered only if KQ=0 and (KSWAVE=1 or 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 KQ=0, 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 HDERR,HDERC C Subroutines HDERR,HDERC... This file. EXTERNAL ERROR,RSEP3I,RSEP3R C ERROR.. File 'error.for'. C RSEP3I,RSEP3R... File 'sep.for'. C Called by HDERC: C Subroutine CHR0 and its entries CHR1,CHR11,CHR2... This file. C SQRTC.. This file. C CUBIC.. File 'cubic.for'. C Called by HDERR: C Subroutine GAA and its entries GABX,GABP,GAAXX,GAAXP,GAAPP... C This file. C EIGEN.. File 'eigen.for'. C C Reference: C Klimes, L. (2006): C Common-ray tracing and dynamic ray tracing for S waves C in a smooth elastic anisotropic medium. C Stud. geophys. geod., 50, 449-461. C C Date: 2022, December 5 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Temporary storage locations: REAL DH,WH,WH1,WH2 C C Imaginary parts of complex-valued arguments: REAL HPI,HSI,H1I,H2I,H3I,H4I,H5I,H6I REAL H11I,H12I,H22I,H13I,H23I,H33I,H14I,H24I,H34I,H44I REAL H15I,H25I,H35I,H45I,H55I,H16I,H26I,H36I,H46I,H56I,H66I C C Degree of homogeneous Hamiltonian function and other input data: INTEGER KQ REAL DEGREE SAVE DEGREE,KQ DATA DEGREE/0./ C C....................................................................... C C Input SEP parameters IF(DEGREE.EQ.0.) THEN CALL RSEP3I('KQ',KQ,0) CALL RSEP3R('DEGREE',DEGREE,-1.) END IF C IF(KQ.EQ.0) THEN C C Real-valued stiffness tensor and Hamiltonian function: CALL HDERR(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) HI=0. C ELSE C Complex-valued stiffness tensor and Hamiltonian function: CALL HDERC(ICB,A,Q,P1,P2,P3,HP,HPI,HS,HSI,H,HI, * H1,H1I,H2,H2I,H3,H3I,H4,H4I,H5,H5I,H6,H6I, * H11,H11I,H12,H12I,H22,H22I,H13,H13I, * H23,H23I,H33,H33I,H14,H14I,H24,H24I, * H34,H34I,H44,H44I,H15,H15I,H25,H25I, * H35,H35I,H45,H45I,H55,H55I,H16,H16I, * H26,H26I,H36,H36I,H46,H46I,H56,H56I, * H66,H66I,E) IERR=0 C The reference real-valued Hamiltonian function equals C the real part of the complex-valued Hamiltonian function. C C Second-degree Hamiltonian function (Klimes, 2006, eqs.36-38) C and the renormalized imaginary part of the Hamiltonian function IF(DEGREE.EQ.2.) THEN HP=HP+HP HS=HS+HS H=H+H HI=HI/H ELSE HP=DEGREE*HP HS=DEGREE*HS DH=DEGREE*H IF(DEGREE.EQ.1.) THEN HP=HP*HP HS=HS*HS WH=DH*DH ELSE IF(DEGREE.EQ.-1.) THEN HP=1./(HP*HP) HS=1./(HS*HS) WH=1./(DH*DH) ELSE IF(DEGREE.EQ.-2.) THEN HP=1./HP HS=1./HS WH=1./DH ELSE C 5A6 CALL ERROR('5A6 in HDER: Wrong degree of Hamiltonian') C Only homogeneous Hamiltonian functions of degrees -2, -1, 1 C and 2 are now allowed. END IF HI=HI/DH WH1=WH/DH WH2=(2.-DEGREE)*WH1/DH H11=WH2*H1*H1+WH1*H11 H12=WH2*H1*H2+WH1*H12 H22=WH2*H2*H2+WH1*H22 H13=WH2*H1*H3+WH1*H13 H23=WH2*H2*H3+WH1*H23 H33=WH2*H3*H3+WH1*H33 H14=WH2*H1*H4+WH1*H14 H24=WH2*H2*H4+WH1*H24 H34=WH2*H3*H4+WH1*H34 H44=WH2*H4*H4+WH1*H44 H15=WH2*H1*H5+WH1*H15 H25=WH2*H2*H5+WH1*H25 H35=WH2*H3*H5+WH1*H35 H45=WH2*H4*H5+WH1*H45 H55=WH2*H5*H5+WH1*H55 H16=WH2*H1*H6+WH1*H16 H26=WH2*H2*H6+WH1*H26 H36=WH2*H3*H6+WH1*H36 H46=WH2*H4*H6+WH1*H46 H56=WH2*H5*H6+WH1*H56 H66=WH2*H6*H6+WH1*H66 H1=WH1*H1 H2=WH1*H2 H3=WH1*H3 H4=WH1*H4 H5=WH1*H5 H6=WH1*H6 H=WH END IF C C Renormalizing the derivatives of the Hamiltonian function WH=H H1=H1/WH H2=H2/WH H3=H3/WH H11=H11/WH H12=H12/WH H22=H22/WH H13=H13/WH H23=H23/WH H33=H33/WH WH=SQRT(WH) H4=H4/WH H5=H5/WH H6=H6/WH H14=H14/WH H24=H24/WH H34=H34/WH H15=H15/WH H25=H25/WH H35=H35/WH H16=H16/WH H26=H26/WH H36=H36/WH END IF C RETURN END C C======================================================================= C SUBROUTINE HDERC(ICB,A,Q,P1,P2,P3,HPR,HPI,HSR,HSI,HR,HI, * H1R,H1I,H2R,H2I,H3R,H3I,H4R,H4I,H5R,H5I,H6R,H6I, * H11R,H11I,H12R,H12I,H22R,H22I,H13R,H13I, * H23R,H23I,H33R,H33I,H14R,H14I,H24R,H24I, * H34R,H34I,H44R,H44I,H15R,H15I,H25R,H25I, * H35R,H35I,H45R,H45I,H55R,H55I,H16R,H16I, * H26R,H26I,H36R,H36I,H46R,H46I,H56R,H56I, * H66R,H66I,E) INTEGER ICB REAL A(10,21),Q(10,21),P1,P2,P3,HPR,HPI,HSR,HSI,HR,HI REAL H1R,H1I,H2R,H2I,H3R,H3I,H4R,H4I,H5R,H5I,H6R,H6I REAL H11R,H11I,H12R,H12I,H22R,H22I,H13R,H13I REAL H23R,H23I,H33R,H33I,H14R,H14I,H24R,H24I REAL H34R,H34I,H44R,H44I,H15R,H15I,H25R,H25I REAL H35R,H35I,H45R,H45I,H55R,H55I,H16R,H16I REAL H26R,H26I,H36R,H36I,H46R,H46I,H56R,H56I REAL H66R,H66I,E(9) C C Subroutine designed to calculate the phase-space derivatives of the C complex-valued reference Hamiltonian function in anisotropic C viscoelastic media. C C Input: C ICB... Wave type: C ICB.GT.0: P wave, C ICB.LT.0: S wave. C A... Real parts of the density-reduced viscoelastic moduli C and their first and second 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 Q... Imaginary parts of the density-reduced viscoelastic moduli C and their first and second spatial derivatives, C ordered analogoustly to the real parts. C P1,P2,P3... Real-valued reference slowness vector. C C Output: C HPR,HPI,HSR,HSI... Complex values of the Hamiltonian function C corresponding to P and S waves, respectively. C HR,HI.. Real and imaginary parts of the value of the Hamiltonian C function corresponding to the given wave. C H1R,H1I,H2R,H2I,H3R,H3I,H4R,H4I,H5R,H5I,H6R,H6I... First-order C phase-space derivatives of the Hamiltonian function, C calculated for the given slowness vector. C H11R,H11I,H12R,H12I,H22R,H22I,H13R,H13I,H23R,H23I,H33R,H33I, C H14R,H14I,H24R,H24I,H34R,H34I,H44R,H44I,H15R,H15I,H25R,H25I, C H35R,H35I,H45R,H45I,H55R,H55I,H16R,H16I,H26R,H26I,H36R,H36I, C H46R,H46I,H56R,H56I,H66R,H66I... Second-order partial phase-space C derivatives of the Hamiltonian function, C calculated for the given slowness vector. C E... Components of the reference polarization vectors: C E(1:3),E(4:6): S-wave polarization vectors. C E(7:9): P-wave polarization vector. C C Input SEP parameters: C DEGREE=real... Degree of the homogeneous Hamiltonian function to C be arithmetically averaged for common S-wave ray tracing. C DEGREE must equal 2, 1, -1 or -2. 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 If specified, there must be KSWAVE=0 (common S-wave rays). C Default: KSWAVE=0 C C Subroutines and external functions required: EXTERNAL CHR0,CHR1,CHR11,CHR2,SQRTC,CUBIC,ERROR,RSEP3I,RSEP3R C Subroutine CHR0 and its entries CHR1,CHR11,CHR2... This file. C SQRTC.. This file. C CUBIC.. File 'cubic.for'. C ERROR.. File 'error.for'. C RSEP3I,RSEP3R... File 'sep.for'. C C References: C Klimes, L. (2022a): C Tracing real-valued reference rays in anisotropic viscoelastic C media. Stud. geophys. geod., 66, 124-144. C Klimes, L. (2022b): C S-wave polarization vectors in anisotropic viscoelastic media. C Seismic Waves in Complex 3-D Structures, 31, 69-87. C C Date: 2022, November 26 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Temporary storage locations: C C Christoffel matrix REAL CH11R,CH12R,CH22R,CH13R,CH23R,CH33R,TRCHR REAL CH11I,CH12I,CH22I,CH13I,CH23I,CH33I,TRCHI C Traceless S-wave matrix (Klimes, 2022b, eq.11) REAL DCH11R,DCH12R,DCH22R,DCH13R,DCH23R,DCH33R REAL DCH11I,DCH12I,DCH22I,DCH13I,DCH23I,DCH33I C C Characteristic equation REAL AR,AI,BR,BI,CR,CI C C Eigenvalues of the Christoffel matrix REAL G1R,G1I,G2R,G2I,G3R,G3I,G1RNEW,G1INEW,G2RNEW,G2INEW C C Numerator and denominator in fractions REAL ENUMR,ENUMI,DENOR,DENOI,DD C C Tensor square of the P-wave eigenvector REAL G11R,G11I,G12R,G12I,G22R,G22I,G13R,G13I,G23R,G23I,G33R,G33I C C Projection matrices in the direction of the P-wave eigenvector REAL E11R,E11I,E12R,E12I,E22R,E22I,E13R,E13I,E23R,E23I,E33R,E33I REAL F11R,F11I,F12R,F12I,F22R,F22I,F13R,F13I,F23R,F23I,F33R,F33I C Temporary factors of these matrices REAL ER,EI,FR,FI C C P-wave eigenvector REAL R13,R23,R33,Y13,Y23,Y33 C S-wave eigenvectors REAL R11,R21,R31,Y11,Y21,Y31,R12,R22,R32,Y12,Y22,Y32 REAL R1R1,Y1Y1,R2R2,Y2Y2,RR11,RR12,RR22,SR11,SR12,SR22,DETSR C Eigenvector REAL R1,R2,R3,Y1,Y2,Y3 LOGICAL LONE,LTWO C LONE=.TRUE. ... Singularity with one nearly real-valued S-wave C eigenvector. Its imaginary part cannot be determined. C LTWO=.TRUE. ... Singularity with two S-wave eigenvectors. C The eigenvectors cannot be determined. C C Auxiliary storage locations REAL AUXR,AUXI,AUXXR,AUXXI,AUX1,AUX2,AUX3 C C Calculating the derivatives of the Hamiltonian function REAL C11R,C11I,C12R,C12I,C22R,C22I,C13R,C13I,C23R,C23I,C33R,C33I REAL D11R,D11I,D12R,D12I,D22R,D22I,D13R,D13I,D23R,D23I,D33R,D33I C Linear combinations of eigenvalues REAL G2G1R,G2G1I,G3G1R,G3G1I,G3G2R,G3G2I,G1G2R,G1G2I C C Factors C0, C1, C2, C3, C4 and C REAL C0R,C0I,C1R,C1I,C2R,C2I,C3R,C3I,C4R,C4I C Algebraic expressions containg S-wave eigenvalues REAL SUMR,SUMI,PRDR,PRDI,G1G1R,G1G1I,G2G2R,G2G2I,SUMGGR,SUMGGI C Algebraic expressions containg square roots of S-wave eigenvalues REAL G1SR,G1SI,G2SR,G2SI REAL SUMSR,SUMSI,PRDSR,PRDSI,PSUSR,PSUSI REAL G1S3R,G1S3I,G2S3R,G2S3I,SUMS3R,SUMS3I C Square roots of P-wave eigenvalue REAL G3SR,G3SI C Algebraic expressions containg fourth roots of S-wave eigenvalues REAL G1SSR,G1SSI,G2SSR,G2SSI REAL SUMSSR,SUMSSI,PRDSSR,PRDSSI,PSUSSR,PSUSSI REAL G1SS3R,G1SS3I,G2SS3R,G2SS3I,SUSS3R,SUSS3I C C Tensor square of the first derivatives of the Hamiltonian function REAL H1H1R,H1H1I,H1H2R,H1H2I,H2H2R,H2H2I,H1H3R,H1H3I REAL H2H3R,H2H3I,H3H3R,H3H3I,H1H4R,H1H4I,H2H4R,H2H4I REAL H3H4R,H3H4I,H4H4R,H4H4I,H1H5R,H1H5I,H2H5R,H2H5I REAL H3H5R,H3H5I,H4H5R,H4H5I,H5H5R,H5H5I,H1H6R,H1H6I REAL H2H6R,H2H6I,H3H6R,H3H6I,H4H6R,H4H6I,H5H6R,H5H6I REAL H6H6R,H6H6I C C Degree of homogeneous Hamiltonian function and other input data: INTEGER KSWAVE REAL DEGREE SAVE DEGREE DATA DEGREE/0./ C C....................................................................... C C Input SEP parameters IF(DEGREE.EQ.0.) THEN CALL RSEP3R('DEGREE',DEGREE,-1.) CALL RSEP3I('KSWAVE',KSWAVE,0) IF(KSWAVE.NE.0) THEN C 5A9 CALL ERROR('5A9 in HDERC: Common S-wave rays allowed only') END IF END IF C C Christoffel matrix (Klimes, 2022a, eq.5): CALL CHR0(A,Q,P1,P2,P3, * CH11R,CH11I,CH12R,CH12I,CH22R,CH22I, * CH13R,CH13I,CH23R,CH23I,CH33R,CH33I) C C Characteristic equation (Klimes, 2022a, eq.6): C Part of the matrix of cofactors of the Christoffel matrix C13R=CH12R*CH23R-CH13R*CH22R-CH12I*CH23I+CH13I*CH22I C13I=CH12R*CH23I-CH13R*CH22I+CH12I*CH23R-CH13I*CH22R C23R=CH12R*CH13R-CH11R*CH23R-CH12I*CH13I+CH11I*CH23I C23I=CH12R*CH13I-CH11R*CH23I+CH12I*CH13R-CH11I*CH23R C33R=CH11R*CH22R-CH12R*CH12R-CH11I*CH22I+CH12I*CH12I C33I=CH11R*CH22I-CH12R*CH12I+CH11I*CH22R-CH12I*CH12R C Trace of the Christoffel matrix TRCHR=CH11R+CH22R+CH33R TRCHI=CH11I+CH22I+CH33I IF(TRCHR.LE.0.) THEN C 5A11 CALL ERROR('5A11 in HDERC: Indefinite Christoffel matrix') C Real part of the trace of the Christoffel matrix is not C positive. C Check the stability conditions of the stiffness matrix. END IF C Minus trace of the Christoffel matrix AR=-TRCHR AI=-TRCHI C Trace of the matrix of cofactors of the Christoffel matrix C (Klimes, 2022a, eq.7) BR=C33R * +CH11R*CH33R-CH13R*CH13R-CH11I*CH33I+CH13I*CH13I * +CH22R*CH33R-CH23R*CH23R-CH22I*CH33I+CH23I*CH23I BI=C33I * +CH11R*CH33I-2.*CH13R*CH13I+CH11I*CH33R * +CH22R*CH33I-2.*CH23R*CH23I+CH22I*CH33R C Minus determinant of the Christoffel matrix CR=-CH13R*C13R-CH23R*C23R-CH33R*C33R * +CH13I*C13I+CH23I*C23I+CH33I*C33I CI=-CH13R*C13I-CH23R*C23I-CH33R*C33I * -CH13I*C13R-CH23I*C23R-CH33I*C33R C C Solving the characteristic equation: AUXR=TRCHR AR=AR/AUXR AI=AI/AUXR AUXR=AUXR*TRCHR BR=BR/AUXR BI=BI/AUXR AUXR=AUXR*TRCHR CR=CR/AUXR CI=CI/AUXR CALL CUBIC(AR,AI,BR,BI,CR,CI,G1R,G1I,G2R,G2I,G3R,G3I) C The numerical error of the P-wave eigenvalue corresponds C to the machine precission, whereas the numerical error C of the S-wave eigenvalues corresponds to the square root C of the machine precission (Klimes, 2022a, sec.2.2). G1R=G1R*TRCHR G1I=G1I*TRCHR G2R=G2R*TRCHR G2I=G2I*TRCHR G3R=G3R*TRCHR G3I=G3I*TRCHR IF(G3R.LE.0.) THEN C 5A12 CALL ERROR('5A12 in HDERC: Negative P-wave eigenvalue') C Real part of the P-wave eigenvalue of the Christoffel matrix is C not positive. C Check the stability conditions of the stiffness matrix. END IF IF(G3I.GT.0.) THEN IF(G3I.LT.0.000001*G3R) THEN G3I=0. ELSE C 5A13 CALL ERROR('5A13 in HDERC: Negative P-wave attenuation') C Imaginary part of the P-wave eigenvalue of the Christoffel C matrix is positive. C Check the stability conditions of the stiffness matrix. END IF END IF C C Tensor square of the P-wave eigenvector: C Modifying the Christoffel matrix for the P-wave C11R=CH11R-G3R C11I=CH11I-G3I C12R=CH12R C12I=CH12I C13R=CH13R C13I=CH13I C22R=CH22R-G3R C22I=CH22I-G3I C23R=CH23R C23I=CH23I C33R=CH33R-G3R C33I=CH33I-G3I C Matrix of cofactors (Klimes, 2022a, eq.35; Klimes, 2022b, eq.7) D11R=C22R*C33R-C23R*C23R-C22I*C33I+C23I*C23I D11I=C22R*C33I-C23R*C23I+C22I*C33R-C23I*C23R D12R=C13R*C23R-C12R*C33R-C13I*C23I+C12I*C33I D12I=C13R*C23I-C12R*C33I+C13I*C23R-C12I*C33R D22R=C11R*C33R-C13R*C13R-C11I*C33I+C13I*C13I D22I=C11R*C33I-C13R*C13I+C11I*C33R-C13I*C13R D13R=C12R*C23R-C13R*C22R-C12I*C23I+C13I*C22I D13I=C12R*C23I-C13R*C22I+C12I*C23R-C13I*C22R D23R=C12R*C13R-C11R*C23R-C12I*C13I+C11I*C23I D23I=C12R*C13I-C11R*C23I+C12I*C13R-C11I*C23R D33R=C11R*C22R-C12R*C12R-C11I*C22I+C12I*C12I D33I=C11R*C22I-C12R*C12I+C11I*C22R-C12I*C12R C Trace of the matrix of cofactors DENOR=D11R+D22R+D33R DENOI=D11I+D22I+D33I DD=DENOR*DENOR+DENOI*DENOI DENOR= DENOR/DD DENOI=-DENOI/DD C Tensor square of the P-wave eigenvector C (Klimes, 2022a, eq.36; 2022b, eq.8) G11R=D11R*DENOR-D11I*DENOI G11I=D11R*DENOI+D11I*DENOR G12R=D12R*DENOR-D12I*DENOI G12I=D12R*DENOI+D12I*DENOR G13R=D13R*DENOR-D13I*DENOI G13I=D13R*DENOI+D13I*DENOR G22R=D22R*DENOR-D22I*DENOI G22I=D22R*DENOI+D22I*DENOR G23R=D23R*DENOR-D23I*DENOI G23I=D23R*DENOI+D23I*DENOR G33R=D33R*DENOR-D33I*DENOI G33I=D33R*DENOI+D33I*DENOR C C Traceless S-wave matrix (Klimes, 2022b, eq.11): AUXR=0.5*(3.*G3R-TRCHR) AUXI=0.5*(3.*G3I-TRCHI) DCH11R=CH11R-G11R*AUXR+G11I*AUXI DCH11I=CH11I-G11R*AUXI-G11I*AUXR DCH12R=CH12R-G12R*AUXR+G12I*AUXI DCH12I=CH12I-G12R*AUXI-G12I*AUXR DCH13R=CH13R-G13R*AUXR+G13I*AUXI DCH13I=CH13I-G13R*AUXI-G13I*AUXR DCH22R=CH22R-G22R*AUXR+G22I*AUXI DCH22I=CH22I-G22R*AUXI-G22I*AUXR DCH23R=CH23R-G23R*AUXR+G23I*AUXI DCH23I=CH23I-G23R*AUXI-G23I*AUXR DCH33R=CH33R-G33R*AUXR+G33I*AUXI DCH33I=CH33I-G33R*AUXI-G33I*AUXR AUXR=(DCH11R+DCH22R+DCH33R)/3. AUXI=(DCH11I+DCH22I+DCH33I)/3. DCH11R=DCH11R-AUXR DCH11I=DCH11I-AUXI DCH22R=DCH22R-AUXR DCH22I=DCH22I-AUXI DCH33R=DCH33R-AUXR DCH33I=DCH33I-AUXI C Trace of the matrix of cofactors (equal to -(G2-G1)**2/4) AUXR=DCH11R*DCH22R-DCH12R*DCH12R-DCH11I*DCH22I+DCH12I*DCH12I * +DCH11R*DCH33R-DCH13R*DCH13R-DCH11I*DCH33I+DCH13I*DCH13I * +DCH22R*DCH33R-DCH23R*DCH23R-DCH22I*DCH33I+DCH23I*DCH23I AUXI=DCH11R*DCH22I-DCH12R*DCH12I+DCH11I*DCH22R-DCH12I*DCH12R * +DCH11R*DCH33I-2.*DCH13R*DCH13I+DCH11I*DCH33R * +DCH22R*DCH33I-2.*DCH23R*DCH23I+DCH22I*DCH33R C Difference (G2-G1)/2 of S-wave eigenvalues CALL SQRTC(-AUXR,-AUXI,G2G1R,G2G1I) C Average S-wave eigenvalue AUXR=0.5*(TRCHR-G3R) AUXI=0.5*(TRCHI-G3I) C Improving S-wave eigenvalues G1RNEW=AUXR-G2G1R G1INEW=AUXI-G2G1I G2RNEW=AUXR+G2G1R G2INEW=AUXI+G2G1I IF((G1RNEW-G1R)**2+(G1INEW-G1I)**2 * +(G2RNEW-G2R)**2+(G2INEW-G2I)**2.GT.2.*(0.001*G3R)**2) THEN C 5A10 CALL ERROR('5A10 in HDERC: Inaccurate S-wave eigenvalues') C This error should not appear, contact the authors. END IF G1R=G1RNEW G1I=G1INEW G2R=G2RNEW G2I=G2INEW C The numerical error of the S-wave eigenvalues G1 and G2 should now C correspond to the machine precission. C C P-wave and common S-wave Hamiltonian functions C (Klimes, 2022a, eqs.23,26) IF(DEGREE.EQ.2.) THEN HPR=0.5*G3R HPI=0.5*G3I HSR=0.25*(G1R+G2R) HSI=0.25*(G1I+G2I) ELSE IF(DEGREE.EQ.1.) THEN CALL SQRTC(G3R,G3I,HPR,HPI) CALL SQRTC(G1R,G1I,G1SR,G1SI) CALL SQRTC(G2R,G2I,G2SR,G2SI) HSR=0.5*(G1SR+G2SR) HSI=0.5*(G1SI+G2SI) ELSE IF(DEGREE.EQ.-1.) THEN CALL SQRTC(G3R,G3I,G3SR,G3SI) CALL SQRTC(G1R,G1I,G1SR,G1SI) CALL SQRTC(G2R,G2I,G2SR,G2SI) DD=G3SR*G3SR+G3SI*G3SI HPR=-G3SR/DD HPI= G3SI/DD DD=G1SR*G1SR+G1SI*G1SI HSR= G1SR/DD HSI=-G1SI/DD DD=G2SR*G2SR+G2SI*G2SI HSR=-0.5*(HSR+G2SR/DD) HSI=-0.5*(HSI-G2SI/DD) ELSE IF(DEGREE.EQ.-2.) THEN DD=-2.*(G3R*G3R+G3I*G3I) HPR= G3R/DD HPI=-G3I/DD DD=G1R*G1R+G1I*G1I HSR= G1R/DD HSI=-G1I/DD DD=G2R*G2R+G2I*G2I HSR=-0.25*(HSR+G2R/DD) HSI=-0.25*(HSI-G2I/DD) END IF C C Projection matrix in the direction of the P-wave eigenvector C (Klimes, 2022a, eq.38) E11R=1.-G11R E12R= -G12R E22R=1.-G22R E13R= -G13R E23R= -G23R E33R=1.-G33R E11I= -G11I E12I= -G12I E22I= -G22I E13I= -G13I E23I= -G23I E33I= -G33I C C Projection of the Christoffel matrix (Klimes, 2022a, eq.39) F11R=CH11R-G11R*G3R+G11I*G3I F12R=CH12R-G12R*G3R+G12I*G3I F22R=CH22R-G22R*G3R+G22I*G3I F13R=CH13R-G13R*G3R+G13I*G3I F23R=CH23R-G23R*G3R+G23I*G3I F33R=CH33R-G33R*G3R+G33I*G3I F11I=CH11I-G11R*G3I-G11I*G3R F12I=CH12I-G12R*G3I-G12I*G3R F22I=CH22I-G22R*G3I-G22I*G3R F13I=CH13I-G13R*G3I-G13I*G3R F23I=CH23I-G23R*G3I-G23I*G3R F33I=CH33I-G33R*G3I-G33I*G3R C IF(ICB.GT.0) THEN C Reference P-wave rays C C Hamiltonian function HR=HPR HI=HPI C C Converting the eigenvalue derivatives to Hamiltonian derivatives C (Klimes, 2022a, eqs.24-25) AUXXR=0.5*DEGREE/(G3R*G3R+G3I*G3I) AUXR=( HR*G3R+HI*G3I)*AUXXR AUXI=(-HR*G3I+HI*G3R)*AUXXR C11R=G11R*AUXR-G11I*AUXI C12R=G12R*AUXR-G12I*AUXI C22R=G22R*AUXR-G22I*AUXI C13R=G13R*AUXR-G13I*AUXI C23R=G23R*AUXR-G23I*AUXI C33R=G33R*AUXR-G33I*AUXI C11I=G11R*AUXI+G11I*AUXR C12I=G12R*AUXI+G12I*AUXR C22I=G22R*AUXI+G22I*AUXR C13I=G13R*AUXI+G13I*AUXR C23I=G23R*AUXI+G23I*AUXR C33I=G33R*AUXI+G33I*AUXR C C First-order phase-space derivatives of the Hamiltonian function C (Klimes, 2022a, eqs.24,40) CALL CHR1(A,Q,P1,P2,P3, * C11R,C11I,C12R,C12I,C22R,C22I,C13R,C13I,C23R,C23I,C33R,C33I, * H1R,H1I,H2R,H2I,H3R,H3I,H4R,H4I,H5R,H5I,H6R,H6I) C C Second-order phase-space derivatives of the Hamiltonian function C (Klimes, 2022a, eqs.25,41) C C Terms with the second-order derivatives of Christoffel matrix CALL CHR2(A,Q,P1,P2,P3, * C11R,C11I,C12R,C12I,C22R,C22I,C13R,C13I,C23R,C23I,C33R,C33I, * H11R,H11I,H12R,H12I,H22R,H22I,H13R,H13I,H23R,H23I,H33R,H33I, * H14R,H14I,H24R,H24I,H34R,H34I,H44R,H44I, * H15R,H15I,H25R,H25I,H35R,H35I,H45R,H45I,H55R,H55I, * H16R,H16I,H26R,H26I,H36R,H36I,H46R,H46I,H56R,H56I,H66R,H66I) C C Terms with the first-order derivatives of Christoffel matrix C Matrix Dij=2*Pij=2*(Eij*(G3-G1-G2)+Fij)/((G3-G1)*(G3-G2)) C (Klimes, 2022a, eq.42) G3G1R=G3R-G1R G3G1I=G3I-G1I G3G2R=G3R-G2R G3G2I=G3I-G2I G1G2R=G3G1R-G2R G1G2I=G3G1I-G2I DENOR=G3G1R*G3G2R-G3G1I*G3G2I DENOI=G3G1R*G3G2I+G3G1I*G3G2R DD=(DENOR*DENOR+DENOI*DENOI)/2. FR= DENOR/DD FI=-DENOI/DD ER=G1G2R*FR-G1G2I*FI EI=G1G2R*FI+G1G2I*FR D11R=E11R*ER-E11I*EI+F11R*FR-F11I*FI D12R=E12R*ER-E12I*EI+F12R*FR-F12I*FI D22R=E22R*ER-E22I*EI+F22R*FR-F22I*FI D13R=E13R*ER-E13I*EI+F13R*FR-F13I*FI D23R=E23R*ER-E23I*EI+F23R*FR-F23I*FI D33R=E33R*ER-E33I*EI+F33R*FR-F33I*FI D11I=E11R*EI+E11I*ER+F11R*FI+F11I*FR D12I=E12R*EI+E12I*ER+F12R*FI+F12I*FR D22I=E22R*EI+E22I*ER+F22R*FI+F22I*FR D13I=E13R*EI+E13I*ER+F13R*FI+F13I*FR D23I=E23R*EI+E23I*ER+F23R*FI+F23I*FR D33I=E33R*EI+E33I*ER+F33R*FI+F33I*FR CALL CHR11( * C11R,C11I,C12R,C12I,C22R,C22I,C13R,C13I,C23R,C23I,C33R,C33I, * D11R,D11I,D12R,D12I,D22R,D22I,D13R,D13I,D23R,D23I,D33R,D33I, * H11R,H11I,H12R,H12I,H22R,H22I,H13R,H13I,H23R,H23I,H33R,H33I, * H14R,H14I,H24R,H24I,H34R,H34I,H44R,H44I, * H15R,H15I,H25R,H25I,H35R,H35I,H45R,H45I,H55R,H55I, * H16R,H16I,H26R,H26I,H36R,H36I,H46R,H46I,H56R,H56I,H66R,H66I) C IF(DEGREE.NE.2.) THEN C Tensor square of the first derivatives of Hamiltonian function H1H1R=H1R*H1R-H1I*H1I H1H2R=H1R*H2R-H1I*H2I H2H2R=H2R*H2R-H2I*H2I H1H3R=H1R*H3R-H1I*H3I H2H3R=H2R*H3R-H2I*H3I H3H3R=H3R*H3R-H3I*H3I H1H4R=H1R*H4R-H1I*H4I H2H4R=H2R*H4R-H2I*H4I H3H4R=H3R*H4R-H3I*H4I H4H4R=H4R*H4R-H4I*H4I H1H5R=H1R*H5R-H1I*H5I H2H5R=H2R*H5R-H2I*H5I H3H5R=H3R*H5R-H3I*H5I H4H5R=H4R*H5R-H4I*H5I H5H5R=H5R*H5R-H5I*H5I H1H6R=H1R*H6R-H1I*H6I H2H6R=H2R*H6R-H2I*H6I H3H6R=H3R*H6R-H3I*H6I H4H6R=H4R*H6R-H4I*H6I H5H6R=H5R*H6R-H5I*H6I H6H6R=H6R*H6R-H6I*H6I H1H1I=H1R*H1I*2. H1H2I=H1R*H2I+H1I*H2R H2H2I=H2R*H2I*2. H1H3I=H1R*H3I+H1I*H3R H2H3I=H2R*H3I+H2I*H3R H3H3I=H3R*H3I*2. H1H4I=H1R*H4I+H1I*H4R H2H4I=H2R*H4I+H2I*H4R H3H4I=H3R*H4I+H3I*H4R H4H4I=H4R*H4I*2. H1H5I=H1R*H5I+H1I*H5R H2H5I=H2R*H5I+H2I*H5R H3H5I=H3R*H5I+H3I*H5R H4H5I=H4R*H5I+H4I*H5R H5H5I=H5R*H5I*2. H1H6I=H1R*H6I+H1I*H6R H2H6I=H2R*H6I+H2I*H6R H3H6I=H3R*H6I+H3I*H6R H4H6I=H4R*H6I+H4I*H6R H5H6I=H5R*H6I+H5I*H6R H6H6I=H6R*H6I*2. C C Terms with the first-order derivatives of Hamiltonian function C (Klimes, 2022a, eq.25 with eqs.23-24) AUXXR=(DEGREE-2.)/(DEGREE*(HR*HR+HI*HI)) AUXR= HR*AUXXR AUXI=-HI*AUXXR H11R=H11R+H1H1R*AUXR-H1H1I*AUXI H12R=H12R+H1H2R*AUXR-H1H2I*AUXI H22R=H22R+H2H2R*AUXR-H2H2I*AUXI H13R=H13R+H1H3R*AUXR-H1H3I*AUXI H23R=H23R+H2H3R*AUXR-H2H3I*AUXI H33R=H33R+H3H3R*AUXR-H3H3I*AUXI H14R=H14R+H1H4R*AUXR-H1H4I*AUXI H24R=H24R+H2H4R*AUXR-H2H4I*AUXI H34R=H34R+H3H4R*AUXR-H3H4I*AUXI H44R=H44R+H4H4R*AUXR-H4H4I*AUXI H15R=H15R+H1H5R*AUXR-H1H5I*AUXI H25R=H25R+H2H5R*AUXR-H2H5I*AUXI H35R=H35R+H3H5R*AUXR-H3H5I*AUXI H45R=H45R+H4H5R*AUXR-H4H5I*AUXI H55R=H55R+H5H5R*AUXR-H5H5I*AUXI H16R=H16R+H1H6R*AUXR-H1H6I*AUXI H26R=H26R+H2H6R*AUXR-H2H6I*AUXI H36R=H36R+H3H6R*AUXR-H3H6I*AUXI H46R=H46R+H4H6R*AUXR-H4H6I*AUXI H56R=H56R+H5H6R*AUXR-H5H6I*AUXI H66R=H66R+H6H6R*AUXR-H6H6I*AUXI H11I=H11I+H1H1R*AUXI+H1H1I*AUXR H12I=H12I+H1H2R*AUXI+H1H2I*AUXR H22I=H22I+H2H2R*AUXI+H2H2I*AUXR H13I=H13I+H1H3R*AUXI+H1H3I*AUXR H23I=H23I+H2H3R*AUXI+H2H3I*AUXR H33I=H33I+H3H3R*AUXI+H3H3I*AUXR H14I=H14I+H1H4R*AUXI+H1H4I*AUXR H24I=H24I+H2H4R*AUXI+H2H4I*AUXR H34I=H34I+H3H4R*AUXI+H3H4I*AUXR H44I=H44I+H4H4R*AUXI+H4H4I*AUXR H15I=H15I+H1H5R*AUXI+H1H5I*AUXR H25I=H25I+H2H5R*AUXI+H2H5I*AUXR H35I=H35I+H3H5R*AUXI+H3H5I*AUXR H45I=H45I+H4H5R*AUXI+H4H5I*AUXR H55I=H55I+H5H5R*AUXI+H5H5I*AUXR H16I=H16I+H1H6R*AUXI+H1H6I*AUXR H26I=H26I+H2H6R*AUXI+H2H6I*AUXR H36I=H36I+H3H6R*AUXI+H3H6I*AUXR H46I=H46I+H4H6R*AUXI+H4H6I*AUXR H56I=H56I+H5H6R*AUXI+H5H6I*AUXR H66I=H66I+H6H6R*AUXI+H6H6I*AUXR END IF C ELSE IF(ICB.LT.0) THEN C Reference common S-wave rays C C Hamiltonian function HR=HSR HI=HSI C C Factors C0, C1, C2, C3, C4 and C IF(DEGREE.EQ.2.) THEN C (Klimes, 2022a, sec.4.2.1) C0R= 0. C0I= 0. C1R=-1. C1I= 0. C2R=-G1R-G2R C2I=-G1I-G2I C3R=0. C3I=0. C4R=0. C4I=0. CR=0. CI=0. ELSE IF(DEGREE.EQ.1.) THEN C (Klimes, 2022a, sec.4.2.2) SUMR=G1R+G2R SUMI=G1I+G2I PRDR=G1R*G2R-G1I*G2I PRDI=G1R*G2I+G1I*G2R CALL SQRTC(G1R,G1I,G1SR,G1SI) CALL SQRTC(G2R,G2I,G2SR,G2SI) SUMSR=G1SR+G2SR SUMSI=G1SI+G2SI PRDSR=G1SR*G2SR-G1SI*G2SI PRDSI=G1SR*G2SI+G1SI*G2SR PSUSR=PRDSR*SUMSR-PRDSI*SUMSI PSUSI=PRDSR*SUMSI+PRDSI*SUMSR G1S3R=G1SR*G1R-G1SI*G1I G1S3I=G1SR*G1I+G1SI*G1R G2S3R=G2SR*G2R-G2SI*G2I G2S3I=G2SR*G2I+G2SI*G2R SUMS3R=G1S3R+G2S3R SUMS3I=G1S3I+G2S3I CALL SQRTC(G1SR,G1SI,G1SSR,G1SSI) CALL SQRTC(G2SR,G2SI,G2SSR,G2SSI) SUMSSR=G1SSR+G2SSR SUMSSI=G1SSI+G2SSI PRDSSR=G1SSR*G2SSR-G1SSI*G2SSI PRDSSI=G1SSR*G2SSI+G1SSI*G2SSR PSUSSR=PRDSSR*SUMSSR-PRDSSI*SUMSSI PSUSSI=PRDSSR*SUMSSI+PRDSSI*SUMSSR DENOR=PSUSR DENOI=PSUSI DD=DENOR*DENOR+DENOI*DENOI DENOR= DENOR/DD DENOI=-DENOI/DD C0R=-DENOR C0I=-DENOI ENUMR=G1R+PRDSR+G2R ENUMI=G1I+PRDSI+G2I C1R=-ENUMR*DENOR+ENUMI*DENOI C1I=-ENUMR*DENOI-ENUMI*DENOR ENUMR=SUMSR*SUMS3R-SUMSI*SUMS3I+PRDR ENUMI=SUMSR*SUMS3I+SUMSI*SUMS3R+PRDI C2R=-ENUMR*DENOR+ENUMI*DENOI C2I=-ENUMR*DENOI-ENUMI*DENOR ENUMR=SUMSR+PRDSSR ENUMI=SUMSI+PRDSSI DENOR=PSUSSR*PSUSR-PSUSSI*PSUSI DENOI=PSUSSR*PSUSI+PSUSSI*PSUSR DD=DENOR*DENOR+DENOI*DENOI DENOR= DENOR/DD DENOI=-DENOI/DD C3R=-ENUMR*DENOR+ENUMI*DENOI C3I=-ENUMR*DENOI-ENUMI*DENOR AUXR=ENUMR*SUMR-ENUMI*SUMI AUXI=ENUMR*SUMI+ENUMI*SUMR ENUMR=AUXR+PRDSR*PRDSSR-PRDSI*PRDSSI ENUMI=AUXI+PRDSR*PRDSSI+PRDSI*PRDSSR C4R=-ENUMR*DENOR+ENUMI*DENOI C4I=-ENUMR*DENOI-ENUMI*DENOR AUXR=SUMSR*SUMSR-SUMSI*SUMSI AUXI=SUMSR*SUMSI+SUMSI*SUMSR AUXXR=AUXR*SUMSSR-AUXI*SUMSSI AUXXI=AUXR*SUMSSI+AUXI*SUMSSR DD=(AUXXR*AUXXR+AUXXI*AUXXI)*2. AUXXR=AUXXR/DD AUXXI=AUXXI/DD CR=DENOR*AUXXR-DENOI*AUXXI CI=DENOR*AUXXI+DENOI*AUXXR ELSE IF(DEGREE.EQ.-1.) THEN C (Klimes, 2022a, sec.4.2.3) SUMR=G1R+G2R SUMI=G1I+G2I PRDR=G1R*G2R-G1I*G2I PRDI=G1R*G2I+G1I*G2R G1G1R=(G1R+G1I)*(G1R-G1I) G1G1I=G1R*G1I*2. G2G2R=(G2R+G2I)*(G2R-G2I) G2G2I=G2R*G2I*2. SUMGGR=G1G1R+G2G2R SUMGGI=G1G1I+G2G2I CALL SQRTC(G1R,G1I,G1SR,G1SI) CALL SQRTC(G2R,G2I,G2SR,G2SI) SUMSR=G1SR+G2SR SUMSI=G1SI+G2SI PRDSR=G1SR*G2SR-G1SI*G2SI PRDSI=G1SR*G2SI+G1SI*G2SR PSUSR=PRDSR*SUMSR-PRDSI*SUMSI PSUSI=PRDSR*SUMSI+PRDSI*SUMSR G1S3R=G1SR*G1R-G1SI*G1I G1S3I=G1SR*G1I+G1SI*G1R G2S3R=G2SR*G2R-G2SI*G2I G2S3I=G2SR*G2I+G2SI*G2R SUMS3R=G1S3R+G2S3R SUMS3I=G1S3I+G2S3I CALL SQRTC(G1SR,G1SI,G1SSR,G1SSI) CALL SQRTC(G2SR,G2SI,G2SSR,G2SSI) SUMSSR=G1SSR+G2SSR SUMSSI=G1SSI+G2SSI PRDSSR=G1SSR*G2SSR-G1SSI*G2SSI PRDSSI=G1SSR*G2SSI+G1SSI*G2SSR PSUSSR=PRDSSR*SUMSSR-PRDSSI*SUMSSI PSUSSI=PRDSSR*SUMSSI+PRDSSI*SUMSSR G1SS3R=G1SSR*G1SR-G1SSI*G1SI G1SS3I=G1SSR*G1SI+G1SSI*G1SR G2SS3R=G2SSR*G2SR-G2SSI*G2SI G2SS3I=G2SSR*G2SI+G2SSI*G2SR SUSS3R=G1SS3R+G2SS3R SUSS3I=G1SS3I+G2SS3I DENOR=PRDR*PSUSR-PRDI*PSUSI DENOI=PRDR*PSUSI+PRDI*PSUSR DD=DENOR*DENOR+DENOI*DENOI DENOR= DENOR/DD DENOI=-DENOI/DD ENUMR=SUMR+PRDSR ENUMI=SUMI+PRDSI C0R=-ENUMR*DENOR+ENUMI*DENOI C0I=-ENUMR*DENOI-ENUMI*DENOR ENUMR=SUMSR*SUMS3R-SUMSI*SUMS3I+PRDR ENUMI=SUMSR*SUMS3I+SUMSI*SUMS3R+PRDI C1R=-ENUMR*DENOR+ENUMI*DENOI C1I=-ENUMR*DENOI-ENUMI*DENOR AUXR=SUMR+PRDSR AUXI=SUMI+PRDSI ENUMR=AUXR*SUMGGR-AUXI*SUMGGI+G1S3R*G2S3R-G1S3I*G2S3I ENUMI=AUXR*SUMGGI+AUXI*SUMGGR+G1S3R*G2S3I+G1S3I*G2S3R C2R=-ENUMR*DENOR+ENUMI*DENOI C2I=-ENUMR*DENOI-ENUMI*DENOR AUXR=PSUSSR*SUMSR-PSUSSI*SUMSI AUXI=PSUSSR*SUMSI+PSUSSI*SUMSR DENOR=PRDR*AUXR-PRDI*AUXI DENOI=PRDR*AUXI+PRDI*AUXR DD=DENOR*DENOR+DENOI*DENOI DENOR= DENOR/DD DENOI=-DENOI/DD ENUMR=SUMSSR*SUSS3R-SUMSSI*SUSS3I+PRDSR ENUMI=SUMSSR*SUSS3I+SUMSSI*SUSS3R+PRDSI C3R=-ENUMR*DENOR+ENUMI*DENOI C3I=-ENUMR*DENOI-ENUMI*DENOR DD=AUXR*AUXR+AUXI*AUXI DENOR= AUXR/DD DENOI=-AUXI/DD ENUMR=G1R*G1SSR-G1I*G1SSI+G2R*G2SSR-G2I*G2SSI ENUMI=G1R*G1SSI+G1I*G1SSR+G2R*G2SSI+G2I*G2SSR AUXR=PRDR*PRDSSR-PRDI*PRDSSI AUXI=PRDR*PRDSSI+PRDI*PRDSSR DD=AUXR*AUXR+AUXI*AUXI AUXR= AUXR/DD AUXI=-AUXI/DD C4R=-ENUMR*AUXR+ENUMI*AUXI-DENOR C4I=-ENUMR*AUXI-ENUMI*AUXR-DENOI AUXR=PRDR*SUMSR-PRDI*SUMSI AUXI=PRDR*SUMSI+PRDI*SUMSR DD=AUXR*AUXR+AUXI*AUXI AUXR= AUXR/DD AUXI=-AUXI/DD AUXXR=(DENOR+DENOI)*(DENOR-DENOI) AUXXI=DENOR*DENOI*2. DENOR=AUXR*AUXXR-AUXI*AUXXI DENOI=AUXR*AUXXI+AUXI*AUXXR ENUMR=SUMSR+0.5*PRDSSR ENUMI=SUMSI+0.5*PRDSSI CR=ENUMR*DENOR-ENUMI*DENOI CI=ENUMR*DENOI+ENUMI*DENOR ELSE IF(DEGREE.EQ.-2.) THEN C (Klimes, 2022a, sec.4.2.4) SUMR=G1R+G2R SUMI=G1I+G2I G1G1R=(G1R+G1I)*(G1R-G1I) G1G1I=G1R*G1I*2. G2G2R=(G2R+G2I)*(G2R-G2I) G2G2I=G2R*G2I*2. SUMGGR=G1G1R+G2G2R SUMGGI=G1G1I+G2G2I PRDR=G1R*G2R-G1I*G2I PRDI=G1R*G2I+G1I*G2R CALL SQRTC(G1R,G1I,G1SR,G1SI) CALL SQRTC(G2R,G2I,G2SR,G2SI) SUMSR=G1SR+G2SR SUMSI=G1SI+G2SI PRDSR=G1SR*G2SR-G1SI*G2SI PRDSI=G1SR*G2SI+G1SI*G2SR PSUSR=PRDSR*SUMSR-PRDSI*SUMSI PSUSI=PRDSR*SUMSI+PRDSI*SUMSR G1S3R=G1SR*G1R-G1SI*G1I G1S3I=G1SR*G1I+G1SI*G1R G2S3R=G2SR*G2R-G2SI*G2I G2S3I=G2SR*G2I+G2SI*G2R SUMS3R=G1S3R+G2S3R SUMS3I=G1S3I+G2S3I DENOR=(PRDR+PRDI)*(PRDR-PRDI) DENOI=PRDR*PRDI*2. DD=DENOR*DENOR+DENOI*DENOI DENOR= DENOR/DD DENOI=-DENOI/DD ENUMR=SUMR ENUMI=SUMI C0R=-ENUMR*DENOR+ENUMI*DENOI C0I=-ENUMR*DENOI-ENUMI*DENOR ENUMR=G1G1R+PRDR+G2G2R ENUMI=G1G1I+PRDI+G2G2I C1R=-ENUMR*DENOR+ENUMI*DENOI C1I=-ENUMR*DENOI-ENUMI*DENOR ENUMR=SUMR*SUMGGR-SUMI*SUMGGI ENUMI=SUMR*SUMGGI+SUMI*SUMGGR C2R=-ENUMR*DENOR+ENUMI*DENOI C2I=-ENUMR*DENOI-ENUMI*DENOR DENOR=PRDR*PSUSR-PRDI*PSUSI DENOI=PRDR*PSUSI+PRDI*PSUSR DD=DENOR*DENOR+DENOI*DENOI DENOR= DENOR/DD DENOI=-DENOI/DD ENUMR=SUMR+PRDSR ENUMI=SUMI+PRDSI C3R=-ENUMR*DENOR+ENUMI*DENOI C3I=-ENUMR*DENOI-ENUMI*DENOR ENUMR=SUMSR*SUMS3R-SUMSI*SUMS3I+PRDR ENUMI=SUMSR*SUMS3I+SUMSI*SUMS3R+PRDI C4R=-ENUMR*DENOR+ENUMI*DENOI C4I=-ENUMR*DENOI-ENUMI*DENOR AUXR=PRDR*SUMSR-PRDI*SUMSI AUXI=PRDR*SUMSI+PRDI*SUMSR DENOR=(AUXR+AUXI)*(AUXR-AUXI) DENOI=AUXR*AUXI*2. DD=DENOR*DENOR+DENOI*DENOI CR= DENOR/DD CI=-DENOI/DD ELSE C 5A7 CALL ERROR('5A7 in HDERC: Wrong degree of Hamiltonian') C Only homogeneous Hamiltonian functions of degrees -2, -1, 1 C and 2 are now allowed. END IF C C One quarter of matrix Cij0 C (Klimes, 2022a, eq.52) C0R=0.25*C0R C0I=0.25*C0I AUXR=0.25*C1R AUXI=0.25*C1I C11R=F11R*C0R-F11I*C0I-E11R*AUXR+E11I*AUXI C11I=F11R*C0I+F11I*C0R-E11R*AUXI-E11I*AUXR C12R=F12R*C0R-F12I*C0I-E12R*AUXR+E12I*AUXI C12I=F12R*C0I+F12I*C0R-E12R*AUXI-E12I*AUXR C22R=F22R*C0R-F22I*C0I-E22R*AUXR+E22I*AUXI C22I=F22R*C0I+F22I*C0R-E22R*AUXI-E22I*AUXR C13R=F13R*C0R-F13I*C0I-E13R*AUXR+E13I*AUXI C13I=F13R*C0I+F13I*C0R-E13R*AUXI-E13I*AUXR C23R=F23R*C0R-F23I*C0I-E23R*AUXR+E23I*AUXI C23I=F23R*C0I+F23I*C0R-E23R*AUXI-E23I*AUXR C33R=F33R*C0R-F33I*C0I-E33R*AUXR+E33I*AUXI C33I=F33R*C0I+F33I*C0R-E33R*AUXI-E33I*AUXR C C First-order phase-space derivatives of the Hamiltonian function C (Klimes, 2022a, eq.63) CALL CHR1(A,Q,P1,P2,P3, * C11R,C11I,C12R,C12I,C22R,C22I,C13R,C13I,C23R,C23I,C33R,C33I, * H1R,H1I,H2R,H2I,H3R,H3I,H4R,H4I,H5R,H5I,H6R,H6I) C C Second-order phase-space derivatives of the Hamiltonian function C (Klimes, 2022a, eq.64) C C Terms with the second-order derivatives of Christoffel matrix CALL CHR2(A,Q,P1,P2,P3, * C11R,C11I,C12R,C12I,C22R,C22I,C13R,C13I,C23R,C23I,C33R,C33I, * H11R,H11I,H12R,H12I,H22R,H22I,H13R,H13I,H23R,H23I,H33R,H33I, * H14R,H14I,H24R,H24I,H34R,H34I,H44R,H44I, * H15R,H15I,H25R,H25I,H35R,H35I,H45R,H45I,H55R,H55I, * H16R,H16I,H26R,H26I,H36R,H36I,H46R,H46I,H56R,H56I,H66R,H66I) C C Terms with the first-order derivatives of Christoffel matrix C Matrix Dij=0.5*(Cij1-G3*Cij0)/((G3-G1)*(G3-G1)) G3G1R=G3R-G1R G3G1I=G3I-G1I G3G2R=G3R-G2R G3G2I=G3I-G2I DENOR=G3G1R*G3G2R-G3G1I*G3G2I DENOI=G3G1R*G3G2I+G3G1I*G3G2R DD=DENOR*DENOR+DENOI*DENOI DENOR= DENOR/DD DENOI=-DENOI/DD AUXR=0.5*DENOR AUXI=0.5*DENOI FR=C1R*AUXR-C1I*AUXI FI=C1R*AUXI+C1I*AUXR ER=C2R*AUXR-C2I*AUXI EI=C2R*AUXI+C2I*AUXR AUXR=2.*(DENOR*G3R-DENOI*G3I) AUXI=2.*(DENOR*G3I+DENOI*G3R) D11R=F11R*FR-F11I*FI-E11R*ER+E11I*EI-C11R*AUXR+C11I*AUXI D11I=F11R*FI+F11I*FR-E11R*EI-E11I*ER-C11R*AUXI-C11I*AUXR D12R=F12R*FR-F12I*FI-E12R*ER+E12I*EI-C12R*AUXR+C12I*AUXI D12I=F12R*FI+F12I*FR-E12R*EI-E12I*ER-C12R*AUXI-C12I*AUXR D22R=F22R*FR-F22I*FI-E22R*ER+E22I*EI-C22R*AUXR+C22I*AUXI D22I=F22R*FI+F22I*FR-E22R*EI-E22I*ER-C22R*AUXI-C22I*AUXR D13R=F13R*FR-F13I*FI-E13R*ER+E13I*EI-C13R*AUXR+C13I*AUXI D13I=F13R*FI+F13I*FR-E13R*EI-E13I*ER-C13R*AUXI-C13I*AUXR D23R=F23R*FR-F23I*FI-E23R*ER+E23I*EI-C23R*AUXR+C23I*AUXI D23I=F23R*FI+F23I*FR-E23R*EI-E23I*ER-C23R*AUXI-C23I*AUXR D33R=F33R*FR-F33I*FI-E33R*ER+E33I*EI-C33R*AUXR+C33I*AUXI D33I=F33R*FI+F33I*FR-E33R*EI-E33I*ER-C33R*AUXI-C33I*AUXR CALL CHR11( * D11R,D11I,D12R,D12I,D22R,D22I,D13R,D13I,D23R,D23I,D33R,D33I, * G11R,G11I,G12R,G12I,G22R,G22I,G13R,G13I,G23R,G23I,G33R,G33I, * H11R,H11I,H12R,H12I,H22R,H22I,H13R,H13I,H23R,H23I,H33R,H33I, * H14R,H14I,H24R,H24I,H34R,H34I,H44R,H44I, * H15R,H15I,H25R,H25I,H35R,H35I,H45R,H45I,H55R,H55I, * H16R,H16I,H26R,H26I,H36R,H36I,H46R,H46I,H56R,H56I,H66R,H66I) IF(DEGREE.NE.2.) THEN C Matrix Cij=0.5*C*(Fij-G1*Eij) FR=0.5*CR FI=0.5*CI ER=FR*G1R-FI*G1I EI=FR*G1I+FI*G1R C11R=F11R*FR-F11I*FI-E11R*ER+E11I*EI C11I=F11R*FI+F11I*FR-E11R*EI-E11I*ER C12R=F12R*FR-F12I*FI-E12R*ER+E12I*EI C12I=F12R*FI+F12I*FR-E12R*EI-E12I*ER C22R=F22R*FR-F22I*FI-E22R*ER+E22I*EI C22I=F22R*FI+F22I*FR-E22R*EI-E22I*ER C13R=F13R*FR-F13I*FI-E13R*ER+E13I*EI C13I=F13R*FI+F13I*FR-E13R*EI-E13I*ER C23R=F23R*FR-F23I*FI-E23R*ER+E23I*EI C23I=F23R*FI+F23I*FR-E23R*EI-E23I*ER C33R=F33R*FR-F33I*FI-E33R*ER+E33I*EI C33I=F33R*FI+F33I*FR-E33R*EI-E33I*ER C Matrix Dij=Fij-G2*Eij D11R=F11R-E11R*G2R+E11I*G2I D11I=F11I-E11R*G2I-E11I*G2R D12R=F12R-E12R*G2R+E12I*G2I D12I=F12I-E12R*G2I-E12I*G2R D22R=F22R-E22R*G2R+E22I*G2I D22I=F22I-E22R*G2I-E22I*G2R D13R=F13R-E13R*G2R+E13I*G2I D13I=F13I-E13R*G2I-E13I*G2R D23R=F23R-E23R*G2R+E23I*G2I D23I=F23I-E23R*G2I-E23I*G2R D33R=F33R-E33R*G2R+E33I*G2I D33I=F33I-E33R*G2I-E33I*G2R CALL CHR11( * C11R,C11I,C12R,C12I,C22R,C22I,C13R,C13I,C23R,C23I,C33R,C33I, * D11R,D11I,D12R,D12I,D22R,D22I,D13R,D13I,D23R,D23I,D33R,D33I, * H11R,H11I,H12R,H12I,H22R,H22I,H13R,H13I,H23R,H23I,H33R,H33I, * H14R,H14I,H24R,H24I,H34R,H34I,H44R,H44I, * H15R,H15I,H25R,H25I,H35R,H35I,H45R,H45I,H55R,H55I, * H16R,H16I,H26R,H26I,H36R,H36I,H46R,H46I,H56R,H56I,H66R,H66I) C Matrix Cij=Cij3 C11R=F11R*C3R-F11I*C3I-E11R*C4R+E11I*C4I C11I=F11R*C3I+F11I*C3R-E11R*C4I-E11I*C4R C12R=F12R*C3R-F12I*C3I-E12R*C4R+E12I*C4I C12I=F12R*C3I+F12I*C3R-E12R*C4I-E12I*C4R C22R=F22R*C3R-F22I*C3I-E22R*C4R+E22I*C4I C22I=F22R*C3I+F22I*C3R-E22R*C4I-E22I*C4R C13R=F13R*C3R-F13I*C3I-E13R*C4R+E13I*C4I C13I=F13R*C3I+F13I*C3R-E13R*C4I-E13I*C4R C23R=F23R*C3R-F23I*C3I-E23R*C4R+E23I*C4I C23I=F23R*C3I+F23I*C3R-E23R*C4I-E23I*C4R C33R=F33R*C3R-F33I*C3I-E33R*C4R+E33I*C4I C33I=F33R*C3I+F33I*C3R-E33R*C4I-E33I*C4R C Matrix Dij=((N-2)/8)*Cij3 AUXR=0.125*DEGREE-0.25 D11R=C11R*AUXR D11I=C11I*AUXR D12R=C12R*AUXR D12I=C12I*AUXR D22R=C22R*AUXR D22I=C22I*AUXR D13R=C13R*AUXR D13I=C13I*AUXR D23R=C23R*AUXR D23I=C23I*AUXR D33R=C33R*AUXR D33I=C33I*AUXR CALL CHR11( * C11R,C11I,C12R,C12I,C22R,C22I,C13R,C13I,C23R,C23I,C33R,C33I, * D11R,D11I,D12R,D12I,D22R,D22I,D13R,D13I,D23R,D23I,D33R,D33I, * H11R,H11I,H12R,H12I,H22R,H22I,H13R,H13I,H23R,H23I,H33R,H33I, * H14R,H14I,H24R,H24I,H34R,H34I,H44R,H44I, * H15R,H15I,H25R,H25I,H35R,H35I,H45R,H45I,H55R,H55I, * H16R,H16I,H26R,H26I,H36R,H36I,H46R,H46I,H56R,H56I,H66R,H66I) END IF C ELSE C 5A8 CALL ERROR('5A8 in HDERC: Zero wave type') C This error should not appear, contact the authors. END IF C C....................................................................... C C Reference real-valued polarization vectors C C S1-wave eigenvector of Christoffel matrix (Klimes, 2022a, eq.47) R11=F11R-E11R*G2R+E11I*G2I R21=F22R-E22R*G2R+E22I*G2I R31=F33R-E33R*G2R+E33I*G2I Y11=F11I-E11R*G2I-E11I*G2R Y21=F22I-E22R*G2I-E22I*G2R Y31=F33I-E33R*G2I-E33I*G2R AUX1=R11*R11+Y11*Y11 AUX2=R21*R21+Y21*Y21 AUX3=R31*R31+Y31*Y31 IF(AUX1.GE.AUX2.AND.AUX1.GE.AUX3) THEN R21=F12R-E12R*G2R+E12I*G2I R31=F13R-E13R*G2R+E13I*G2I Y21=F12I-E12R*G2I-E12I*G2R Y31=F13I-E13R*G2I-E13I*G2R ELSE IF(AUX2.GE.AUX3) THEN R11=F12R-E12R*G2R+E12I*G2I R31=F23R-E23R*G2R+E23I*G2I Y11=F12I-E12R*G2I-E12I*G2R Y31=F23I-E23R*G2I-E23I*G2R ELSE R11=F13R-E13R*G2R+E13I*G2I R21=F23R-E23R*G2R+E23I*G2I Y11=F13I-E13R*G2I-E13I*G2R Y21=F23I-E23R*G2I-E23I*G2R END IF C C S2-wave eigenvector of Christoffel matrix (Klimes, 2022a, eq.48) R12=F11R-E11R*G1R+E11I*G1I R22=F22R-E22R*G1R+E22I*G1I R32=F33R-E33R*G1R+E33I*G1I Y12=F11I-E11R*G1I-E11I*G1R Y22=F22I-E22R*G1I-E22I*G1R Y32=F33I-E33R*G1I-E33I*G1R AUX1=R12*R12+Y12*Y12 AUX2=R22*R22+Y22*Y22 AUX3=R32*R32+Y32*Y32 IF(AUX1.GE.AUX2.AND.AUX1.GE.AUX3) THEN R22=F12R-E12R*G1R+E12I*G1I R32=F13R-E13R*G1R+E13I*G1I Y22=F12I-E12R*G1I-E12I*G1R Y32=F13I-E13R*G1I-E13I*G1R ELSE IF(AUX2.GE.AUX3) THEN R12=F12R-E12R*G1R+E12I*G1I R32=F23R-E23R*G1R+E23I*G1I Y12=F12I-E12R*G1I-E12I*G1R Y32=F23I-E23R*G1I-E23I*G1R ELSE R12=F13R-E13R*G1R+E13I*G1I R22=F23R-E23R*G1R+E23I*G1I Y12=F13I-E13R*G1I-E13I*G1R Y22=F23I-E23R*G1I-E23I*G1R END IF C C Squared norms of the above non-normalized S-wave eigenvectors R1R1=R11*R11+R21*R21+R31*R31 Y1Y1=Y11*Y11+Y21*Y21+Y31*Y31 R2R2=R12*R12+R22*R22+R32*R32 Y2Y2=Y12*Y12+Y22*Y22+Y32*Y32 AUX3=1.E-10*(G3R*G3R+G3I*G3I) LONE=.FALSE. LTWO=.FALSE. IF(R1R1+Y1Y1.LT.AUX3.OR.R2R2+Y2Y2.LT.AUX3) THEN C Singularity with two S-wave eigenvectors of Christoffel matrix C (Klimes, 2022b, sec.5.3) LTWO=.TRUE. GO TO 20 ELSE C Squared complex-valued pseudonorms of the S-wave eigenvectors R1R1=R1R1-Y1Y1 Y1Y1=2.*(R11*Y11+R21*Y21+R31*Y31) R2R2=R2R2-Y2Y2 Y2Y2=2.*(R12*Y12+R22*Y22+R32*Y32) IF(R1R1*R1R1+Y1Y1*Y1Y1.GT.AUX3*AUX3.AND. * R2R2*R2R2+Y2Y2*Y2Y2.GT.AUX3*AUX3) THEN C S-wave eigenvectors have non-zero pseudonorms C (Klimes, 2022b, sec.5.1) C and are pseudonormalized (Klimes, 2022b, eq.112) CALL SQRTC(R1R1,Y1Y1,AUXR,AUXI) AUXXR=AUXR*AUXR+AUXI*AUXI AUXR=AUXR/AUXXR AUXI=-AUXI/AUXXR R11=R11*AUXR-Y11*AUXI R21=R21*AUXR-Y21*AUXI R31=R31*AUXR-Y31*AUXI CALL SQRTC(R2R2,Y2Y2,AUXR,AUXI) AUXXR=AUXR*AUXR+AUXI*AUXI AUXR=AUXR/AUXXR AUXI=-AUXI/AUXXR R12=R12*AUXR-Y12*AUXI R22=R22*AUXR-Y22*AUXI R32=R32*AUXR-Y32*AUXI ELSE C S-wave eigenvectors cannot be pseudonormalized C at a singularity with a single S-wave eigenvector C (Klimes, 2022b, sec.5.2). C Vector Gj0=Rj+i*Yj (Klimes, 2022b, eq.11): AUX1=DCH11R*DCH11R+DCH11I*DCH11I AUX2=DCH22R*DCH22R+DCH22I*DCH22I AUX3=DCH33R*DCH33R+DCH33I*DCH33I IF(AUX1.GE.AUX2.AND.AUX1.GE.AUX3) THEN R1=DCH11R Y1=DCH11I R2=DCH12R Y2=DCH12I R3=DCH13R Y3=DCH13I ELSE IF(AUX2.GE.AUX3) THEN R1=DCH12R Y1=DCH12I R2=DCH22R Y2=DCH22I R3=DCH23R Y3=DCH23I ELSE R1=DCH13R Y1=DCH13I R2=DCH23R Y2=DCH23I R3=DCH33R Y3=DCH33I END IF C Checking the norm of vector Gj0=Rj+i*Yj AUXR=R1*R1+R2*R2+R3*R3 AUXI=Y1*Y1+Y2*Y2+Y3*Y3 IF(AUXR+AUXI.LT.1.E-10*(G3R*G3R+G3I*G3I)) THEN C Vector Gj0=Rj+i*Yj is nearly zero. The singularity thus C should be treated as one with two S-wave eigenvectors C (Klimes, 2022b, sec.5.3) LTWO=.TRUE. GO TO 20 END IF C Vector Ej2 (Klimes, 2022b, eqs.129,130): AUXR=SQRT(AUXR) R12=R1/AUXR R22=R2/AUXR R32=R3/AUXR E(4)=R12 E(5)=R22 E(6)=R32 AUXI=SQRT(AUXI) IF(AUXI.GE.1.E-10**AUXR) THEN C Vector Ej1 (Klimes, 2022b, eqs.128,130): R11=Y1/AUXI R21=Y2/AUXI R31=Y3/AUXI E(1)=R11 E(2)=R21 E(3)=R31 GO TO 80 ELSE C Too small vector Im(Gj0)=Yj LONE=.TRUE. END IF END IF END IF C C P-wave eigenvector of Christoffel matrix (Klimes, 2022b, sec.2.3) 20 CONTINUE IF(ICB.GT.0.OR.LONE.OR.LTWO) THEN AUX1=G11R*G11R+G11I*G11I AUX2=G22R*G22R+G22I*G22I AUX3=G33R*G33R+G33I*G33I IF(AUX1.GE.AUX2.AND.AUX1.GE.AUX3) THEN R13=G11R R23=G12R R33=G13R Y13=G11I Y23=G12I Y33=G13I ELSE IF(AUX2.GE.AUX3) THEN R13=G12R R23=G22R R33=G23R Y13=G12I Y23=G22I Y33=G23I ELSE R13=G13R R23=G23R R33=G33R Y13=G13I Y23=G23I Y33=G33I END IF C C Improving the P-wave eigenvector of the Christoffel matrix C (Klimes, 2022b, eq.9) R1=G11R*R13+G12R*R23+G13R*R33-G11I*Y13-G12I*Y23-G13I*Y33 R2=G12R*R13+G22R*R23+G23R*R33-G12I*Y13-G22I*Y23-G23I*Y33 R3=G13R*R13+G23R*R23+G33R*R33-G13I*Y13-G23I*Y23-G33I*Y33 Y1=G11R*Y13+G12R*Y23+G13R*Y33+G11I*R13+G12I*R23+G13I*R33 Y2=G12R*Y13+G22R*Y23+G23R*Y33+G12I*R13+G22I*R23+G23I*R33 Y3=G13R*Y13+G23R*Y23+G33R*Y33+G13I*R13+G23I*R23+G33I*R33 C C Normalizing the P-wave eigenvector of the Christoffel matrix C (Klimes, 2022b, eq.10) AUXXR=R1*R1+R2*R2+R3*R3-Y1*Y1-Y2*Y2-Y3*Y3 AUXXI=2.*(R1*Y1+R2*Y2+R3*Y3) CALL SQRTC(AUXXR,AUXXI,AUXR,AUXI) AUXXR=AUXR*AUXR+AUXI*AUXI AUXR=AUXR/AUXXR AUXI=-AUXI/AUXXR R13=R1*AUXR-Y1*AUXI R23=R2*AUXR-Y2*AUXI R33=R3*AUXR-Y3*AUXI Y13=R1*AUXI+Y1*AUXR Y23=R2*AUXI+Y2*AUXR Y33=R3*AUXI+Y3*AUXR C C Reference real-valued P-wave polarization vector C (Klimes, 2022b, eqs.17,112,115,136) AUXR=SQRT(R13*R13+R23*R23+R33*R33) R13=R13/AUXR R23=R23/AUXR R33=R33/AUXR E(7)=R13 E(8)=R23 E(9)=R33 END IF C C LTWO=.TRUE. ... Reference vector Ej3 is defined only. C LONE=.TRUE. ... Reference vectors Ej2 and Ej3 are defined. C else C ICB.GT.0 ... Reference vectors Ej1, Ej2 and Ej3 are defined. C ICB.LT.0 ... Reference vectors Ej1 and Ej2 are defined. C Continuating at label 80. C C Singularity with two S-wave eigenvectors of the Christoffel matrix C (Klimes, 2022b, sec.5.3) IF(LTWO) THEN AUXR=R13*R13+R23*R23+R33*R33 AUXI=Y13*Y13+Y23*Y23+Y33*Y33 IF(AUXI.GE.1.E-10**AUXR) THEN C (Klimes, 2022b, eq.136) AUX2=R13*Y13+R23*Y23+R33*Y33 R12=Y13-R13*AUX2 R22=Y23-R23*AUX2 R32=Y33-R33*AUX2 ELSE C Too small vector Im(Gj3)=Yj3. C Replacing Yj3 by the basis vector most perpendicular to Rj3 IF(R13*R13.LE.R23*R23.AND.R13*R13.LE.R33*R33) THEN R12=1.0-R13*R13 R22= -R23*R13 R32= -R33*R13 ELSE IF(R23*R23.LE.R33*R33) THEN R12= -R13*R23 R22=1.0-R23*R23 R32= -R33*R23 ELSE R12= -R13*R33 R22= -R23*R33 R32=1.0-R33*R33 END IF END IF C Normalizing Rj2 (Klimes, 2022b, eq.18) AUX2=SQRT(R12*R12+R22*R22+R32*R32) R12=R12/AUX2 R22=R22/AUX2 R32=R32/AUX2 E(4)=R12 E(5)=R22 E(6)=R32 END IF C C S-wave singularity IF(LONE.OR.LTWO) THEN C (Klimes, 2022b, eq.19) E(1)=R22*R33-R32*R23 E(2)=R32*R13-R12*R33 E(3)=R12*R23-R22*R13 GO TO 90 END IF C C P wave IF(ICB.GT.0) THEN C Projection of the real parts of the S-wave eigenvectors C (Klimes, 2022b, eq.116) AUX1=R13*R11+R23*R21+R33*R31 R11=R11-R13*AUX1 R21=R21-R23*AUX1 R31=R31-R33*AUX1 AUX2=R13*R12+R23*R22+R33*R32 R12=R12-R13*AUX2 R22=R22-R23*AUX2 R32=R32-R33*AUX2 END IF C C Normalization of the reference S-wave polarization vectors C (Klimes, 2022b, eq.117) AUX1=SQRT(R11*R11+R21*R21+R31*R31) R11=R11/AUX1 R21=R21/AUX1 R31=R31/AUX1 AUX2=SQRT(R12*R12+R22*R22+R32*R32) R12=R12/AUX2 R22=R22/AUX2 R32=R32/AUX2 C Matrix of scalar products (Klimes, 2022b, eqs.118,120) RR11=R11*R11+R21*R21+R31*R31 RR12=R11*R12+R21*R22+R31*R32 RR22=R12*R12+R22*R22+R32*R32 DETSR=SQRT(RR11*RR22-RR12*RR12) C Square root of the matrix of scalar products AUX1=SQRT(RR11+RR22+DETSR+DETSR) SR11=(RR11+DETSR)/AUX1 SR12=RR12/AUX1 SR22=(RR22+DETSR)/AUX1 C Inverse square root of the matrix of scalar products RR11= SR22/DETSR RR12=-SR12/DETSR RR22= SR11/DETSR C Reference real-valued S-wave polarization vectors C (Klimes, 2022b, eqs.119,121) E(1)=R11*RR11+R12*RR12 E(2)=R21*RR11+R22*RR12 E(3)=R31*RR11+R32*RR12 E(4)=R11*RR12+R12*RR22 E(5)=R21*RR12+R22*RR22 E(6)=R31*RR12+R32*RR22 C C Reference real-valued P-wave polarization vector for S wave C (Klimes, 2022b, eqs.122,131) 80 CONTINUE IF(ICB.LT.0) THEN E(7)=E(2)*E(6)-E(3)*E(5) E(8)=E(3)*E(4)-E(1)*E(6) E(9)=E(1)*E(5)-E(2)*E(4) END IF C 90 CONTINUE RETURN END C C======================================================================= C SUBROUTINE CHR0(A,Q,P1,P2,P3, * CH11R,CH11I,CH12R,CH12I,CH22R,CH22I, * CH13R,CH13I,CH23R,CH23I,CH33R,CH33I) C C Subroutine calculating the complex-valued Christoffel matrix. C C ENTRY CHR1(A,Q,P1,P2,P3, C * E11R,E11I,E12R,E12I,E22R,E22I,E13R,E13I,E23R,E23I,E33R,E33I, C * G1R,G1I,G2R,G2I,G3R,G3I,G4R,G4I,G5R,G5I,G6R,G6I) C C Subroutine calculating and storing the first-order phase-space C derivatives Gij' of the Christoffel matrix Gij, C and returning the products G'=Gij'*Eij of these first-order C phase-space derivatives with given symmetric matrix Eik. C CHR0 must be invoked before CHR1. C C ENTRY CHR11( C * E11R,E11I,E12R,E12I,E22R,E22I,E13R,E13I,E23R,E23I,E33R,E33I, C * F11R,F11I,F12R,F12I,F22R,F22I,F13R,F13I,F23R,F23I,F33R,F33I, C * G11R,G11I,G12R,G12I,G22R,G22I,G13R,G13I,G23R,G23I,G33R,G33I, C * G14R,G14I,G24R,G24I,G34R,G34I,G44R,G44I, C * G15R,G15I,G25R,G25I,G35R,G35I,G45R,G45I,G55R,G55I, C * G16R,G16I,G26R,G26I,G36R,G36I,G46R,G46I,G56R,G56I,G66R,G66I) C C Subroutine calculating products C GG'"=Gij'*Eik*Fjl*Gkl" C where Gij' and Gij" are the first-order phase-space derivatives of C symmetric Christoffel matrix Gij stored by subroutine CHR1, C and Eik and Fik are given symmetric matrices. C Subroutine CHR11 increases the input values G'" by products GG'". C CHR0 and CHR1 must be invoked before CHR11. C C ENTRY CHR2(A,Q,P1,P2,P3, C * E11R,E11I,E12R,E12I,E22R,E22I,E13R,E13I,E23R,E23I,E33R,E33I, C * G11R,G11I,G12R,G12I,G22R,G22I,G13R,G13I,G23R,G23I,G33R,G33I, C * G14R,G14I,G24R,G24I,G34R,G34I,G44R,G44I, C * G15R,G15I,G25R,G25I,G35R,G35I,G45R,G45I,G55R,G55I, C * G16R,G16I,G26R,G26I,G36R,G36I,G46R,G46I,G56R,G56I,G66R,G66I) C C Subroutine calculating the second-order phase-space derivatives C Gij'" of the Christoffel matrix Gij, C and returning the products G'"=Gij'"*Eij of these second-order C phase-space derivatives with given symmetric matrix Eik. C CHR0 must be invoked before CHR2. C REAL A(10,21),Q(10,21),P1,P2,P3 REAL CH11R,CH11I,CH12R,CH12I,CH22R,CH22I REAL CH13R,CH13I,CH23R,CH23I,CH33R,CH33I REAL E11R,E11I,E12R,E12I,E22R,E22I,E13R,E13I,E23R,E23I,E33R,E33I REAL F11R,F11I,F12R,F12I,F22R,F22I,F13R,F13I,F23R,F23I,F33R,F33I REAL G1R,G1I,G2R,G2I,G3R,G3I,G4R,G4I,G5R,G5I,G6R,G6I REAL G11R,G11I,G12R,G12I,G22R,G22I,G13R,G13I,G23R,G23I,G33R,G33I REAL G14R,G14I,G24R,G24I,G34R,G34I,G44R,G44I REAL G15R,G15I,G25R,G25I,G35R,G35I,G45R,G45I,G55R,G55I REAL G16R,G16I,G26R,G26I,G36R,G36I,G46R,G46I,G56R,G56I,G66R,G66I C C Input: C A... Real parts of the elastic moduli and their first and C second spatial derivatives. C Q... Imaginary parts of the elastic moduli and their first and C second spatial derivatives. C P1,P2,P3... Slowness vector. C E11R,E11I,E12R,E12I,E22R,E22I,E13R,E13I,E23R,E23I,E33R,E33I... C Symmetric complex-valued 3*3 matrix. C F11R,F11I,F12R,F12I,F22R,F22I,F13R,F13I,F23R,F23I,F33R,F33I... C Additional symmetric complex-valued 3*3 matrix. C C Output: C CH11R,CH11I,CH12R,CH12I,CH22R,CH22I,CH13R,CH13I,CH23R,CH23I,CH33R, C CH33I... Symmetric complex-valued Christoffel matrix. C G1R,G1I,G2R,G2I,G3R,G3I,G4R,G4I,G5R,G5I,G6R,G6I... First-order C partial phase-space derivatives of the Christoffel matrix, C multiplied by given matrix Eij in CHR1. C G11R,G11I,G12R,G12I,G22R,G22I,G13R,G13I,G23R,G23I,G33R,G33I, C G14R,G14I,G24R,G24I,G34R,G34I,G44R,G44I,G15R,G15I,G25R,G25I, C G35R,G35I,G45R,G45I,G55R,G55I,G16R,G16I,G26R,G26I,G36R,G36I, C G46R,G46I,G56R,G56I,G66R,G66I... Input values increased by the C product of the first-order partial phase-space derivatives C of the Christoffel matrix, multiplied by given matrices C Eij and Fij in CRH11. C Second-order partial phase-space derivatives of the C Christoffel matrix, multiplied by given matrix Eij in C CRH2. C C No subroutines and external functions required: C C Date: 2022, November 9 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Quantities to be saved: C C Real parts of elastic moduli A(i,j,k,l): REAL AIJKL(21) SAVE AIJKL 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),A1111) EQUIVALENCE (AIJKL( 2),A1122,A2211) EQUIVALENCE (AIJKL( 3),A2222) EQUIVALENCE (AIJKL( 4),A1133,A3311) EQUIVALENCE (AIJKL( 5),A2233,A3322) EQUIVALENCE (AIJKL( 6),A3333) EQUIVALENCE (AIJKL( 7),A1123,A1132,A2311,A3211) EQUIVALENCE (AIJKL( 8),A2223,A2232,A2322,A3222) EQUIVALENCE (AIJKL( 9),A3323,A3332,A2333,A3233) EQUIVALENCE (AIJKL(10),A2323,A3223,A2332,A3232) EQUIVALENCE (AIJKL(11),A1113,A1131,A1311,A3111) EQUIVALENCE (AIJKL(12),A2213,A2231,A1322,A3122) EQUIVALENCE (AIJKL(13),A3313,A3331,A1333,A3133) EQUIVALENCE (AIJKL(14),A2313,A3213,A2331,A3231, * A1323,A1332,A3123,A3132) EQUIVALENCE (AIJKL(15),A1313,A3113,A1331,A3131) EQUIVALENCE (AIJKL(16),A1112,A1121,A1211,A2111) EQUIVALENCE (AIJKL(17),A2212,A2221,A1222,A2122) EQUIVALENCE (AIJKL(18),A3312,A3321,A1233,A2133) EQUIVALENCE (AIJKL(19),A2312,A3212,A2321,A3221, * A1223,A1232,A2123,A2132) EQUIVALENCE (AIJKL(20),A1312,A3112,A1321,A3121, * A1213,A1231,A2113,A2131) EQUIVALENCE (AIJKL(21),A1212,A2112,A1221,A2121) C C Imaginary parts of elastic moduli Q(i,j,k,l): REAL QIJKL(21) SAVE QIJKL REAL Q1111 REAL Q1122,Q2211 REAL Q2222 REAL Q1133,Q3311 REAL Q2233,Q3322 REAL Q3333 REAL Q1123,Q1132,Q2311,Q3211 REAL Q2223,Q2232,Q2322,Q3222 REAL Q3323,Q3332,Q2333,Q3233 REAL Q2323,Q3223,Q2332,Q3232 REAL Q1113,Q1131,Q1311,Q3111 REAL Q2213,Q2231,Q1322,Q3122 REAL Q3313,Q3331,Q1333,Q3133 REAL Q2313,Q3213,Q2331,Q3231,Q1323,Q3123,Q1332,Q3132 REAL Q1313,Q3113,Q1331,Q3131 REAL Q1112,Q1121,Q1211,Q2111 REAL Q2212,Q2221,Q1222,Q2122 REAL Q3312,Q3321,Q1233,Q2133 REAL Q2312,Q3212,Q2321,Q3221,Q1223,Q2123,Q1232,Q2132 REAL Q1312,Q3112,Q1321,Q3121,Q1213,Q2113,Q1231,Q2131 REAL Q1212,Q2112,Q1221,Q2121 EQUIVALENCE (QIJKL( 1),Q1111) EQUIVALENCE (QIJKL( 2),Q1122,Q2211) EQUIVALENCE (QIJKL( 3),Q2222) EQUIVALENCE (QIJKL( 4),Q1133,Q3311) EQUIVALENCE (QIJKL( 5),Q2233,Q3322) EQUIVALENCE (QIJKL( 6),Q3333) EQUIVALENCE (QIJKL( 7),Q1123,Q1132,Q2311,Q3211) EQUIVALENCE (QIJKL( 8),Q2223,Q2232,Q2322,Q3222) EQUIVALENCE (QIJKL( 9),Q3323,Q3332,Q2333,Q3233) EQUIVALENCE (QIJKL(10),Q2323,Q3223,Q2332,Q3232) EQUIVALENCE (QIJKL(11),Q1113,Q1131,Q1311,Q3111) EQUIVALENCE (QIJKL(12),Q2213,Q2231,Q1322,Q3122) EQUIVALENCE (QIJKL(13),Q3313,Q3331,Q1333,Q3133) EQUIVALENCE (QIJKL(14),Q2313,Q3213,Q2331,Q3231, * Q1323,Q1332,Q3123,Q3132) EQUIVALENCE (QIJKL(15),Q1313,Q3113,Q1331,Q3131) EQUIVALENCE (QIJKL(16),Q1112,Q1121,Q1211,Q2111) EQUIVALENCE (QIJKL(17),Q2212,Q2221,Q1222,Q2122) EQUIVALENCE (QIJKL(18),Q3312,Q3321,Q1233,Q2133) EQUIVALENCE (QIJKL(19),Q2312,Q3212,Q2321,Q3221, * Q1223,Q1232,Q2123,Q2132) EQUIVALENCE (QIJKL(20),Q1312,Q3112,Q1321,Q3121, * Q1213,Q1231,Q2113,Q2131) EQUIVALENCE (QIJKL(21),Q1212,Q2112,Q1221,Q2121) 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 Q(i,j,k)=Q(i,j,k,l)*P(l): REAL Q111,Q211,Q311,Q121,Q221,Q321,Q131,Q231,Q331 REAL Q112,Q212,Q312,Q122,Q222,Q322,Q132,Q232,Q332 REAL Q113,Q213,Q313,Q123,Q223,Q323,Q133,Q233,Q333 SAVE Q111,Q121,Q221,Q131,Q231,Q331 SAVE Q112,Q122,Q222,Q132,Q232,Q332 SAVE Q113,Q123,Q223,Q133,Q233,Q333 EQUIVALENCE (Q121,Q211),(Q131,Q311),(Q231,Q321) EQUIVALENCE (Q122,Q212),(Q132,Q312),(Q232,Q322) EQUIVALENCE (Q123,Q213),(Q133,Q313),(Q233,Q323) C C Tensor product of the slowness vector C calculated by CHR1 and then also used by CHR2: REAL P11,P12,P22,P13,P23,P33,D12,D13,D23 SAVE P11,P12,P22,P13,P23,P33,D12,D13,D23 C C First-order phase-space derivatives of the Christoffel matrix C calculated by CHR1 and then also used by CHR11: REAL G111R,G121R,G221R,G131R,G231R,G331R REAL G112R,G122R,G222R,G132R,G232R,G332R REAL G113R,G123R,G223R,G133R,G233R,G333R REAL G114R,G124R,G224R,G134R,G234R,G334R REAL G115R,G125R,G225R,G135R,G235R,G335R REAL G116R,G126R,G226R,G136R,G236R,G336R REAL G111I,G121I,G221I,G131I,G231I,G331I REAL G112I,G122I,G222I,G132I,G232I,G332I REAL G113I,G123I,G223I,G133I,G233I,G333I REAL G114I,G124I,G224I,G134I,G234I,G334I REAL G115I,G125I,G225I,G135I,G235I,G335I REAL G116I,G126I,G226I,G136I,G236I,G336I SAVE G111R,G121R,G221R,G131R,G231R,G331R SAVE G112R,G122R,G222R,G132R,G232R,G332R SAVE G113R,G123R,G223R,G133R,G233R,G333R SAVE G114R,G124R,G224R,G134R,G234R,G334R SAVE G115R,G125R,G225R,G135R,G235R,G335R SAVE G116R,G126R,G226R,G136R,G236R,G336R SAVE G111I,G121I,G221I,G131I,G231I,G331I SAVE G112I,G122I,G222I,G132I,G232I,G332I SAVE G113I,G123I,G223I,G133I,G233I,G333I SAVE G114I,G124I,G224I,G134I,G234I,G334I SAVE G115I,G125I,G225I,G135I,G235I,G335I SAVE G116I,G126I,G226I,G136I,G236I,G336I C C----------------------------------------------------------------------- C C Temporary storage locations: C C Loop variable INTEGER I C C Spatial derivatives of elastic moduli A(i,j,k,l) or Q(i,j,k,l): REAL BIJKL(21) REAL B1111 REAL B1122,B2211 REAL B2222 REAL B1133,B3311 REAL B2233,B3322 REAL B3333 REAL B1123,B1132,B2311,B3211 REAL B2223,B2232,B2322,B3222 REAL B3323,B3332,B2333,B3233 REAL B2323,B3223,B2332,B3232 REAL B1113,B1131,B1311,B3111 REAL B2213,B2231,B1322,B3122 REAL B3313,B3331,B1333,B3133 REAL B2313,B3213,B2331,B3231,B1323,B3123,B1332,B3132 REAL B1313,B3113,B1331,B3131 REAL B1112,B1121,B1211,B2111 REAL B2212,B2221,B1222,B2122 REAL B3312,B3321,B1233,B2133 REAL B2312,B3212,B2321,B3221,B1223,B2123,B1232,B2132 REAL B1312,B3112,B1321,B3121,B1213,B2113,B1231,B2131 REAL B1212,B2112,B1221,B2121 EQUIVALENCE (BIJKL( 1),B1111) EQUIVALENCE (BIJKL( 2),B1122,B2211) EQUIVALENCE (BIJKL( 3),B2222) EQUIVALENCE (BIJKL( 4),B1133,B3311) EQUIVALENCE (BIJKL( 5),B2233,B3322) EQUIVALENCE (BIJKL( 6),B3333) EQUIVALENCE (BIJKL( 7),B1123,B1132,B2311,B3211) EQUIVALENCE (BIJKL( 8),B2223,B2232,B2322,B3222) EQUIVALENCE (BIJKL( 9),B3323,B3332,B2333,B3233) EQUIVALENCE (BIJKL(10),B2323,B3223,B2332,B3232) EQUIVALENCE (BIJKL(11),B1113,B1131,B1311,B3111) EQUIVALENCE (BIJKL(12),B2213,B2231,B1322,B3122) EQUIVALENCE (BIJKL(13),B3313,B3331,B1333,B3133) EQUIVALENCE (BIJKL(14),B2313,B3213,B2331,B3231, * B1323,B1332,B3123,B3132) EQUIVALENCE (BIJKL(15),B1313,B3113,B1331,B3131) EQUIVALENCE (BIJKL(16),B1112,B1121,B1211,B2111) EQUIVALENCE (BIJKL(17),B2212,B2221,B1222,B2122) EQUIVALENCE (BIJKL(18),B3312,B3321,B1233,B2133) EQUIVALENCE (BIJKL(19),B2312,B3212,B2321,B3221, * B1223,B1232,B2123,B2132) EQUIVALENCE (BIJKL(20),B1312,B3112,B1321,B3121, * B1213,B1231,B2113,B2131) EQUIVALENCE (BIJKL(21),B1212,B2112,B1221,B2121) C C Temporary storage locations in CHR1: REAL C1122,C1133,C2233 REAL C1132,C1123 REAL C1232,C2123 REAL C1233,C2133 EQUIVALENCE (C1123,C1132) EQUIVALENCE (C1232,C2123) EQUIVALENCE (C1233,C2133) C C Weighting factors Wijkl=Eik*Fjl in CHR11: REAL W1111R,W1112R,W1122R,W1113R,W1123R,W1133R REAL W1211R,W1212R,W1222R,W1213R,W1223R,W1233R REAL W2211R,W2212R,W2222R,W2213R,W2223R,W2233R REAL W1311R,W1312R,W1322R,W1313R,W1323R,W1333R REAL W2311R,W2312R,W2322R,W2313R,W2323R,W2333R REAL W3311R,W3312R,W3322R,W3313R,W3323R,W3333R REAL W1111I,W1112I,W1122I,W1113I,W1123I,W1133I REAL W1211I,W1212I,W1222I,W1213I,W1223I,W1233I REAL W2211I,W2212I,W2222I,W2213I,W2223I,W2233I REAL W1311I,W1312I,W1322I,W1313I,W1323I,W1333I REAL W2311I,W2312I,W2322I,W2313I,W2323I,W2333I REAL W3311I,W3312I,W3322I,W3313I,W3323I,W3333I C Wklij=Wijkl EQUIVALENCE (W1112R,W1211R),(W1122R,W2211R),(W1222R,W2212R) EQUIVALENCE (W1113R,W1311R),(W1213R,W1312R),(W2213R,W1322R) EQUIVALENCE (W1123R,W2311R),(W1223R,W2312R),(W2223R,W2322R) EQUIVALENCE (W1323R,W2313R),(W1133R,W3311R),(W1233R,W3312R) EQUIVALENCE (W2233R,W3322R),(W1333R,W3313R),(W2333R,W3323R) EQUIVALENCE (W1112I,W1211I),(W1122I,W2211I),(W1222I,W2212I) EQUIVALENCE (W1113I,W1311I),(W1213I,W1312I),(W2213I,W1322I) EQUIVALENCE (W1123I,W2311I),(W1223I,W2312I),(W2223I,W2322I) EQUIVALENCE (W1323I,W2313I),(W1133I,W3311I),(W1233I,W3312I) EQUIVALENCE (W2233I,W3322I),(W1333I,W3313I),(W2333I,W3323I) C C Temporary storage locations in CHR11: REAL X111R,X121R,X221R,X131R,X231R,X331R REAL X112R,X122R,X222R,X132R,X232R,X332R REAL X113R,X123R,X223R,X133R,X233R,X333R REAL X114R,X124R,X224R,X134R,X234R,X334R REAL X115R,X125R,X225R,X135R,X235R,X335R REAL X116R,X126R,X226R,X136R,X236R,X336R REAL X111I,X121I,X221I,X131I,X231I,X331I REAL X112I,X122I,X222I,X132I,X232I,X332I REAL X113I,X123I,X223I,X133I,X233I,X333I REAL X114I,X124I,X224I,X134I,X234I,X334I REAL X115I,X125I,X225I,X135I,X235I,X335I REAL X116I,X126I,X226I,X136I,X236I,X336I C C Weighting factors in CHR2: REAL W11R,W12R,W13R,W14R,W15R,W16R REAL W22R,W23R,W24R,W25R,W26R REAL W33R,W34R,W35R,W36R REAL W44R,W45R,W46R REAL W55R,W56R REAL W66R REAL W11I,W12I,W13I,W14I,W15I,W16I REAL W22I,W23I,W24I,W25I,W26I REAL W33I,W34I,W35I,W36I REAL W44I,W45I,W46I REAL W55I,W56I REAL W66I C Temporary storage locations in CHR2: REAL G(10) C C----------------------------------------------------------------------- C C Storing the elastic moduli: DO 10 I=1,21 AIJKL(I)=A(1,I) QIJKL(I)=Q(1,I) 10 CONTINUE C C Multiplying the elastic tensor by the slowness vector: 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 C Q111=Q1111*P1+Q1112*P2+Q1113*P3 Q121=Q1211*P1+Q1212*P2+Q1213*P3 Q221=Q2211*P1+Q2212*P2+Q2213*P3 Q131=Q1311*P1+Q1312*P2+Q1313*P3 Q231=Q2311*P1+Q2312*P2+Q2313*P3 Q331=Q3311*P1+Q3312*P2+Q3313*P3 Q112=Q1121*P1+Q1122*P2+Q1123*P3 Q122=Q1221*P1+Q1222*P2+Q1223*P3 Q222=Q2221*P1+Q2222*P2+Q2223*P3 Q132=Q1321*P1+Q1322*P2+Q1323*P3 Q232=Q2321*P1+Q2322*P2+Q2323*P3 Q332=Q3321*P1+Q3322*P2+Q3323*P3 Q113=Q1131*P1+Q1132*P2+Q1133*P3 Q123=Q1231*P1+Q1232*P2+Q1233*P3 Q223=Q2231*P1+Q2232*P2+Q2233*P3 Q133=Q1331*P1+Q1332*P2+Q1333*P3 Q233=Q2331*P1+Q2332*P2+Q2333*P3 Q333=Q3331*P1+Q3332*P2+Q3333*P3 C C Christoffel matrix: CH11R=P1*A111+P2*A211+P3*A311 CH12R=P1*A112+P2*A212+P3*A312 CH22R=P1*A122+P2*A222+P3*A322 CH13R=P1*A113+P2*A213+P3*A313 CH23R=P1*A123+P2*A223+P3*A323 CH33R=P1*A133+P2*A233+P3*A333 C CH11I=P1*Q111+P2*Q211+P3*Q311 CH12I=P1*Q112+P2*Q212+P3*Q312 CH22I=P1*Q122+P2*Q222+P3*Q322 CH13I=P1*Q113+P2*Q213+P3*Q313 CH23I=P1*Q123+P2*Q223+P3*Q323 CH33I=P1*Q133+P2*Q233+P3*Q333 C RETURN C C----------------------------------------------------------------------- C ENTRY CHR1(A,Q,P1,P2,P3, * E11R,E11I,E12R,E12I,E22R,E22I,E13R,E13I,E23R,E23I,E33R,E33I, * G1R,G1I,G2R,G2I,G3R,G3I,G4R,G4I,G5R,G5I,G6R,G6I) C C Entry calculating and storing the first-order phase-space C derivatives Gij' of the Christoffel matrix Gij, C and returning the products G'=Gij'*Eij of these first-order C phase-space derivatives with given symmetric matrix Eik. C C Tensor square of the slowness vector P11=P1*P1 P12=P1*P2 P22=P2*P2 P13=P1*P3 P23=P2*P3 P33=P3*P3 D12=P12+P12 D13=P13+P13 D23=P23+P23 C C Real parts of the spatial derivatives of the Christoffel matrix C (Klimes, 2022a, eq.8) DO 21 I=1,21 BIJKL(I)=A(2,I) 21 CONTINUE C1122=B1122+B2121 C1133=B1133+B3131 C2233=B2233+B3232 C1132=B1132+B2131 C1232=B1232+B2231 C1233=B1233+B3231 G111R=B1111*P11+B1112*D12+B2112*P22+B1113*D13+B2113*D23+B3113*P33 G121R=B1121*P11+C1122*P12+B2122*P22+C1123*P13+C2123*P23+B3123*P33 G221R=B1221*P11+B1222*D12+B2222*P22+B1223*D13+B2223*D23+B3223*P33 G131R=B1131*P11+C1132*P12+B2132*P22+C1133*P13+C2133*P23+B3133*P33 G231R=B1231*P11+C1232*P12+B2232*P22+C1233*P13+C2233*P23+B3233*P33 G331R=B1331*P11+B1332*D12+B2332*P22+B1333*D13+B2333*D23+B3333*P33 DO 22 I=1,21 BIJKL(I)=A(3,I) 22 CONTINUE C1122=B1122+B2121 C1133=B1133+B3131 C2233=B2233+B3232 C1132=B1132+B2131 C1232=B1232+B2231 C1233=B1233+B3231 G112R=B1111*P11+B1112*D12+B2112*P22+B1113*D13+B2113*D23+B3113*P33 G122R=B1121*P11+C1122*P12+B2122*P22+C1123*P13+C2123*P23+B3123*P33 G222R=B1221*P11+B1222*D12+B2222*P22+B1223*D13+B2223*D23+B3223*P33 G132R=B1131*P11+C1132*P12+B2132*P22+C1133*P13+C2133*P23+B3133*P33 G232R=B1231*P11+C1232*P12+B2232*P22+C1233*P13+C2233*P23+B3233*P33 G332R=B1331*P11+B1332*D12+B2332*P22+B1333*D13+B2333*D23+B3333*P33 DO 23 I=1,21 BIJKL(I)=A(4,I) 23 CONTINUE C1122=B1122+B2121 C1133=B1133+B3131 C2233=B2233+B3232 C1132=B1132+B2131 C1232=B1232+B2231 C1233=B1233+B3231 G113R=B1111*P11+B1112*D12+B2112*P22+B1113*D13+B2113*D23+B3113*P33 G123R=B1121*P11+C1122*P12+B2122*P22+C1123*P13+C2123*P23+B3123*P33 G223R=B1221*P11+B1222*D12+B2222*P22+B1223*D13+B2223*D23+B3223*P33 G133R=B1131*P11+C1132*P12+B2132*P22+C1133*P13+C2133*P23+B3133*P33 G233R=B1231*P11+C1232*P12+B2232*P22+C1233*P13+C2233*P23+B3233*P33 G333R=B1331*P11+B1332*D12+B2332*P22+B1333*D13+B2333*D23+B3333*P33 C C Imaginary parts of spatial derivatives of the Christoffel matrix DO 24 I=1,21 BIJKL(I)=Q(2,I) 24 CONTINUE C1122=B1122+B2121 C1133=B1133+B3131 C2233=B2233+B3232 C1132=B1132+B2131 C1232=B1232+B2231 C1233=B1233+B3231 G111I=B1111*P11+B1112*D12+B2112*P22+B1113*D13+B2113*D23+B3113*P33 G121I=B1121*P11+C1122*P12+B2122*P22+C1123*P13+C2123*P23+B3123*P33 G221I=B1221*P11+B1222*D12+B2222*P22+B1223*D13+B2223*D23+B3223*P33 G131I=B1131*P11+C1132*P12+B2132*P22+C1133*P13+C2133*P23+B3133*P33 G231I=B1231*P11+C1232*P12+B2232*P22+C1233*P13+C2233*P23+B3233*P33 G331I=B1331*P11+B1332*D12+B2332*P22+B1333*D13+B2333*D23+B3333*P33 DO 25 I=1,21 BIJKL(I)=Q(3,I) 25 CONTINUE C1122=B1122+B2121 C1133=B1133+B3131 C2233=B2233+B3232 C1132=B1132+B2131 C1232=B1232+B2231 C1233=B1233+B3231 G112I=B1111*P11+B1112*D12+B2112*P22+B1113*D13+B2113*D23+B3113*P33 G122I=B1121*P11+C1122*P12+B2122*P22+C1123*P13+C2123*P23+B3123*P33 G222I=B1221*P11+B1222*D12+B2222*P22+B1223*D13+B2223*D23+B3223*P33 G132I=B1131*P11+C1132*P12+B2132*P22+C1133*P13+C2133*P23+B3133*P33 G232I=B1231*P11+C1232*P12+B2232*P22+C1233*P13+C2233*P23+B3233*P33 G332I=B1331*P11+B1332*D12+B2332*P22+B1333*D13+B2333*D23+B3333*P33 DO 26 I=1,21 BIJKL(I)=Q(4,I) 26 CONTINUE C1122=B1122+B2121 C1133=B1133+B3131 C2233=B2233+B3232 C1132=B1132+B2131 C1232=B1232+B2231 C1233=B1233+B3231 G113I=B1111*P11+B1112*D12+B2112*P22+B1113*D13+B2113*D23+B3113*P33 G123I=B1121*P11+C1122*P12+B2122*P22+C1123*P13+C2123*P23+B3123*P33 G223I=B1221*P11+B1222*D12+B2222*P22+B1223*D13+B2223*D23+B3223*P33 G133I=B1131*P11+C1132*P12+B2132*P22+C1133*P13+C2133*P23+B3133*P33 G233I=B1231*P11+C1232*P12+B2232*P22+C1233*P13+C2233*P23+B3233*P33 G333I=B1331*P11+B1332*D12+B2332*P22+B1333*D13+B2333*D23+B3333*P33 C C Real parts of the slowness derivatives of the Christoffel matrix C (Klimes, 2022a, eq.9) G114R=A111+A111 G124R=A112+A121 G224R=A122+A122 G134R=A113+A131 G234R=A123+A132 G334R=A133+A133 G115R=A211+A211 G125R=A212+A221 G225R=A222+A222 G135R=A213+A231 G235R=A223+A232 G335R=A233+A233 G116R=A311+A311 G126R=A312+A321 G226R=A322+A322 G136R=A313+A331 G236R=A323+A332 G336R=A333+A333 C C Imaginary parts of slowness derivatives of the Christoffel matrix G114I=Q111+Q111 G124I=Q112+Q121 G224I=Q122+Q122 G134I=Q113+Q131 G234I=Q123+Q132 G334I=Q133+Q133 G115I=Q211+Q211 G125I=Q212+Q221 G225I=Q222+Q222 G135I=Q213+Q231 G235I=Q223+Q232 G335I=Q233+Q233 G116I=Q311+Q311 G126I=Q312+Q321 G226I=Q322+Q322 G136I=Q313+Q331 G236I=Q323+Q332 G336I=Q333+Q333 C C Calculating the trace of the matrix product C (Klimes, 2022a, eq.40) G1R=G111R*E11R+G221R*E22R+G331R*E33R * -G111I*E11I-G221I*E22I-G331I*E33I * +(G121R*E12R+G131R*E13R+G231R*E23R * -G121I*E12I-G131I*E13I-G231I*E23I)*2. G1I=G111R*E11I+G221R*E22I+G331R*E33I * +G111I*E11R+G221I*E22R+G331I*E33R * +(G121R*E12I+G131R*E13I+G231R*E23I * +G121I*E12R+G131I*E13R+G231I*E23R)*2. G2R=G112R*E11R+G222R*E22R+G332R*E33R * -G112I*E11I-G222I*E22I-G332I*E33I * +(G122R*E12R+G132R*E13R+G232R*E23R * -G122I*E12I-G132I*E13I-G232I*E23I)*2. G2I=G112R*E11I+G222R*E22I+G332R*E33I * +G112I*E11R+G222I*E22R+G332I*E33R * +(G122R*E12I+G132R*E13I+G232R*E23I * +G122I*E12R+G132I*E13R+G232I*E23R)*2. G3R=G113R*E11R+G223R*E22R+G333R*E33R * -G113I*E11I-G223I*E22I-G333I*E33I * +(G123R*E12R+G133R*E13R+G233R*E23R * -G123I*E12I-G133I*E13I-G233I*E23I)*2. G3I=G113R*E11I+G223R*E22I+G333R*E33I * +G113I*E11R+G223I*E22R+G333I*E33R * +(G123R*E12I+G133R*E13I+G233R*E23I * +G123I*E12R+G133I*E13R+G233I*E23R)*2. G4R=G114R*E11R+G224R*E22R+G334R*E33R * -G114I*E11I-G224I*E22I-G334I*E33I * +(G124R*E12R+G134R*E13R+G234R*E23R * -G124I*E12I-G134I*E13I-G234I*E23I)*2. G4I=G114R*E11I+G224R*E22I+G334R*E33I * +G114I*E11R+G224I*E22R+G334I*E33R * +(G124R*E12I+G134R*E13I+G234R*E23I * +G124I*E12R+G134I*E13R+G234I*E23R)*2. G5R=G115R*E11R+G225R*E22R+G335R*E33R * -G115I*E11I-G225I*E22I-G335I*E33I * +(G125R*E12R+G135R*E13R+G235R*E23R * -G125I*E12I-G135I*E13I-G235I*E23I)*2. G5I=G115R*E11I+G225R*E22I+G335R*E33I * +G115I*E11R+G225I*E22R+G335I*E33R * +(G125R*E12I+G135R*E13I+G235R*E23I * +G125I*E12R+G135I*E13R+G235I*E23R)*2. G6R=G116R*E11R+G226R*E22R+G336R*E33R * -G116I*E11I-G226I*E22I-G336I*E33I * +(G126R*E12R+G136R*E13R+G236R*E23R * -G126I*E12I-G136I*E13I-G236I*E23I)*2. G6I=G116R*E11I+G226R*E22I+G336R*E33I * +G116I*E11R+G226I*E22R+G336I*E33R * +(G126R*E12I+G136R*E13I+G236R*E23I * +G126I*E12R+G136I*E13R+G236I*E23R)*2. C RETURN C C----------------------------------------------------------------------- C ENTRY CHR11( * E11R,E11I,E12R,E12I,E22R,E22I,E13R,E13I,E23R,E23I,E33R,E33I, * F11R,F11I,F12R,F12I,F22R,F22I,F13R,F13I,F23R,F23I,F33R,F33I, * G11R,G11I,G12R,G12I,G22R,G22I,G13R,G13I,G23R,G23I,G33R,G33I, * G14R,G14I,G24R,G24I,G34R,G34I,G44R,G44I, * G15R,G15I,G25R,G25I,G35R,G35I,G45R,G45I,G55R,G55I, * G16R,G16I,G26R,G26I,G36R,G36I,G46R,G46I,G56R,G56I,G66R,G66I) C C Entry calculating products C GG'"=Gij'*Eik*Fjl*Gkl" C where Gij' and Gij" are the first-order phase-space derivatives of C symmetric Christoffel matrix Gij stored by subroutine CHR1, C and Eik and Fik are given symmetric matrices. C Subroutine CHR11 increases the input values G'" by products GG'". C C Wijkl=Eik*Fjl+Eil*Fjk+Ejk*Fil+Ejl*Fik, C Wiikl=Eik*Fil+Eil*Fik, Wijkk=Eik*Fjk+Ejk*Fik, Wiikk=Eik*Fik, W1111R=E11R*F11R * -E11I*F11I W1112R=E11R*F12R+E12R*F11R * -E11I*F12I-E12I*F11I W1212R=E11R*F22R+E22R*F11R+E12R*F12R*2. * -E11I*F22I-E22I*F11I-E12I*F12I*2. W1122R=E12R*F12R * -E12I*F12I W1222R=E12R*F22R+E22R*F12R * -E12I*F22I-E22I*F12I W2222R=E22R*F22R * -E22I*F22I W1113R=E11R*F13R+E13R*F11R * -E11I*F13I-E13I*F11I W1213R=E11R*F23R+E23R*F11R+E12R*F13R+E13R*F12R * -E11I*F23I-E23I*F11I-E12I*F13I-E13I*F12I W2213R=E12R*F23R+E23R*F12R * -E12I*F23I-E23I*F12I W1313R=E11R*F33R+E33R*F11R+E13R*F13R*2. * -E11I*F33I-E33I*F11I-E13I*F13I*2. W1123R=E12R*F13R+E13R*F12R * -E12I*F13I-E13I*F12I W1223R=E12R*F23R+E23R*F12R+E13R*F22R+E22R*F13R * -E12I*F23I-E23I*F12I-E13I*F22I-E22I*F13I W2223R=E22R*F23R+E23R*F22R * -E22I*F23I-E23I*F22I W1323R=E12R*F33R+E33R*F12R+E13R*F23R+E23R*F13R * -E12I*F33I-E33I*F12I-E13I*F23I-E23I*F13I W2323R=E22R*F33R+E33R*F22R+E23R*F23R*2. * -E22I*F33I-E33I*F22I-E23I*F23I*2. W1133R=E13R*F13R * -E13I*F13I W1233R=E13R*F23R+E23R*F13R * -E13I*F23I-E23I*F13I W2233R=E23R*F23R * -E23I*F23I W1333R=E13R*F33R+E33R*F13R * -E13I*F33I-E33I*F13I W2333R=E23R*F33R+E33R*F23R * -E23I*F33I-E33I*F23I W3333R=E33R*F33R * -E33I*F33I W1111I=E11R*F11I * +E11I*F11R W1112I=E11R*F12I+E12R*F11I * +E11I*F12R+E12I*F11R W1212I=E11R*F22I+E22R*F11I+E12R*F12I*2. * +E11I*F22R+E22I*F11R+E12I*F12R*2. W1122I=E12R*F12I * +E12I*F12R W1222I=E12R*F22I+E22R*F12I * +E12I*F22R+E22I*F12R W2222I=E22R*F22I * +E22I*F22R W1113I=E11R*F13I+E13R*F11I * +E11I*F13R+E13I*F11R W1213I=E11R*F23I+E23R*F11I+E12R*F13I+E13R*F12I * +E11I*F23R+E23I*F11R+E12I*F13R+E13I*F12R W2213I=E12R*F23I+E23R*F12I * +E12I*F23R+E23I*F12R W1313I=E11R*F33I+E33R*F11I+E13R*F13I*2. * +E11I*F33R+E33I*F11R+E13I*F13R*2. W1123I=E12R*F13I+E13R*F12I * +E12I*F13R+E13I*F12R W1223I=E12R*F23I+E23R*F12I+E13R*F22I+E22R*F13I * +E12I*F23R+E23I*F12R+E13I*F22R+E22I*F13R W2223I=E22R*F23I+E23R*F22I * +E22I*F23R+E23I*F22R W1323I=E12R*F33I+E33R*F12I+E13R*F23I+E23R*F13I * +E12I*F33R+E33I*F12R+E13I*F23R+E23I*F13R W2323I=E22R*F33I+E33R*F22I+E23R*F23I*2. * +E22I*F33R+E33I*F22R+E23I*F23R*2. W1133I=E13R*F13I * +E13I*F13R W1233I=E13R*F23I+E23R*F13I * +E13I*F23R+E23I*F13R W2233I=E23R*F23I * +E23I*F23R W1333I=E13R*F33I+E33R*F13I * +E13I*F33R+E33I*F13R W2333I=E23R*F33I+E33R*F23I * +E23I*F33R+E33I*F23R W3333I=E33R*F33I * +E33I*F33R C C Xij"=Wijkl*Gkl" X111R=W1111R*G111R+W1112R*G121R+W1122R*G221R * +W1113R*G131R+W1123R*G231R+W1133R*G331R * -W1111I*G111I-W1112I*G121I-W1122I*G221I * -W1113I*G131I-W1123I*G231I-W1133I*G331I X121R=W1211R*G111R+W1212R*G121R+W1222R*G221R * +W1213R*G131R+W1223R*G231R+W1233R*G331R * -W1211I*G111I-W1212I*G121I-W1222I*G221I * -W1213I*G131I-W1223I*G231I-W1233I*G331I X221R=W2211R*G111R+W2212R*G121R+W2222R*G221R * +W2213R*G131R+W2223R*G231R+W2233R*G331R * -W2211I*G111I-W2212I*G121I-W2222I*G221I * -W2213I*G131I-W2223I*G231I-W2233I*G331I X131R=W1311R*G111R+W1312R*G121R+W1322R*G221R * +W1313R*G131R+W1323R*G231R+W1333R*G331R * -W1311I*G111I-W1312I*G121I-W1322I*G221I * -W1313I*G131I-W1323I*G231I-W1333I*G331I X231R=W2311R*G111R+W2312R*G121R+W2322R*G221R * +W2313R*G131R+W2323R*G231R+W2333R*G331R * -W2311I*G111I-W2312I*G121I-W2322I*G221I * -W2313I*G131I-W2323I*G231I-W2333I*G331I X331R=W3311R*G111R+W3312R*G121R+W3322R*G221R * +W3313R*G131R+W3323R*G231R+W3333R*G331R * -W3311I*G111I-W3312I*G121I-W3322I*G221I * -W3313I*G131I-W3323I*G231I-W3333I*G331I X112R=W1111R*G112R+W1112R*G122R+W1122R*G222R * +W1113R*G132R+W1123R*G232R+W1133R*G332R * -W1111I*G112I-W1112I*G122I-W1122I*G222I * -W1113I*G132I-W1123I*G232I-W1133I*G332I X122R=W1211R*G112R+W1212R*G122R+W1222R*G222R * +W1213R*G132R+W1223R*G232R+W1233R*G332R * -W1211I*G112I-W1212I*G122I-W1222I*G222I * -W1213I*G132I-W1223I*G232I-W1233I*G332I X222R=W2211R*G112R+W2212R*G122R+W2222R*G222R * +W2213R*G132R+W2223R*G232R+W2233R*G332R * -W2211I*G112I-W2212I*G122I-W2222I*G222I * -W2213I*G132I-W2223I*G232I-W2233I*G332I X132R=W1311R*G112R+W1312R*G122R+W1322R*G222R * +W1313R*G132R+W1323R*G232R+W1333R*G332R * -W1311I*G112I-W1312I*G122I-W1322I*G222I * -W1313I*G132I-W1323I*G232I-W1333I*G332I X232R=W2311R*G112R+W2312R*G122R+W2322R*G222R * +W2313R*G132R+W2323R*G232R+W2333R*G332R * -W2311I*G112I-W2312I*G122I-W2322I*G222I * -W2313I*G132I-W2323I*G232I-W2333I*G332I X332R=W3311R*G112R+W3312R*G122R+W3322R*G222R * +W3313R*G132R+W3323R*G232R+W3333R*G332R * -W3311I*G112I-W3312I*G122I-W3322I*G222I * -W3313I*G132I-W3323I*G232I-W3333I*G332I X113R=W1111R*G113R+W1112R*G123R+W1122R*G223R * +W1113R*G133R+W1123R*G233R+W1133R*G333R * -W1111I*G113I-W1112I*G123I-W1122I*G223I * -W1113I*G133I-W1123I*G233I-W1133I*G333I X123R=W1211R*G113R+W1212R*G123R+W1222R*G223R * +W1213R*G133R+W1223R*G233R+W1233R*G333R * -W1211I*G113I-W1212I*G123I-W1222I*G223I * -W1213I*G133I-W1223I*G233I-W1233I*G333I X223R=W2211R*G113R+W2212R*G123R+W2222R*G223R * +W2213R*G133R+W2223R*G233R+W2233R*G333R * -W2211I*G113I-W2212I*G123I-W2222I*G223I * -W2213I*G133I-W2223I*G233I-W2233I*G333I X133R=W1311R*G113R+W1312R*G123R+W1322R*G223R * +W1313R*G133R+W1323R*G233R+W1333R*G333R * -W1311I*G113I-W1312I*G123I-W1322I*G223I * -W1313I*G133I-W1323I*G233I-W1333I*G333I X233R=W2311R*G113R+W2312R*G123R+W2322R*G223R * +W2313R*G133R+W2323R*G233R+W2333R*G333R * -W2311I*G113I-W2312I*G123I-W2322I*G223I * -W2313I*G133I-W2323I*G233I-W2333I*G333I X333R=W3311R*G113R+W3312R*G123R+W3322R*G223R * +W3313R*G133R+W3323R*G233R+W3333R*G333R * -W3311I*G113I-W3312I*G123I-W3322I*G223I * -W3313I*G133I-W3323I*G233I-W3333I*G333I X114R=W1111R*G114R+W1112R*G124R+W1122R*G224R * +W1113R*G134R+W1123R*G234R+W1133R*G334R * -W1111I*G114I-W1112I*G124I-W1122I*G224I * -W1113I*G134I-W1123I*G234I-W1133I*G334I X124R=W1211R*G114R+W1212R*G124R+W1222R*G224R * +W1213R*G134R+W1223R*G234R+W1233R*G334R * -W1211I*G114I-W1212I*G124I-W1222I*G224I * -W1213I*G134I-W1223I*G234I-W1233I*G334I X224R=W2211R*G114R+W2212R*G124R+W2222R*G224R * +W2213R*G134R+W2223R*G234R+W2233R*G334R * -W2211I*G114I-W2212I*G124I-W2222I*G224I * -W2213I*G134I-W2223I*G234I-W2233I*G334I X134R=W1311R*G114R+W1312R*G124R+W1322R*G224R * +W1313R*G134R+W1323R*G234R+W1333R*G334R * -W1311I*G114I-W1312I*G124I-W1322I*G224I * -W1313I*G134I-W1323I*G234I-W1333I*G334I X234R=W2311R*G114R+W2312R*G124R+W2322R*G224R * +W2313R*G134R+W2323R*G234R+W2333R*G334R * -W2311I*G114I-W2312I*G124I-W2322I*G224I * -W2313I*G134I-W2323I*G234I-W2333I*G334I X334R=W3311R*G114R+W3312R*G124R+W3322R*G224R * +W3313R*G134R+W3323R*G234R+W3333R*G334R * -W3311I*G114I-W3312I*G124I-W3322I*G224I * -W3313I*G134I-W3323I*G234I-W3333I*G334I X115R=W1111R*G115R+W1112R*G125R+W1122R*G225R * +W1113R*G135R+W1123R*G235R+W1133R*G335R * -W1111I*G115I-W1112I*G125I-W1122I*G225I * -W1113I*G135I-W1123I*G235I-W1133I*G335I X125R=W1211R*G115R+W1212R*G125R+W1222R*G225R * +W1213R*G135R+W1223R*G235R+W1233R*G335R * -W1211I*G115I-W1212I*G125I-W1222I*G225I * -W1213I*G135I-W1223I*G235I-W1233I*G335I X225R=W2211R*G115R+W2212R*G125R+W2222R*G225R * +W2213R*G135R+W2223R*G235R+W2233R*G335R * -W2211I*G115I-W2212I*G125I-W2222I*G225I * -W2213I*G135I-W2223I*G235I-W2233I*G335I X135R=W1311R*G115R+W1312R*G125R+W1322R*G225R * +W1313R*G135R+W1323R*G235R+W1333R*G335R * -W1311I*G115I-W1312I*G125I-W1322I*G225I * -W1313I*G135I-W1323I*G235I-W1333I*G335I X235R=W2311R*G115R+W2312R*G125R+W2322R*G225R * +W2313R*G135R+W2323R*G235R+W2333R*G335R * -W2311I*G115I-W2312I*G125I-W2322I*G225I * -W2313I*G135I-W2323I*G235I-W2333I*G335I X335R=W3311R*G115R+W3312R*G125R+W3322R*G225R * +W3313R*G135R+W3323R*G235R+W3333R*G335R * -W3311I*G115I-W3312I*G125I-W3322I*G225I * -W3313I*G135I-W3323I*G235I-W3333I*G335I X116R=W1111R*G116R+W1112R*G126R+W1122R*G226R * +W1113R*G136R+W1123R*G236R+W1133R*G336R * -W1111I*G116I-W1112I*G126I-W1122I*G226I * -W1113I*G136I-W1123I*G236I-W1133I*G336I X126R=W1211R*G116R+W1212R*G126R+W1222R*G226R * +W1213R*G136R+W1223R*G236R+W1233R*G336R * -W1211I*G116I-W1212I*G126I-W1222I*G226I * -W1213I*G136I-W1223I*G236I-W1233I*G336I X226R=W2211R*G116R+W2212R*G126R+W2222R*G226R * +W2213R*G136R+W2223R*G236R+W2233R*G336R * -W2211I*G116I-W2212I*G126I-W2222I*G226I * -W2213I*G136I-W2223I*G236I-W2233I*G336I X136R=W1311R*G116R+W1312R*G126R+W1322R*G226R * +W1313R*G136R+W1323R*G236R+W1333R*G336R * -W1311I*G116I-W1312I*G126I-W1322I*G226I * -W1313I*G136I-W1323I*G236I-W1333I*G336I X236R=W2311R*G116R+W2312R*G126R+W2322R*G226R * +W2313R*G136R+W2323R*G236R+W2333R*G336R * -W2311I*G116I-W2312I*G126I-W2322I*G226I * -W2313I*G136I-W2323I*G236I-W2333I*G336I X336R=W3311R*G116R+W3312R*G126R+W3322R*G226R * +W3313R*G136R+W3323R*G236R+W3333R*G336R * -W3311I*G116I-W3312I*G126I-W3322I*G226I * -W3313I*G136I-W3323I*G236I-W3333I*G336I X111I=W1111R*G111I+W1112R*G121I+W1122R*G221I * +W1113R*G131I+W1123R*G231I+W1133R*G331I * +W1111I*G111R+W1112I*G121R+W1122I*G221R * +W1113I*G131R+W1123I*G231R+W1133I*G331R X121I=W1211R*G111I+W1212R*G121I+W1222R*G221I * +W1213R*G131I+W1223R*G231I+W1233R*G331I * +W1211I*G111R+W1212I*G121R+W1222I*G221R * +W1213I*G131R+W1223I*G231R+W1233I*G331R X221I=W2211R*G111I+W2212R*G121I+W2222R*G221I * +W2213R*G131I+W2223R*G231I+W2233R*G331I * +W2211I*G111R+W2212I*G121R+W2222I*G221R * +W2213I*G131R+W2223I*G231R+W2233I*G331R X131I=W1311R*G111I+W1312R*G121I+W1322R*G221I * +W1313R*G131I+W1323R*G231I+W1333R*G331I * +W1311I*G111R+W1312I*G121R+W1322I*G221R * +W1313I*G131R+W1323I*G231R+W1333I*G331R X231I=W2311R*G111I+W2312R*G121I+W2322R*G221I * +W2313R*G131I+W2323R*G231I+W2333R*G331I * +W2311I*G111R+W2312I*G121R+W2322I*G221R * +W2313I*G131R+W2323I*G231R+W2333I*G331R X331I=W3311R*G111I+W3312R*G121I+W3322R*G221I * +W3313R*G131I+W3323R*G231I+W3333R*G331I * +W3311I*G111R+W3312I*G121R+W3322I*G221R * +W3313I*G131R+W3323I*G231R+W3333I*G331R X112I=W1111R*G112I+W1112R*G122I+W1122R*G222I * +W1113R*G132I+W1123R*G232I+W1133R*G332I * +W1111I*G112R+W1112I*G122R+W1122I*G222R * +W1113I*G132R+W1123I*G232R+W1133I*G332R X122I=W1211R*G112I+W1212R*G122I+W1222R*G222I * +W1213R*G132I+W1223R*G232I+W1233R*G332I * +W1211I*G112R+W1212I*G122R+W1222I*G222R * +W1213I*G132R+W1223I*G232R+W1233I*G332R X222I=W2211R*G112I+W2212R*G122I+W2222R*G222I * +W2213R*G132I+W2223R*G232I+W2233R*G332I * +W2211I*G112R+W2212I*G122R+W2222I*G222R * +W2213I*G132R+W2223I*G232R+W2233I*G332R X132I=W1311R*G112I+W1312R*G122I+W1322R*G222I * +W1313R*G132I+W1323R*G232I+W1333R*G332I * +W1311I*G112R+W1312I*G122R+W1322I*G222R * +W1313I*G132R+W1323I*G232R+W1333I*G332R X232I=W2311R*G112I+W2312R*G122I+W2322R*G222I * +W2313R*G132I+W2323R*G232I+W2333R*G332I * +W2311I*G112R+W2312I*G122R+W2322I*G222R * +W2313I*G132R+W2323I*G232R+W2333I*G332R X332I=W3311R*G112I+W3312R*G122I+W3322R*G222I * +W3313R*G132I+W3323R*G232I+W3333R*G332I * +W3311I*G112R+W3312I*G122R+W3322I*G222R * +W3313I*G132R+W3323I*G232R+W3333I*G332R X113I=W1111R*G113I+W1112R*G123I+W1122R*G223I * +W1113R*G133I+W1123R*G233I+W1133R*G333I * +W1111I*G113R+W1112I*G123R+W1122I*G223R * +W1113I*G133R+W1123I*G233R+W1133I*G333R X123I=W1211R*G113I+W1212R*G123I+W1222R*G223I * +W1213R*G133I+W1223R*G233I+W1233R*G333I * +W1211I*G113R+W1212I*G123R+W1222I*G223R * +W1213I*G133R+W1223I*G233R+W1233I*G333R X223I=W2211R*G113I+W2212R*G123I+W2222R*G223I * +W2213R*G133I+W2223R*G233I+W2233R*G333I * +W2211I*G113R+W2212I*G123R+W2222I*G223R * +W2213I*G133R+W2223I*G233R+W2233I*G333R X133I=W1311R*G113I+W1312R*G123I+W1322R*G223I * +W1313R*G133I+W1323R*G233I+W1333R*G333I * +W1311I*G113R+W1312I*G123R+W1322I*G223R * +W1313I*G133R+W1323I*G233R+W1333I*G333R X233I=W2311R*G113I+W2312R*G123I+W2322R*G223I * +W2313R*G133I+W2323R*G233I+W2333R*G333I * +W2311I*G113R+W2312I*G123R+W2322I*G223R * +W2313I*G133R+W2323I*G233R+W2333I*G333R X333I=W3311R*G113I+W3312R*G123I+W3322R*G223I * +W3313R*G133I+W3323R*G233I+W3333R*G333I * +W3311I*G113R+W3312I*G123R+W3322I*G223R * +W3313I*G133R+W3323I*G233R+W3333I*G333R X114I=W1111R*G114I+W1112R*G124I+W1122R*G224I * +W1113R*G134I+W1123R*G234I+W1133R*G334I * +W1111I*G114R+W1112I*G124R+W1122I*G224R * +W1113I*G134R+W1123I*G234R+W1133I*G334R X124I=W1211R*G114I+W1212R*G124I+W1222R*G224I * +W1213R*G134I+W1223R*G234I+W1233R*G334I * +W1211I*G114R+W1212I*G124R+W1222I*G224R * +W1213I*G134R+W1223I*G234R+W1233I*G334R X224I=W2211R*G114I+W2212R*G124I+W2222R*G224I * +W2213R*G134I+W2223R*G234I+W2233R*G334I * +W2211I*G114R+W2212I*G124R+W2222I*G224R * +W2213I*G134R+W2223I*G234R+W2233I*G334R X134I=W1311R*G114I+W1312R*G124I+W1322R*G224I * +W1313R*G134I+W1323R*G234I+W1333R*G334I * +W1311I*G114R+W1312I*G124R+W1322I*G224R * +W1313I*G134R+W1323I*G234R+W1333I*G334R X234I=W2311R*G114I+W2312R*G124I+W2322R*G224I * +W2313R*G134I+W2323R*G234I+W2333R*G334I * +W2311I*G114R+W2312I*G124R+W2322I*G224R * +W2313I*G134R+W2323I*G234R+W2333I*G334R X334I=W3311R*G114I+W3312R*G124I+W3322R*G224I * +W3313R*G134I+W3323R*G234I+W3333R*G334I * +W3311I*G114R+W3312I*G124R+W3322I*G224R * +W3313I*G134R+W3323I*G234R+W3333I*G334R X115I=W1111R*G115I+W1112R*G125I+W1122R*G225I * +W1113R*G135I+W1123R*G235I+W1133R*G335I * +W1111I*G115R+W1112I*G125R+W1122I*G225R * +W1113I*G135R+W1123I*G235R+W1133I*G335R X125I=W1211R*G115I+W1212R*G125I+W1222R*G225I * +W1213R*G135I+W1223R*G235I+W1233R*G335I * +W1211I*G115R+W1212I*G125R+W1222I*G225R * +W1213I*G135R+W1223I*G235R+W1233I*G335R X225I=W2211R*G115I+W2212R*G125I+W2222R*G225I * +W2213R*G135I+W2223R*G235I+W2233R*G335I * +W2211I*G115R+W2212I*G125R+W2222I*G225R * +W2213I*G135R+W2223I*G235R+W2233I*G335R X135I=W1311R*G115I+W1312R*G125I+W1322R*G225I * +W1313R*G135I+W1323R*G235I+W1333R*G335I * +W1311I*G115R+W1312I*G125R+W1322I*G225R * +W1313I*G135R+W1323I*G235R+W1333I*G335R X235I=W2311R*G115I+W2312R*G125I+W2322R*G225I * +W2313R*G135I+W2323R*G235I+W2333R*G335I * +W2311I*G115R+W2312I*G125R+W2322I*G225R * +W2313I*G135R+W2323I*G235R+W2333I*G335R X335I=W3311R*G115I+W3312R*G125I+W3322R*G225I * +W3313R*G135I+W3323R*G235I+W3333R*G335I * +W3311I*G115R+W3312I*G125R+W3322I*G225R * +W3313I*G135R+W3323I*G235R+W3333I*G335R X116I=W1111R*G116I+W1112R*G126I+W1122R*G226I * +W1113R*G136I+W1123R*G236I+W1133R*G336I * +W1111I*G116R+W1112I*G126R+W1122I*G226R * +W1113I*G136R+W1123I*G236R+W1133I*G336R X126I=W1211R*G116I+W1212R*G126I+W1222R*G226I * +W1213R*G136I+W1223R*G236I+W1233R*G336I * +W1211I*G116R+W1212I*G126R+W1222I*G226R * +W1213I*G136R+W1223I*G236R+W1233I*G336R X226I=W2211R*G116I+W2212R*G126I+W2222R*G226I * +W2213R*G136I+W2223R*G236I+W2233R*G336I * +W2211I*G116R+W2212I*G126R+W2222I*G226R * +W2213I*G136R+W2223I*G236R+W2233I*G336R X136I=W1311R*G116I+W1312R*G126I+W1322R*G226I * +W1313R*G136I+W1323R*G236I+W1333R*G336I * +W1311I*G116R+W1312I*G126R+W1322I*G226R * +W1313I*G136R+W1323I*G236R+W1333I*G336R X236I=W2311R*G116I+W2312R*G126I+W2322R*G226I * +W2313R*G136I+W2323R*G236I+W2333R*G336I * +W2311I*G116R+W2312I*G126R+W2322I*G226R * +W2313I*G136R+W2323I*G236R+W2333I*G336R X336I=W3311R*G116I+W3312R*G126I+W3322R*G226I * +W3313R*G136I+W3323R*G236I+W3333R*G336I * +W3311I*G116R+W3312I*G126R+W3322I*G226R * +W3313I*G136R+W3323I*G236R+W3333I*G336R C C G'"=Gab+Gij'*Xij"=Gab+Gij'*Wijkl*Gkl" G11R=G11R+G111R*X111R+G121R*X121R+G221R*X221R * +G131R*X131R+G231R*X231R+G331R*X331R * -G111I*X111I-G121I*X121I-G221I*X221I * -G131I*X131I-G231I*X231I-G331I*X331I G12R=G12R+G111R*X112R+G121R*X122R+G221R*X222R * +G131R*X132R+G231R*X232R+G331R*X332R * -G111I*X112I-G121I*X122I-G221I*X222I * -G131I*X132I-G231I*X232I-G331I*X332I G22R=G22R+G112R*X112R+G122R*X122R+G222R*X222R * +G132R*X132R+G232R*X232R+G332R*X332R * -G112I*X112I-G122I*X122I-G222I*X222I * -G132I*X132I-G232I*X232I-G332I*X332I G13R=G13R+G111R*X113R+G121R*X123R+G221R*X223R * +G131R*X133R+G231R*X233R+G331R*X333R * -G111I*X113I-G121I*X123I-G221I*X223I * -G131I*X133I-G231I*X233I-G331I*X333I G23R=G23R+G112R*X113R+G122R*X123R+G222R*X223R * +G132R*X133R+G232R*X233R+G332R*X333R * -G112I*X113I-G122I*X123I-G222I*X223I * -G132I*X133I-G232I*X233I-G332I*X333I G33R=G33R+G113R*X113R+G123R*X123R+G223R*X223R * +G133R*X133R+G233R*X233R+G333R*X333R * -G113I*X113I-G123I*X123I-G223I*X223I * -G133I*X133I-G233I*X233I-G333I*X333I G14R=G14R+G111R*X114R+G121R*X124R+G221R*X224R * +G131R*X134R+G231R*X234R+G331R*X334R * -G111I*X114I-G121I*X124I-G221I*X224I * -G131I*X134I-G231I*X234I-G331I*X334I G24R=G24R+G112R*X114R+G122R*X124R+G222R*X224R * +G132R*X134R+G232R*X234R+G332R*X334R * -G112I*X114I-G122I*X124I-G222I*X224I * -G132I*X134I-G232I*X234I-G332I*X334I G34R=G34R+G113R*X114R+G123R*X124R+G223R*X224R * +G133R*X134R+G233R*X234R+G333R*X334R * -G113I*X114I-G123I*X124I-G223I*X224I * -G133I*X134I-G233I*X234I-G333I*X334I G44R=G44R+G114R*X114R+G124R*X124R+G224R*X224R * +G134R*X134R+G234R*X234R+G334R*X334R * -G114I*X114I-G124I*X124I-G224I*X224I * -G134I*X134I-G234I*X234I-G334I*X334I G15R=G15R+G111R*X115R+G121R*X125R+G221R*X225R * +G131R*X135R+G231R*X235R+G331R*X335R * -G111I*X115I-G121I*X125I-G221I*X225I * -G131I*X135I-G231I*X235I-G331I*X335I G25R=G25R+G112R*X115R+G122R*X125R+G222R*X225R * +G132R*X135R+G232R*X235R+G332R*X335R * -G112I*X115I-G122I*X125I-G222I*X225I * -G132I*X135I-G232I*X235I-G332I*X335I G35R=G35R+G113R*X115R+G123R*X125R+G223R*X225R * +G133R*X135R+G233R*X235R+G333R*X335R * -G113I*X115I-G123I*X125I-G223I*X225I * -G133I*X135I-G233I*X235I-G333I*X335I G45R=G45R+G114R*X115R+G124R*X125R+G224R*X225R * +G134R*X135R+G234R*X235R+G334R*X335R * -G114I*X115I-G124I*X125I-G224I*X225I * -G134I*X135I-G234I*X235I-G334I*X335I G55R=G55R+G115R*X115R+G125R*X125R+G225R*X225R * +G135R*X135R+G235R*X235R+G335R*X335R * -G115I*X115I-G125I*X125I-G225I*X225I * -G135I*X135I-G235I*X235I-G335I*X335I G16R=G16R+G111R*X116R+G121R*X126R+G221R*X226R * +G131R*X136R+G231R*X236R+G331R*X336R * -G111I*X116I-G121I*X126I-G221I*X226I * -G131I*X136I-G231I*X236I-G331I*X336I G26R=G26R+G112R*X116R+G122R*X126R+G222R*X226R * +G132R*X136R+G232R*X236R+G332R*X336R * -G112I*X116I-G122I*X126I-G222I*X226I * -G132I*X136I-G232I*X236I-G332I*X336I G36R=G36R+G113R*X116R+G123R*X126R+G223R*X226R * +G133R*X136R+G233R*X236R+G333R*X336R * -G113I*X116I-G123I*X126I-G223I*X226I * -G133I*X136I-G233I*X236I-G333I*X336I G46R=G46R+G114R*X116R+G124R*X126R+G224R*X226R * +G134R*X136R+G234R*X236R+G334R*X336R * -G114I*X116I-G124I*X126I-G224I*X226I * -G134I*X136I-G234I*X236I-G334I*X336I G56R=G56R+G115R*X116R+G125R*X126R+G225R*X226R * +G135R*X136R+G235R*X236R+G335R*X336R * -G115I*X116I-G125I*X126I-G225I*X226I * -G135I*X136I-G235I*X236I-G335I*X336I G66R=G66R+G116R*X116R+G126R*X126R+G226R*X226R * +G136R*X136R+G236R*X236R+G336R*X336R * -G116I*X116I-G126I*X126I-G226I*X226I * -G136I*X136I-G236I*X236I-G336I*X336I G11I=G11I+G111R*X111I+G121R*X121I+G221R*X221I * +G131R*X131I+G231R*X231I+G331R*X331I * +G111I*X111R+G121I*X121R+G221I*X221R * +G131I*X131R+G231I*X231R+G331I*X331R G12I=G12I+G111R*X112I+G121R*X122I+G221R*X222I * +G131R*X132I+G231R*X232I+G331R*X332I * +G111I*X112R+G121I*X122R+G221I*X222R * +G131I*X132R+G231I*X232R+G331I*X332R G22I=G22I+G112R*X112I+G122R*X122I+G222R*X222I * +G132R*X132I+G232R*X232I+G332R*X332I * +G112I*X112R+G122I*X122R+G222I*X222R * +G132I*X132R+G232I*X232R+G332I*X332R G13I=G13I+G111R*X113I+G121R*X123I+G221R*X223I * +G131R*X133I+G231R*X233I+G331R*X333I * +G111I*X113R+G121I*X123R+G221I*X223R * +G131I*X133R+G231I*X233R+G331I*X333R G23I=G23I+G112R*X113I+G122R*X123I+G222R*X223I * +G132R*X133I+G232R*X233I+G332R*X333I * +G112I*X113R+G122I*X123R+G222I*X223R * +G132I*X133R+G232I*X233R+G332I*X333R G33I=G33I+G113R*X113I+G123R*X123I+G223R*X223I * +G133R*X133I+G233R*X233I+G333R*X333I * +G113I*X113R+G123I*X123R+G223I*X223R * +G133I*X133R+G233I*X233R+G333I*X333R G14I=G14I+G111R*X114I+G121R*X124I+G221R*X224I * +G131R*X134I+G231R*X234I+G331R*X334I * +G111I*X114R+G121I*X124R+G221I*X224R * +G131I*X134R+G231I*X234R+G331I*X334R G24I=G24I+G112R*X114I+G122R*X124I+G222R*X224I * +G132R*X134I+G232R*X234I+G332R*X334I * +G112I*X114R+G122I*X124R+G222I*X224R * +G132I*X134R+G232I*X234R+G332I*X334R G34I=G34I+G113R*X114I+G123R*X124I+G223R*X224I * +G133R*X134I+G233R*X234I+G333R*X334I * +G113I*X114R+G123I*X124R+G223I*X224R * +G133I*X134R+G233I*X234R+G333I*X334R G44I=G44I+G114R*X114I+G124R*X124I+G224R*X224I * +G134R*X134I+G234R*X234I+G334R*X334I * +G114I*X114R+G124I*X124R+G224I*X224R * +G134I*X134R+G234I*X234R+G334I*X334R G15I=G15I+G111R*X115I+G121R*X125I+G221R*X225I * +G131R*X135I+G231R*X235I+G331R*X335I * +G111I*X115R+G121I*X125R+G221I*X225R * +G131I*X135R+G231I*X235R+G331I*X335R G25I=G25I+G112R*X115I+G122R*X125I+G222R*X225I * +G132R*X135I+G232R*X235I+G332R*X335I * +G112I*X115R+G122I*X125R+G222I*X225R * +G132I*X135R+G232I*X235R+G332I*X335R G35I=G35I+G113R*X115I+G123R*X125I+G223R*X225I * +G133R*X135I+G233R*X235I+G333R*X335I * +G113I*X115R+G123I*X125R+G223I*X225R * +G133I*X135R+G233I*X235R+G333I*X335R G45I=G45I+G114R*X115I+G124R*X125I+G224R*X225I * +G134R*X135I+G234R*X235I+G334R*X335I * +G114I*X115R+G124I*X125R+G224I*X225R * +G134I*X135R+G234I*X235R+G334I*X335R G55I=G55I+G115R*X115I+G125R*X125I+G225R*X225I * +G135R*X135I+G235R*X235I+G335R*X335I * +G115I*X115R+G125I*X125R+G225I*X225R * +G135I*X135R+G235I*X235R+G335I*X335R G16I=G16I+G111R*X116I+G121R*X126I+G221R*X226I * +G131R*X136I+G231R*X236I+G331R*X336I * +G111I*X116R+G121I*X126R+G221I*X226R * +G131I*X136R+G231I*X236R+G331I*X336R G26I=G26I+G112R*X116I+G122R*X126I+G222R*X226I * +G132R*X136I+G232R*X236I+G332R*X336I * +G112I*X116R+G122I*X126R+G222I*X226R * +G132I*X136R+G232I*X236R+G332I*X336R G36I=G36I+G113R*X116I+G123R*X126I+G223R*X226I * +G133R*X136I+G233R*X236I+G333R*X336I * +G113I*X116R+G123I*X126R+G223I*X226R * +G133I*X136R+G233I*X236R+G333I*X336R G46I=G46I+G114R*X116I+G124R*X126I+G224R*X226I * +G134R*X136I+G234R*X236I+G334R*X336I * +G114I*X116R+G124I*X126R+G224I*X226R * +G134I*X136R+G234I*X236R+G334I*X336R G56I=G56I+G115R*X116I+G125R*X126I+G225R*X226I * +G135R*X136I+G235R*X236I+G335R*X336I * +G115I*X116R+G125I*X126R+G225I*X226R * +G135I*X136R+G235I*X236R+G335I*X336R G66I=G66I+G116R*X116I+G126R*X126I+G226R*X226I * +G136R*X136I+G236R*X236I+G336R*X336I * +G116I*X116R+G126I*X126R+G226I*X226R * +G136I*X136R+G236I*X236R+G336I*X336R C RETURN C C----------------------------------------------------------------------- C ENTRY CHR2(A,Q,P1,P2,P3, * E11R,E11I,E12R,E12I,E22R,E22I,E13R,E13I,E23R,E23I,E33R,E33I, * G11R,G11I,G12R,G12I,G22R,G22I,G13R,G13I,G23R,G23I,G33R,G33I, * G14R,G14I,G24R,G24I,G34R,G34I,G44R,G44I, * G15R,G15I,G25R,G25I,G35R,G35I,G45R,G45I,G55R,G55I, * G16R,G16I,G26R,G26I,G36R,G36I,G46R,G46I,G56R,G56I,G66R,G66I) C C Entry calculating the second-order phase-space derivatives C G'"ij of the Christoffel matrix Gij, C and returning the products G'"=Gij'"*Eij of these second-order C phase-space derivatives with given symmetric matrix Eik. C C Spatial derivatives of the Christoffel matrix. C W11R=E11R*P11 W12R=E12R*P12*2. W22R=E22R*P22 W13R=E13R*P13*2. W23R=E23R*P23*2. W33R=E33R*P33 W14R=(E12R*P13+E13R*P12)*2. W24R=(E22R*P23+E23R*P22)*2. W34R=(E23R*P33+E33R*P23)*2. W44R=E22R*P33+E33R*P22+E23R*P23*2. W15R=(E11R*P13+E13R*P11)*2. W25R=(E12R*P23+E23R*P12)*2. W35R=(E13R*P33+E33R*P13)*2. W45R=(E12R*P33+E13R*P23+E23R*P13+E33R*P12)*2. W55R=E11R*P33+E33R*P11+E13R*P13*2. W16R=(E11R*P12+E12R*P11)*2. W26R=(E12R*P22+E22R*P12)*2. W36R=(E13R*P23+E23R*P13)*2. W46R=(E12R*P23+E13R*P22+E22R*P13+E23R*P12)*2. W56R=(E11R*P23+E13R*P12+E12R*P13+E23R*P11)*2. W66R=E11R*P22+E22R*P11+E12R*P12*2. W11I=E11I*P11 W12I=E12I*P12*2. W22I=E22I*P22 W13I=E13I*P13*2. W23I=E23I*P23*2. W33I=E33I*P33 W14I=(E12I*P13+E13I*P12)*2. W24I=(E22I*P23+E23I*P22)*2. W34I=(E23I*P33+E33I*P23)*2. W44I=E22I*P33+E33I*P22+E23I*P23*2. W15I=(E11I*P13+E13I*P11)*2. W25I=(E12I*P23+E23I*P12)*2. W35I=(E13I*P33+E33I*P13)*2. W45I=(E12I*P33+E13I*P23+E23I*P13+E33I*P12)*2. W55I=E11I*P33+E33I*P11+E13I*P13*2. W16I=(E11I*P12+E12I*P11)*2. W26I=(E12I*P22+E22I*P12)*2. W36I=(E13I*P23+E23I*P13)*2. W46I=(E12I*P23+E13I*P22+E22I*P13+E23I*P12)*2. W56I=(E11I*P23+E13I*P12+E12I*P13+E23I*P11)*2. W66I=E11I*P22+E22I*P11+E12I*P12*2. DO 41 I=5,10 G(I)=W11R*A(I, 1)+W12R*A(I, 2)+W22R*A(I, 3)+W13R*A(I, 4) * +W23R*A(I, 5)+W33R*A(I, 6)+W14R*A(I, 7)+W24R*A(I, 8) * +W34R*A(I, 9)+W44R*A(I,10)+W15R*A(I,11)+W25R*A(I,12) * +W35R*A(I,13)+W45R*A(I,14)+W55R*A(I,15)+W16R*A(I,16) * +W26R*A(I,17)+W36R*A(I,18)+W46R*A(I,19)+W56R*A(I,20) * +W66R*A(I,21) * -W11I*Q(I, 1)-W12I*Q(I, 2)-W22I*Q(I, 3)-W13I*Q(I, 4) * -W23I*Q(I, 5)-W33I*Q(I, 6)-W14I*Q(I, 7)-W24I*Q(I, 8) * -W34I*Q(I, 9)-W44I*Q(I,10)-W15I*Q(I,11)-W25I*Q(I,12) * -W35I*Q(I,13)-W45I*Q(I,14)-W55I*Q(I,15)-W16I*Q(I,16) * -W26I*Q(I,17)-W36I*Q(I,18)-W46I*Q(I,19)-W56I*Q(I,20) * -W66I*Q(I,21) 41 CONTINUE G11R=G(5) G12R=G(6) G22R=G(7) G13R=G(8) G23R=G(9) G33R=G(10) DO 42 I=5,10 G(I)=W11R*Q(I, 1)+W12R*Q(I, 2)+W22R*Q(I, 3)+W13R*Q(I, 4) * +W23R*Q(I, 5)+W33R*Q(I, 6)+W14R*Q(I, 7)+W24R*Q(I, 8) * +W34R*Q(I, 9)+W44R*Q(I,10)+W15R*Q(I,11)+W25R*Q(I,12) * +W35R*Q(I,13)+W45R*Q(I,14)+W55R*Q(I,15)+W16R*Q(I,16) * +W26R*Q(I,17)+W36R*Q(I,18)+W46R*Q(I,19)+W56R*Q(I,20) * +W66R*Q(I,21) * +W11I*A(I, 1)+W12I*A(I, 2)+W22I*A(I, 3)+W13I*A(I, 4) * +W23I*A(I, 5)+W33I*A(I, 6)+W14I*A(I, 7)+W24I*A(I, 8) * +W34I*A(I, 9)+W44I*A(I,10)+W15I*A(I,11)+W25I*A(I,12) * +W35I*A(I,13)+W45I*A(I,14)+W55I*A(I,15)+W16I*A(I,16) * +W26I*A(I,17)+W36I*A(I,18)+W46I*A(I,19)+W56I*A(I,20) * +W66I*A(I,21) 42 CONTINUE G11I=G(5) G12I=G(6) G22I=G(7) G13I=G(8) G23I=G(9) G33I=G(10) C C....................................................................... C C Mixed derivatives of the Christoffel matrix. C C Derivative with respect to P1: W11R=E11R*P1 W12R=E12R*P2 W13R=E13R*P3 W14R=E12R*P3+E13R*P2 W15R=E11R*P3+E13R*P1*2. W25R=E23R*P2 W35R=E33R*P3 W45R=E23R*P3+E33R*P2 W55R=E33R*P1+E13R*P3 W16R=E11R*P2+E12R*P1*2. W26R=E22R*P2 W36R=E23R*P3 W46R=E22R*P3+E23R*P2 W56R=E13R*P2+E12R*P3+E23R*P1*2. W66R=E22R*P1+E12R*P2 W11I=E11I*P1 W12I=E12I*P2 W13I=E13I*P3 W14I=E12I*P3+E13I*P2 W15I=E11I*P3+E13I*P1*2. W25I=E23I*P2 W35I=E33I*P3 W45I=E23I*P3+E33I*P2 W55I=E33I*P1+E13I*P3 W16I=E11I*P2+E12I*P1*2. W26I=E22I*P2 W36I=E23I*P3 W46I=E22I*P3+E23I*P2 W56I=E13I*P2+E12I*P3+E23I*P1*2. W66I=E22I*P1+E12I*P2 DO 51 I=2,4 G(I)=W11R*A(I, 1)+W12R*A(I, 2)+W13R*A(I, 4)+W14R*A(I, 7) * +W15R*A(I,11)+W25R*A(I,12)+W35R*A(I,13)+W45R*A(I,14) * +W55R*A(I,15)+W16R*A(I,16)+W26R*A(I,17)+W36R*A(I,18) * +W46R*A(I,19)+W56R*A(I,20)+W66R*A(I,21) * -W11I*Q(I, 1)-W12I*Q(I, 2)-W13I*Q(I, 4)-W14I*Q(I, 7) * -W15I*Q(I,11)-W25I*Q(I,12)-W35I*Q(I,13)-W45I*Q(I,14) * -W55I*Q(I,15)-W16I*Q(I,16)-W26I*Q(I,17)-W36I*Q(I,18) * -W46I*Q(I,19)-W56I*Q(I,20)-W66I*Q(I,21) 51 CONTINUE G14R=2.*G(2) G24R=2.*G(3) G34R=2.*G(4) DO 52 I=2,4 G(I)=W11R*Q(I, 1)+W12R*Q(I, 2)+W13R*Q(I, 4)+W14R*Q(I, 7) * +W15R*Q(I,11)+W25R*Q(I,12)+W35R*Q(I,13)+W45R*Q(I,14) * +W55R*Q(I,15)+W16R*Q(I,16)+W26R*Q(I,17)+W36R*Q(I,18) * +W46R*Q(I,19)+W56R*Q(I,20)+W66R*Q(I,21) * +W11I*A(I, 1)+W12I*A(I, 2)+W13I*A(I, 4)+W14I*A(I, 7) * +W15I*A(I,11)+W25I*A(I,12)+W35I*A(I,13)+W45I*A(I,14) * +W55I*A(I,15)+W16I*A(I,16)+W26I*A(I,17)+W36I*A(I,18) * +W46I*A(I,19)+W56I*A(I,20)+W66I*A(I,21) 52 CONTINUE G14I=2.*G(2) G24I=2.*G(3) G34I=2.*G(4) C C Derivative with respect to P2: W12R=E12R*P1 W22R=E22R*P2 W23R=E23R*P3 W14R=E13R*P1 W24R=E22R*P3+E23R*P2*2. W34R=E33R*P3 W44R=E33R*P2+E23R*P3 W25R=E12R*P3+E23R*P1 W45R=E13R*P3+E33R*P1 W16R=E11R*P1 W26R=E12R*P2*2.+E22R*P1 W36R=E13R*P3 W46R=E12R*P3+E13R*P2*2.+E23R*P1 W56R=E11R*P3+E13R*P1 W66R=E11R*P2+E12R*P1 W12I=E12I*P1 W22I=E22I*P2 W23I=E23I*P3 W14I=E13I*P1 W24I=E22I*P3+E23I*P2*2. W34I=E33I*P3 W44I=E33I*P2+E23I*P3 W25I=E12I*P3+E23I*P1 W45I=E13I*P3+E33I*P1 W16I=E11I*P1 W26I=E12I*P2*2.+E22I*P1 W36I=E13I*P3 W46I=E12I*P3+E13I*P2*2.+E23I*P1 W56I=E11I*P3+E13I*P1 W66I=E11I*P2+E12I*P1 DO 53 I=2,4 G(I)=W12R*A(I, 2)+W22R*A(I, 3)+W23R*A(I, 5)+W14R*A(I, 7) * +W24R*A(I, 8)+W34R*A(I, 9)+W44R*A(I,10)+W25R*A(I,12) * +W45R*A(I,14)+W16R*A(I,16)+W26R*A(I,17)+W36R*A(I,18) * +W46R*A(I,19)+W56R*A(I,20)+W66R*A(I,21) * -W12I*Q(I, 2)-W22I*Q(I, 3)-W23I*Q(I, 5)-W14I*Q(I, 7) * -W24I*Q(I, 8)-W34I*Q(I, 9)-W44I*Q(I,10)-W25I*Q(I,12) * -W45I*Q(I,14)-W16I*Q(I,16)-W26I*Q(I,17)-W36I*Q(I,18) * -W46I*Q(I,19)-W56I*Q(I,20)-W66I*Q(I,21) 53 CONTINUE G15R=2.*G(2) G25R=2.*G(3) G35R=2.*G(4) DO 54 I=2,4 G(I)=W12R*Q(I, 2)+W22R*Q(I, 3)+W23R*Q(I, 5)+W14R*Q(I, 7) * +W24R*Q(I, 8)+W34R*Q(I, 9)+W44R*Q(I,10)+W25R*Q(I,12) * +W45R*Q(I,14)+W16R*Q(I,16)+W26R*Q(I,17)+W36R*Q(I,18) * +W46R*Q(I,19)+W56R*Q(I,20)+W66R*Q(I,21) * +W12I*A(I, 2)+W22I*A(I, 3)+W23I*A(I, 5)+W14I*A(I, 7) * +W24I*A(I, 8)+W34I*A(I, 9)+W44I*A(I,10)+W25I*A(I,12) * +W45I*A(I,14)+W16I*A(I,16)+W26I*A(I,17)+W36I*A(I,18) * +W46I*A(I,19)+W56I*A(I,20)+W66I*A(I,21) 54 CONTINUE G15I=2.*G(2) G25I=2.*G(3) G35I=2.*G(4) C C Derivative with respect to P3: W13R=E13R*P1 W23R=E23R*P2 W33R=E33R*P3 W14R=E12R*P1 W24R=E22R*P2 W34R=E23R*P3*2.+E33R*P2 W44R=E22R*P3+E23R*P2 W15R=E11R*P1 W25R=E12R*P2 W35R=E13R*P3*2.+E33R*P1 W45R=E12R*P3*2.+E13R*P2+E23R*P1 W55R=E11R*P3+E13R*P1 W36R=E13R*P2+E23R*P1 W46R=E12R*P2+E22R*P1 W56R=E11R*P2+E12R*P1 W13I=E13I*P1 W23I=E23I*P2 W33I=E33I*P3 W14I=E12I*P1 W24I=E22I*P2 W34I=E23I*P3*2.+E33I*P2 W44I=E22I*P3+E23I*P2 W15I=E11I*P1 W25I=E12I*P2 W35I=E13I*P3*2.+E33I*P1 W45I=E12I*P3*2.+E13I*P2+E23I*P1 W55I=E11I*P3+E13I*P1 W36I=E13I*P2+E23I*P1 W46I=E12I*P2+E22I*P1 W56I=E11I*P2+E12I*P1 DO 55 I=2,4 G(I)=W13R*A(I, 4)+W23R*A(I, 5)+W33R*A(I, 6)+W14R*A(I, 7) * +W24R*A(I, 8)+W34R*A(I, 9)+W44R*A(I,10)+W15R*A(I,11) * +W25R*A(I,12)+W35R*A(I,13)+W45R*A(I,14)+W55R*A(I,15) * +W36R*A(I,18)+W46R*A(I,19)+W56R*A(I,20) * -W13I*Q(I, 4)-W23I*Q(I, 5)-W33I*Q(I, 6)-W14I*Q(I, 7) * -W24I*Q(I, 8)-W34I*Q(I, 9)-W44I*Q(I,10)-W15I*Q(I,11) * -W25I*Q(I,12)-W35I*Q(I,13)-W45I*Q(I,14)-W55I*Q(I,15) * -W36I*Q(I,18)-W46I*Q(I,19)-W56I*Q(I,20) 55 CONTINUE G16R=2.*G(2) G26R=2.*G(3) G36R=2.*G(4) DO 56 I=2,4 G(I)=W13R*Q(I, 4)+W23R*Q(I, 5)+W33R*Q(I, 6)+W14R*Q(I, 7) * +W24R*Q(I, 8)+W34R*Q(I, 9)+W44R*Q(I,10)+W15R*Q(I,11) * +W25R*Q(I,12)+W35R*Q(I,13)+W45R*Q(I,14)+W55R*Q(I,15) * +W36R*Q(I,18)+W46R*Q(I,19)+W56R*Q(I,20) * +W13I*A(I, 4)+W23I*A(I, 5)+W33I*A(I, 6)+W14I*A(I, 7) * +W24I*A(I, 8)+W34I*A(I, 9)+W44I*A(I,10)+W15I*A(I,11) * +W25I*A(I,12)+W35I*A(I,13)+W45I*A(I,14)+W55I*A(I,15) * +W36I*A(I,18)+W46I*A(I,19)+W56I*A(I,20) 56 CONTINUE G16I=2.*G(2) G26I=2.*G(3) G36I=2.*G(4) C C....................................................................... C C Slowness derivatives of the Christoffel matrix. C W11R=E11R*2. W12R=E12R*2. W22R=E22R*2. W13R=E13R*2. W23R=E23R*2. W33R=E33R*2. W11I=E11I*2. W12I=E12I*2. W22I=E22I*2. W13I=E13I*2. W23I=E23I*2. W33I=E33I*2. C GAB(i+3,l+3)=[A(i,j,k,l)+A(i,k,j,l)]*EA(j)*EA(k): G44R= A1111*W11R+(A1121+A1211)*W12R+(A1131+A1311)*W13R * +A1221*W22R+(A1231+A1321)*W23R+A1331*W33R G45R= A1112*W11R+(A1122+A1212)*W12R+(A1132+A1312)*W13R * +A1222*W22R+(A1232+A1322)*W23R+A1332*W33R G55R= A2112*W11R+(A2122+A2212)*W12R+(A2132+A2312)*W13R * +A2222*W22R+(A2232+A2322)*W23R+A2332*W33R G46R= A1113*W11R+(A1123+A1213)*W12R+(A1133+A1313)*W13R * +A1223*W22R+(A1233+A1323)*W23R+A1333*W33R G56R= A2113*W11R+(A2123+A2213)*W12R+(A2133+A2313)*W13R * +A2223*W22R+(A2233+A2323)*W23R+A2333*W33R G66R= A3113*W11R+(A3123+A3213)*W12R+(A3133+A3313)*W13R * +A3223*W22R+(A3233+A3323)*W23R+A3333*W33R G44R=G44R-Q1111*W11I-(Q1121+Q1211)*W12I-(Q1131+Q1311)*W13I * -Q1221*W22I-(Q1231+Q1321)*W23I-Q1331*W33I G45R=G45R-Q1112*W11I-(Q1122+Q1212)*W12I-(Q1132+Q1312)*W13I * -Q1222*W22I-(Q1232+Q1322)*W23I-Q1332*W33I G55R=G55R-Q2112*W11I-(Q2122+Q2212)*W12I-(Q2132+Q2312)*W13I * -Q2222*W22I-(Q2232+Q2322)*W23I-Q2332*W33I G46R=G46R-Q1113*W11I-(Q1123+Q1213)*W12I-(Q1133+Q1313)*W13I * -Q1223*W22I-(Q1233+Q1323)*W23I-Q1333*W33I G56R=G56R-Q2113*W11I-(Q2123+Q2213)*W12I-(Q2133+Q2313)*W13I * -Q2223*W22I-(Q2233+Q2323)*W23I-Q2333*W33I G66R=G66R-Q3113*W11I-(Q3123+Q3213)*W12I-(Q3133+Q3313)*W13I * -Q3223*W22I-(Q3233+Q3323)*W23I-Q3333*W33I G44I= A1111*W11I+(A1121+A1211)*W12I+(A1131+A1311)*W13I * +A1221*W22I+(A1231+A1321)*W23I+A1331*W33I G45I= A1112*W11I+(A1122+A1212)*W12I+(A1132+A1312)*W13I * +A1222*W22I+(A1232+A1322)*W23I+A1332*W33I G55I= A2112*W11I+(A2122+A2212)*W12I+(A2132+A2312)*W13I * +A2222*W22I+(A2232+A2322)*W23I+A2332*W33I G46I= A1113*W11I+(A1123+A1213)*W12I+(A1133+A1313)*W13I * +A1223*W22I+(A1233+A1323)*W23I+A1333*W33I G56I= A2113*W11I+(A2123+A2213)*W12I+(A2133+A2313)*W13I * +A2223*W22I+(A2233+A2323)*W23I+A2333*W33I G66I= A3113*W11I+(A3123+A3213)*W12I+(A3133+A3313)*W13I * +A3223*W22I+(A3233+A3323)*W23I+A3333*W33I G44I=G44I+Q1111*W11R+(Q1121+Q1211)*W12R+(Q1131+Q1311)*W13R * +Q1221*W22R+(Q1231+Q1321)*W23R+Q1331*W33R G45I=G45I+Q1112*W11R+(Q1122+Q1212)*W12R+(Q1132+Q1312)*W13R * +Q1222*W22R+(Q1232+Q1322)*W23R+Q1332*W33R G55I=G55I+Q2112*W11R+(Q2122+Q2212)*W12R+(Q2132+Q2312)*W13R * +Q2222*W22R+(Q2232+Q2322)*W23R+Q2332*W33R G46I=G46I+Q1113*W11R+(Q1123+Q1213)*W12R+(Q1133+Q1313)*W13R * +Q1223*W22R+(Q1233+Q1323)*W23R+Q1333*W33R G56I=G56I+Q2113*W11R+(Q2123+Q2213)*W12R+(Q2133+Q2313)*W13R * +Q2223*W22R+(Q2233+Q2323)*W23R+Q2333*W33R G66I=G66I+Q3113*W11R+(Q3123+Q3213)*W12R+(Q3133+Q3313)*W13R * +Q3223*W22R+(Q3233+Q3323)*W23R+Q3333*W33R C RETURN END C C======================================================================= C SUBROUTINE SQRTC(CR,CI,SR,SI) REAL CR,CI,SR,SI C C Subroutine designed to calculate the square root of a complex number. C C Input: C CR,CI.. Real and imaginary parts of a complex number. C C Output: C SR,SI.. Real and imaginary parts of the square root of the complex C number. The square root branch with the nonnegative real C part is selected. C C Date: 2022, September 15 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C IF(CR.GE.0.) THEN SR=SQRT(0.5*(SQRT(CR*CR+CI*CI)+CR)) IF(SR.EQ.0.) THEN SI=0. ELSE SI=0.5*CI/SR END IF ELSE SI=SQRT(0.5*(SQRT(CR*CR+CI*CI)-CR)) SR=0.5*CI/SI IF(SR.LT.0.) THEN SR=-SR SI=-SI END IF END IF RETURN END C C======================================================================= C SUBROUTINE HDERR(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 real-valued Hamiltonian function in anisotropic elastic media. 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 (first 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: 2022, October 30 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 HDERR: 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 HDERR: 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 HDERR: 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 HDERR: 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 HDERR: 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 Auxiliary subroutine to HDERR calculating the eigenvalues and C eigenvectors of the real-valued Christoffel matrix. C C ENTRY GABX(A,P1,P2,P3,E1A,E2A,E3A,E1B,E2B,E3B,GAB1,GAB2,GAB3) C Subroutine calculating the first-order spatial derivatives C of the Christoffel matrix, transformed into the given vectors. C C ENTRY GABP(H,E1A,E2A,E3A,E1B,E2B,E3B,GAB4,GAB5,GAB6) C Subroutine calculating the first-order partial derivatives C of the Christoffel matrix with respect to the slowness vector, C transformed into the given vectors. C C ENTRY GAAXX(A,P1,P2,P3,E1A,E2A,E3A, C * GAB11,GAB12,GAB22,GAB13,GAB23,GAB33) C Subroutine calculating the second-order spatial derivatives C of the diagonal elements of the Christoffel matrix, transformed C into the given vectors. 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-order mixed derivatives C of the diagonal elements of the Christoffel matrix transformed C into the given vectors. C C ENTRY GAAPP(E1A,E2A,E3A,GAB44,GAB45,GAB55,GAB46,GAB56,GAB66) C Subroutine calculating the second-order slowness derivatives C of the diagonal elements of the Christoffel matrix transformed C into the given vectors. 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 transformed into C the given vectors. 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 transformed into C the given vectors. 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 transformed into C the given vectors. 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