C Subroutine file 'trans.for' to transform the computed quantities at an C interface C C By Vlastislav Cerveny, Ludek Klimes, Ivan Psencik C C This file consists of the following subroutines: C TRANS...Subroutine transforming the computed quantities across a C curved interface, see C.R.T.5.9. C CONV... Subroutine designed to replace the amplitudes Y(28) to C Y(IY(1)) expressed in the ray-centred coordinate system by C the amplitudes involving appropriate conversion C coefficients, see C.R.T.5.5.4. C C======================================================================= C SUBROUTINE TRANS(YL,Y,YY,IY,KODNEW,ICBNEW,IEND) REAL YL(6),Y(35),YY(5) INTEGER IY(12),KODNEW,ICBNEW,IEND C C This subroutine transforms the computed quantities across a curved C interface, see C.R.T.5.9. C C Input: C YL... Array containing local quantities at the point of C incidence. C Y... Array containing basic quantities computed along the ray, C corresponding to the incident wave. C YY... Array containing real auxiliary quantities computed along C the ray, corresponding to the incident wave. C IY... Array containing integer auxiliary quantities computed C along the ray, corresponding to the incident wave. C KODNEW..Position in the code (index in the array KODE) C corresponding to the next element of the ray. Output from C subroutine CODE. C ICBNEW..The index of the complex block in which the next element C of the ray is to be situated, supplemented by the sign '+' C for P wave or '-' for S wave. Output from subroutine CODE. C C Output: C YL... Array containing local quantities at the R/T point, C corresponding to the generated wave. C Y... Array containing basic quantities computed along the ray, C corresponding to the generated wave. C YY... Array containing real auxiliary quantities computed along C the ray, corresponding to the generated wave. Unchanged C input values. C IY... Array containing integer auxiliary quantities computed C along the ray, corresponding to the generated wave. C IEND... Indication of the termination of the ray computation (see C C.R.T.5.4). C IEND.EQ.0... Computation of the ray continues. C IEND.NE.0... The computation of the ray terminates, C different values of IEND may specify the reason of the C termination: C 24... Transmission is required by the code at a free C surface. C 25... Ray of the required reflected or transmitted wave C is not real-valued (overcritical reflection or C transmission). C 26... S wave in a liquid block is required by the code. C 32... Reflection or type conversion at the fictitious C part of the interface. Here, the interface is C considered to be fictitious if P and S wave velocities C and the density are the same at both sides of the C interface at the point of incidence. C 33... Zero reflection/transmission coefficient. C 37... The next element of the ray is required by the code C to be situated in a complex block that does not touch C the point of incidence. C C Subroutines and external functions required: EXTERNAL ISIDE,VELOC,KOOR,METRIC,SRFC2,PARM2,SMVPRD,COEF50,COEFSH INTEGER ISIDE,KOOR C ISIDE,VELOC... File 'model.for' of the package 'MODEL'. C KOOR,METRIC... File 'metric.for' of the package 'MODEL'. C SRFC2 and subsequent routines... File 'srfc.for' and subsequent C files (like 'val.for' and 'fit.for') of the package C 'MODEL'. C PARM2 and subsequent routines... File 'parm.for' and subsequent C files (like 'val.for' and 'fit.for') of the package C 'MODEL'. C SMVPRD... File 'means.for' of the package 'MODEL'. C COEF50,COEFSH... File 'coef.for'. C C Error messages: C 591... Wrong interface: C This error should not appear. Contact the authors. C C Date: 1996, September 30 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations for local model parameters: FAUX(10), C G(12),GAMMA(18),GSQRD, UP(10),US(10),RO,QP,QS, VP,VS,VD(10),QL: INCLUDE 'auxmod.inc' C C....................................................................... C C Auxiliary storage locations: C INTEGER NAMPL1,NAMPL2,ISB1,ICB1,ICB2,I,J,NCODE,ND REAL VP1,VS1,RO1,VP2,VS2,RO2,V1,VD1(4),V2,VD2(10) REAL H11,H13,H43,HHH,Z11,Z12,Z13,Z41,Z42,Z43,ZZZ REAL H21,H23,H53, Z21,Z22,Z23,Z51,Z52,Z53 REAL H31,H33,H63, Z31,Z32,Z33,Z61,Z62,Z63 REAL E11,E21,E31, C1,S1,C2,S2,C,S REAL V11,V21,V31,V12,V22,V32,D11,D12,D22,AUX1,AUX2,AUX3 REAL P,AUX,E12,CFC11,CFC12,CFC22,Q1,Q2,P1,P2 REAL COEF,RMOD,RPH,RMODSH,RPHSH,RR,YR REAL R1A1,Y1A1,R2A1,Y2A1,R1B1,Y1B1,R2B1,Y2B1 REAL R1A2,Y1A2,R2A2,Y2A2,R1B2,Y1B2,R2B2,Y2B2 C C NAMPL1..Number of the amplitude coefficients of the incident wave C (5.9.1a). C NAMPL2..Number of the amplitude coefficients of the R/T wave C (5.9.1b). C ISB1... Index of a simple block of the incident wave (5.9.1e), C (5.9.3). C ICB1... Index of a complex block of the incident wave (5.9.1e). C ICB2... Index of a complex block at the other side (5.9.1e). C I,J... Temporary storage locations (5.9.3), (5.9.5). C NCODE,ND... Type of the R/T coefficient (5.9.6). C VP1,VS1,RO1... Parameters of the complex block IABS(ICB1) (5.9.2). C VP2,VS2,RO2... Parameters of the complex block ICB2 (5.9.2). C V1,VD1(2:4)... Incident wave velocity and its derivatives (5.9.2). C V2,VD2(2:4)...R/T wave velocity and its derivatives (5.9.2). C VD2(1),VD2(5:10)... Temporary storage locations (5.9.2). C H11,H21,H31,H13,H23,H33... Covariant components of the basis C vectors of ray-centred coordinate system of the incident C wave (5.9.3). C H43,H53,H63... Contravariant components of the slowness vector of C the incident wave (5.9.3). Not used in Cartesian C coordinates. C HHH... Norm of the slowness vector (5.9.3). C Z11,Z21,Z31,Z12,Z22,Z32,Z13,Z23,Z33... Covariant components of the C local Cartesian coordinate system (5.9.3). C Z41,Z51,Z61,Z42,Z52,Z62,Z43,Z53,Z63... Contravariant components of C the local Cartesian coordinate system (5.9.3). C ZZZ... Norm of the gradient of the function describing the C interface (5.9.3,5.9.4). C E11,E21,E31... Covariant components of the unit vector C perpendicular to the R/T ray in the plane of incidence C (5.9.3). C E11 is also a temporary storage location in (5.9.5). C E13,E23,E33... Covariant components of the unit vector tangent to C the R/T ray (5.9.3), not used in this version. C C1,S1...Cosine and sine of the angle of incidence (5.9.3,5.9.5). C C2,S2...Cosine and sine of the R/T angle (5.9.3,5.9.5). C C,S... Cosine and sine of the revolution angle C (5.9.3,5.9.5,5.9.6). C V11,V21,V31... Velocity gradient in local Cartesian coordinates C corresponding to the incident wave (5.9.4). C V12,V22,V32... Velocity gradient in local Cartesian coordinates C corresponding to the incident wave (5.9.4). C D11,D12,D22... Matrix of the curvature of the interface (5.9.4). C AUX1,AUX2,AUX3... Temporary storage locations (5.9.4). C P... Snell ray parameter(5.9.5,5.9.6) C AUX,E12,CFC11,CFC12,CFC22,Q1,Q2,P1,P2... Temporary storage C locations (5.9.5). C COEF... (5.9.6). C RMOD,RPH... P-SV R/T coefficient (5.9.6). C RMODSH,RPHSH... SH R/T coefficient (5.9.6). C RR,YR...Real and imaginary part of a complex-valued R/T C coefficient (5.9.6). C ( R1A1,Y1A1,R2A1,Y2A1 ) C ( R1B1,Y1B1,R2B1,Y2B1 )... Amplidude coefficients (5.9.6): C ( R1A2,Y1A2,R2A2,Y2A2 ) C ( R1B2,Y1B2,R2B2,Y2B2 ) Initial C Part: Component: polarisation: Wave: C R...Real 1...P-SV A...First 1...Incident C Y...Imaginary 2...SH B...Second 2...R/T C C....................................................................... C C C (5.9.1) Transformation of auxiliary quantities, travel time and C coordinates: C (a) Number of real-valued quantities describing the reduced C amplitudes of the incident wave: NAMPL1=IY(1)-27 C (b) Number of real-valued quantities describing the reduced C amplitudes of the generated wave: IF(IY(5).GT.0.AND.ICBNEW.LT.0) THEN NAMPL2=NAMPL1*2 ELSE IF(IY(5).LT.0.AND.ICBNEW.GT.0) THEN NAMPL2=NAMPL1/2 ELSE NAMPL2=NAMPL1 END IF C (c) Output number of basic quantities: IY(1)=27+NAMPL2 C (d) New position in the code: IY(2)=KODNEW C (e) Indices of simple and complex blocks: ISB1=IY(4) ICB1=IY(5) ICB2=IABS(IY(8)) IF(IABS(ICBNEW).EQ.ICB2) THEN C Transmission: IY(3)=IABS(IY(5)) IY(4)=IY(7) IY(5)=ICBNEW IY(6)=-IY(6) ELSE IF(IABS(ICBNEW).EQ.IABS(IY(5))) THEN C Reflection: IY(5)=ICBNEW ELSE C Required complex block is not attainable: IEND=37 RETURN END IF C (f) IY(7), IY(8) may be undefined IY(7)=0 IY(8)=0 C (g) Number of transformations at an interface: IY(11)=IY(11)+1 C (h) IY(9), IY(10), IY(12), YY(1:5) are unchanged C (i) Travel time and coordinates Y(1:5) remain unchanged C C (5.9.2) Velocities (metric tensor is determined later on) C (b) Material parameters corresponding to the incident wave: VP1=YL(1) VS1=YL(2) RO1=YL(3) VD1(2)=YL(4) VD1(3)=YL(5) VD1(4)=YL(6) IF(ICB1.GE.0) THEN V1=VP1 ELSE V1=VS1 END IF C (c) Material parameters on the other side of the interface: IF(ICB2.NE.0) THEN CALL PARM2(ICB2,Y(3),UP,US,RO2,QP,QS) CALL VELOC(ISIGN(ICB2,ICBNEW),UP,US,QP,QS,VP2,VS2,VD2,QL) ELSE C Free space: RO2=0. END IF IF(VP1.EQ.VP2) THEN IF(ICBNEW.NE.ISIGN(ICB2,ICB1)) THEN C Transmission without conversion is excluded IF(VS1.EQ.VS2.AND.RO1.EQ.RO2) THEN C Reflection or wave type conversion at a fictitious interface IEND=32 RETURN END IF END IF END IF C (d) Material parameters corresponding to the generated wave: IF(ICBNEW.EQ.ICB1) THEN C Reflection without conversion: V2=V1 VD2(2)=VD1(2) VD2(3)=VD1(3) VD2(4)=VD1(4) ELSE IF(ICBNEW.EQ.-ICB1) THEN C Reflection with conversion: CALL PARM2(IABS(ICBNEW),Y(3),UP,US,AUX,QP,QS) CALL VELOC(ICBNEW,UP,US,QP,QS,AUX,AUX,VD2,QL) V2=VD2(1) ELSE C Transmission: V2=VD2(1) IF(RO2.EQ.0.) THEN C Transmission into free space: IEND=24 RETURN END IF IF(V2.EQ.0.) THEN C Zero velocity (S-wave in liquid): IEND=26 RETURN END IF YL(1)=VP2 YL(2)=VS2 YL(3)=RO2 END IF YL(4)=VD2(2) YL(5)=VD2(3) YL(6)=VD2(4) C C (5.9.3) Transformation of the slowness vector and of the C basis of the ray-centred coordinate system: CALL SRFC2(IABS(IY(6)),Y(3),FAUX) C Non-unit vectors - covariant components H13=Y(6) H23=Y(7) H33=Y(8) H11=Y(9) H21=Y(10) H31=Y(11) Z13=FAUX(2) Z23=FAUX(3) Z33=FAUX(4) C (c) Non-unit vectors - contravariant components IF(KOOR().EQ.0) THEN C Cartesian coordinates: C Norm of the slowness vector HHH=SQRT(H13*H13+H23*H23+H33*H33) C Norm of the gradient of the function describing the interface ZZZ=SQRT(Z13*Z13+Z23*Z23+Z33*Z33) C Unit normal to the interface: C Covariant components Z13=Z13/ZZZ Z23=Z23/ZZZ Z33=Z33/ZZZ C Contravariant components Z43=Z13 Z53=Z23 Z63=Z33 ELSE C Curvilinear coordinates: CALL METRIC(Y(3),GSQRD,G,GAMMA) C Second covariant derivatives of the interface function FAUX(5)=FAUX(5) * -GAMMA(1)*FAUX(2)-GAMMA(7)*FAUX(3)-GAMMA(13)*FAUX(4) FAUX(6)=FAUX(6) * -GAMMA(2)*FAUX(2)-GAMMA(8)*FAUX(3)-GAMMA(14)*FAUX(4) FAUX(7)=FAUX(7) * -GAMMA(3)*FAUX(2)-GAMMA(9)*FAUX(3)-GAMMA(15)*FAUX(4) FAUX(8)=FAUX(8) * -GAMMA(4)*FAUX(2)-GAMMA(10)*FAUX(3)-GAMMA(16)*FAUX(4) FAUX(9)=FAUX(9) * -GAMMA(5)*FAUX(2)-GAMMA(11)*FAUX(3)-GAMMA(17)*FAUX(4) FAUX(10)=FAUX(10) * -GAMMA(6)*FAUX(2)-GAMMA(12)*FAUX(3)-GAMMA(18)*FAUX(4) C Slowness vector - contravariant components CALL SMVPRD(G(7),H13,H23,H33,H43,H53,H63) C Norm of the slowness vector HHH=SQRT(H13*H43+H23*H53+H33*H63) C Contravariant gradient of the function describing the interface CALL SMVPRD(G(7),Z13,Z23,Z33,Z43,Z53,Z63) C Norm of the gradient of the function describing the interface ZZZ=SQRT(Z13*Z43+Z23*Z53+Z33*Z63) C Unit normal to the interface: C Covariant components Z13=Z13/ZZZ Z23=Z23/ZZZ Z33=Z33/ZZZ C Contravariant components Z43=Z43/ZZZ Z53=Z53/ZZZ Z63=Z63/ZZZ END IF C Unit vector tangent to the incident ray H13=H13/HHH H23=H23/HHH H33=H33/HHH C Cosine of the angle of incidence C1=H13*Z43+H23*Z53+H33*Z63 C Vector tangent to the interface in the plane of incidence Z11=H13-Z13*C1 Z21=H23-Z23*C1 Z31=H33-Z33*C1 IF(KOOR().EQ.0) THEN C Cartesian coordinates: AUX=Z11*Z11+Z21*Z21+Z31*Z31 ELSE C Curvilinear coordinates: Z41=H43/HHH-Z43*C1 Z51=H53/HHH-Z53*C1 Z61=H63/HHH-Z63*C1 AUX=Z11*Z41+Z21*Z51+Z31*Z61 END IF C Sine of the angle of incidence IF(AUX.LT.0.5) THEN S1=SQRT(AUX) ELSE S1=SQRT(1.-C1*C1) END IF IF(S1.LT.0.0001) THEN C Correction for nearly normal incidence AUX1=Z43*H11+Z53*H21+Z63*H31 AUX2=SQRT(1.-AUX1*AUX1) Z11=(H11-Z13*AUX1)/AUX2 Z21=(H21-Z23*AUX1)/AUX2 Z31=(H31-Z33*AUX1)/AUX2 ELSE Z11=Z11/S1 Z21=Z21/S1 Z31=Z31/S1 END IF C Angle of reflection/transmission S2=S1*V2/V1 C2=1.-S2*S2 IF(C2.LE.0.) THEN C Overcritical ray IEND=25 RETURN ELSE C2=SIGN(SQRT(C2),C1) END IF IF(IABS(ICBNEW).EQ.IABS(ICB1)) C2=-C2 C Correction for nearly grazing rays (not included in the paper) I=ISIDE(IY(6),ISB1) IF(I.EQ.2) THEN PAUSE 'Error 591 in TRANS: Wrong interface' STOP ELSE IF(FLOAT(I)*C1.GT.0.) THEN C Incident ray goes from the interface - interchanging R/T C2=-C2 END IF IF(KOOR().EQ.0) THEN C Cartesian coordinates: C Vectors tangent to the interface: C Vector in the plane of incidence - contravariant components Z41=Z11 Z51=Z21 Z61=Z31 C Vector perpendicular to the plane of incidence - covariant comp. Z12=Z23*Z31-Z33*Z21 Z22=Z33*Z11-Z13*Z31 Z32=Z13*Z21-Z23*Z11 C Vector perpendicular to the plane of incidence - contravariant Z42=Z12 Z52=Z22 Z62=Z32 ELSE C Curvilinear coordinates: C Vectors tangent to the interface: C Vector in the plane of incidence - contravariant components CALL SMVPRD(G(7),Z11,Z21,Z31,Z41,Z51,Z61) C Vector perpendicular to the plane of incidence - covariant comp. Z12=(Z53*Z61-Z63*Z51)*GSQRD Z22=(Z63*Z41-Z43*Z61)*GSQRD Z32=(Z43*Z51-Z53*Z41)*GSQRD C Vector perpendicular to the plane of incidence - contravariant Z42=(Z23*Z31-Z33*Z21)/GSQRD Z52=(Z33*Z11-Z13*Z31)/GSQRD Z62=(Z13*Z21-Z23*Z11)/GSQRD END IF C Revolution angle of the ray centred coordinate system C=(H11*Z41+H21*Z51+H31*Z61)/C1 S=-H11*Z42-H21*Z52-H31*Z62 C Unit vector perpendicular to the R/T ray in the plane of incidence E11=Z11*C2-Z13*S2 E21=Z21*C2-Z23*S2 E31=Z31*C2-Z33*S2 C Unit vector tangent to the R/T ray C1 E13=Z11*S2+Z13*C2 C2 E23=Z21*S2+Z23*C2 C3 E33=Z31*S2+Z33*C2 C Slowness vector C4 Y(6)=E13/V2 C5 Y(7)=E23/V2 C6 Y(8)=E33/V2 C For a nearly normal incidence, the vector (Z11,Z21,Z31) is not C situated in the plane of incidence and the above six equations C cannot be used. Thus, they are replaced by the following three C ones Y(6)=Y(6)+Z13*(C2/V2-C1/V1) Y(7)=Y(7)+Z23*(C2/V2-C1/V1) Y(8)=Y(8)+Z33*(C2/V2-C1/V1) C First basis vector of the ray centred coordinate system Y(9) =E11*C-Z12*S Y(10)=E21*C-Z22*S Y(11)=E31*C-Z32*S C C (5.9.4) Curvature of the interface and velocity gradients in the C local Cartesian coordinate system: CALL SMVPRD(FAUX(5),Z41,Z51,Z61,AUX1,AUX2,AUX3) C Note: In the original paper on C.R.T., the curvature matrix D is C erroneously defined with opposite sign. D11=-(AUX1*Z41+AUX2*Z51+AUX3*Z61)/ZZZ D12=-(AUX1*Z42+AUX2*Z52+AUX3*Z62)/ZZZ CALL SMVPRD(FAUX(5),Z42,Z52,Z62,AUX1,AUX2,AUX3) D22=-(AUX1*Z42+AUX2*Z52+AUX3*Z62)/ZZZ V11=VD1(2)*Z41+VD1(3)*Z51+VD1(4)*Z61 V21=VD1(2)*Z42+VD1(3)*Z52+VD1(4)*Z62 V31=VD1(2)*Z43+VD1(3)*Z53+VD1(4)*Z63 V12=VD2(2)*Z41+VD2(3)*Z51+VD2(4)*Z61 V22=VD2(2)*Z42+VD2(3)*Z52+VD2(4)*Z62 V32=VD2(2)*Z43+VD2(3)*Z53+VD2(4)*Z63 C C (5.9.5) Dynamic ray tracing across a curved interface P=S1/V1 AUX=C1/V1-C2/V2 E11=(S1*C1*V31-(1.+C1*C1)*V11)/V1-(S2*C2*V32-(1.+C2*C2)*V12)/V2 E12=V22/V2-V21/V1 CFC11=(AUX*D11+P*E11)/C2/C2 CFC12=(AUX*D12+P*E12)/C2 CFC22= AUX*D22 AUX=C1/C2 DO 51 I=1,4 J=4*I Q1= C*Y(8+J)+S*Y(9+J) Q2=-S*Y(8+J)+C*Y(9+J) P1= C*Y(10+J)+S*Y(11+J) P2=-S*Y(10+J)+C*Y(11+J) Q1=Q1/AUX P1=P1*AUX+CFC11*Q1+CFC12*Q2 P2=P2 +CFC12*Q1+CFC22*Q2 Y(8+J)=C*Q1-S*Q2 Y(9+J)=S*Q1+C*Q2 Y(10+J)=C*P1-S*P2 Y(11+J)=S*P1+C*P2 51 CONTINUE C C (5.9.6) Transformation of reduced amplitudes IF(IABS(ICBNEW).EQ.IABS(ICB1)) THEN C Reflection: NCODE=1 COEF=1. ELSE C Transmission: NCODE=3 COEF=RO2/RO1 END IF COEF=SQRT(COEF*ABS(C2/C1)*V2/V1) IF(C1.GE.0) THEN C Epsilon=+1. (see C.R.T.5.9.7) ND=0 ELSE C Epsilon=-1. (see C.R.T.5.9.7) ND=1 END IF RMODSH=0. R2A2=0. Y2A2=0. R2B2=0. Y2B2=0. IF(ICB1.GE.0) THEN C Incident P wave: IF(ICBNEW.LT.0) THEN C Outgoing S wave: NCODE=NCODE+1 END IF R1A1=Y(28) Y1A1=Y(29) IF(NAMPL1.GT.2) THEN R1B1= Y(30) Y1B1= Y(31) END IF ELSE C Incident S wave: IF(ICB2.EQ.0) THEN NCODE=NCODE+2 ELSE NCODE=NCODE+4 END IF R1A1= C*Y(28)+S*Y(30) Y1A1= C*Y(29)+S*Y(31) IF(NAMPL1.GT.4) THEN R1B1= C*Y(32)+S*Y(34) Y1B1= C*Y(33)+S*Y(35) END IF IF(ICBNEW.LT.0) THEN C Outgoing S wave: NCODE=NCODE+1 IF(ICBNEW.EQ.ICB1) THEN CALL COEFSH(P,VS1,RO1,VS2,RO2,1,RMODSH,RPHSH) ELSE CALL COEFSH(P,VS1,RO1,VS2,RO2,2,RMODSH,RPHSH) END IF RMODSH=COEF*RMODSH RR=RMODSH*COS(RPHSH) YR=RMODSH*SIN(RPHSH) R2A1=-S*Y(28)+C*Y(30) Y2A1=-S*Y(29)+C*Y(31) R2A2=RR*R2A1-YR*Y2A1 Y2A2=YR*R2A1+RR*Y2A1 IF(NAMPL1.GT.4) THEN R2B1=-S*Y(32)+C*Y(34) Y2B1=-S*Y(33)+C*Y(35) R2B2=RR*R2B1-YR*Y2B1 Y2B2=YR*R2B1+RR*Y2B1 END IF END IF END IF CALL COEF50(P,VP1,VS1,RO1,VP2,VS2,RO2,NCODE,ND,RMOD,RPH) RMOD=COEF*RMOD IF(RMOD.EQ.0..AND.RMODSH.EQ.0.) THEN C Zero reflection/transmission coefficient IEND=33 RETURN END IF RR=RMOD*COS(RPH) YR=RMOD*SIN(RPH) R1A2=RR*R1A1-YR*Y1A1 Y1A2=YR*R1A1+RR*Y1A1 IF(ICBNEW.GE.0) THEN C Outgoing P wave: Y(28)=R1A2 Y(29)=Y1A2 IF(NAMPL2.GT.2) THEN R1B2=RR*R1B1-YR*Y1B1 Y1B2=YR*R1B1+RR*Y1B1 Y(30)=R1B2 Y(31)=Y1B2 END IF ELSE C Outgoing S wave: Y(28)= C*R1A2-S*R2A2 Y(29)= C*Y1A2-S*Y2A2 Y(30)= S*R1A2+C*R2A2 Y(31)= S*Y1A2+C*Y2A2 IF(NAMPL2.GT.4) THEN R1B2=RR*R1B1-YR*Y1B1 Y1B2=YR*R1B1+RR*Y1B1 Y(32)= C*R1B2-S*R2B2 Y(33)= C*Y1B2-S*Y2B2 Y(34)= S*R1B2+C*R2B2 Y(35)= S*Y1B2+C*Y2B2 END IF END IF C IEND=0 RETURN END C C======================================================================= C SUBROUTINE CONV(KSTORE,YL,Y,IY,NAMPL,YC) INTEGER KSTORE,IY(12),NAMPL REAL YL(6),Y(35),YC(12) C C This subroutine replaces the amplitudes Y(28) to Y(IY(1)) expressed in C the ray-centred coordinate system by the amplitudes involving C appropriate conversion coefficients, see C.R.T.5.5.4. C C Input: C KSTORE..Specifies whether the conversion coefficients are to be C considered (see C.R.T.5.5.4 and 5.6e): C KSTORE.LE.1... No amplitude conversion coefficients are C applied at the point of intersection with an interface. C KSTORE.GE.2... Amplitude conversion coefficients are C applied at the point of intersection with an interface. C YL... Array containing local quantities at the point of C incidence. C Y... Array of the dimension at least 39. The locations Y(1) to C Y(IY(1)) contain basic quantities computed along the ray, C corresponding to the incident wave. C IY... Array containing integer auxiliary quantities computed C along the ray, corresponding to the incident wave. C YC... Array of the dimension at least 12. C C Output: C NAMPL...Number of real quantities representing complex-valued C vectorial reduced amplitudes. If no conversion C coefficients are applied NAMPL=IY(1)-27, otherwise C NAMPL=6 or 12 (see C.R.T.5.5.4). C YC... Array containing real quantities representing complex- C -valued vectorial reduced amplitudes. If no conversion C coefficients are applied, YC is a copy of Y(28) to C Y(IY(1)). Otherwise, yc represents the vectorial reduced C amplitudes involving appropriate conversion coefficients, C expressed in ray-centred coordinate system (see C C.R.T.5.5.4): C P wave at the initial point of the ray: C NAMPL=6, C YC(1)=REAL(A13), YC(2)=AIMAG(A13), C YC(3)=REAL(A23), YC(4)=AIMAG(A23), C YC(5)=REAL(A33), YC(6)=AIMAG(A33). C S wave at the initial point of the ray: C NAMPL=12, C YC(1)=REAL(A11), YC(2)=AIMAG(A11), C YC(3)=REAL(A21), YC(4)=AIMAG(A21), C YC(5)=REAL(A31), YC(6)=AIMAG(A31), C YC(7)=REAL(A12), YC(8)=AIMAG(A12), C YC(9)=REAL(A22), YC(10)=AIMAG(A22), C YC(11)=REAL(A32), YC(12)=AIMAG(A32). C C Subroutines and external functions required: EXTERNAL VELOC,KOOR,METRIC,SRFC2,PARM2,SMVPRD,COEF50 INTEGER KOOR C VELOC... File 'model.for' of the package 'MODEL'. C KOOR,METRIC... File 'metric.for' of the package 'MODEL'. C SRFC2 and subsequent routines... File 'srfc.for' and subsequent C files (like 'val.for' and 'fit.for') of the package C 'MODEL'. C PARM2 and subsequent routines... File 'parm.for' and subsequent C files (like 'val.for' and 'fit.for') of the package C 'MODEL'. C SMVPRD... File 'means.for' of the package 'MODEL'. C COEF50... File 'coef.for'. C C Date: 1990, October 1 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations for local model parameters: FAUX(10), C G(12),GAMMA(18),GSQRD, UP(10),US(10),RO,QP,QS, VP,VS,VD(10),QL: INCLUDE 'auxmod.inc' C C....................................................................... C C Auxiliary storage locations: C INTEGER ICB2,NCODE,ND,NW,IW,I REAL VP1,VS1,RO1,V1,VP2,VS2,RO2,V2 REAL H11,H13,H43,HHH,Z13,Z42,Z43,ZZZ REAL H21,H23,H53, Z23,Z52,Z53 REAL H31,H33,H63, Z33,Z62,Z63 REAL P,C1,S1,C2,S2,C,S,CR,CI,SR,SI REAL RMOD,RPH,RR,YR REAL R1A1,Y1A1,R2A1,Y2A1,R1B1,Y1B1,R2B1,Y2B1 REAL R1A2,Y1A2,R2A2,Y2A2,R1B2,Y1B2,R2B2,Y2B2,RA2,YA2,RB2,YB2 C C ICB2... Index of a complex block at the other side of the C interface. Auxiliary variable. C NCODE,ND... Type of the R/T coefficient. C NW... Total number of R/T waves. C IW... Loop variable. Index of an R/T wave. C I... Auxiliary loop variable. C VP1,VS1,RO1... Parameters of the complex block IABS(IY(5)). C V1... Incident wave velocity. C VP2,VS2,RO2... Parameters of the complex block ICB2. C V2... R/T wave velocity. C H11,H21,H31,H13,H23,H33... Covariant components of the basis C vectors of ray-centred coordinate system of the incident C wave. C H43,H53,H63... Contravariant components of the slowness vector of C the incident wave. Not used in Cartesian coordinates. C HHH... Norm of the slowness vector. C Z13,Z23,Z33... Covariant components of the local Cartesian C coordinate system. C Z42,Z52,Z62,Z43,Z53,Z63... Contravariant components of the local C Cartesian coordinate system. C ZZZ... Norm of the gradient of the function describing the C interface. C C1,S1...Cosine and sine of the angle of incidence. C C2,S2...Cosine and sine of the R/T angle. Auxiliary variables. C C,S... Cosine and sine of the revolution angle. C P... Snell ray parameter. C CR,CI,SR,SI...Real and imaginary parts of the cosine and sine of C the rotation angle in the plane of incidence. C RMOD,RPH... Any R/T coefficient. C RR,YR...Real and imaginary part of a complex-valued R/T C coefficient. C ( R1A1,Y1A1,R2A1,Y2A1 )... Amplidude coefficients in ray-centred C ( R1B1,Y1B1,R2B1,Y2B1 ) coordinate system of the incident wave C ( R1A2,Y1A2,R2A2,Y2A2 ) C ( R1B2,Y1B2,R2B2,Y2B2 ) Initial C Part: Component: polarisation: Wave: C R...Real 1...P-SV, SV A...First 1...Incident C Y...Imaginary 2...SH I I B...Second 2...R/T C I I...R/T wave C I...Incident wave C ( RA2,YA2 )... P-sv amplidude coefficients in ray-centred C ( RB2,YB2 ) coordinate system of the R/T wave C Part: Polarisation: Wave: C R...Real A...First 2...R/T C Y...Imaginary B...Second C C....................................................................... C C Check whether the amplitude conversion coefficients are applied IF(KSTORE.LE.1) THEN NAMPL=IY(1)-27 DO 10 I=1,NAMPL YC(I)=Y(27+I) 10 CONTINUE RETURN END IF C C Material parameters corresponding to the incident wave: VP1=YL(1) VS1=YL(2) RO1=YL(3) IF(IY(5).GE.0) THEN V1=VP1 ELSE V1=VS1 END IF C Material parameters on the other side of the interface: ICB2=IABS(IY(8)) IF(ICB2.NE.0) THEN CALL PARM2(ICB2,Y(3),UP,US,RO2,QP,QS) CALL VELOC(ICB2,UP,US,QP,QS,VP2,VS2,VD,QL) NW=4 ELSE C Free space: RO2=0. NW=2 END IF C VP1, VS1, RO1, V1, VP2, VS2, RO2 and NW are defined. C C Basis of the ray-centred coordinate system and the local Cartesian C coordinate system connected with the interface: CALL SRFC2(IABS(IY(6)),Y(3),FAUX) C Non-unit vectors - covariant components H13=Y(6) H23=Y(7) H33=Y(8) H11=Y(9) H21=Y(10) H31=Y(11) Z13=FAUX(2) Z23=FAUX(3) Z33=FAUX(4) C (c) Non-unit vectors - contravariant components IF(KOOR().EQ.0) THEN C Cartesian coordinates: GSQRD=1. C Norm of the slowness vector HHH=SQRT(H13*H13+H23*H23+H33*H33) C Norm of the gradient of the function describing the interface ZZZ=SQRT(Z13*Z13+Z23*Z23+Z33*Z33) C Unit normal to the interface: C Covariant components Z13=Z13/ZZZ Z23=Z23/ZZZ Z33=Z33/ZZZ C Contravariant components Z43=Z13 Z53=Z23 Z63=Z33 ELSE C Curvilinear coordinates: CALL METRIC(Y(3),GSQRD,G,GAMMA) C Slowness vector - contravariant components CALL SMVPRD(G(7),H13,H23,H33,H43,H53,H63) C Norm of the slowness vector HHH=SQRT(H13*H43+H23*H53+H33*H63) C Contravariant gradient of the function describing the interface CALL SMVPRD(G(7),Z13,Z23,Z33,Z43,Z53,Z63) C Norm of the gradient of the function describing the interface ZZZ=SQRT(Z13*Z43+Z23*Z53+Z33*Z63) C Unit normal to the interface: C Covariant components Z13=Z13/ZZZ Z23=Z23/ZZZ Z33=Z33/ZZZ C Contravariant components Z43=Z43/ZZZ Z53=Z53/ZZZ Z63=Z63/ZZZ END IF C Unit vector tangent to the incident ray H13=H13/HHH H23=H23/HHH H33=H33/HHH C Cosine of the angle of incidence C1=H13*Z43+H23*Z53+H33*Z63 C Sine of the angle of incidence S1=SQRT(1.-C1*C1) C C1 and S1 are the cosine and sine of the angle of incidence. C C Revolution angle of the basis of the ray-centred coordinate system C around ray: IF(S1.LT.0.0001) THEN C Nearly normal incidence C=1. S=0. ELSE C Vector perpendicular to the plane of incidence - contravariant Z42=(Z23*H33-Z33*H23) Z52=(Z33*H13-Z13*H33) Z62=(Z13*H23-Z23*H13) C Revolution angle of the ray centred coordinate system C=-(H11*Z43+H23*Z53+H31*Z63)/S1 S=-(H11*Z42+H21*Z52+H31*Z62)/S1/GSQRD END IF C C and S are the cosine and sine of the revolution angle. C C Amplitudes of the incident wave in the rotated ray-centred C coordinate system of the incident wave: NAMPL=6 IF(IY(5).GE.0) THEN C Incident P wave: NCODE=0 R1A1=Y(28) Y1A1=Y(29) R1A2=0. Y1A2=0. R2A2=0. Y2A2=0. YC(5)=R1A1 YC(6)=Y1A1 IF(IY(1).GT.29) THEN NAMPL=12 R1B1= Y(30) Y1B1= Y(31) R1B2=0. Y1B2=0. R2B2=0. Y2B2=0. YC(11)=R1B1 YC(12)=Y1B1 END IF ELSE C Incident S wave: NCODE=NW R1A1= C*Y(28)+S*Y(30) Y1A1= C*Y(29)+S*Y(31) R2A1=-S*Y(28)+C*Y(30) Y2A1=-S*Y(29)+C*Y(31) R1A2=R1A1 Y1A2=Y1A1 R2A2=R2A1 Y2A2=Y2A1 YC(5)=0. YC(6)=0. IF(IY(1).GT.31) THEN NAMPL=12 R1B1= C*Y(32)+S*Y(34) Y1B1= C*Y(33)+S*Y(35) R2B1=-S*Y(32)+C*Y(34) Y2B1=-S*Y(33)+C*Y(35) R1B2=R1B1 Y1B2=Y1B1 R2B2=R2B1 Y2B2=Y2B1 YC(11)=0. YC(12)=0. END IF END IF C C Snell parameter P=S1/V1 C IF(C1.GE.0) THEN C Epsilon=+1. (see C.R.T.5.9.7) ND=0 ELSE C Epsilon=-1. (see C.R.T.5.9.7) ND=1 END IF C C Loop over the generated waves: DO 90 IW=1,2 NCODE=NCODE+1 C C Propagation velocity of the generated wave IF(IW.EQ.1) THEN V2=VP1 ELSE IF(IW.EQ.2) THEN V2=VS1 ELSE IF(IW.EQ.3) THEN V2=VP2 ELSE V2=VS2 END IF C C Rotation angle of the basis of the ray-centred coordinate system C in the plane of incidence: S2=P*V2 C2=1.-S2*S2 CR=S1*S2 SR=C1*S2 IF(C2.GE.0.) THEN C2=SIGN(SQRT(C2),C1) IF(IW.LE.2) C2=-C2 CR=CR+C1*C2 SR=SR-S1*C2 CI=0. SI=0. ELSE C2=SIGN(SQRT(-C2),C1) IF(IW.LE.2) C2=-C2 CI= C1*C2 SI=-S1*C2 END IF C Cosine of the rotation angle is CR+i*CI, where i=SQRT(-1), C sine of the rotation angle is SR+i*SI. C C Transformation of reduced amplitudes (5.9.6) CALL COEF50(P,VP1,VS1,RO1,VP2,VS2,RO2,NCODE,ND,RMOD,RPH) RR=RMOD*COS(RPH) YR=RMOD*SIN(RPH) RA2=RR*R1A1-YR*Y1A1 YA2=YR*R1A1+RR*Y1A1 IF(NAMPL.GT.6) THEN RB2=RR*R1B1-YR*Y1B1 YB2=YR*R1B1+RR*Y1B1 END IF IF(IW.EQ.1.OR.IW.EQ.3) THEN C Outgoing P wave: R1A2=R1A2+SR*RA2-SI*YA2 Y1A2=Y1A2+SI*RA2+SR*YA2 YC(5)=YC(5)+CR*RA2-CI*YA2 YC(6)=YC(6)+CI*RA2+CR*YA2 IF(NAMPL.GT.6) THEN R1B2=R1B2+SR*RB2-SI*YB2 Y1B2=Y1B2+SI*RB2+SR*YB2 YC(11)=YC(11)+CR*RB2-CI*YB2 YC(12)=YC(12)+CI*RB2+CR*YB2 END IF ELSE C Outgoing S wave: R1A2=R1A2+CR*RA2-CI*YA2 Y1A2=Y1A2+CI*RA2+CR*YA2 YC(5)=YC(5)-SR*RA2+SI*YA2 YC(6)=YC(6)-SI*RA2-SR*YA2 IF(NAMPL.GT.6) THEN R1B2=R1B2+CR*RB2-CI*YB2 Y1B2=Y1B2+CI*RB2+CR*YB2 YC(11)=YC(11)-SR*RB2+SI*YB2 YC(12)=YC(12)-SI*RB2-SR*YB2 END IF IF(IY(5).LT.0) THEN C SH wave coefficients: IF(IW.EQ.2) THEN CALL COEFSH(P,VS1,RO1,VS2,RO2,1,RMOD,RPH) ELSE CALL COEFSH(P,VS1,RO1,VS2,RO2,2,RMOD,RPH) END IF RR=RMOD*COS(RPH) YR=RMOD*SIN(RPH) R2A2=R2A2+RR*R2A1-YR*Y2A1 Y2A2=Y2A2+YR*R2A1+RR*Y2A1 IF(NAMPL.GT.6) THEN R2B2=R2B2+RR*R2B1-YR*Y2B1 Y2B2=Y2B2+YR*R2B1+RR*Y2B1 END IF END IF END IF 90 CONTINUE C C Rotation of the basis of the ray-centred coordinate system around C ray: YC(1)= C*R1A2-S*R2A2 YC(2)= C*Y1A2-S*Y2A2 YC(3)= S*R1A2+C*R2A2 YC(4)= S*Y1A2+C*Y2A2 IF(NAMPL.GT.6) THEN YC(7)= C*R1B2-S*R2B2 YC(8)= C*Y1B2-S*Y2B2 YC(9)= S*R1B2+C*R2B2 YC(10)=S*Y1B2+C*Y2B2 END IF C RETURN END C C======================================================================= C