C
C Subroutine file 'rm.for' to support program 'greenss.for'. C C Version: 5.20 C Date: 1998, September 6 C C======================================================================= C SUBROUTINE RM(NLAY,IFREE,ITER,OMEGA,TR,P1,P2,P3, * AR1,AI1,AR2,AI2,AR3,AI3,BR1,BI1,BR2,BI2,BR3,BI3) INTEGER NLAY,IFREE,ITER REAL OMEGA,TR,P1,P2,P3 REAL AR1,AI1,AR2,AI2,AR3,AI3,BR1,BI1,BR2,BI2,BR3,BI3 C C Subroutine to calculate the response of the transmission through the C stack of fine horizontal layers at the receiver site. The parameters C of the layer stack should have already been stored in the memory by C the invocation of subroutine INPUT. The frequency-dependent C transmission coefficiens are calculated by subroutines of package C RMATRIX C written by C. J. Thomson (Copyright (c) 1996 C. J. Thomson). C C Input: C NLAY... Number of layers, including the bottom and top halfspaces. C IFREE...IFREE=1: The thickness THICKN(1) of the first (top) layer C is used for the calculation of the travel time C correction at the geophone situated in the top halfspace C at the elevation THICKN(1) above the uppermost C interface. C IFREE=0: Has the same effect as THICKN(1)=0. C ITER... ITER=1: Polarization vectors in the layers are calculated. C ITER.GT.1: Polarization vectors from the previous C invocation are used. C OMEGA...Angular frequency. C TR... Travel time at the receiver. C P1,P2,P3... Slowness vector at the receiver. C AR1,AI1,AR2,AI2,AR3,AI3... Complex-valued displacement at the C bottom of the layer stack. C C Output: C TR... Travel time at the bottom of the layer stack. C BR1,BI1,BR2,BI2,BR3,BI3... Complex-valued displacement at the C geophone (IFREE=1) or at the top of the layer stack C (IFREE=0). C C Date: 1998, September 4 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Common blocks for RMATRIX: INTEGER NLAYMX PARAMETER (NLAYMX=30) INTEGER ISOFLG DOUBLE PRECISION THICKN,A,RHO COMMON /ELAST/ISOFLG(NLAYMX),THICKN(NLAYMX),A(3,3,3,3,NLAYMX), * RHO(NLAYMX) COMPLEX*16 EVALS,EVECS COMMON /EVPROB/EVALS(6,NLAYMX),EVECS(6,6,NLAYMX) C C Input and output arguments for RMATRIX: COMPLEX*16 PX,PY,COMEGA COMPLEX*16 TUS(3,3),RUS(3,3),TDS(3,3),RDS(3,3) C C Temporary storage locations: COMPLEX*16 E11,E21,E31,E12,E22,E32,E13,E23,E33 COMPLEX F11,F21,F31,F12,F22,F32,F13,F23,F33,E,A1,A2,A3,B1,B2,B3 C C----------------------------------------------------------------------- C C C Travel-time correction corresponding to the elevation of the C bottom of the layer stack with respect to the receiver: TR=TR+THICKN(NLAY)*ABS(P3) C C Calculation of the transmission matrices: COMEGA=DCMPLX(DBLE(OMEGA),0.0D0) PX=DCMPLX(DBLE(P1),0.0D0) PY=DCMPLX(DBLE(P2),0.0D0) CALL RMATRIX(PX,PY,COMEGA,NLAY,IFREE,TUS,RUS,TDS,RDS,ITER) C C Complex-valued amplitude at the bottom of the layer stack A1=CMPLX(AR1,AI1) A2=CMPLX(AR2,AI2) A3=CMPLX(AR3,AI3) C C Transforming the amplitude from Cartesian coordinates to C the eigenvectors Eij E11=EVECS(1,4,NLAY) E21=EVECS(2,4,NLAY) E31=EVECS(3,4,NLAY) E12=EVECS(1,5,NLAY) E22=EVECS(2,5,NLAY) E32=EVECS(3,5,NLAY) E13=EVECS(1,6,NLAY) E23=EVECS(2,6,NLAY) E33=EVECS(3,6,NLAY) C Fij=inv(Eij)*det(E) F11=E22*E33-E23*E32 F12=E32*E13-E33*E12 F13=E12*E23-E13*E22 F21=E23*E31-E21*E33 F22=E33*E11-E31*E13 F23=E13*E21-E11*E23 F31=E21*E32-E22*E31 F32=E31*E12-E32*E11 F33=E11*E22-E12*E21 C E=det(E) E=E11*F11+E12*F21+E13*F31 C Bi=inv(Eij)*Ai B1=(F11*A1+F12*A2+F13*A3)/E B2=(F21*A1+F22*A2+F23*A3)/E B3=(F31*A1+F32*A2+F33*A3)/E C C Transition through the layer stack A1=TUS(1,1)*B1+TUS(1,2)*B2+TUS(1,3)*B3 A2=TUS(2,1)*B1+TUS(2,2)*B2+TUS(2,3)*B3 A3=TUS(3,1)*B1+TUS(3,2)*B2+TUS(3,3)*B3 C C Transforming the amplitude from the eigenvectors to the C Cartesian coordinates B1=EVECS(1,4,1)*A1+EVECS(1,5,1)*A2+EVECS(1,6,1)*A3 B2=EVECS(2,4,1)*A1+EVECS(2,5,1)*A2+EVECS(2,6,1)*A3 B3=EVECS(3,4,1)*A1+EVECS(3,5,1)*A2+EVECS(3,6,1)*A3 C C Complex-valued amplitude at the top of the layer stack BR1=REAL (B1) BI1=AIMAG(B1) BR2=REAL (B2) BI2=AIMAG(B2) BR3=REAL (B3) BI3=AIMAG(B3) C RETURN END C C======================================================================= C