C
C Subroutine file 'coef.for' - reflection/transmission coefficients C C Date: 1989, December 8 C Coded by: Ludek Klimes C C....................................................................... C C This file consists of: C COEFSH..Subroutine designed to evaluate SH R/T coefficients C (see C.R.T.5.9.7). C COEFSH C COEF50..Subroutine designed to evaluate P-SV R/T coefficients C (see C.R.T.5.9.7). C COEF50 C C C======================================================================= C C C SUBROUTINE COEFSH(P,VS1,RO1,VS2,RO2,NCODE,RMOD,RPH) C C ****************************************************************** C C The routine COEFSH is designed for the computation of reflection C and transmission coefficients at a plane interface between two C homogeneous solid halfspaces or at a free surface of a homogeneous C solid halfspace. C C I n p u t p a r a m e t e r s C P...Ray parameter C VS1,RO1...Parameters of the first halfspace C VS2,RO2...Parameters of second halfspace. For the free C surface take RO2=0. And arbitrary C value of VS2 C NCODE...Code of the computed coefficient C S1S1...NCODE=1 C S1S2...NCODE=2 C C O u t p u t p a r a m e t e r s C RMOD,RPH...Modul and argument of the coefficient C C N o t e s C 1/ Time factor of incident wave...EXP(-i*OMEGA*T) C 2/ Formulae are taken from Cerveny,Molotkov and Psencik, C Ray Method in Seismology, page 34 C C ****************************************************************** C COMPLEX A,B X= 1.-P*P*VS1*VS1 Y= 1.-P*P*VS2*VS2 C= RO1*VS1*SQRT(ABS(X)) D= RO2*VS2*SQRT(ABS(Y)) A= CMPLX(C,0.) IF(X.LT.0.) A= CMPLX(0.,C) B= CMPLX(D,0.) IF(Y.LT.0.) B= CMPLX(0.,D) GO TO (1,2), NCODE 1 A= (A-B)/(A+B) GO TO 3 2 A= (A+A)/(A+B) 3 RMOD= SQRT(REAL(A)*REAL(A)+AIMAG(A)*AIMAG(A)) RPH= ATAN2(AIMAG(A),REAL(A)) RETURN END C C======================================================================= C C C SUBROUTINE COEF50(P,VP1,VS1,RO1,VP2,VS2,RO2,NCODE,ND,RMOD,RPH) C C ****************************************************************** C C The routine COEF50 is designed to evaluate the reflection and C transmission coefficients at a plane interface between two C homogeneous solid halfspaces or at a free surface of a homogeneous C solid halfspace. C C The kinds of individual coefficients are specified by the C following numbers C a/ Interface between two solid halfspaces C P1P1...1 P1S1...2 P1P2...3 P1S2...4 C S1P1...5 S1S1...6 S1P2...7 S1S2...8 C b/ Free surface (for RO2.LT.0.00001) C PP.....1 PX.....5 PX,PZ...X- and Z- components of the C PS.....2 PZ.....6 coef.of conversion,incident P wave C SP.....3 SX.....7 SX,SZ...X- and Z- components of the C SS.....4 SZ.....8 coef.of conversion,incident S wave C C I n p u t p a r a m e t e r s C P...Ray parameter C VP1,VS1,RO1...Parameters of the first halfspace C VP2,VS2,RO2...Parameters of second halfspace. For the free C surface take RO2.LT.0.000001,eg.RO2=0., and C arbitrary values of VP2 and VS2 C NCODE...Code of the computed coefficient C ND...=0 when the interface is situated at the right-hand C side of the incident ray (X against P) C =1 when the interface is situated at the left-hand C side of the incident ray (X along P) C C O u t p u t p a r a m e t e r s C RMOD,RPH...Modul and argument of the coefficient C C N o t e s C 1/ Positive P...In the direction of propagation C 2/ Positive S...To the left from P C 3/ Positive X...To the right from Z (to the right from P) C 4/ Positive Z...In the direction of P C 5/ Time factor of incident wave...EXP(-i*OMEGA*T) C 6/ Formulae are taken from Cerveny and Ravindra, Theory of Seismic C Head Waves,Pages 63-67. Due to the note 2, the signs at certain C coefficients are opposite (and time factor is changed, L.K.) C (signs of conversion coefficients are opposite to Cerveny, C Molotkov and Psencik, Ray Method in Seismology, page 35) C C ****************************************************************** C COMPLEX B(4),RR,C1,C2,C3,C4,H1,H2,H3,H4,H5,H6,H,HB,HC DIMENSION PRMT(4),D(4),DD(4) PRMT(1)=VP1 PRMT(2)=VS1 PRMT(3)=VP2 PRMT(4)=VS2 IF(RO2.LT.0.000001)GO TO 150 A1=VP1*VS1 A2=VP2*VS2 A3=VP1*RO1 A4=VP2*RO2 A5=VS1*RO1 A6=VS2*RO2 Q=2.*(A6*VS2-A5*VS1) PP=P*P QP=Q*PP X=RO2-QP Y=RO1+QP Z=RO2-RO1-QP G1=A1*A2*PP*Z*Z G2=A2*X*X G3=A1*Y*Y G4=A4*A5 G5=A3*A6 G6=Q*Q*PP DO 21 I=1,4 DD(I)=P*PRMT(I) 21 D(I)=SQRT(ABS(1.-DD(I)*DD(I))) IF(DD(1).LE.1..AND.DD(2).LE.1..AND.DD(3).LE.1..AND.DD(4).LE.1.) 1GO TO 100 C C Complex coefficients DO 22 I=1,4 IF(DD(I).GT.1.)GO TO 23 B(I)=CMPLX(D(I),0.) GO TO 22 23 B(I)= CMPLX(0.,D(I)) 22 CONTINUE C1=B(1)*B(2) C2=B(3)*B(4) C3=B(1)*B(4) C4=B(2)*B(3) H1=CMPLX(G1,0.) H2=G2*C1 H3=G3*C2 H4=G4*C3 H5=G5*C4 H6=G6*C1*C2 H=1./(H1+H2+H3+H4+H5+H6) HB=2.*H HC=HB*P GO TO (1,2,3,4,5,6,7,8),NCODE 1 RR=H*(H2+H4+H6-H1-H3-H5) GO TO 26 2 RR=VP1*B(1)*HC*(Q*Y*C2+A2*X*Z) IF(ND.NE.0)RR=-RR GO TO 26 3 RR=A3*B(1)*HB*(VS2*B(2)*X+VS1*B(4)*Y) GO TO 26 4 RR=-A3*B(1)*HC*(Q*C4-VS1*VP2*Z) IF(ND.NE.0)RR=-RR GO TO 26 5 RR=-VS1*B(2)*HC*(Q*Y*C2+A2*X*Z) IF(ND.NE.0)RR=-RR GO TO 26 6 RR=H*(H2+H5+H6-H1-H3-H4) GO TO 26 7 RR=A5*B(2)*HC*(Q*C3-VP1*VS2*Z) IF(ND.NE.0)RR=-RR GO TO 26 8 RR=A5*B(2)*HB*(VP1*B(3)*Y+VP2*B(1)*X) GO TO 26 C Real coefficients 100 E1=D(1)*D(2) E2=D(3)*D(4) E3=D(1)*D(4) E4=D(2)*D(3) S1=G1 S2=G2*E1 S3=G3*E2 S4=G4*E3 S5=G5*E4 S6=G6*E1*E2 S=1./(S1+S2+S3+S4+S5+S6) SB=2.*S SC=SB*P GO TO (101,102,103,104,105,106,107,108),NCODE 101 R=S*(S2+S4+S6-S1-S3-S5) GO TO 250 102 R=VP1*D(1)*SC*(Q*Y*E2+A2*X*Z) IF(ND.NE.0)R=-R GO TO 250 103 R=A3*D(1)*SB*(VS2*D(2)*X+VS1*D(4)*Y) GO TO 250 104 R=-A3*D(1)*SC*(Q*E4-VS1*VP2*Z) IF(ND.NE.0)R=-R GO TO 250 105 R=-VS1*D(2)*SC*(Q*Y*E2+A2*X*Z) IF(ND.NE.0)R=-R GO TO 250 106 R=S*(S2+S5+S6-S1-S3-S4) GO TO 250 107 R=A5*D(2)*SC*(Q*E3-VP1*VS2*Z) IF(ND.NE.0)R=-R GO TO 250 108 R=A5*D(2)*SB*(VP1*D(3)*Y+VP2*D(1)*X) GO TO 250 C C Earths surface,complex coefficients and coefficients of conversion 150 A1=VS1*P A2=A1*A1 A3=2.*A2 A4=2.*A1 A5=A4+A4 A6=1.-A3 A7=2.*A6 A8=2.*A3*VS1/VP1 A9=A6*A6 DD(1)=P*VP1 DD(2)=P*VS1 DO 151 I=1,2 151 D(I)=SQRT(ABS(1.-DD(I)*DD(I))) IF(DD(1).LE.1..AND.DD(2).LE.1.)GO TO 200 DO 154 I=1,2 IF(DD(I).GT.1.)GO TO 155 B(I)=CMPLX(D(I),0.) GO TO 154 155 B(I)= CMPLX(0.,D(I)) 154 CONTINUE H1=B(1)*B(2) H2=H1*A8 H=1./(A9+H2) GO TO (161,162,163,164,165,166,167,168),NCODE 161 RR=(-A9+H2)*H GO TO 26 162 RR=-A5*B(1)*H*A6 IF(ND.NE.0)RR=-RR GO TO 26 163 RR=A5*B(2)*H*A6*VS1/VP1 IF(ND.NE.0)RR=-RR GO TO 26 164 RR=-(A9-H2)*H GO TO 26 165 RR=-A5*H1*H IF(ND.NE.0)RR=-RR GO TO 26 166 RR=A7*B(1)*H GO TO 26 167 RR=-A7*B(2)*H GO TO 26 168 RR=-A5*VS1*H1*H/VP1 IF(ND.NE.0)RR=-RR 26 Z2=REAL(RR) Z3=AIMAG(RR) IF(Z2.EQ.0..AND.Z3.EQ.0.)GO TO 157 RMOD=SQRT(Z2*Z2+Z3*Z3) RPH=ATAN2(Z3,Z2) RETURN 157 RMOD=0. RPH=0. RETURN C C Earths surface,real coefficients and coefficients of conversion 200 S1=D(1)*D(2) S2=A8*S1 S=1./(A9+S2) GO TO (201,202,203,204,205,206,207,208),NCODE 201 R=(-A9+S2)*S GO TO 250 202 R=-A5*D(1)*S*A6 IF(ND.NE.0)R=-R GO TO 250 203 R=A5*D(2)*S*A6*VS1/VP1 IF(ND.NE.0)R=-R GO TO 250 204 R=(S2-A9)*S GO TO 250 205 R=-A5*S1*S IF(ND.NE.0)R=-R GO TO 250 206 R=A7*D(1)*S GO TO 250 207 R=-A7*D(2)*S GO TO 250 208 R=-A5*VS1*S1*S/VP1 IF(ND.NE.0)R=-R 250 IF(R.LT.0.)GO TO 251 RMOD=R RPH=0. RETURN 251 RMOD=-R RPH=3.14159 RETURN END C C======================================================================= C