C
C Routine COEF52 C ************** C C Routine COEF52 is designed for the computation of the C displacement reflection/transmission coefficients (R/T C coefficients) of inhomogeneous P, SV and SH plane waves at a C stack of homogeneous isotropic dissipative layers between two C homogeneous isotropic dissipative halfspaces. The stack of C layers may be reduced to a single interface. C As special cases, it is also possible to compute the R/T C coefficients of pressure waves at a stack of liquid layers C between two liquid halfspaces, the R/T coefficients at a C surficial stack of layers, and the conversion coefficients C (receiver functions) at a free surface of the surficial C stack of layers. In all these cases, analogous R/T coefficients C for non-dissipative media may be also computed. C The routine is written in Fortran 77 programming language. C C Calling: C ******** C CALL COEF52(NC,NH,NQ,NUM,LIN,LOU,ANGLE,GAMMA,FREQ,RMOD,RPHASE, C RMOD0,RPH0) C Input parameters: C NC ..... Type of the R/T coefficient C NH ..... NH=0 ... the second halfspace is solid or liquid C NH=1 ... the second halfspace is a vacuum C NQ ..... NQ=0 ... model is non-dissipative C NQ=1 ... model is dissipative C NUM .... NUM=0 ... reading the model, no computations C NUM=1 ... computations of R/T coefficients C LIN .... Number of input file logical unit C LOU .... Number of output file logical unit C ANGLE .. Angle of incidence (in degrees) C GAMMA .. Attenuation angle of the incident inhomogeneous C plane wave (in degrees) C FREQ ... Frequency (in Hz) C Output parameters: C RMOD ... Modul of computed R/T coefficient C RPHASE.. Phase of computed R/T coefficient C RMOD0 .. Modul of analogous R/T coefficient for a single interface C between two non-dissipative halfspaces C RPH0 ... Phase of analogous R/T coefficient. C For NUM=0, all other parameters, with the exception of LIN and C LOU, may be arbitrary. C C Specification of the model: C *************************** C The number of layers in the model (including both C halfspaces) is specified by NZ. For NZ=2, the two halfspaces C are in contact and the transition layer reduces to a single C interface. The layer number 1 corresponds to the incidence C (first) halfspace, the layer number NZ to the second halfspace. C For each layer (including halfspaces), the following C real-valued parameters should be specified: C VP ... P velocity (in km/s), C VS ... S velocity (in km/s), C RHO .. density (in g/cm**3), C QP ... quality factor of P waves, C QS ... quality factor of S waves, C D ... thickness of the layer (in km). C The quality factors QP and QS of P and S waves are assumed C to be independent of frequency FREQ (constant-Q model). C The real-valued velocities VP and VS of P and S waves are C material properties of the model for the reference frequency C FREF. They correspond to the real parts of the Lame's elastic C moduli for the reference frequency FREF. For known VP(FREF), C VS(FREF), QP and QS, the frequency-dependent complex valued C velocities CVP(FREQ) and CVS(FREQ) of P and S waves for the C constant Q model are given by relations C CVP(FREQ) = VP(FREF)*(1.+ln(FREQ/FREF)/(pi*QP) - i/(2.*QP)), C CVS(FREQ) = VS(FREF)*(1.+ln(FREQ/FREF)/(pi*QS) - i/(2.*QS)). C The frequency-dependent real-valued velocities VP(FREQ) and C VS(FREQ) equal the real parts of CVP(FREQ) and CVS(FREQ). C They specify the velocity dispersion. C In both halfspaces, the thicknesses (D(1) and D(NZ)) may C be taken arbitrarily; they are automatically adjusted to zero. C For pressure waves in liquid media, VS and QS are not used C in the computations. The values of VS and QS may be specified C arbitrarily in this case. C For SH waves, VP and QP are not used in the computations. C The values of VP and QP may be specified arbitrarily in this case. C For non-dissipative media (NQ=0), QP and QS are not used C in computations. The values of QP and QS may be specified C arbitrarily in this case. C If the second halfspace represents a vacuum (NH=1), C the density in the second halfspace is automatically adjusted C to zero. Similarly, the velocities VP and VS are adjusted to 0.001. C For NUM=0, routine COEF52 reads the input data for the C model (including the reference frequency FREF), and does not C perform any computations. For NUM=1, the computation of the R/T C coefficients is performed, for the known model. C C Type of R/T coefficients. C ************************* C The type of the R/T coefficient we wish to compute is C specified by NC (input parameter of COEF52). The following C table shows the values of NC for individual R/T coefficients: C a) P/SV waves, solids: C P1P1...1 P1S1...2 P1P2...3 P1S2...4 C S1P1...5 S1S1...6 S1P2...7 S1S2...8 C b) SH waves, solids: C S1S1...9 S1S2...10 C c) Pressure waves, liquids: C P1P1...13 P1P2...14 C For example, P1S1 (NC=2) corresponds to the reflection C coefficient, for incident P wave and reflected S wave. C Similarly, S1P2 (NC=7) corresponds to the transmission C coefficient, for incident S wave and transmitted P wave. C The second halfspace may be solid or liquid (NH=0) or C a vacuum (NH=1). For NH=1, the interface between the (NZ-1)-th C layer and the NZ-th layer (second halfspace) represents a free C surface, e.g. the surface of the Earth. For NH=1, the reflection C coefficients P1P1, P1S1, S1P1 and S1S1 have a standard meaning. C The physical meaning of the transmission coefficients P1P2, C P1S2, S1P2 and S1S2 for NH=1 is, however, different. They C give the so-called conversion coefficients, also called receiver C functions: C For NH=1: P1P2=PZ, P1S2=PX, S1P2=SZ, S1S2=SX. C The conversion coefficients (receiver functions) represent C horizontal or vertical displacement components of the free surface C (top of the stack of layers) due to a plane wave incident on the C bottom of the stack of layers from below. In the notation PZ, PX, C SZ and SX for the conversion coefficients, the first letter C specifies the incident wave (P or S), and the second letter the C Cartesian component of the displacement vector (X or Z). The C Cartesian axis X is horizontal, tangential to the free surface, C situated in the plane of incidence, oriented "against" the C direction of propagation. (For the X-axis oriented "along" the C direction of propagation, we must take opposite signs of PX and C SX.) The Cartesian axis Z is vertical, perpendicular to the C interface, and oriented away from the transition layer. For SH C waves and NH=1, the R/T coefficient S1S2 represents the conversion C coefficient SY. The Cartesian axis Y is horizontal, perpendicular C to the plane of incidence. C C Incident wave: C ************** C In dissipative media, the slowness vector of the incident C wave is complex-valued. The real part of the slowness vector C is called the propagation vector, and its imaginary part the C attenuation vector. The plane of incidence is specified by the C normal vector to the interface(s) and by the propagation vector. C The attenuation vector of the incident wave is assumed to be C situated in the plane of incidence. The direction of the C propagation vector is specified by the real-valued angle of C incidence ANGLE, and the direction of the attenuation vector is C specified by the real-valued attenuation angle GAMMA. The angles C ANGLE and GAMMA are input parameters of COEF52. ANGLE represents C the acute angle between the propagation vector of the incident C wave and normal to the interface. It is specified in degrees, and C should be in the range <0.,90.>. The attenuation angle is the C angle between the propagation and attenuation vectors of the C incident wave, is specified in degrees, and should be in the C range (-90.,90.). It is positive if the attenuation vector is C pointing to the left from the propagation vector. C The complex-valued polarisation vector of the incident P C wave is proportional to the slowness vector. The complex-valued C polarisation vector of the incident SV wave is perpendicular to C the slowness vector (in the complex sense), and its real part is C oriented to the left from the propagation vector. The polarisation C vector of the incident SH wave is perpendicular to the plane of C incidence. Analogous orientation of polarisation vectors is used C even for R/T waves. Consequently, the R/T coefficients are C independent of the used coordinate system. C For given NC, NH, NQ, ANGLE, GAMMA, and FREQ, and for a C given model, the routine COEF52 returns the complex-valued C R/T coefficient RCOEF. RCOEF is expressed in terms of its modulus C RMOD and phase RPHASE. Both RMOD and RPHASE are real-valued. C RPHASE is specified in the range <-180.,180.>. C The time factor of the incident wave is exp(-i*omega*t), C where t is running time and omega=2.*pi*FREQ. For the time factor C exp(i*omega*t), it would be necessary to change the sign of C RPHASE. C For RMOD=0., RPHASE cannot be computed. Consequently, C COEF52 returns 'the warning' RPHASE=999 for RMOD<0.00001. In C plotting RPHASE, the points with RPHASE=999 should be ignored. C In each computation, also quantities RMOD0 and RPH0 are C computed, in addition to RMOD and RPHASE. They correspond to the C same R/T coefficient as RMOD and RPHASE, but to the plane C interface between the two non-dissipative halfspaces (NZ=2, NQ=0). C C Routines used C ************* C For matrix computations of R/T coefficients for P and SV C waves, routine RTMATQ is used in dissipative media, and routine C RTMAT in non-dissipative media. They further use routines RTQ, C RT, etc. C For matrix computation of R/T coefficients of SH waves and C for pressure waves in liquids, routine RTMQ2 is used in C dissipative media, and routine RTM2 in non-dissipative media. C For comparisons, R/T coefficients at a single interface C between two non-dissipative media (NZ=2, NQ=0), are computed C in routine COEF8, using explicit (non-matrix) expressions for C R/T coefficients. C Routine CRITAN tests a possible existence of nonelastic C discrete critical angles, and computes them if they exist. C Nonelastic discrete critical angles correspond to angles of C incidence for which both the real and imaginary parts of the C vertical component of the slowness vector of any R/T wave in the C dissipative model under consideration vanish. The R/T coefficients C may display anomalous jumps at the discrete critical angles. The C application of CRITAN removes such jumps, see Brokesova and C Cerveny (1998). In COEF52, routine CRITAN is applied only if NZ=2, C not for NZ.GT.2. Note that the anomalous behaviour of R/T C coefficients was also observed in the so-called Rayleigh window. C These cases will be subjects of further investigation. C Routines RTMAT and RTMATQ are modifications of analogous C routines, written by G. Muller, used in his reflectivity packages, C and described in the paper Muller (1985). The recursive algorithm C is used to compute the R/T coefficients. The potential R/T C coefficients, computed by Muller, are adjusted to displacement C R/T coefficients, considered here. Although G.Muller C considers the time convention exp(i*omega*t), the resulting phase C is adjusted to our time convention exp(-i*omega*t). Moreover, C inhomogeneous plane incident waves with arbitrary attenuation C angles are introduced. The routines are used in COEF52 with the C author's permission. C Routine COEF8 is a modification of an analogous routine, C used in the program package SEIS88 by V. Cerveny and I. Psencik. C See explicit analytical expressions for the coefficients in C Cerveny (2001). The book also displays pictures of many R/T C coefficients at a single interface between two non-dissipative C halfspaces. C Routine COEF52 itself is a modification of the routine FDEP, C used in the program package BEAM87 by V.Cerveny to study the C propagation of Gaussian beams in complex, 2-D, laterally varying C layered structures, containing thin stacks of layers. See C Cerveny (1989). C The theoretical treatment of reflection and transmission of C inhomogeneous plane waves at a single interface between two C dissipative elastic halfspaces can be found in Brokesova and C Cerveny (1998). The reference also offers a more detailed C description of routine COEF52 and presents numerous examples C of R/T coefficients. See also Brokesova (2001). C C Demonstration program for testing C ********************************* C For testing purposes, a brief main program C RTCOEF is included. C C References C ********** C Brokesova,J. (2000). Reflection/transmission coefficients at a C plane interface in dissipative and non-dissipative media: C A comparison. J.Comput.Acoustics, 9,623 -641. C Brokesova,J., and Cerveny,V. (1998). Inhomogeneous plane waves C in dissipative, isotropic and anisotropic media. Reflection/ C transmission coefficients. In Seismic waves in complex 3-D C structures, Report No. 7, pp. 57 - 146. Department of C Geophysics, Charles University, Prague. C Cerveny,V. (1989). Synthetic body wave seismograms for laterally C varying media containing thin transition layers. Geophys. J. C Int., 99, 331-349, C Cerveny,V. (2001). Seismic ray theory. Cambridge Univ. Press, C Cambridge. C Muller,G. (1985). The reflectivity method. A tutorial. J.Geophys., C 58, 153-174. C C C Consortium Project "Seismic Waves C in Complex 3-D Structures" C Praha, April 2003 C J.Brokesova, V.Cerveny C ********************************************************************* c subroutine coef52(nc,nh,nq,num,lin,lou,angle,gamma,freq,rmod, 1rphase,rmod0,rph0) complex rpp,rps,rss,rsp,tpp,tps,tss,tsp,rcoef,rp,tp,rs,ts,cunit, 1u common/model/nz,vp(50),vs(50),rho(50),d(50),qp(50),qs(50),fref, 1rho2,vp2,vs2,qp2,qs2 c cunit=(0.,1.) pi=3.141593 ang=angle*pi/180. gam=gamma*pi/180. inc=0 if(nc.gt.4.and.nc.lt.11)inc=1 ntype=2 if(nc.ge.13)ntype=0 if(nc.gt.8.and.nc.lt.13)ntype=1 if(num.eq.1)go to 20 c c reading the model read(lin,100)nz write(lou,100)nz do 1 i=1,nz read(lin,101)vp(i),vs(i),rho(i),qp(i),qs(i),d(i) write(lou,101)vp(i),vs(i),rho(i),qp(i),qs(i),d(i) 1 continue read(lin,101)fref write(lou,101)fref rho2=rho(nz) vp2=vp(nz) vs2=vs(nz) d(1)=0. d(2)=0. return c 20 cvr=0. if(nq.eq.1)cvr=alog(freq/fref)/pi if(nh.eq.1)go to 21 vp(nz)=vp2 vs(nz)=vs2 rho(nz)=rho2 go to 24 21 vp(nz)=0.001 vs(nz)=0.001 rho(nz)=0. 24 continue c c computation of the tangential component of the slowness vector c of the incident wave if(inc.eq.1)go to 22 v1=vp(1) q1=qp(1) urr=1./v1 if(nq.eq.0)go to 23 v1=vp(1)*(1.+cvr/qp(1)) go to 23 22 v1=vs(1) q1=qs(1) urr=1./v1 if(nq.eq.0)go to 23 v1=vs(1)*(1.+cvr/qs(1)) 23 vv1=v1*v1 par=1/v1 ur=par*sin(ang) urr=urr*sin(ang) if(nq.eq.0)go to 26 gaux=1./(q1*q1*cos(gam)*cos(gam)) gg=sqrt(gaux+1.) g1=sqrt(gg+1.) g2=sqrt(gg-1.) u=(g1*sin(ang)-cunit*g2*sin(ang-gam))/sqrt(2.*vv1*(1.+ *1./(q1*q1))) c c computation of R/T coefficients 26 continue call coef8(urr,vp(1),vs(1),rho(1),vp(nz),vs(nz),rho(nz),nc, 1rmod0,rph0) rph0=rph0*180./pi 31 if(ntype.eq.2)go to 33 if(nq.eq.1)go to 32 call rtmat2(nz,vp,vs,rho,d,ur,freq,rp,tp,rs,ts,ntype,nc,rcoef) go to 40 32 call rtmq2(fref,nz,vp,vs,rho,qp,qs,d,u,freq,rp,tp,rs,ts,ntype, 1nc,rcoef,v1,q1,ang,gam) go to 40 33 if(nq.eq.1)go to 34 call rtmat(nz,vp,vs,rho,d,ur,freq,rpp,rps,rss,rsp,tpp,tps,tss, 1tsp,nc,rcoef) go to 40 34 call rtmatq(fref,nz,vp,vs,rho,qp,qs,d,u,freq,rpp,rps,rss,rsp,tpp, 1tps,tss,tsp,nc,rcoef,v1,q1,ang,gam) c 40 b=real(rcoef) c=-aimag(rcoef) rmod=sqrt(b*b+c*c) if(rmod.lt.0.00001)go to 41 rphase=atan2(c,b)*180./pi go to 42 41 rphase=999. c 100 format(8i5) 101 format(8f10.3) 42 return end C C********************************************************************* c SUBROUTINE RTMAT2(N,A,B,RHO,D,U,FREQ,RP,TP,RS,TS,NTYPE,NC, 1 RCOEF) C C R/T COEFFICIENTS AT A STACK OF LAYERS FOR PRESSURE WAVES IN C LIQUIDS AND FOR SH WAVES IN NON-DISSIPATIVE MEDIA C C NTYPE=0....PRESSURE WAVES, LIQUID C NTYPE=1....SH WAVES C C DIMENSION A(N),B(N),RHO(N),D(N) COMPLEX RPD,RPU,TPD,TPU,RP,TP,RSD,RSU,TSD,TSU,RS,TS, 1 A1,A2,B1,B2,AA,EP,ES,DD,FF,RCOEF C C UQ=U*U W=6.283185*FREQ M=N-1 X=-W*D(M) IF(NTYPE.EQ.1)GO TO 5 C C PRESSURE WAVES, LIQUID C AQ1=1./(A(M)*A(M)) RHO1=RHO(M) AQ2=1./(A(N)*A(N)) RHO2=RHO(N) V=SQRT(ABS(AQ1-UQ)) ARG=X*V A1=CMPLX(V,0.) if(uq.gt.aq1)a1=cmplx(0.,-v) IF(UQ.LE.AQ1)EP=CMPLX(COS(ARG),SIN(ARG)) IF(UQ.GT.AQ1)EP=CMPLX(EXP(ARG),0.) V=SQRT(ABS(AQ2-UQ)) A2=CMPLX(V,0.) IF(UQ.GT.AQ2)A2=CMPLX(0.,-V) AA=RHO2*A1+RHO1*A2 RPD=(RHO2*A1-RHO1*A2)/AA RPU=-RPD TPD=2.*RHO2*A1/AA TPU=2.*RHO1*A2/AA RP=RPD*EP*EP TP=TPD*EP GO TO 10 C C SH WAVES C 5 BQ1=1./(B(M)*B(M)) Z1=RHO(M)*B(M)*B(M) BQ2=1./(B(N)*B(N)) Z2=RHO(N)*B(N)*B(N) V=SQRT(ABS(BQ1-UQ)) ARG=X*V B1=CMPLX(V,0.) IF(UQ.LE.BQ1)ES=CMPLX(COS(ARG),SIN(ARG)) IF(UQ.GT.BQ1)B1=CMPLX(0.,-V) IF(UQ.GT.BQ1)ES=CMPLX(EXP(ARG),0.) V=SQRT(ABS(BQ2-UQ)) B2=CMPLX(V,0.) IF(UQ.GT.BQ2)B2=CMPLX(0.,-V) AA=Z1*B1+Z2*B2 RSD=(Z1*B1-Z2*B2)/AA RSU=-RSD TSD=2.*Z1*B1/AA TSU=2.*Z2*B2/AA RS=RSD*ES*ES TS=TSD*ES C C MATRIX MULTIPLICATION FOR LAYERED MEDIUM C 10 IF(N.EQ.2)GO TO 1001 II=N-2 DO 1000 I=II,1,-1 M=I+1 X=-W*D(I) IF(NTYPE.EQ.1)GO TO 15 AQ1=1./(A(I)*A(I)) RHO1=RHO(I) AQ2=1./(A(M)*A(M)) RHO2=RHO(M) V=SQRT(ABS(AQ1-UQ)) ARG=X*V IF(UQ.LE.AQ1)EP=CMPLX(COS(ARG),SIN(ARG)) IF(UQ.GT.AQ1)EP=CMPLX(EXP(ARG),0.) A1=CMPLX(V,0.) IF(UQ.GT.AQ1)A1=CMPLX(0.,-V) V=SQRT(ABS(AQ2-UQ)) A2=CMPLX(V,0.) IF(UQ.GT.AQ2)A2=CMPLX(0.,-V) AA=RHO2*A1+RHO1*A2 RPD=(RHO2*A1-RHO1*A2)/AA RPU=-RPD TPD=2.*RHO2*A1/AA TPU=2.*RHO1*A2/AA DD=1.-RPU*RP RP=(RPD+TPU*TPD*RP/DD)*EP*EP IF(I.EQ.1)FF=TPD/DD IF(I.GT.1)FF=TPD*EP/DD TP=FF*TP GO TO 1000 15 BQ1=1./(B(I)*B(I)) Z1=RHO(I)*B(I)*B(I) BQ2=1./(B(M)*B(M)) Z2=RHO(M)*B(M)*B(M) V=SQRT(ABS(BQ1-UQ)) ARG=X*V IF(UQ.LE.BQ1)ES=CMPLX(COS(ARG),SIN(ARG)) IF(UQ.GT.BQ1)ES=CMPLX(EXP(ARG),0.) B1=CMPLX(V,0.) IF(UQ.GT.BQ1)B1=CMPLX(0.,-V) V=SQRT(ABS(BQ2-UQ)) B2=CMPLX(V,0.) IF(UQ.GT.BQ2)B2=CMPLX(0.,-V) AA=Z1*B1+Z2*B2 RSD=(Z1*B1-Z2*B2)/AA RSU=-RSD TSD=2.*Z1*B1/AA TSU=2.*Z2*B2/AA DD=1.-RSU*RS RS=(RSD+TSU*TSD*RS/DD)*ES*ES IF(I.EQ.1)FF=TSD/DD IF(I.GT.1)FF=TSD*ES/DD TS=FF*TS C 1000 CONTINUE 1001 CONTINUE if(nc.eq.9)rcoef=rs if(nc.eq.10)rcoef=ts if(nc.eq.13)rcoef=rp if(nc.eq.14)rcoef=tp C RETURN END c c******************************************************************** c SUBROUTINE rtmq2(FR,N,A,B,RHO,QA,QB,D,U,FREQ,RP,TP,RS,TS,NTYPE, 1 nc,rcoef,v1,q1,aincr,gamma) C C R/T COEFFICIENTS AT A TRANSITION LAYER FOR PRESSURE WAVES C IN LIQUIDS AND FOR SH WAVES IN DISSIPATIVE MEDIA C C NTYPE=0....PRESSURE WAVES, LIQUID C NTYPE=1....SH WAVES C DIMENSION A(N),B(N),RHO(N),D(N),QA(N),QB(N) COMPLEX RPD,RPU,TPD,TPU,RP,TP,RSD,RSU,TSD,TSU,RS,TS,EIN, 1 A1,A2,B1,B2,AA,EP,ES,Z1,Z2,CW,CV,aq1,aq2,bq1,bq2,rcoef, 2 DD,FF,U,UQ C UQ=U*U WR=6.283185*FR EIN=CMPLX(0.,1.) W=6.283185*FREQ CW=W CV=ALOG(W/WR)/3.141593+0.5*EIN M=N-1 IF(NTYPE.EQ.1)GO TO 5 C C PRESSURE WAVES, DISSIPATIVE LIQUIDS C RHO1=RHO(M) AQ1=A(M)*(1.+CV/QA(M)) AQ1=1./(AQ1*AQ1) A1=CSQRT(AQ1-UQ) RHO2=RHO(N) AQ2=A(N)*(1.+CV/QA(N)) AQ2=1./(AQ2*AQ2) A2=CSQRT(AQ2-UQ) c if(n.gt.2)go to 1 cvr=alog(w/wr)/3.14159 va1=a(m)*(1.+cvr/qa(m)) qa1=qa(m) va2=a(n)*(1.+cvr/qa(n)) qa2=qa(n) call critan(v1,va1,q1,qa1,gamma,ncrit,acrit1,acrit2) if(ncrit.ne.0)then if(aincr.gt.acrit1.and.aincr.lt.acrit2)a1=-a1 endif call critan(v1,va2,q1,qa2,gamma,ncrit,acrit1,acrit2) if(ncrit.ne.0)then if(aincr.gt.acrit1.and.aincr.lt.acrit2)a2=-a2 endif 1 continue c AA=RHO2*A1+RHO1*A2 RPD=(RHO2*A1-RHO1*A2)/AA RPU=-RPD TPD=2.*RHO2*A1/AA TPU=2.*RHO1*A2/AA EP=CEXP(-D(M)*CW*EIN*A1) RP=RPD*EP*EP TP=TPD*EP GO TO 10 C C SH WAVES, DISSIPATIVE MEDIA C 5 continue BQ1=B(M)*(1.+CV/QB(M)) BQ1=1./(BQ1*BQ1) B1=CSQRT(BQ1-UQ) Z1=RHO(M)/BQ1 BQ2=B(N)*(1.+CV/QB(N)) BQ2=1./(BQ2*BQ2) B2=CSQRT(BQ2-UQ) Z2=RHO(N)/BQ2 c if (n.gt.2)go to 2 cvr=alog(w/wr)/3.141593 vb1=b(m)*(1.+cvr/qb(m)) qb1=qb(m) vb2=b(n)*(1.+cvr/qb(n)) qb2=qb(n) call critan(v1,vb1,q1,qb1,gamma,ncrit,acrit1,acrit2) if(ncrit.ne.0)then if(aincr.gt.acrit1.and.aincr.lt.acrit2)b1=-b1 endif call critan(v1,vb2,q1,qb2,gamma,ncrit,acrit1,acrit2) if(ncrit.ne.0)then if(aincr.gt.acrit1.and.aincr.lt.acrit2)b2=-b2 endif 2 continue c AA=Z1*B1+Z2*B2 RSD=(Z1*B1-Z2*B2)/AA RSU=-RSD TSD=2.*Z1*B1/AA TSU=2.*Z2*B2/AA ES=CEXP(-D(M)*CW*EIN*B1) RS=RSD*ES*ES TS=TSD*ES C C MATRIX MULTIPLICATION FOR LAYERED MEDIUM C 10 IF(N.EQ.2)GO TO 1001 II=N-2 DO 1000 I=II,1,-1 M=I+1 IF(NTYPE.EQ.1)GO TO 15 AQ1=A(I)*(1.+CV/QA(I)) AQ1=1./(AQ1*AQ1) RHO1=RHO(I) AQ2=A(M)*(1.+CV/QA(M)) AQ2=1./(AQ2*AQ2) RHO2=RHO(M) A1=CSQRT(AQ1-UQ) A2=CSQRT(AQ2-UQ) AA=RHO2*A1+RHO1*A2 RPD=(RHO2*A1-RHO1*A2)/AA RPU=-RPD TPD=2.*RHO2*A1/AA TPU=2.*RHO1*A2/AA EP=CEXP(-D(I)*CW*EIN*A1) DD=1.-RPU*RP RP=(RPD+TPU*TPD*RP/DD)*EP*EP IF(I.EQ.1)FF=TPD/DD IF(I.GT.1)FF=TPD*EP/DD TP=FF*TP GO TO 1000 15 BQ1=B(I)*(1.+CV/QB(I)) BQ1=1./(BQ1*BQ1) Z1=RHO(I)/BQ1 BQ2=B(M)*(1.+CV/QB(M)) BQ2=1./(BQ2*BQ2) Z2=RHO(M)/BQ2 B1=CSQRT(BQ1-UQ) B2=CSQRT(BQ2-UQ) AA=Z1*B1+Z2*B2 RSD=(Z1*B1-Z2*B2)/AA RSU=-RSD TSD=2.*Z1*B1/AA TSU=2.*Z2*B2/AA ES=CEXP(-D(I)*CW*EIN*B1) DD=1.-RSU*RS RS=(RSD+TSU*TSD*RS/DD) IF(I.EQ.1)FF=TSD/DD IF(I.GT.1)FF=TSD*ES/DD TS=FF*TS C 1000 CONTINUE 1001 CONTINUE if(nc.eq.9)rcoef=rs if(nc.eq.10)rcoef=ts if(nc.eq.13)rcoef=rp if(nc.eq.14)rcoef=tp C RETURN END c c********************************************************************* c SUBROUTINE RTMAT(N,A,B,RHO,D,U,FRQ,RPP,RPS,RSS,RSP, * TPP,TPS,TSS,TSP,NC,RCOEF) C C P/SV REFLECTION AND TRANSMISSION COEFFICIENTS FOR A C NON-DISSIPATIVE TRANSITION LAYER, USING RECURSIVE FORMALISM C C N = NUMBER OF LAYERS (HALFSPACES INCLUDED) C A,B,RHO,D = MODEL C U = TANGENTIAL SLOWNESS C FRQ = FREQUENCY C RCOEF = COMPLEX-VALUED R/T COEFFICIENT OF THE TYPE NC C DIMENSION A(1),B(1),RHO(1),D(1) COMPLEX RPP,RPS,RSS,RSP,TPP,TPS,TSS,TSP,MU11,MU12,MU21,MU22, * RPPD,RSPD,RPSD,RSSD,TPPD,TSPD,TPSD,TSSD,RPPU,RSPU,RPSU,RSSU, * TPPU,TSPU,TPSU,TSSU,F11,F12,F21,F22,G11,G12,G21,G22,H11,H12, * H21,H22,I11,I12,I21,I22,EP,ES,EX,rcoef M=N-1 AQ1=1./(A(M)*A(M)) BQ1=1./(B(M)*B(M)) RHO1=RHO(M) AQ2=1./(A(N)*A(N)) BQ2=1./(B(N)*B(N)) RHO2=RHO(N) C=2.*(RHO1/BQ1-RHO2/BQ2) CALL RT(AQ1,BQ1,RHO1,AQ2,BQ2,RHO2,C,U,MU11,MU12,MU21,MU22, * TPPD,TSPD,TPSD,TSSD,RPPU,RSPU,RPSU,RSSU,TPPU,TSPU,TPSU,TSSU) W=6.283185*FRQ UQ=U*U CALL PHASE(UQ,W,AQ1,BQ1,D(M),EP,ES) EX=EP*ES RPP=MU11*EP*EP RSP=MU12*EX RPS=MU21*EX RSS=MU22*ES*ES TPP=TPPD*EP TSP=TSPD*ES TPS=TPSD*EP TSS=TSSD*ES IF(N.EQ.2) go to 1001 II=N-2 DO 1000 I=II,1,-1 M=I+1 AQ1=1./(A(I)*A(I)) BQ1=1./(B(I)*B(I)) RHO1=RHO(I) AQ2=1./(A(M)*A(M)) BQ2=1./(B(M)*B(M)) RHO2=RHO(M) C=2.*(RHO1/BQ1-RHO2/BQ2) CALL RT(AQ1,BQ1,RHO1,AQ2,BQ2,RHO2,C,U,RPPD,RSPD,RPSD,RSSD, * TPPD,TSPD,TPSD,TSSD,RPPU,RSPU,RPSU,RSSU,TPPU,TSPU,TPSU,TSSU) CALL MULMAT(RPP,RSP,RPS,RSS,TPPD,TSPD,TPSD,TSSD,G11,G12,G21, * G22) CALL MULMAT(RPP,RSP,RPS,RSS,RPPU,RSPU,RPSU,RSSU,H11,H12,H21, * H22) H11=1.-H11 H12=-H12 H21=-H21 H22=1.-H22 CALL MATINV(H11,H12,H21,H22,I11,I12,I21,I22) CALL MULMAT(I11,I12,I21,I22,G11,G12,G21,G22,H11,H12,H21,H22) CALL MULMAT(TPPU,TSPU,TPSU,TSSU,H11,H12,H21,H22,G11,G12,G21, * G22) MU11=RPPD+G11 MU12=RSPD+G12 MU21=RPSD+G21 MU22=RSSD+G22 CALL MULMAT(RPPU,RSPU,RPSU,RSSU,RPP,RSP,RPS,RSS,G11,G12,G21, * G22) G11=1.-G11 G12=-G12 G21=-G21 G22=1.-G22 CALL MATINV(G11,G12,G21,G22,H11,H12,H21,H22) CALL MULMAT(H11,H12,H21,H22,TPPD,TSPD,TPSD,TSSD,G11,G12,G21, * G22) CALL PHASE(UQ,W,AQ1,BQ1,D(I),EP,ES) F11=G11*EP F12=G12*ES F21=G21*EP F22=G22*ES CALL MULMAT(TPP,TSP,TPS,TSS,F11,F12,F21,F22,G11,G12,G21,G22) TPP=G11 TSP=G12 TPS=G21 TSS=G22 EX=EP*ES RPP=MU11*EP*EP RSP=MU12*EX RPS=MU21*EX 1000 RSS=MU22*ES*ES 1001 CONTINUE if(nc.eq.1)rcoef=rpp if(nc.eq.2)rcoef=-a(1)*rps/b(1) if(nc.eq.3)rcoef=a(1)*tpp/a(n) if(nc.eq.4)rcoef=-a(1)*tps/b(n) if(nc.eq.5)rcoef=-b(1)*rsp/a(1) if(nc.eq.6)rcoef=rss if(nc.eq.7)rcoef=-b(1)*tsp/a(n) if(nc.eq.8)rcoef=b(1)*tss/b(n) return END C C****************************************************************** C SUBROUTINE RT(AQ1,BQ1,RHO1,AQ2,BQ2,RHO2,C,U,RPPD,RSPD,RPSD, * RSSD,TPPD,TSPD,TPSD,TSSD,RPPU,RSPU,RPSU,RSSU,TPPU,TSPU, * TPSU,TSSU) C C A ROUTINE TO RTMAT. C MATRICES OF R/T COEFFICIENTS AT A SINGLE INTERFACE BETWEEN C TWO HALFSPACES C COMPLEX RPPD,RSPD,RPSD,RSSD,TPPD,TSPD,TPSD,TSSD,RPPU,RSPU, * RPSU,RSSU,TPPU,TSPU,TPSU,TSSU,A1,B1,A2,B2,D1D,D2D,D1U,D2U, * D,A1B1,A2B2,A1B2,A2B1,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10 UQ=U*U V=SQRT(ABS(AQ1-UQ)) A1=CMPLX(V,0.) IF(UQ.GT.AQ1) A1=CMPLX(0.,-V) V=SQRT(ABS(BQ1-UQ)) B1=CMPLX(V,0.) IF(UQ.GT.BQ1) B1=CMPLX(0.,-V) V=SQRT(ABS(AQ2-UQ)) A2=CMPLX(V,0.) IF(UQ.GT.AQ2) A2=CMPLX(0.,-V) V=SQRT(ABS(BQ2-UQ)) B2=CMPLX(V,0.) IF(UQ.GT.BQ2) B2=CMPLX(0.,-V) C0=C*UQ C1=C0-RHO1 C2=C0+RHO2 C3=C1+RHO2 A1B1=A1*B1 A2B2=A2*B2 A1B2=A1*B2 A2B1=A2*B1 RHO12=RHO1*RHO2 T1=C1*C1*A2B2+RHO12*A2B1 T2=C2*C2*A1B1+RHO12*A1B2 T3=C3*C3*UQ T4=C*C0*A1B1*A2B2 D1D=T3+T1 D2D=T4+T2 D=D1D+D2D D1U=T3+T2 D2U=T4+T1 T5=2./D T6=A1*T5 T7=B1*T5 T8=A2*T5 T9=B2*T5 RPPD=(D2D-D1D)/D RPPU=(D2U-D1U)/D T10=U*(C3*C2+C*C1*A2B2) RPSD=-T6*T10 RSPD=T7*T10 T10=RHO12*(A1B2-A2B1)*T5 RSSD=RPPD-T10 RSSU=RPPU+T10 T10=U*(C3*C1+C*C2*A1B1) RPSU=T8*T10 RSPU=-T9*T10 T10=C2*B1-C1*B2 TPPD=RHO1*T6*T10 TPPU=RHO2*T8*T10 T10=U*(C3+C*A2B1) TPSD=-RHO1*T6*T10 TSPU=RHO2*T9*T10 T10=C2*A1-C1*A2 TSSD=RHO1*T7*T10 TSSU=RHO2*T9*T10 T10=U*(C3+C*A1B2) TSPD=RHO1*T7*T10 TPSU=-RHO2*T8*T10 RETURN END C C****************************************************************** C SUBROUTINE MULMAT(A11,A12,A21,A22,B11,B12,B21,B22,C11,C12, * C21,C22) C C A ROUTINE TO RTMAT AND RTMATQ. C MATRIX MULTIPLICATION C COMPLEX A11,A12,A21,A22,B11,B12,B21,B22,C11,C12,C21,C22 C11=A11*B11+A12*B21 C12=A11*B12+A12*B22 C21=A21*B11+A22*B21 C22=A21*B12+A22*B22 RETURN END C C****************************************************************** C SUBROUTINE MATINV(A11,A12,A21,A22,B11,B12,B21,B22) C C A ROUTINE TO RTMAT AND RTMATQ. C MATRIX INVERSION C COMPLEX A11,A12,A21,A22,B11,B12,B21,B22,D D=1./(A11*A22-A12*A21) B11=A22*D B12=-A12*D B21=-A21*D B22=A11*D RETURN END C C****************************************************************** C SUBROUTINE PHASE(UQ,W,AQ,BQ,D,EP,ES) C C A ROUTINE TO RTMAT AND RTMATQ C COMPLEX EP,ES X=-W*D V=SQRT(ABS(AQ-UQ)) A=X*V IF(UQ.GT.AQ) GO TO 100 EP=CMPLX(COS(A),SIN(A)) GO TO 200 100 EP=CMPLX(EXP(A),0.) 200 V=SQRT(ABS(BQ-UQ)) A=X*V IF(UQ.GT.BQ) GO TO 300 ES=CMPLX(COS(A),SIN(A)) GO TO 400 300 ES=CMPLX(EXP(A),0.) 400 RETURN END C C *********************************************************** C SUBROUTINE COEF8(P,VP1,VS1,RO1,VP2,VS2,RO2,NCODE,RMOD,RPH) C C ROUTINE COEF8 IS DESIGNED FOR THE COMPUTATION OF REFLECTION C AND TRANSMISSION COEFFICIENTS OF P, SV AND SH PLANE WAVES C AT A SINGLE PLANE INTERFACE BETWEEN TWO HOMOGENEOUS C NON-DISSIPATIVE SOLID HALFSPACES OR AT A FREE SURFACE OF A C NON-DISSIPATIVE HOMOGENEOUS SOLID HALFSPACE. EXPLICIT EQUATIONS C ARE USED. ALSO PRESSURE WAVES IN LIQUIDS ARE CONSIDERED. 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) C AUX=999.*3.14159/180. IF(NCODE.GE.9)GO TO 300 IF(RO2.LT.0.0001)GO TO 150 PRMT(1)=VP1 PRMT(2)=VS1 PRMT(3)=VP2 PRMT(4)=VS2 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. 11.)GO 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=G1 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) 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) GO TO 26 5 RR=-VS1*B(2)*HC*(Q*Y*C2+A2*X*Z) 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) 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) 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) GO TO 250 105 R=-VS1*D(2)*SC*(Q*Y*E2+A2*X*Z) 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) 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 CONVERSION COEFFICIENTS 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,166,165,163,164,168,167),NCODE 161 RR=(-A9+H2)*H GO TO 26 162 RR=-A5*B(1)*H*A6 GO TO 26 163 RR=A5*B(2)*H*A6*VS1/VP1 GO TO 26 164 RR=-(A9-H2)*H GO TO 26 165 RR=A5*H1*H 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 C 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) IF(RMOD.LT.0.00001)RPH=AUX RETURN 157 RMOD=0. RPH=AUX RETURN C C EARTHS SURFACE,REAL COEFFICIENTS AND CONVERSION COEFFICIENTS 200 S1=D(1)*D(2) S2=A8*S1 S=1./(A9+S2) GO TO (201,202,206,205,203,204,208,207),NCODE 201 R=(-A9+S2)*S GO TO 250 202 R=-A5*D(1)*S*A6 GO TO 250 203 R=A5*D(2)*S*A6*VS1/VP1 GO TO 250 204 R=(S2-A9)*S GO TO 250 205 R=A5*S1*S 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 250 IF(R.LT.0.)GO TO 251 RMOD=R RPH=0. IF(RMOD.LT.0.00001)RPH=AUX RR=CMPLX(R,0.) RETURN 251 RMOD=-R RPH=-3.14159 IF(RMOD.LT.0.00001)RPH=AUX RR=CMPLX(R,0.) RETURN C C SH COEFFICIENTS OF REFLECTION, TRANSMISSION AND CONVERSION C 300 IF(NCODE.GE.13)GO TO 400 IF(RO2.LT.0.0001) GO TO 302 A1=P*VS1 A2=P*VS2 A3=RO1*VS1 A4=RO2*VS2 A5=SQRT(ABS(1.-A1*A1)) A6=SQRT(ABS(1.-A2*A2)) C1=A5 IF(A2.LE.1.)C2=A6 IF(A2.GT.1.)C2=CMPLX(0.,A6) C1=C1*A3 C2=C2*A4 H=1./(C1+C2) IF(NCODE.EQ.10)GO TO 301 RR=H*(C1-C2) GO TO 26 301 RR=2.*C1*H GO TO 26 302 RMOD=1. RPH=0. IF(NCODE.EQ.10)GO TO 303 RETURN 303 RMOD=2. RPH=0. RETURN C C PRESSURE R/T COEFFICIENTS, LIQUIDS C 400 IF(RO2.LT.0.0001) GO TO 402 A1=P*VP1 A2=P*VP2 A3=RO1*VP1 A4=RO2*VP2 A5=SQRT(ABS(1.-A1*A1)) A6=SQRT(ABS(1.-A2*A2)) C1=A5 IF(A2.LT.1.)C2=A6 IF(A2.GE.1.)C2=CMPLX(0.,A6) C1=A4*C1 C2=A3*C2 H=1./(C1+C2) IF(NCODE.EQ.14)GO TO 401 RR=H*(C1-C2) GO TO 26 401 RR=2.*C1*H GO TO 26 402 RMOD=1. RPH=3.141593 IF(NCODE.EQ.14)GO TO 403 RETURN 403 RMOD=0. RPH=AUX RETURN END C C ************************************************************* C SUBROUTINE RTMATQ(FR,N,A,B,RHO,QA,QB,D,U,FRQ, * RPP,RPS,RSS,RSP,TPP,TPS,TSS,TSP,NC,RCOEF,v1,q1,aincr,gamma) C C P/SV REFLECTION AND TRANSMISSION COEFFICIENTS FOR A DISSIPATIVE C TRANSITION LAYER, USING RECURSIVE ALGORITHM. C C FR = REFERENCE FREQUENCY C N = NUMBER OF LAYERS (HALFSPACES INCLUDED) C A,B,RHO,QA,QB,D = MODEL (QA AND QB INDEPENDENT OF FREQUENCY) C U = TANGENTIAL SLOWNESS (COMPLEX-VALUED) C FRQ = FREQUENCY C RCOEF = COMPLEX-VALUED R/T COEFFICIENT OF THE TYPE NC C C LAYERS 1 AND N ARE HALFSPACES C C DIMENSION A(1),B(1),RHO(1),D(1),QA(1),QB(1),nc(1) COMPLEX RPP,RPS,RSS,RSP,TPP,TPS,TSS,TSP,MU11,MU12,MU21,MU22, * RPPD,RSPD,RPSD,RSSD,TPPD,TSPD,TPSD,TSSD,RPPU,RSPU,RPSU,RSSU, * TPPU,TSPU,TPSU,TSSU,F11,F12,F21,F22,G11,G12,G21,G22,H11,H12, * H21,H22,I11,I12,I21,I22,EP,ES,EX,AQ1,BQ1,AQ2,BQ2,C,CW,CV, * EIN,rcoef,apom,bpom,apom2,bpom2,u,uq C WR=6.28319*FR EIN=CMPLX(0.,1.) W=6.28319*FRQ CW=W CV=ALOG(W/WR)/3.141593+0.5*EIN C M=N-1 RHO1=RHO(M) AQ1=A(M)*(1.+CV/QA(M)) APOM=AQ1 AQ1=1./(AQ1*AQ1) BQ1=B(M)*(1.+CV/QB(M)) BPOM=BQ1 BQ1=1./(BQ1*BQ1) RHO2=RHO(N) AQ2=A(N)*(1.+CV/QA(N)) APOM2=AQ2 AQ2=1./(AQ2*AQ2) BQ2=B(N)*(1.+CV/QB(N)) BPOM2=BQ2 BQ2=1./(BQ2*BQ2) c c determination of possible critical angles in critan (for n=2) if(n.gt.2) go to 1 cvr=alog(w/wr)/3.141593 va1=a(m)*(1.+cvr/qa(m)) vb1=b(m)*(1.+cvr/qb(m)) qa1=qa(m) qb1=qb(m) va2=a(n)*(1.+cvr/qa(n)) vb2=b(n)*(1.+cvr/qb(n)) qa2=qa(n) qb2=qb(n) 1 continue c C=2.*(RHO1/BQ1-RHO2/BQ2) CALL RTQ(AQ1,BQ1,RHO1,AQ2,BQ2,RHO2,C,U,MU11,MU12,MU21,MU22, * TPPD,TSPD,TPSD,TSSD,RPPU,RSPU,RPSU,RSSU,TPPU,TSPU,TPSU,TSSU, * aincr,gamma,v1,q1,va1,vb1,va2,vb2,qa1,qb1,qa2,qb2,n) UQ=U*U EP=CEXP(-D(M)*CW*EIN*CSQRT(AQ1-UQ)) ES=CEXP(-D(M)*CW*EIN*CSQRT(BQ1-UQ)) EX=EP*ES RPP=MU11*EP*EP RSP=MU12*EX RPS=MU21*EX RSS=MU22*ES*ES TPP=TPPD*EP TSP=TSPD*ES TPS=TPSD*EP TSS=TSSD*ES IF(N.EQ.2) go to 1001 C C LOOP FOR LAYERS C II=N-2 DO 1000 I=II,1,-1 M=I+1 AQ1=A(I)*(1.+CV/QA(I)) apom=aq1 AQ1=1./(AQ1*AQ1) BQ1=B(I)*(1.+CV/QB(I)) bpom=bq1 BQ1=1./(BQ1*BQ1) RHO1=RHO(I) AQ2=A(M)*(1.+CV/QA(M)) AQ2=1./(AQ2*AQ2) BQ2=B(M)*(1.+CV/QB(M)) BQ2=1./(BQ2*BQ2) RHO2=RHO(M) C=2.*(RHO1/BQ1-RHO2/BQ2) CALL RTQ(AQ1,BQ1,RHO1,AQ2,BQ2,RHO2,C,U,RPPD,RSPD,RPSD,RSSD, * TPPD,TSPD,TPSD,TSSD,RPPU,RSPU,RPSU,RSSU,TPPU,TSPU,TPSU,TSSU, * aincr,gamma,v1,q1,va1,vb1,va2,vb2,qa1,qb1,qa2,qb2,n) CALL MULMAT(RPP,RSP,RPS,RSS,TPPD,TSPD,TPSD,TSSD,G11,G12,G21, * G22) CALL MULMAT(RPP,RSP,RPS,RSS,RPPU,RSPU,RPSU,RSSU,H11,H12,H21, * H22) H11=1.-H11 H12=-H12 H21=-H21 H22=1.-H22 CALL MATINV(H11,H12,H21,H22,I11,I12,I21,I22) CALL MULMAT(I11,I12,I21,I22,G11,G12,G21,G22,H11,H12,H21,H22) CALL MULMAT(TPPU,TSPU,TPSU,TSSU,H11,H12,H21,H22,G11,G12,G21, * G22) MU11=RPPD+G11 MU12=RSPD+G12 MU21=RPSD+G21 MU22=RSSD+G22 CALL MULMAT(RPPU,RSPU,RPSU,RSSU,RPP,RSP,RPS,RSS,G11,G12,G21, * G22) G11=1.-G11 G12=-G12 G21=-G21 G22=1.-G22 CALL MATINV(G11,G12,G21,G22,H11,H12,H21,H22) CALL MULMAT(H11,H12,H21,H22,TPPD,TSPD,TPSD,TSSD,G11,G12,G21, * G22) EP=CEXP(-D(I)*CW*EIN*CSQRT(AQ1-UQ)) ES=CEXP(-D(I)*CW*EIN*CSQRT(BQ1-UQ)) F11=G11*EP F12=G12*ES F21=G21*EP F22=G22*ES CALL MULMAT(TPP,TSP,TPS,TSS,F11,F12,F21,F22,G11,G12,G21,G22) TPP=G11 TSP=G12 TPS=G21 TSS=G22 EX=EP*ES RPP=MU11*EP*EP RSP=MU12*EX RPS=MU21*EX 1000 RSS=MU22*ES*ES 1001 continue if (nc(1).eq.1)rcoef=rpp if(nc(1).eq.2)rcoef=-apom*rps/bpom if(nc(1).eq.3)rcoef=apom*tpp/apom2 if(nc(1).eq.4)rcoef=-apom*tps/bpom2 if(nc(1).eq.5)rcoef=-bpom*rsp/apom if(nc(1).eq.6)rcoef=rss if(nc(1).eq.7)rcoef=-bpom*tsp/apom2 if(nc(1).eq.8)rcoef=bpom*tss/bpom2 RETURN END C C ************************************************************ C SUBROUTINE RTQ(AQ1,BQ1,RHO1,AQ2,BQ2,RHO2,C,U,RPPD,RSPD,RPSD, * RSSD,TPPD,TSPD,TPSD,TSSD,RPPU,RSPU,RPSU,RSSU,TPPU,TSPU, * TPSU,TSSU,aincr,gamma, * v1,q1,va1,vb1,va2,vb2,qa1,qb1,qa2,qb2,n) C C A ROUTINE TO RTMATQ C R/T COEFFICIENTS AT A SINGLE INTERFACE BETWEEN C TWO HOMOGENEOUS DISSIPATIVE HALFSPACES C COMPLEX RPPD,RSPD,RPSD,RSSD,TPPD,TSPD,TPSD,TSSD,RPPU,RSPU, * RPSU,RSSU,TPPU,TSPU,TPSU,TSSU,A1,B1,A2,B2,D1D,D2D,D1U,D2U, * D,A1B1,A2B2,A1B2,A2B1,T1,T2,T3,T4,T5,T6,T7,T8,T9,T10, * AQ1,BQ1,AQ2,BQ2,C,C1,C2,C3,C0,U,UQ C uq=u*u A1=CSQRT(AQ1-UQ) B1=CSQRT(BQ1-UQ) A2=CSQRT(AQ2-UQ) B2=CSQRT(BQ2-UQ) c c determination of possible critical angles in critan, for n=2 if(n.gt.2)go to 1 call critan(v1,va1,q1,qa1,gamma,ncrit,acrit1,acrit2) if(ncrit.ne.0)then if(aincr.gt.acrit1.and.aincr.lt.acrit2)a1=-a1 endif call critan(v1,vb1,q1,qb1,gamma,ncrit,acrit1,acrit2) if(ncrit.ne.0)then if(aincr.gt.acrit1.and.aincr.lt.acrit2)b1=-b1 endif call critan(v1,va2,q1,qa2,gamma,ncrit,acrit1,acrit2) if(ncrit.ne.0)then if(aincr.gt.acrit1.and.aincr.lt.acrit2)a2=-a2 endif call critan(v1,vb2,q1,qb2,gamma,ncrit,acrit1,acrit2) if(ncrit.ne.0)then if(aincr.gt.acrit1.and.aincr.lt.acrit2)b2=-b2 endif c 1 C0=C*UQ C1=C0-RHO1 C2=C0+RHO2 C3=C1+RHO2 A1B1=A1*B1 A2B2=A2*B2 A1B2=A1*B2 A2B1=A2*B1 RHO12=RHO1*RHO2 T1=C1*C1*A2B2+RHO12*A2B1 T2=C2*C2*A1B1+RHO12*A1B2 T3=C3*C3*UQ T4=C*C0*A1B1*A2B2 D1D=T3+T1 D2D=T4+T2 D=D1D+D2D D1U=T3+T2 D2U=T4+T1 T5=2./D T6=A1*T5 T7=B1*T5 T8=A2*T5 T9=B2*T5 RPPD=(D2D-D1D)/D RPPU=(D2U-D1U)/D T10=U*(C3*C2+C*C1*A2B2) RPSD=-T6*T10 RSPD=T7*T10 T10=RHO12*(A1B2-A2B1)*T5 RSSD=RPPD-T10 RSSU=RPPU+T10 T10=U*(C3*C1+C*C2*A1B1) RPSU=T8*T10 RSPU=-T9*T10 T10=C2*B1-C1*B2 TPPD=RHO1*T6*T10 TPPU=RHO2*T8*T10 T10=U*(C3+C*A2B1) TPSD=-RHO1*T6*T10 TSPU=RHO2*T9*T10 T10=C2*A1-C1*A2 TSSD=RHO1*T7*T10 TSSU=RHO2*T9*T10 T10=U*(C3+C*A1B2) TSPD=RHO1*T7*T10 TPSU=-RHO2*T8*T10 RETURN END C C ************************************************************* C COMPLEX FUNCTION CSQ(X) C C A ROUTINE TO RTMATQ. C COMPLEX X A=CABS(X) B=REAL(X) RE=0.7071*SQRT(A+B) AI=-0.7071*SQRT(A-B) CSQ=CMPLX(RE,AI) RETURN END c c********************************************************************** c subroutine critan(v1,v2,q1,q2,gamma,n,acrit1,acrit2) c c A routine to compute the discrete critical angles c c auxw1=sqrt(1.+1./(q1*q1)) c w1=2./(v1*v1*(1.+auxw1)) c auxw2=sqrt(1.+1./(q2*q2)) c w2=2./(v2*v2*(1.+auxw2)) auxw1=1.+1./(q1*q1) auxw2=1.+1./(q2*q2) w1=1./(v1*v1*auxw1) w2=1./(v2*v2*auxw2) c=(2.*w2*q1)/(w1*q2)-1. a=(w2*q1*cos(gamma))/(w1*q2) c n = 0 acrit1=1000. acrit2=1000. c caux=1./(c*c) cos2g=cos(gamma)*cos(gamma) if(cos2g.le.caux)then aux=sqrt(1.-c*c*cos2g) omega1=(1.+c*cos2g+sin(gamma)*aux)/2. omega2=(1.+c*cos2g-sin(gamma)*aux)/2. if(0.lt.omega1.and.omega1.lt.1)then acrit1=asin(sqrt(omega1)) if(acrit1.lt.gamma)then acrit1=1000. go to 2 endif if(abs(acrit1-gamma).lt.0.000001)then acrit1=1000. go to 2 endif aim1 = a - sin(acrit1)*sin(acrit1-gamma) if(abs(aim1).gt.0.0001)then acrit1 = 1000. go to 2 endif c n=1 g=sqrt(1.+1./(q1*q1*cos2g)) rea1=w2-(w1*((1.+g)*sin(acrit1)*sin(acrit1)- * (g-1.)*sin(acrit1-gamma)*sin(acrit1-gamma)))/2. if(rea1.gt.0)then acrit1=1000. c n=0 else n=1 endif endif c 2 if(0.lt.omega2.and.omega2.lt.1)then acrit2=asin(sqrt(omega2)) if(acrit2.lt.gamma)then acrit2=1000. go to 3 endif if(abs(acrit2-gamma).lt.0.000001)then acrit2=1000. go to 3 endif aim2=a - sin(acrit2)*sin(acrit2-gamma) if(abs(aim2).gt.0.0001)then acrit2=1000. go to 3 endif c n=n+1 g=sqrt(1.+1./(q1*q1*cos2g)) rea2=w2-(w1*((1.+g)*sin(acrit2)*sin(acrit2)- * (g-1.)*sin(acrit2-gamma)*sin(acrit2-gamma)))/2. if(rea2.gt.0)then acrit2=1000. c n=n-1 else n=n+1 endif endif endif c 3 if(acrit2.lt.acrit1)then acr11=acrit2 acrit2=acrit1 acrit1=acr11 endif c if(abs(acrit1-acrit2).lt.0.00005)then n=1 acrit2=1000. endif return end c c*********************************************************************** c