! LAYERXT calculates dx and dt for a ray in a layer with a linear ! velocity gradient. This is a highly modified version of a ! subroutine in Chris Chapman's WKBJ program. ! ! Inputs: p = horizontal slowness ! h = layer thickness ! utop = slowness at top of layer ! ubot = slowness at bottom of layer ! Returns: dx = range offset ! dt = travel time ! irtr = return code ! = -1, zero thickness layer ! = 0, ray turned above layer ! = 1, ray passed through layer ! = 2, ray turned in layer, 1 leg counted in dx,dt ! subroutine LAYERXT(p,h,utop,ubot,dx,dt,irtr) if (p.ge.utop) then !ray turned above layer dx=0. dt=0. irtr=0 return else if (h.eq.0.) then !zero thickness layer dx=0. dt=0. irtr=-1 return end if u1=utop u2=ubot v1=1./u1 v2=1./u2 b=(v2-v1)/h !slope of velocity gradient eta1=sqrt(u1**2-p**2) if (b.eq.0.) then !constant velocity layer dx=h*p/eta1 dt=h*u1**2/eta1 irtr=1 return end if x1=eta1/(u1*b*p) tau1=(alog((u1+eta1)/p)-eta1/u1)/b if (p.ge.ubot) then !ray turned within layer, dx=x1 !no contribution to integral dtau=tau1 !from bottom point dt=dtau+p*dx irtr=2 return end if irtr=1 eta2=sqrt(u2**2-p**2) x2=eta2/(u2*b*p) tau2=(alog((u2+eta2)/p)-eta2/u2)/b dx=x1-x2 dtau=tau1-tau2 dt=dtau+p*dx return end ! RTCOEF calculates P/SV reflection/transmission coefficients ! for an interface between two solid layers, based on the ! equations on p. 149-150 of Aki and Richards. This version ! is modified from an older routine provided by Tom Sereno. ! ! Inputs: vp1 = P-wave velocity of layer 1 (top layer) ! (real) vs1 = S-wave velocity of layer 1 ! den1 = density of layer 1 ! vp2 = P-wave velocity of layer 2 (bottom layer) ! vs2 = S-wave velocity of layer 2 ! den2 = density of layer 2 ! hslow = horizontal slowness (ray parameter) ! Returns: rt(1) = down P to P up (refl) ! (complex) rt(2) = down P to S up (refl) ! rt(3) = down P to P down (tran) ! rt(4) = down P to S down (tran) ! rt(5) = down S to P up (refl) ! rt(6) = down S to S up (refl) ! rt(7) = down S to P down (tran) ! rt(8) = down S to S down (tran) ! rt(9) = up P to P up (tran) ! rt(10) = up P to S up (tran) ! rt(11) = up P to P down (refl) ! rt(12) = up P to S down (refl) ! rt(13) = up S to P up (tran) ! rt(14) = up S to S up (tran) ! rt(15) = up S to P down (refl) ! rt(16) = up S to S down (refl) ! ! NOTE: All input variables are real. ! All output variables are complex! ! Coefficients are not energy normalized. ! SUBROUTINE RTCOEF(vp1,vs1,den1,vp2,vs2,den2,hslow,rt) implicit complex (a-h,o-z) complex rt(16) real vp1,vs1,den1,vp2,vs2,den2,hslow alpha1=cmplx(vp1,0.) beta1=cmplx(vs1,0.) rho1=cmplx(den1,0.) alpha2=cmplx(vp2,0.) beta2=cmplx(vs2,0.) rho2=cmplx(den2,0.) p=cmplx(hslow,0.) cone=cmplx(1.,0.) ctwo=cmplx(2.,0.) term1=(cone-ctwo*beta1**2*p**2) term2=(cone-ctwo*beta2**2*p**2) a=rho2*term2-rho1*term1 b=rho2*term2+ctwo*rho1*beta1**2*p**2 c=rho1*term1+ctwo*rho2*beta2**2*p**2 d=ctwo*(rho2*beta2**2-rho1*beta1**2) ! compute signs and cosines, allowing for complex incidence angles si1=alpha1*p si2=alpha2*p sj1=beta1*p sj2=beta2*p ci1=csqrt(cone-si1**2) ci2=csqrt(cone-si2**2) cj1=csqrt(cone-sj1**2) cj2=csqrt(cone-sj2**2) E=b*ci1/alpha1+c*ci2/alpha2 F=b*cj1/beta1+c*cj2/beta2 G=a-d*ci1*cj2/(alpha1*beta2) H=a-d*ci2*cj1/(alpha2*beta1) DEN=E*F+G*H*p**2 trm1=b*ci1/alpha1-c*ci2/alpha2 trm2=a+d*ci1*cj2/(alpha1*beta2) rt(1)=(trm1*F-trm2*H*p**2)/DEN !refl down P to P up trm1=a*b+c*d*ci2*cj2/(alpha2*beta2) rt(2)=(-ctwo*ci1*trm1*p)/(beta1*DEN) !refl down P to S up rt(3)=ctwo*rho1*ci1*F/(alpha2*DEN) !trans down P to P down rt(4)=ctwo*rho1*ci1*H*p/(beta2*DEN) !trans down P to S down trm1=a*b+c*d*ci2*cj2/(alpha2*beta2) rt(5)=(-ctwo*cj1*trm1*p)/(alpha1*DEN) !refl down S to P up trm1=b*cj1/beta1-c*cj2/beta2 trm2=a+d*ci2*cj1/(alpha2*beta1) rt(6)=-(trm1*E-trm2*G*p**2)/DEN !refl down S to S up rt(7)=-ctwo*rho1*cj1*G*p/(alpha2*DEN) !trans down S to P down rt(8)=ctwo*rho1*cj1*E/(beta2*DEN) !trans down S to S down rt(9)=ctwo*rho2*ci2*F/(alpha1*DEN) !trans up P to P up rt(10)=-ctwo*rho2*ci2*G*p/(beta1*DEN) !trans up P to S up trm1=b*ci1/alpha1-c*ci2/alpha2 trm2=a+d*ci2*cj1/(alpha2*beta1) rt(11)=-(trm1*F+trm2*G*p**2)/DEN !refl up P to P down trm1=a*c+b*d*ci1*cj1/(alpha1*beta1) rt(12)=(ctwo*ci2*trm1*p)/(beta2*DEN) !refl up P to S down rt(13)=ctwo*rho2*cj2*H*p/(alpha1*DEN) !trans up S to P up rt(14)=ctwo*rho2*cj2*E/(beta1*DEN) !trans up S to S up trm1=a*c+b*d*ci1*cj1/(alpha1*beta1) rt(15)=(ctwo*cj2*trm1*p)/(alpha2*DEN) !refl up S to P down trm1=b*cj1/beta1-c*cj2/beta2 trm2=a+d*ci1*cj2/(alpha1*beta2) rt(16)=(trm1*E+trm2*H*p**2)/DEN !refl up S to S down return end ! GETAUX returns auxiliary fault plane, given strike,dip,rake ! of main fault plane. ! ! Inputs: strike1, dip1, rake1 (degrees, primary fault plane) ! Returns: strike2, dip2, rake2 (degrees, auxiliary fault plane) ! subroutine GETAUX(strike1,dip1,rake1,strike2,dip2,rake2) degrad=180./3.1415927 s1=strike1/degrad d1=dip1/degrad r1=rake1/degrad d2=acos(sin(r1)*sin(d1)) sr2=cos(d1)/sin(d2) cr2=-sin(d1)*cos(r1)/sin(d2) r2=atan2(sr2,cr2) s12=cos(r1)/sin(d2) c12=-1./(tan(d1)*tan(d2)) s2=s1-atan2(s12,c12) strike2=s2*degrad dip2=d2*degrad rake2=r2*degrad if (dip2.gt.90.) then strike2=strike2+180. dip2=180.-dip2 rake2=360.-rake2 end if if (strike2.gt.360.) strike2=strike2-360. return end