C
C Subroutine file 'wanpfa.for' with subroutines used in 'crtpfv.for' and C 'wan.for' C C Version: 7.00 C Date: 2013, June 12 C C Coded by: Petr Bulant C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C http://sw3d.cz/staff/bulant.htm C C======================================================================= C SUBROUTINE WACHRI(P1,P2,P3,B11,B12,B22,B13,B23,B33, * B14,B24,B34,B44,B15,B25,B35,B45,B55, * B16,B26,B36,B46,B56,B66, * G1,G2,G3,EE,DER) C C---------------------------------------------------------------------- C Subroutine to compute the Christoffel matrix, its eigenvalues C and eigenvectors. REAL P1,P2,P3 REAL B11,B12,B22,B13,B23,B33,B14,B24,B34,B44 REAL B15,B25,B35,B45,B55,B16,B26,B36,B46,B56,B66 REAL G1,G2,G3,EE(9),DER(9) C C Input: C P1,P2,P3... Slowness vector. C Bii ... Values of real parts of 21 reduced C (divided by the density) elastic parameters. C Output: C G1,G2,G3 ... Eigenvalues of the Christoffel matrix. C EE ... Eigenvectors of the Christoffel matrix. C DER ... Derivatives dx/dt=dH/dp=Aijkl Ej Ek pl stored columnwise C for qP, qS1 and qS2 waves. C C----------------------------------------------------------------------- C REAL A11,A12,A13,A14,A21,A22,A23,A24,A31,A32,A33,A34,A44 REAL A15,A25,A35,A45,A55,A16,A26,A36,A46,A56,A66 REAL A1111,A2111,A3111,A1211,A2211,A3211,A1311,A2311,A3311 REAL A1121,A2121,A3121,A1221,A2221,A3221,A1321,A2321,A3321 REAL A1131,A2131,A3131,A1231,A2231,A3231,A1331,A2331,A3331 REAL A1112,A2112,A3112,A1212,A2212,A3212,A1312,A2312,A3312 REAL A1122,A2122,A3122,A1222,A2222,A3222,A1322,A2322,A3322 REAL A1132,A2132,A3132,A1232,A2232,A3232,A1332,A2332,A3332 REAL A1113,A2113,A3113,A1213,A2213,A3213,A1313,A2313,A3313 REAL A1123,A2123,A3123,A1223,A2223,A3223,A1323,A2323,A3323 REAL A1133,A2133,A3133,A1233,A2233,A3233,A1333,A2333,A3333 EQUIVALENCE (A11,A1111) EQUIVALENCE (A22,A2222) EQUIVALENCE (A33,A3333) EQUIVALENCE (A16,A1112,A1121,A1211,A2111) EQUIVALENCE (A26,A2221,A2212,A2122,A1222) EQUIVALENCE (A15,A1113,A1131,A1311,A3111) EQUIVALENCE (A35,A3331,A3313,A3133,A1333) EQUIVALENCE (A24,A2223,A2232,A2322,A3222) EQUIVALENCE (A34,A3332,A3323,A3233,A2333) EQUIVALENCE (A23,A2233,A3322) EQUIVALENCE (A13,A1133,A3311) EQUIVALENCE (A12,A1122,A2211) EQUIVALENCE (A44, A2323,A3232,A2332,A3223) EQUIVALENCE (A55, A1313,A3131,A1331,A3113) EQUIVALENCE (A66, A1212,A2121,A1221,A2112) EQUIVALENCE (A14,A1123,A1132,A2311,A3211) EQUIVALENCE (A25,A2213,A2231,A1322,A3122) EQUIVALENCE (A36,A3312,A3321,A1233,A2133) EQUIVALENCE (A56,A1321,A3112,A2113,A1231,A1213,A2131,A1312,A3121) EQUIVALENCE (A46,A2312,A3221,A1223,A2132,A2123,A1232,A2321,A3212) EQUIVALENCE (A45,A3213,A2331,A1332,A3123,A3132,A1323,A3231,A2313) REAL GAMMA(6),E11,E21,E31,E12,E22,E32,E13,E23,E33 C EQUIVALENCE (GAMMA(1),G1),(GAMMA(3),G2),(GAMMA(6),G3) C EQUIVALENCE (EE(1),E11),(EE(4),E12),(EE(7),E13) C EQUIVALENCE (EE(2),E21),(EE(5),E22),(EE(8),E23) C EQUIVALENCE (EE(3),E31),(EE(6),E32),(EE(9),E33) REAL A111,A112,A121,A122,A113,A123,A131,A132,A133 REAL A211,A212,A221,A222,A213,A223,A231,A232,A233 REAL A311,A312,A322,A313,A321,A323,A331,A332,A333 REAL AUX INTEGER KQICHM DATA KQICHM/-1/ C C GAMMA,G1,G2,G3...Christoffel matrix, later its eigenvalues. C (E11,E12,E13) C EE=(E21,E22,E23)... Eigenvectors of the christoffel matrix. C (E31,E32,E33) C A111,A211,A311,A112,A212,A312,A122,A222,A322,A113,A213,A313,A123, C A223,A323,A133,A233,A333... A(I,J,K,L)*P(L) summed over L. C A11,A21,A31,A12,A22,A32,A13,A23,A33 ... Aijk*Ek C C....................................................................... C A11=B11 A22=B22 A33=B33 A16=B16 A26=B26 A15=B15 A35=B35 A24=B24 A34=B34 A23=B23 A13=B13 A12=B12 A44=B44 A55=B55 A66=B66 A14=B14 A25=B25 A36=B36 A56=B56 A46=B46 A45=B45 C Christoffel matrix: A111=A1111*P1+A1112*P2+A1113*P3 A112=A1121*P1+A1122*P2+A1123*P3 A121=A1211*P1+A1212*P2+A1213*P3 A122=A1221*P1+A1222*P2+A1223*P3 A113=A1131*P1+A1132*P2+A1133*P3 A123=A1231*P1+A1232*P2+A1233*P3 A131=A1311*P1+A1312*P2+A1313*P3 A132=A1321*P1+A1322*P2+A1323*P3 A133=A1331*P1+A1332*P2+A1333*P3 A211=A2111*P1+A2112*P2+A2113*P3 A212=A2121*P1+A2122*P2+A2123*P3 A221=A2211*P1+A2212*P2+A2213*P3 A222=A2221*P1+A2222*P2+A2223*P3 A213=A2131*P1+A2132*P2+A2133*P3 A223=A2231*P1+A2232*P2+A2233*P3 A231=A2311*P1+A2312*P2+A2313*P3 A232=A2321*P1+A2322*P2+A2323*P3 A233=A2331*P1+A2332*P2+A2333*P3 A311=A3111*P1+A3112*P2+A3113*P3 A312=A3121*P1+A3122*P2+A3123*P3 A322=A3221*P1+A3222*P2+A3223*P3 A313=A3131*P1+A3132*P2+A3133*P3 A321=A3211*P1+A3212*P2+A3213*P3 A323=A3231*P1+A3232*P2+A3233*P3 A331=A3311*P1+A3312*P2+A3313*P3 A332=A3321*P1+A3322*P2+A3323*P3 A333=A3331*P1+A3332*P2+A3333*P3 GAMMA(1)=P1*A111+P2*A211+P3*A311 GAMMA(2)=P1*A112+P2*A212+P3*A312 GAMMA(3)=P1*A122+P2*A222+P3*A322 GAMMA(4)=P1*A113+P2*A213+P3*A313 GAMMA(5)=P1*A123+P2*A223+P3*A323 GAMMA(6)=P1*A133+P2*A233+P3*A333 C C Quasi-isotropic modification of the Christoffel matrix: IF(KQICHM.LT.0) THEN CALL RSEP3I('QICHM',KQICHM,0) END IF IF(KQICHM.GT.0) THEN AUX=SQRT(P1*P1+P2*P2+P3*P3) E11=P1/AUX E21=P2/AUX E31=P3/AUX E12=GAMMA(1)*E11+GAMMA(2)*E21+GAMMA(4)*E31 E22=GAMMA(2)*E11+GAMMA(3)*E21+GAMMA(5)*E31 E32=GAMMA(4)*E11+GAMMA(5)*E21+GAMMA(6)*E31 AUX=2.*(E11*E12+E21*E22+E31*E32) GAMMA(1)=GAMMA(1)-E11*E12-E12*E11+AUX*E11*E11 GAMMA(2)=GAMMA(2)-E11*E22-E12*E21+AUX*E11*E21 GAMMA(3)=GAMMA(3)-E21*E22-E22*E21+AUX*E21*E21 GAMMA(4)=GAMMA(4)-E11*E32-E12*E31+AUX*E11*E31 GAMMA(5)=GAMMA(5)-E21*E32-E22*E31+AUX*E21*E31 GAMMA(6)=GAMMA(6)-E31*E32-E32*E31+AUX*E31*E31 END IF C C Selecting eigenvalue and eigenvector of the Christoffel matrix: CALL EIGEN(GAMMA,EE,3,0) G1=GAMMA(1) G2=GAMMA(3) G3=GAMMA(6) E11=EE(1) E21=EE(2) E31=EE(3) E12=EE(4) E22=EE(5) E32=EE(6) E13=EE(7) E23=EE(8) E33=EE(9) IF (G3.LT.0.) THEN C WAN-11 CALL ERROR('WAN-11: Negative eigenvalue of Christoffel matrix') C This error should not appear. END IF AUX=E11*E11+E21*E21+E31*E31 IF (ABS(AUX-1.).GT.0.00001) THEN C WAN-12 CALL ERROR('WAN-12: Eigenvector is not normalized') C This error should not appear. ENDIF AUX=E12*E12+E22*E22+E32*E32 IF (ABS(AUX-1.).GT.0.00001) THEN C WAN-13 CALL ERROR('WAN-13: Eigenvector is not normalized') C This error should not appear. ENDIF AUX=E13*E13+E23*E23+E33*E33 IF (ABS(AUX-1.).GT.0.00001) THEN C WAN-14 CALL ERROR('WAN-14: Eigenvector is not normalized') C This error should not appear. ENDIF C C Computation of derivatives dx/dt: A11= A111*E11+A112*E21+A113*E31 A21= A211*E11+A212*E21+A213*E31 A31= A311*E11+A312*E21+A313*E31 A12= A121*E11+A122*E21+A123*E31 A22= A221*E11+A222*E21+A223*E31 A32= A321*E11+A322*E21+A323*E31 A13= A131*E11+A132*E21+A133*E31 A23= A231*E11+A232*E21+A233*E31 A33= A331*E11+A332*E21+A333*E31 DER(1)=A11*E11+ A12*E21+ A13*E31 DER(2)=A21*E11+ A22*E21+ A23*E31 DER(3)=A31*E11+ A32*E21+ A33*E31 A11= A111*E12+A112*E22+A113*E32 A21= A211*E12+A212*E22+A213*E32 A31= A311*E12+A312*E22+A313*E32 A12= A121*E12+A122*E22+A123*E32 A22= A221*E12+A222*E22+A223*E32 A32= A321*E12+A322*E22+A323*E32 A13= A131*E12+A132*E22+A133*E32 A23= A231*E12+A232*E22+A233*E32 A33= A331*E12+A332*E22+A333*E32 DER(4)=A11*E12+ A12*E22+ A13*E32 DER(5)=A21*E12+ A22*E22+ A23*E32 DER(6)=A31*E12+ A32*E22+ A33*E32 A11= A111*E13+A112*E23+A113*E33 A21= A211*E13+A212*E23+A213*E33 A31= A311*E13+A312*E23+A313*E33 A12= A121*E13+A122*E23+A123*E33 A22= A221*E13+A222*E23+A223*E33 A32= A321*E13+A322*E23+A323*E33 A13= A131*E13+A132*E23+A133*E33 A23= A231*E13+A232*E23+A233*E33 A33= A331*E13+A332*E23+A333*E33 DER(7)=A11*E13+ A12*E23+ A13*E33 DER(8)=A21*E13+ A22*E23+ A23*E33 DER(9)=A31*E13+ A32*E23+ A33*E33 C RETURN END C C======================================================================= C SUBROUTINE WAMAT(A,B,C,D) C C---------------------------------------------------------------------- C Subroutine to compute the product of two 2x2 complex matrices. C The second matrix (C+iD) is destroyed in the computation. REAL A(2,2),B(2,2),C(2,2),D(2,2) C Input: C A,B,C,D ... Real and imaginary parts of the two matrices. C Output: C C,D ... Real and imaginary parts of resulting matrix. C C....................................................................... C Auxiliary storage locations: REAL CR11,CR21,CR12,CR22,CI11,CI21,CI12,CI22 C....................................................................... CR11=A(1,1)*C(1,1)-B(1,1)*D(1,1)+A(1,2)*C(2,1)-B(1,2)*D(2,1) CR21=A(2,1)*C(1,1)-B(2,1)*D(1,1)+A(2,2)*C(2,1)-B(2,2)*D(2,1) CR12=A(1,1)*C(1,2)-B(1,1)*D(1,2)+A(1,2)*C(2,2)-B(1,2)*D(2,2) CR22=A(2,1)*C(1,2)-B(2,1)*D(1,2)+A(2,2)*C(2,2)-B(2,2)*D(2,2) C CI11=A(1,1)*D(1,1)+B(1,1)*C(1,1)+A(1,2)*D(2,1)+B(1,2)*C(2,1) CI21=A(2,1)*D(1,1)+B(2,1)*C(1,1)+A(2,2)*D(2,1)+B(2,2)*C(2,1) CI12=A(1,1)*D(1,2)+B(1,1)*C(1,2)+A(1,2)*D(2,2)+B(1,2)*C(2,2) CI22=A(2,1)*D(1,2)+B(2,1)*C(1,2)+A(2,2)*D(2,2)+B(2,2)*C(2,2) C C(1,1)=CR11 C(2,1)=CR21 C(1,2)=CR12 C(2,2)=CR22 D(1,1)=CI11 D(2,1)=CI21 D(1,2)=CI12 D(2,2)=CI22 C RETURN END C C======================================================================= C SUBROUTINE WAVECR(G,H,U,V,RU,RV) C C---------------------------------------------------------------------- C Subroutine to reorganize the vectors U and V in such way, C that the pair U,V is equal to pair G,H rotated with angle smaller C than 45 degree (clockwise or anticlockwise). C The real numbers RU (associated with U) and RV (associated with V), C are reorganized in the same way. C REAL G(3),H(3),U(3),V(3),RU,RV C Input: C G,H ... A pair of orthonormal vectors. C U,V ... A pair of orthonormal vectors. C RU,RV ... Real numbers associated to U and V. C Output: C U,V ... Selection from input vectors U,V,-U,-V in such way, that C the output pair U,V is equal to pair G,H rotated C with angle smaller than 45 degree. C RU,RV ... Real numbers associated to U and V. C C....................................................................... C C Auxiliary storage locations: REAL E1,E2,E3,AUX C C----------------------------------------------------------------------- AUX=G(1)*U(1)+G(2)*U(2)+G(3)*U(3) IF(ABS(AUX).LT.0.707107) THEN E1=U(1) E2=U(2) E3=U(3) U(1)=V(1) U(2)=V(2) U(3)=V(3) V(1)=E1 V(2)=E2 V(3)=E3 AUX =RU RU =RV RV =AUX END IF AUX=G(1)*U(1)+G(2)*U(2)+G(3)*U(3) IF(AUX.LT.0.0) THEN U(1)=-U(1) U(2)=-U(2) U(3)=-U(3) END IF AUX=H(1)*V(1)+H(2)*V(2)+H(3)*V(3) IF(AUX.LT.0.0) THEN V(1)=-V(1) V(2)=-V(2) V(3)=-V(3) END IF RETURN END C C C======================================================================= C REAL FUNCTION WAANGL(G,H,U,V) C C---------------------------------------------------------------------- C Subroutine to compute the rotation angle between the pair of C orthonormal vectors G and H and the pair of orthonormal vectors U and C V according to equation 18 of Bulant & Klimes, 2002, Pageoph. C REAL G(3),H(3),U(3),V(3) C Input: C G,H ... A pair of orthonormal vectors. C U,V ... A pair of orthonormal vectors. C Output: C WAANGL ... The rotation angle. C C....................................................................... C Subroutines and external functions required: EXTERNAL WASUM REAL WASUM C WASUM ... This file. C----------------------------------------------------------------------- WAANGL=ATAN2((WASUM(U,H)-WASUM(V,G)),(WASUM(U,G)+WASUM(V,H))) RETURN END C C C======================================================================= C SUBROUTINE WAPROJ(P1,P2,P3,E1,E2,E3,G1,G2,G3,H1,H2,H3, * GE1,GE2,HE1,HE2) C C---------------------------------------------------------------------- C Subroutine to project vectors G and H to the plane defined by vector E C and vector PxE. C REAL P1,P2,P3,E1,E2,E3,G1,G2,G3,H1,H2,H3,GE1,GE2,HE1,HE2 C Input: C P1,P2,P3,E1,E2,E3 ... Vectors defining the plane. C G1,G2,G3,H1,H2,H3 ... Vectors to be projected. C Output: C GE1,GE2,HE1,HE2 ... Projected vectors. C C....................................................................... C C Auxiliary storage locations: REAL F1,F2,F3,AUX C C----------------------------------------------------------------------- C Second vector defining the plane: F1=P2*E3-E2*P3 F2=P3*E1-E3*P1 F3=P1*E2-E1*P2 AUX=SQRT(F1*F1+F2*F2+F3*F3) F1=F1/AUX F2=F2/AUX F3=F3/AUX C Projecting vectors G and H to the plane defined by E and F: GE1=E1*G1+E2*G2+E3*G3 GE2=F1*G1+F2*G2+F3*G3 HE1=E1*H1+E2*H2+E3*H3 HE2=F1*H1+F2*H2+F3*H3 C Normalizing vectors GE and HE: AUX=SQRT(GE1*GE1+GE2*GE2) GE1=GE1/AUX GE2=GE2/AUX AUX=SQRT(HE1*HE1+HE2*HE2) HE1=HE1/AUX HE2=HE2/AUX C RETURN END C C======================================================================= C SUBROUTINE WAIJKL(A,B,I) C C---------------------------------------------------------------------- REAL A(10,21),B(3,3,3,3) INTEGER I B(1,1,1,1)=A(I,1) B(1,1,2,2)=A(I,2) B(2,2,1,1)=A(I,2) B(2,2,2,2)=A(I,3) B(1,1,3,3)=A(I,4) B(3,3,1,1)=A(I,4) B(2,2,3,3)=A(I,5) B(3,3,2,2)=A(I,5) B(3,3,3,3)=A(I,6) B(1,1,2,3)=A(I,7) B(1,1,3,2)=A(I,7) B(2,3,1,1)=A(I,7) B(3,2,1,1)=A(I,7) B(2,2,2,3)=A(I,8) B(2,2,3,2)=A(I,8) B(2,3,2,2)=A(I,8) B(3,2,2,2)=A(I,8) B(3,3,3,2)=A(I,9) B(3,3,2,3)=A(I,9) B(3,2,3,3)=A(I,9) B(2,3,3,3)=A(I,9) B(2,3,2,3)=A(I,10) B(3,2,3,2)=A(I,10) B(2,3,3,2)=A(I,10) B(3,2,2,3)=A(I,10) B(1,1,1,3)=A(I,11) B(1,1,3,1)=A(I,11) B(1,3,1,1)=A(I,11) B(3,1,1,1)=A(I,11) B(2,2,1,3)=A(I,12) B(2,2,3,1)=A(I,12) B(1,3,2,2)=A(I,12) B(3,1,2,2)=A(I,12) B(3,3,3,1)=A(I,13) B(3,3,1,3)=A(I,13) B(3,1,3,3)=A(I,13) B(1,3,3,3)=A(I,13) B(3,2,1,3)=A(I,14) B(2,3,3,1)=A(I,14) B(1,3,3,2)=A(I,14) B(3,1,2,3)=A(I,14) B(3,1,3,2)=A(I,14) B(1,3,2,3)=A(I,14) B(3,2,3,1)=A(I,14) B(2,3,1,3)=A(I,14) B(1,3,1,3)=A(I,15) B(3,1,3,1)=A(I,15) B(1,3,3,1)=A(I,15) B(3,1,1,3)=A(I,15) B(1,1,1,2)=A(I,16) B(1,1,2,1)=A(I,16) B(1,2,1,1)=A(I,16) B(2,1,1,1)=A(I,16) B(2,2,2,1)=A(I,17) B(2,2,1,2)=A(I,17) B(2,1,2,2)=A(I,17) B(1,2,2,2)=A(I,17) B(3,3,1,2)=A(I,18) B(3,3,2,1)=A(I,18) B(1,2,3,3)=A(I,18) B(2,1,3,3)=A(I,18) B(2,3,1,2)=A(I,19) B(3,2,2,1)=A(I,19) B(1,2,2,3)=A(I,19) B(2,1,3,2)=A(I,19) B(2,1,2,3)=A(I,19) B(1,2,3,2)=A(I,19) B(2,3,2,1)=A(I,19) B(3,2,1,2)=A(I,19) B(1,3,2,1)=A(I,20) B(3,1,1,2)=A(I,20) B(2,1,1,3)=A(I,20) B(1,2,3,1)=A(I,20) B(1,2,1,3)=A(I,20) B(2,1,3,1)=A(I,20) B(1,3,1,2)=A(I,20) B(3,1,2,1)=A(I,20) B(1,2,1,2)=A(I,21) B(2,1,2,1)=A(I,21) B(1,2,2,1)=A(I,21) B(2,1,1,2)=A(I,21) RETURN END C C======================================================================= C REAL FUNCTION WASUM(A,B) C C---------------------------------------------------------------------- REAL A(3),B(3) WASUM=A(1)*B(1)+A(2)*B(2)+A(3)*B(3) RETURN END C C======================================================================= C REAL FUNCTION WASUM3(AA,B,C) C C---------------------------------------------------------------------- REAL AA(3,3),B(3),C(3) INTEGER I,J WASUM3=0. DO 20, I=1,3 DO 10, J=1,3 WASUM3=WASUM3+AA(I,J)*B(I)*C(J) 10 CONTINUE 20 CONTINUE RETURN END C C======================================================================= C REAL FUNCTION WASUM5(A,B,C,D,E) C C---------------------------------------------------------------------- REAL A(3,3,3,3),B(3),C(3),D(3),E(3) INTEGER I,J,K,L WASUM5=0. DO 13, I=1,3 DO 12, J=1,3 DO 11, K=1,3 DO 10, L=1,3 WASUM5=WASUM5+A(I,J,K,L)*B(I)*C(J)*D(K)*E(L) 10 CONTINUE 11 CONTINUE 12 CONTINUE 13 CONTINUE RETURN END C C======================================================================= C REAL FUNCTION WASUM4(A,B,C,D,J) C C---------------------------------------------------------------------- REAL A(3,3,3,3),B(3),C(3),D(3) INTEGER I,J,K,L WASUM4=0. DO 13, I=1,3 DO 11, K=1,3 DO 10, L=1,3 WASUM4=WASUM4+A(I,J,K,L)*B(I)*C(K)*D(L) 10 CONTINUE 11 CONTINUE 13 CONTINUE RETURN END C C======================================================================= C