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),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: 1992, DECEMBER 19 C CODED BY LUDEK KLIMES C C----------------------------------------------------------------------- C C AUXILIARY STORAGE LOCATIONS FOR LOCAL MODEL PARAMETERS: REAL G(12),GAMMA(18),GSQRD,FAUX(10) REAL UP(10),US(10),RO,QP,QS,VP,VS,VD(10),QL C THESE AUXILIARY VARIABLES AND ARRAYS NEED NOT BE LOCATED IN A C COMMON BLOCK. THERE IS NO REASON TO LOCATE THEM IN THE AUXILIARY C COMMON BLOCK /AUXMOD/ BUT TO SHARE THE MEMORY. COMMON/AUXMOD/ G,GAMMA,GSQRD,FAUX,UP,US,RO,QP,QS,VP,VS,VD,QL C C G,GAMMA,GSQRD... SEE SUBROUTINE METRIC OF THE FILE 'METRIC.FOR'. C FAUX... AUXILIARY ARRAY TO STORE A FUNCTIONAL VALUE AND ITS C DERIVATIVES. C UP,US,RO,QP,QS... SEE SUBROUTINE PARM2 OF THE FILE 'PARM.FOR'. C UP,US,QP,QS,VP,VS,VD,QL... SEE SUBROUTINE VELOC OF THE FILE C 'MODEL.FOR'. 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 NAMPL1..NUMBER OF THE APMLITUDE COEFFICIENTS OF THE INCIDENT WAVE C (5.9.1A). C NAMPL2..NUMBER OF THE APMLITUDE 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(NAMPL1.GT.4) THEN R1B2=RR*R1B1-YR*Y1B1 Y1B2=YR*R1B1+RR*Y1B1 END IF IF(ICBNEW.GE.0) THEN C OUTGOING P WAVE: Y(28)=R1A2 Y(29)=Y1A2 IF(NAMPL2.GT.2) THEN 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(NAMPL1.GT.4) THEN 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: REAL G(12),GAMMA(18),GSQRD,FAUX(10) REAL UP(10),US(10),RO,QP,QS,VP,VS,VD(10),QL C THESE AUXILIARY VARIABLES AND ARRAYS NEED NOT BE LOCATED IN A C COMMON BLOCK. THERE IS NO REASON TO LOCATE THEM IN THE AUXILIARY C COMMON BLOCK /AUXMOD/ BUT TO SHARE THE MEMORY. COMMON/AUXMOD/ G,GAMMA,GSQRD,FAUX,UP,US,RO,QP,QS,VP,VS,VD,QL C C G,GAMMA,GSQRD... SEE SUBROUTINE METRIC OF THE PACKAGE 'METRIC'. C FAUX... AUXILIARY ARRAY TO STORE A FUNCTIONAL VALUE AND ITS C DERIVATIVES. C UP,US,RO,QP,QS... SEE SUBROUTINE PARM2 OF THE FILE 'PARM.FOR'. C UP,US,QP,QS,VP,VS,VD,QL... SEE SUBROUTINE VELOC OF THE FILE C 'MODEL.FOR'. 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 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