SUBROUTINE COEF(STU,UC,ITRANS) C C ROUTINE FOR THE COMPUTATION OF REFLECTION/TRANSMISSION COEFFICIENTS C FOR GENERAL ISOTROPIC/ANISOTROPIC MEDIA C C INPUT PARAMETERS: C STU(1-3)... CARTESIAN COMPONENTS OF THE DISPLACEMENT VECTOR C STU(4-6)... CARTESIAN COMPONETS OF THE TRACTION VECTOR CALCULATED C IN THIS ROUTINE C ITRANS... ITRANS=O: REFLECTION C ITRANS=1: TRANSMISSION C ITRANS=999: SURFACE CONVERSION C C OUTPUT PARAMETERS: C UC(1-3)... RAY-CENTERED COMPONENTS OF THE DISPLACEMENT VECTOR C COMPLEX A(6,6),C(6,6),COA(3,3),COC(3,3),AA(6),TAU(3),P(3),U(3), 1CP1,CP2,CP3,DETA,DETC,UC(3),STU(6),cab DIMENSION Y(18),UN(3) COMMON /AUXI/IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOUT, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori INTEGER CODE COMMON /COD/CODE(50,2),KREF,KC,ITYPE COMPLEX PS COMMON /RAY/ AY(28,2000),DS(20,20),KINT(20),HHH(3,3),TMAX, 1 PS(3,7,20),IS(8,20),N,IREF,IND,IND1 COMMON /ZERO/ RNULL C INDI=0 IF(ITRANS.EQ.999) THEN ITRANS=0 INDI=1 END IF I1=1 I2=3 LAY1=IS(1,IREF) LAY=IS(7,IREF) DO 910 K=1,6 IF(K.GT.3) GOTO 900 UN(K)=DS(K,IREF) 900 Y(K)=AY(K+1,N) 910 CONTINUE IF(ICOEF.GT.0) 1WRITE(LOUT,'(A,/,3(3F14.7,/),3I5)')'COEF: Y(1)-Y(6),UN,IREF,LAY, 2ITYPE',(Y(K),K=1,6),UN,IREF,LAY,ITYPE CALL PARDIS (Y,0) C C DISPLACEMENT AND STRESS VECTOR OF INCIDENT WAVE C BSTU=SQRT(REAL(STU(1)*CONJG(STU(1))+STU(2)*CONJG(STU(2)) 1 +STU(3)*CONJG(STU(3)))) BST=STU(1)*HHH(ITYPE,1)+STU(2)*HHH(ITYPE,2) 1 +STU(3)*HHH(ITYPE,3) BSTU=SIGN(1.,BST)*BSTU DO 950 K=1,3 P(K)=PS(K,7,IREF) STU(K)=STU(K)/BSTU U(k)=STU(K) 950 CONTINUE CALL STRESS (P,UN,U,TAU) DO 960 K=1,3 STU(K+3)=TAU(K) 960 CONTINUE IF(ICOEF.GT.0) 1WRITE(LOUT,'(A,/,6(2F14.7,/))') 'COEF: INCIDENT U,TAU',U,TAU 1000 IF(IANI(LAY).EQ.0) GOTO 3000 C C COMPUTATION OF STRESS-DISPLACEMENT-VECTORS FOR GENERATED WAVES C C ANISOTROPIC CASE C 1010 IF(I1.GT.3)CALL PARDIS(Y,0) DO 1020 K=I1,I2 DO 1030 L=1,3 P(L)=PS(L,K,IREF) 1030 CONTINUE C C EVALUATION OF DISPLACEMENT VECTOR (GENERALLY COMPLEX-VALUED) C CALL DISPL(P,U,K) CALL STRESS(P,UN,U,TAU) IF(ICOEF.GT.0)WRITE(LOUT,'(A,/,9(2F14.7,/))')'COEF: P,U,TAU' 1,P,U,TAU DO 1060 L=1,3 A(K,L)=U(L) A(K,L+3)=TAU(L) 1060 CONTINUE 1020 CONTINUE IF(INDI.EQ.1) GOTO 5000 IF(I2.EQ.6) GOTO 4000 I1=4 I2=6 LAY=IS(8,IREF) IF(LAY.EQ.0) GOTO 5000 IF(IANI(LAY).NE.0) GOTO 1010 C C ISOTROPIC CASE C 3000 IF(I1.GT.3) CALL PARDIS(Y,0) DO 3010 L=1,3 P(L)=PS(L,I1,IREF) 3010 CONTINUE CP1=P(1) CP2=P(2) IVERT=0 IF(CABS(CP1).LT..0000001.AND.CABS(CP2).LT..0000001)IVERT=1 IF(IVERT.EQ.1)THEN U(1)=1. U(2)=0. U(3)=0. GO TO 3020 END IF CP3=-(CP1*CP1+CP2*CP2) CP1=CP1*P(3) CP2=CP2*P(3) CAB=SQRT(REAL(CP1*CONJG(CP1)+CP2*CONJG(CP2)+CP3*CONJG(CP3))) U(1)=CP1/CAB U(2)=CP2/CAB U(3)=CP3/CAB 3020 CONTINUE CALL STRESS(P,UN,U,TAU) IF(ICOEF.GT.0)WRITE(LOUT,'(A,/,9(2F14.7/))')'COEF: P,U,TAU' 1,P,U,TAU DO 3030 L=1,3 A(I1,L)=U(L) A(I1,L+3)=TAU(L) 3030 CONTINUE IF(IVERT.EQ.1)THEN U(1)=0. U(2)=1. IF(REAL(P(3)).LT.0.)U(2)=-1. U(3)=0. GO TO 3040 END IF CP1=-P(2) CP2=P(1) CP3=CMPLX(0.,0.) CAB=SQRT(REAL(CP1*CONJG(CP1)+CP2*CONJG(CP2)+CP3*CONJG(CP3))) U(1)=CP1/CAB U(2)=CP2/CAB U(3)=cp3/CAB 3040 CONTINUE CALL STRESS (P,UN,U,TAU) IF(ICOEF.GT.0)WRITE(LOUT,'(A,/,9(2F14.7/))')'COEF: P,U,TAU' 1,P,U,TAU DO 3050 L=1,3 A(I1+1,L)=U(L) A(I1+1,L+3)=TAU(L) 3050 CONTINUE DO 3210 L=1,3 P(L)=PS(L,I2,IREF) 3210 CONTINUE CP1=P(1) CP2=P(2) CP3=P(3) CAB=SQRT(REAL(CP1*CONJG(CP1)+CP2*CONJG(CP2)+CP3*CONJG(CP3))) U(1)=CP1/CAB U(2)=CP2/CAB U(3)=CP3/CAB CALL STRESS (P,UN,U,TAU) IF(ICOEF.GT.0)WRITE(LOUT,'(A,/,9(2F14.7/))')'COEF: P,U,TAU' 1,P,U,TAU DO 3220 L=1,3 A(I2,L)=U(L) A(I2,L+3)=TAU(L) 3220 CONTINUE IF(INDI.EQ.1) GOTO 5000 IF(I2.EQ.6)GOTO 4000 I1=4 I2=6 LAY=IS(8,IREF) IF(LAY.EQ.0) GOTO 5000 GOTO 1000 C C REVERSE OF SIGNS FOR REFLECTIONS C 4000 DO 4010 L=1,3 DO 4010 K=1,6 A(L,K)=-A(L,K) 4010 CONTINUE IF(ICOEF.EQ.0) GOTO 4012 WRITE(LOUT,'(//,A)') 'COEF: MATRIX A' DO 4011 K=1,6 WRITE(LOUT,'(6(12F10.5))') (A(L,K),L=1,6) 4011 CONTINUE C C SOLUTION OF SYSTEM OF LINEAR EQUATIONS ( CRAMERS METHOD ) C 4012 CALL DET6(A,DETA) DO 4020 K=1,3 UC(K)=CMPLX(0.,0.) 4020 CONTINUE IF(CABS(DETA).LT.1.E-20) THEN WRITE(LOUT,'(A)') 1 'COEF: MATRIX FOR R/T COEF. SINGULAR: AMPLITUDE NOT EVALUATED' IND=11 RETURN END IF IF(IANI(LAY1).NE.0.OR.ITYPE.EQ.3)THEN L1=ITYPE IF(ITRANS.NE.0)L1=L1+3 L2=L1 END IF IF(IANI(LAY1).EQ.0.AND.ITYPE.NE.3)THEN L1=1 IF(ITRANS.NE.0)L1=L1+3 L2=L1+1 END IF IF(ICOEF.NE.0) WRITE(LOUT,'(A,2I5)')'COEF: L1,L2',L1,L2 DO 4044 L=L1,L2 DO 4043 J=1,6 DO 4042 K=1,6 C(J,K)=A(J,K) IF(L.EQ.J)C(J,K)=STU(K) 4042 CONTINUE 4043 CONTINUE CALL DET6(C,DETC) AA(L)=DETC/DETA IF(ICOEF.GT.0)WRITE(LOUT,'(A,2F12.6)') 'COEF: AA=',AA(L) 4044 CONTINUE DO 4100 L=L1,L2 LL=L IF(ITRANS.NE.0)LL=L-3 UC(LL)=AA(L)*BSTU 4100 CONTINUE RETURN C C INTERACTION WITH THE FREE SURFACE C 5000 DO 5010 L=1,3 DO 5010 K=1,3 COA(L,K)=-A(L,K+3) 5010 CONTINUE IF(ICOEF.EQ.0)GO TO 5012 WRITE(LOUT,'(//,A)')'COEF: MATRIX COA' WRITE(LOUT,'(3(6F10.5/))')((COA(L,K),L=1,3),K=1,3) C C SOLUTION OF SYSTEM OF 3 LINEAR EQUATIONS (CRAMER'S METHOD) C 5012 CALL DET3(COA,DETA) IF(CABS(DETA).LT.1.E-09)THEN WRITE(LOUT,'(A)')'COEF: MATRIX COA IS SINGULAR, 1 PROGRAM TERMINATES' STOP END IF DO 5044 L=1,3 DO 5043 J=1,3 DO 5042 K=1,3 COC(J,K)=COA(J,K) IF(L.EQ.J)COC(J,K)=STU(K+3) 5042 CONTINUE 5043 CONTINUE CALL DET3(COC,DETC) AA(L)=DETC/DETA IF(ICOEF.GT.0) 1WRITE(LOUT,'(A,2F12.6)')'COEF: AA=',AA(L) 5044 CONTINUE IF(INDI.NE.1)GO TO 5200 C C CONVERSION AT THE FREE SURFACE C DO 5100 K=1,3 U(K)=CMPLX(0.,0.) DO 5050 J=1,3 U(K)=AA(J)*A(J,K)+U(K) 5050 CONTINUE UC(K)=U(K)+STU(K) 5100 CONTINUE IF(ICOEF.GT.0)WRITE(LOUT,'(A,/,3(2F12.5/))') 'COEF: 1 NORMALIZED CONVERSION-VECTOR',UC DO 5150 K=1,3 UC(K)=UC(K)*BSTU 5150 CONTINUE IF(ICOEF.GT.0)WRITE(LOUT,'(A,/,3(2F12.5/))') 'COEF: 1 CONVERSION-VECTOR',UC RETURN C C REFLECTION FROM THE FREE SURFACE C 5200 L1=1 IF(ITYPE.EQ.3)L1=3 L2=2 IF(ITYPE.EQ.3)L2=3 do 5250 k=1,3 UC(K)=CMPLX(0.,0.) 5250 CONTINUE DO 5300 L=L1,L2 UC(L)=AA(L)*BSTU 5300 CONTINUE RETURN END C C ********************************************************* C SUBROUTINE DET6(A,D6) COMPLEX A(6,6),B(5,5),D6,D D6=(0.,0.) DO 10 J=1,6 IF(A(1,J).NE.(0.,0.)) THEN CALL SET(A,B,1,J,6) CALL DET5(B,D) C=(-1)**(J+1) D6=D6+CMPLX(C,0.)*D*A(1,J) ENDIF 10 CONTINUE RETURN END C C ********************************************************* C SUBROUTINE DET5(A,D5) COMPLEX A(5,5),B(4,4),D5,D D5=(0.,0.) DO 10 J=1,5 IF(A(1,J).NE.(0.,0.)) THEN CALL SET(A,B,1,J,5) CALL DET4(B,D) C=(-1)**(J+1) D5=D5+CMPLX(C,0.)*D*A(1,J) ENDIF 10 CONTINUE RETURN END C C *********************************************************** C SUBROUTINE DET4(A,D4) COMPLEX A(4,4),B(3,3),D4,D D4=(0.,0.) DO 10 J=1,4 IF(A(1,J).NE.(0.,0.)) THEN CALL SET(A,B,1,J,4) CALL DET3(B,D) C=(-1)**(J+1) D4=D4+CMPLX(C,0.)*D*A(1,J) ENDIF 10 CONTINUE RETURN END C C ********************************************************** C SUBROUTINE DET3(A,D3) COMPLEX A(3,3),D3 D3=A(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(3,2)) 1+ A(1,2)*(A(2,3)*A(3,1)-A(2,1)*A(3,3)) 2+ A(1,3)*(A(2,1)*A(3,2)-A(2,2)*A(3,1)) RETURN END C C ************************************************************ C SUBROUTINE SET(A,B,I,J,K) COMPLEX A(K,K),B(K-1,K-1) IFL=0 JFL=0 DO 20 I1=1,K IF(I1.EQ.I) THEN IFL=1 GOTO 20 ENDIF DO 10 J1=1,K IF(J1.EQ.J) THEN JFL=1 GOTO 10 ENDIF IF(IFL.EQ.0.AND.JFL.EQ.0) THEN B(I1,J1)=A(I1,J1) ELSE IF(IFL.EQ.0.AND.JFL.EQ.1) THEN B(I1,J1-1)=A(I1,J1) ELSE IF(IFL.EQ.1.AND.JFL.EQ.0) THEN B(I1-1,J1)=A(I1,J1) ELSE B(I1-1,J1-1)=A(I1,J1) ENDIF 10 CONTINUE JFL=0 20 CONTINUE RETURN END C C ********************************************************** C SUBROUTINE DISC(Y,DEP) C C DETERMINATION OF DEPTH OF 3D INTERFACES AND ITS DERIVATIVES C FOR BICUBIC POLYNOMIAL APPROXIMATION C SAVE B,DX,DY,llay,ibb,iu,il,ju,jl DIMENSION Y(18),DEP(6),B(16,2),DX(2),DY(2) C COMMON /AUXI/ IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOUT, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori COMMON /AUXX/ MMX(20),MMY(20),MMXY(20) COMMON /INTRF/ Z(1000),SX(350),SY(350),NX(20),NY(20),BRD(6),NINT, 1 XINTA COMMON/ZCOEF/ A02(1000),A20(1000),A22(1000) C IBB=0 LB=2 IF(INTR.EQ.LAY)LB=1 MX=NX(INTR) MY=NY(INTR) DO 1 I=2,MX K=MMX(INTR)+I-1 IF(Y(1).LT.SX(K))GO TO 2 1 CONTINUE 2 I1=K DO 3 I=2,MY K=MMY(INTR)+I-1 IF(Y(2).LT.SY(K))GO TO 4 3 CONTINUE 4 J1=K IF(MAUX.EQ.0) GOTO 8 IF(LB.EQ.2) GOTO 5 IF(I1.EQ.IU.AND.J1.EQ.JU.AND.LLAY.EQ.LAY) GOTO 10 GOTO 8 5 IF(I1.EQ.IL.AND.J1.EQ.JL.AND.LLAY.EQ.LAY) GOTO 10 IL=I1 JL=J1 GOTO 9 8 IU=I1 JU=J1 9 LLAY=LAY I=I1-MMX(INTR) J=J1-MMY(INTR) DX(LB)=SX(I1-1) DY(LB)=SY(J1-1) MM=MMXY(INTR)-1 K=MM+(I-1)*MY+J B20=A20(K) B02=A02(K) B22=A22(K) B00=Z(K) K=MM+I*MY+J C20=A20(K) C02=A02(K) C22=A22(K) C00=Z(K) K=MM+(I-1)*MY+J+1 D20=A20(K) D02=A02(K) D22=A22(K) D00=Z(K) K=MM+I*MY+J+1 E20=A20(K) E02=A02(K) E22=A22(K) E00=Z(K) HX=SX(I1)-DX(LB) HY=SY(J1)-DY(LB) XA=3.*HX YA=3.*HY XB=HX/3. YB=HY/3. D32=(E22-D22)/XA D30=(E20-D20)/XA B30=(C20-B20)/XA B32=(C22-B22)/XA D12=(E02-D02)/HX-XB*(E22+2.*D22) D10=(E00-D00)/HX-XB*(E20+2.*D20) B10=(C00-B00)/HX-XB*(C20+2.*B20) B12=(C02-B02)/HX-XB*(C22+2.*B22) B03=(D02-B02)/YA B13=(D12-B12)/YA B23=(D22-B22)/YA B33=(D32-B32)/YA B01=(D00-B00)/HY-YB*(D02+2.*B02) B11=(D10-B10)/HY-YB*(D12+2.*B12) B21=(D20-B20)/HY-YB*(D22+2.*B22) B31=(D30-B30)/HY-YB*(D32+2.*B32) MAUX=1 B(1,LB)=B00 B(2,LB)=B01 B(3,LB)=B02 B(4,LB)=B03 B(5,LB)=B10 B(6,LB)=B11 B(7,LB)=B12 B(8,LB)=B13 B(9,LB)=B20 B(10,LB)=B21 B(11,LB)=B22 B(12,LB)=B23 B(13,LB)=B30 B(14,LB)=B31 B(15,LB)=B32 B(16,LB)=B33 IBB=1 10 AX=Y(1)-DX(LB) AZ=Y(2)-DY(LB) IF(IBB.EQ.1) GOTO 11 B00=B(1,LB) B01=B(2,LB) B02=B(3,LB) B03=B(4,LB) B10=B(5,LB) B11=B(6,LB) B12=B(7,LB) B13=B(8,LB) B20=B(9,LB) B21=B(10,LB) B22=B(11,LB) B23=B(12,LB) B30=B(13,LB) B31=B(14,LB) B32=B(15,LB) B33=B(16,LB) 11 AUX1=((B33*AZ+B32)*AZ+B31)*AZ+B30 AUX2=((B23*AZ+B22)*AZ+B21)*AZ+B20 AUX3=((B13*AZ+B12)*AZ+B11)*AZ+B10 AUX4=((B03*AZ+B02)*AZ+B01)*AZ+B00 DEP(1)=((AUX1*AX+AUX2)*AX+AUX3)*AX+AUX4 IF(NDER.EQ.0) RETURN DEP(2)=(3.*AUX1*AX+2.*AUX2)*AX+AUX3 IF(NDER.EQ.1)GO TO 7 DEP(4)=6.*AUX1*AX+2.*AUX2 7 AUX1=(3.*B33*AZ+2.*B32)*AZ+B31 AUX2=(3.*B23*AZ+2.*B22)*AZ+B21 AUX3=(3.*B13*AZ+2.*B12)*AZ+B11 AUX4=(3.*B03*AZ+2.*B02)*AZ+B01 DEP(3)=((AUX1*AX+AUX2)*AX+AUX3)*AX+AUX4 IF(NDER.EQ.1) RETURN DEP(5)=(3.*AUX1*AX+2.*AUX2)*AX+AUX3 AUX1=3.*B33*AZ+B32 AUX2=3.*B23*AZ+B22 AUX3=3.*B13*AZ+B12 AUX4=3.*B03*AZ+B02 DEP(6)=2.*(((AUX1*AX+AUX2)*AX+AUX3)*AX+AUX4) RETURN END C C ********************************************************** C SUBROUTINE DISPL(P,U,I1) C C ROUTINE FOR THE COMPUTATION OF DISPLACEMENT VECTOR OF GENERATED C WAVES AT AN INTERFACE. IT WORKS EVEN FOR A COMPLEX REFRACTION VECTOR. C ROUTINES CHRM AND POLAR WORK ONLY FOR REAL REFRACTION VECTORS C C INPUT PARAMETERS: C P(1-3)... SLOWNESS VECTOR C I1... TYPE OF THE WAVE WITH THE SLOWNESS VECTOR P C I1=1... REFLECTED S WAVE (SV COMPONENT IN ISOTROPIC CASE) C I1=2... REFLECTED S WAVE (SH COMPONENT IN ISOTROPIC CASE) C I1=3... REFLECTED P WAVE C I1=4... TRANSMITTED S WAVE (SV COMPONENT IN ISOTROPIC CASE) C I1=5... TRANSMITTED S WAVE (SH COMPONENT IN ISOTROPIC CASE) C I1=6... TRANSMITTED P WAVE C C OUTPUT PARAMETERS: C U(1-3)... CORRESPONDING DISPLACEMENT VECTOR C complex uxn(3,3) COMPLEX U,U1,U2,U3,P,P1,P2,P3,P2P3,P1P2,P1P3,P1P1,P2P2,P3P3,C11, 1 C12,C13,C22,C23,C33,C11N,C22N,C33N,C23SQ,C13SQ,C12SQ,CD11,CD12, 2 CD13,CD22,CD23,CD33,CZ01,CZ02,CZ03,CDTR DIMENSION U(3),P(3),xx(3) COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33, 1 A34,A35,A36,A44,A45,A46,A55,A56,A66, 1 DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23, 1 DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44, 1 DXA45,DXA46,DXA55,DXA56,DXA66, 1 DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23, 1 DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44, 1 DYA45,DYA46,DYA55,DYA56,DYA66, 1 DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23, 1 DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44, 1 DZA45,DZA46,DZA55,DZA56,DZA66, 1 A2546,A1266,A1355,A1456,A3645,A2344 COMMON /AUXI/ IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOUT, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori COMPLEX PS COMMON /RAY/ AY(28,2000),DS(20,20),KINT(20),HHH(3,3),tmax, 1 PS(3,7,20),IS(8,20),N,IREF,IND,IND1 C C C THIS PART CORRESPONDS TO CHRM ,BUT FOR COMPLEX REFRACTION VECTORS C 900 P1=P(1) P2=P(2) P3=P(3) 910 continue IF(ICOEF.GT.0)WRITE(LOUT,'(A,3(/2F12.6))')'DISPL: START P', 1p1,p2,p3 P2P3=P2*P3 P1P2=P1*P2 P1P3=P1*P3 P1P1=P1*P1 P2P2=P2*P2 P3P3=P3*P3 C11=P1P1*A11+P2P2*A66+P3P3*A55 1+2.*P2P3*A56+2.*P1P3*A15+2.*P1P2*A16 C22=P1P1*A66+P2P2*A22+P3P3*A44 1+2.*P2P3*A24+2.*P1P3*A46+2.*P1P2*A26 C33=P1P1*A55+P2P2*A44+P3P3*A33 1+2.*P2P3*A34+2.*P1P3*A35+2.*P1P2*A45 C23=P1P1*A56+P2P2*A24+P3P3*A34 1 +P2P3*A2344+P1P3*A3645+P1P2*A2546 C13=P1P1*A15+P2P2*A46+P3P3*A35 1 +P2P3*A3645+P1P3*A1355+P1P2*A1456 C12=P1P1*A16+P2P2*A26+P3P3*A45 1 +P2P3*A2546+P1P3*A1456+P1P2*A1266 C11N=C11-1. C22N=C22-1. C33N=C33-1. C23SQ=C23*C23 C13SQ=C13*C13 C12SQ=C12*C12 CD11=C22N*C33N-C23SQ CD22=C11N*C33N-C13SQ CD33=C11N*C22N-C12SQ CD12=C13*C23-C12*C33N CD13=C12*C23-C13*C22N CD23=C12*C13-C23*C11N CDTR=CD11+CD22+CD33 IF(ICOEF.GT.0) 1WRITE(LOUT,'(A,4(/4F12.6))')'DISPL: CDIJ' 1,CD11,CD12,CD13,CD22,CD23,CD33,CDTR if(abs(cdtr).lt..0000001)write(lout,'(A)')'DISPL: SHEAR WAVE 1 SINGULARITY IN CALCULATION OF R/T COEFFICIENT' C C THIS PART CORRESPONDS TO ROUTINE POLAR BUT FOR COMPLEX POLARIZATION C VECTORS C U1=CD11 U2=CD12 U3=CD13 uxn(1,1)=u1 uxn(1,2)=u2 uxn(1,3)=u3 XN=real(U1*CONJG(U1)+U2*CONJG(U2)+U3*CONJG(U3)) xx(1)=xn IF(ICOEF.GT.0)WRITE(LOUT,'(A,F14.7)')'DISPL: XN',XN U1=CD12 U2=CD22 U3=CD23 uxn(2,1)=u1 uxn(2,2)=u2 uxn(2,3)=u3 XN=real(U1*CONJG(U1)+U2*CONJG(U2)+U3*CONJG(U3)) xx(2)=xn IF(ICOEF.GT.0)WRITE(LOUT,'(A,F14.7)')'DISPL: XN',XN U1=CD13 U2=CD23 U3=CD33 uxn(3,1)=u1 uxn(3,2)=u2 uxn(3,3)=u3 XN=real(U1*CONJG(U1)+U2*CONJG(U2)+U3*CONJG(U3)) xx(3)=xn IF(ICOEF.GT.0)WRITE(LOUT,'(A,F14.7)')'DISPL: XN',XN xn=0. do 920 l=1,3 if(xn.ge.xx(l))go to 920 xn=xx(l) nx=l 920 continue xn=1./sqrt(xn) u1=uxn(nx,1)*xn u2=uxn(nx,2)*xn u3=uxn(nx,3)*xn U(1)=U1 U(2)=U2 U(3)=U3 IF(ICOEF.GT.0)WRITE(LOUT,'(A,3(/2F12.6))')'DISPL: U',U C C CHECK OF PRECISION C CZ01=C11N*u1+C12*u2+C13*u3 CZ02=C12*u1+C22N*u2+C23*u3 CZ03=C13*u1+C23*u2+C33N*u3 IF(ICOEF.GT.0) 1WRITE(LOUT,'(A,3(/2F14.7))')'DISPL: PRECISSION OF DISPL' 1,CZ01,CZ02,CZ03 RETURN END C C ********************************************************** C SUBROUTINE DMAT C C EVALUATES ELEMENTS OF THE MATRIX DIJ C INTEGER CODE COMMON /COD/ CODE(50,2),KREF,KC,ITYPE COMMON/DJK/D11,D12,D13,D22,D23,D33,DTR COMMON/GAM/G11,G12,G13,G22,G23,G33 COMMON/AUXI/ IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOUT, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori COMPLEX PS COMMON /RAY/ AY(28,2000),DS(20,20),KINT(20),HHH(3,3),tmax, 1 PS(3,7,20),IS(8,20),N,IREF,IND,IND1 C G1=G11-1. G2=G22-1. G3=G33-1. D11=G2*G3-G23*G23 D22=G1*G3-G13*G13 D33=G1*G2-G12*G12 D12=G13*G23-G3*G12 D13=G12*G23-G2*G13 D23=G12*G13-G1*G23 DTR=D11+D22+D33 IF(ABS(DTR).LT.0.00001)THEN IND=10 END IF RETURN END C C ********************************************************** C SUBROUTINE DDMAT(DG,g) C C EVALUATES DERIVATIVES OF ELEMENTS OF THE MATRIX DIJ C DIMENSION DG(3,3) COMMON/DDJK/DD11,DD12,DD13,DD22,DD23,DD33,DDTR COMMON/GAM/G11,G12,G13,G22,G23,G33 C DG11=DG(1,1) DG22=DG(2,2) DG33=DG(3,3) DG12=DG(1,2) DG13=DG(1,3) DG23=DG(2,3) G1=G11-1. G2=G22-1. G3=G33-1. DD11=G3*DG22+G2*DG33-2.*G23*DG23+g*(2.-g22-g33) DD22=G3*DG11+G1*DG33-2.*G13*DG13+g*(2.-g33-g11) DD33=G2*DG11+G1*DG22-2.*G12*DG12+g*(2.-g11-g22) DD12=G23*DG13+G13*DG23-DG12*G3-G12*DG33+g*g12 DD13=G23*DG12+G12*DG23-DG13*G2-G13*DG22+g*g13 DD23=G13*DG12+G12*DG13-DG23*G1-G23*DG11+g*g23 DDTR=DD11+DD22+DD33 RETURN END C C ********************************************************** C FUNCTION TR(G) C C TRACE OF PRODUCT OF MATRICES G AND D C DIMENSION G(3,3) COMMON/DJK/ D11,D12,D13,D22,D23,D33,DTR C TR=G(1,1)*D11+G(2,2)*D22+G(3,3)*D33+ 1 2.*(G(1,2)*D12+G(1,3)*D13+G(2,3)*D23) RETURN END C C ********************************************************** C FUNCTION TRD(G) C C TRACE OF PRODUCT OF MATRICES G AND DD C DIMENSION G(3,3) COMMON/DDJK/ DD11,DD12,DD13,DD22,DD23,DD33,DDTR C TRD=G(1,1)*DD11+G(2,2)*DD22+G(3,3)*DD33+ 1 2.*(G(1,2)*DD12+G(1,3)*DD13+G(2,3)*DD23) RETURN END C C ********************************************************** C SUBROUTINE FCT(X,Y,DERY) C C COMPUTATION OF THE RIGHT HAND SIDES OF THE RAY TRACING EQUATIONS C save g,gpp,gpx,gxx,gx,gp,gx1,gx2,gx3,gp1,gp2,gp3 DIMENSION DERY(18),Y(18),AUX(2),AUXX(3,2),VX(3) DIMENSION G(3,3),GPP(3,3),GPX(3,3),GXX(3,3),GX(3),GP(3) DIMENSION GX1(3,3),GX2(3,3),GX3(3,3),GP1(3,3),GP2(3,3),GP3(3,3) COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33, 1 A34,A35,A36,A44,A45,A46,A55,A56,A66, 1 DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23, 1 DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44, 1 DXA45,DXA46,DXA55,DXA56,DXA66, 1 DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23, 1 DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44, 1 DYA45,DYA46,DYA55,DYA56,DYA66, 1 DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23, 1 DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44, 1 DZA45,DZA46,DZA55,DZA56,DZA66, 1 A2546,A1266,A1355,A1456,A3645,A2344 COMMON /APROX1/ E(21),EX(21),EY(21),EZ(21),EXX(21),EXY(21), 1 EXZ(21),EYY(21),EYZ(21),EZZ(21) COMMON /AUXI/ IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOUT, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori INTEGER CODE COMMON /COD/ CODE(50,2),KREF,KC,ITYPE COMPLEX PS COMMON /RAY/ AY(28,2000),DS(20,20),KINT(20),HHH(3,3),tmax, 1 PS(3,7,20),IS(8,20),N,IREF,IND,IND1 COMMON /DJK/ D11,D12,D13,D22,D23,D33,DTR COMMON /DDJK/ DD11,DD12,DD13,DD22,DD23,DD33,DDTR COMMON /ZERO/ RNULL C IF(IANI(LAY).NE.0) GOTO 20 CALL PARDIS(Y,0) IF(ITYPE.LT.3)ITT=16 IF(ITYPE.EQ.3)ITT=1 C C COMPUTATION OF RIGHT HAND SIDES OF RAY TRACING EQUATIONS C IN AN ISOTROPIC LAYER C C P-WAVE FOR ITT=1, S-WAVE FOR ITT=16 C V=E(ITT) VX(1)=EX(ITT) VX(2)=EY(ITT) VX(3)=EZ(ITT) DO 1 I=1,3 DERY(I)=V*Y(3+I) DERY(3+I)=-.5*VX(I)/V 1 CONTINUE IF(NDER.EQ.1)RETURN VXX=EXX(ITT) VXY=EXY(ITT) VXZ=EXZ(ITT) VYY=EYY(ITT) VYZ=EYZ(ITT) VZZ=EZZ(ITT) DO 3 J=1,2 AUX(J)=0. DO 2 I=1,3 II=3+I+3*J AUX(J)=AUX(J)+VX(I)*Y(II) 2 CONTINUE JJ=4+3*J AUXX(1,J)=VXX*Y(JJ)+VXY*Y(JJ+1)+VXZ*Y(JJ+2) AUXX(2,J)=VXY*Y(JJ)+VYY*Y(JJ+1)+VYZ*Y(JJ+2) AUXX(3,J)=VXZ*Y(JJ)+VYZ*Y(JJ+1)+VZZ*Y(JJ+2) 3 CONTINUE DO 4 I=1,3 DERY(I+6)=AUX(1)*Y(I+3)+V*Y(I+12) DERY(I+9)=AUX(2)*Y(I+3)+V*Y(I+15) DERY(I+12)=.5*(VX(I)*AUX(1)-V*AUXX(I,1))/V/V DERY(I+15)=.5*(VX(I)*AUX(2)-V*AUXX(I,2))/V/V 4 CONTINUE RETURN C C C COMPUTATION OF RIGHT HAND SIDES OF RAY TRACING EQUATIONS C IN AN ANISOTROPIC LAYER C C DETERMINATION OF PARAMETERS OF THE MEDIUM C 20 CALL PARDIS(Y,0) CALL CHRM2(Y,G,1) CALL DMAT IF(IND.EQ.10)RETURN CALL PCHRM(Y,GP1,1,1) CALL PCHRM(Y,GP2,2,1) CALL PCHRM(Y,GP3,3,1) CALL CHRM2(Y,GX1,2) CALL CHRM2(Y,GX2,3) CALL CHRM2(Y,GX3,4) GP(1)=TR(GP1)/DTR GP(2)=TR(GP2)/DTR GP(3)=TR(GP3)/DTR GX(1)=TR(GX1)/DTR GX(2)=TR(GX2)/DTR GX(3)=TR(GX3)/DTR DO 5 I=1,3 DERY(I)=.5*GP(I) DERY(I+3)=-.5*GX(I) 5 CONTINUE IF(NDER.EQ.1)RETURN CALL CHRM2(Y,G,1) CALL DDMAT(GP1,gp(1)) GPP(1,1)=TRD(GP1)-GP(1)*DDTR GPP(1,2)=TRD(GP2)-GP(2)*DDTR GPP(1,3)=TRD(GP3)-GP(3)*DDTR GPX(1,1)=TRD(GX1)-GX(1)*DDTR GPX(1,2)=TRD(GX2)-GX(2)*DDTR GPX(1,3)=TRD(GX3)-GX(3)*DDTR CALL DDMAT(GP2,gp(2)) GPP(2,1)=TRD(GP1)-GP(1)*DDTR GPP(2,2)=TRD(GP2)-GP(2)*DDTR GPP(2,3)=TRD(GP3)-GP(3)*DDTR GPX(2,1)=TRD(GX1)-GX(1)*DDTR GPX(2,2)=TRD(GX2)-GX(2)*DDTR GPX(2,3)=TRD(GX3)-GX(3)*DDTR CALL DDMAT(GP3,gp(3)) GPP(3,1)=TRD(GP1)-GP(1)*DDTR GPP(3,2)=TRD(GP2)-GP(2)*DDTR GPP(3,3)=TRD(GP3)-GP(3)*DDTR GPX(3,1)=TRD(GX1)-GX(1)*DDTR GPX(3,2)=TRD(GX2)-GX(2)*DDTR GPX(3,3)=TRD(GX3)-GX(3)*DDTR CALL DDMAT(GX1,gx(1)) GXX(1,1)=TRD(GX1)-GX(1)*DDTR GXX(1,2)=TRD(GX2)-GX(2)*DDTR GXX(1,3)=TRD(GX3)-GX(3)*DDTR CALL DDMAT(GX2,gx(2)) GXX(2,1)=TRD(GX1)-GX(1)*DDTR GXX(2,2)=TRD(GX2)-GX(2)*DDTR GXX(2,3)=TRD(GX3)-GX(3)*DDTR CALL DDMAT(GX3,gx(3)) GXX(3,1)=TRD(GX1)-GX(1)*DDTR GXX(3,2)=TRD(GX2)-GX(2)*DDTR GXX(3,3)=TRD(GX3)-GX(3)*DDTR C DO 11 L=1,3 DO 11 M=L,3 CALL PPCHRM(G,L,M,1) AUX1=TR(G) GPP(M,L)=(GPP(M,L)+AUX1)/DTR IF(L.NE.M)GPP(L,M)=(GPP(L,M)+AUX1)/DTR 11 CONTINUE DO 12 L=1,3 CALL PCHRM(Y,G,L,2) GPX(L,1)=(GPX(L,1)+TR(G))/DTR CALL PCHRM(Y,G,L,3) GPX(L,2)=(GPX(L,2)+TR(G))/DTR CALL PCHRM(Y,G,L,4) GPX(L,3)=(GPX(L,3)+TR(G))/DTR 12 CONTINUE CALL CHRM2(Y,G,5) GXX(1,1)=(GXX(1,1)+TR(G))/DTR CALL CHRM2(Y,G,8) GXX(2,2)=(GXX(2,2)+TR(G))/DTR CALL CHRM2(Y,G,10) GXX(3,3)=(GXX(3,3)+TR(G))/DTR CALL CHRM2(Y,G,6) AUX1=(GXX(1,2)+TR(G))/DTR GXX(1,2)=AUX1 GXX(2,1)=AUX1 CALL CHRM2(Y,G,7) AUX1=(GXX(1,3)+TR(G))/DTR GXX(1,3)=AUX1 GXX(3,1)=AUX1 CALL CHRM2(Y,G,9) AUX1=(GXX(2,3)+TR(G))/DTR GXX(2,3)=AUX1 GXX(3,2)=AUX1 CALL CHRM2(Y,G,1) C DO 13 I=1,3 DERY(I+6)=0. DERY(I+9)=0. DERY(I+12)=0. DERY(I+15)=0. DO 13 K=1,3 DERY(I+6)=DERY(I+6)+.5*(GPX(I,K)*Y(K+6)+GPP(I,K)*Y(K+12)) DERY(I+9)=DERY(I+9)+.5*(GPX(I,K)*Y(K+9)+GPP(I,K)*Y(K+15)) DERY(I+12)=DERY(I+12)-.5*(GXX(I,K)*Y(K+6)+GPX(K,I)*Y(K+12)) DERY(I+15)=DERY(I+15)-.5*(GXX(I,K)*Y(K+9)+GPX(K,I)*Y(K+15)) 13 CONTINUE RETURN END