SUBROUTINE RAYB(W,T0,DT) C C INITIAL-VALUE RAY TRACING BY THE RUNGE-KUTTA METHOD C COMPLEX AY,W(2) COMMON /FORCE/ F(3) COMMON /RAYW/ AY(3,2000),E(3,3,2000),OMEGA,N,NTOT,IND,IND1 COMMON /ISOTR/T(2000),XX(2000),Y(2000),Z(2000),VP(2000),VS(2000), 1RHO(2000) C X=T0 PI=3.14159632 RHO0=RHO(1) VS0=VS(1) AUX=4*PI*SQRT(RHO0)*RHO0*VS0*VS0 W(1)=(E(1,1,1)*F(1)+E(1,2,1)*F(2)+E(1,3,1)*F(3))/AUX W(2)=(E(2,1,1)*F(1)+E(2,2,1)*F(2)+E(2,3,1)*F(3))/AUX AY(1,1)=T0 AY(2,1)=W(1) AY(3,1)=W(2) C C SOLVING OF THE SYSTEM OF COUPLED EQUATIONS C CALL RK(X,W) RETURN END C C *********************************************************** C SUBROUTINE RK(X,Y) C C MODIFIED EULER'S METHOD TO SOLVE A SYSTEM OF LINEAR C ORDINARY DIFFERENTIAL EQUATIONS OF FIRST ORDER C COMPLEX AY,Y(2),DERY(2),Y1(2),Y2(2) COMMON /RAYW/ AY(3,2000),HHH(3,3,2000),OMEGA,N,NTOT,IND,IND1 COMMON /ISOTR/T(2000),XX(2000),YY(2000),ZZ(2000),VP(2000), 1VS(2000),RHO(2000) C NDIM=2 N=1 1 H=T(N+1)-T(N) X=X+H CALL FCTA(Y,DERY,NDIM) DO 2 I=1,NDIM Y1(I)=DERY(I) 2 Y2(I)=Y(I)+H*Y1(I) N=N+1 CALL FCTA(Y2,DERY,NDIM) DO 3 I=1,NDIM 3 Y(I)=Y(I)+.5*H*(Y1(I)+DERY(I)) CALL OUT(X,Y,DERY,NDIM) IF(N.EQ.NTOT) GO TO 4 GO TO 1 4 RETURN END C C******************************************************** C SUBROUTINE FCTA(W,DERY,NDIM) C C COMPUTATION OF THE RIGHT HAND SIDES OF THE QI EQUATIONS C SAVE VS2 COMPLEX AY,DERY(NDIM),W(NDIM),IMAG DIMENSION YY(3) COMMON /RAYW/ AY(3,2000),E(3,3,2000),OMEGA,N,NTOT,IND,IND1 COMMON /ISOTR/T(2000),X(2000),Y(2000),Z(2000),VP(2000),VS(2000), 1RHO(2000) C IMAG=(0.,1.) VS2=VS(N)*VS(N) YY(1)=X(N) YY(2)=Y(N) YY(3)=Z(N) CALL PARD(YY) CALL BIJ(B11,1,1) CALL BIJ(B12,1,2) CALL BIJ(B22,2,2) C C MODIFICATION OF WEAK ANISOTROPY MATRIX BY ROUTINE NEWB C IF(IND1.EQ.1)THEN VSS=VS(N) CALL NEWB(VSS,B11,B12,B22) END IF C DERY(1)=-0.5*IMAG*(OMEGA/VS2)*(W(1)*B11+W(2)*B12) DERY(2)=-0.5*IMAG*(OMEGA/VS2)*(W(1)*B12+W(2)*B22) RETURN END C C**************************************************** C SUBROUTINE OUT(X,Y,DERY,NDIM) C COMPLEX AY,Y(2),DERY(2) COMMON /RAYW/ AY(3,2000),E(3,3,2000),OMEGA,N,NTOT,IND,IND1 C AY(1,N)=X AY(2,N)=Y(1) AY(3,N)=Y(2) RETURN END C C********************************************************* C SUBROUTINE BIJ(B,M,N) C COMPLEX AY COMMON /RAYW/ AY(3,2000),E(3,3,2000),OMEGA,NN,NTOT,IND,IND1 COMMON /ELPAR/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33, 1 A34,A35,A36,A44,A45,A46,A55,A56,A66 COMMON /ISOTR/T(2000),X(2000),Y(2000),Z(2000),VP(2000),VS(2000), 1RHO(2000) C A2344=A23+A44 A3645=A36+A45 A2546=A25+A46 A1355=A13+A55 A1456=A14+A56 A1266=A12+A66 T1=E(3,1,NN) T2=E(3,2,NN) T3=E(3,3,NN) EM1=E(M,1,NN) EM2=E(M,2,NN) EM3=E(M,3,NN) EN1=E(N,1,NN) EN2=E(N,2,NN) EN3=E(N,3,NN) T2T3=T2*T3 T1T2=T1*T2 T1T3=T1*T3 T1T1=T1*T1 T2T2=T2*T2 T3T3=T3*T3 C11=T1T1*A11+T2T2*A66+T3T3*A55 1+2.*(T2T3*A56+T1T3*A15+T1T2*A16) C22=T1T1*A66+T2T2*A22+T3T3*A44 1+2.*(T2T3*A24+T1T3*A46+T1T2*A26) C33=T1T1*A55+T2T2*A44+T3T3*A33 1+2.*(T2T3*A34+T1T3*A35+T1T2*A45) C23=T1T1*A56+T2T2*A24+T3T3*A34 1 +T2T3*A2344+T1T3*A3645+T1T2*A2546 C13=T1T1*A15+T2T2*A46+T3T3*A35 1 +T2T3*A3645+T1T3*A1355+T1T2*A1456 C12=T1T1*A16+T2T2*A26+T3T3*A45 1 +T2T3*A2546+T1T3*A1456+T1T2*A1266 B=C11*EM1*EN1+C22*EM2*EN2+C33*EM3*EN3+C12*(EM1*EN2+EM2*EN1) 1 +C13*(EM1*EN3+EM3*EN1)+C23*(EM2*EN3+EM3*EN2) IF(M.EQ.N)THEN V2=VS(NN)*VS(NN) IF(M.EQ.3)V2=VP(NN)*VP(NN) B=B-V2 END IF RETURN END C C********************************************************* C SUBROUTINE BIAP(MX1,MX,MY1,MY,MXY1) C DIMENSION X(200),FX(200),V(1000) COMMON/ZCOEF/ A02(1000),A20(1000),A22(1000) COMMON /INTRF/ Z(1000),SX(350),SY(350),NX(20),NY(20),BRD(6),NINT, 1 XINTA EQUIVALENCE(Z(1),V(1)) C C ROUTINE DETERMINING THE COEFFICIENTS C OF BICUBIC SPLINE INTERPOLATION C DO 1 J=1,MX L=MX1+J-1 1 X(J)=SX(L) DO 3 I=1,MY DO 2 J=1,MX K=MXY1+(J-1)*MY+I-1 2 FX(J)=V(K) CALL SPLIN(X,FX,1,MX) DO 3 J=1,MX K=MXY1+(J-1)*MY+I-1 3 A20(K)=FX(J) C DO 4 I=1,MY L=MY1+I-1 4 X(I)=SY(L) DO 6 J=1,MX DO 5 I=1,MY K=MXY1+(J-1)*MY+I-1 5 FX(I)=V(K) CALL SPLIN(X,FX,1,MY) DO 6 I=1,MY K=MXY1+(J-1)*MY+I-1 6 A02(K)=FX(I) C DO 7 J=1,MX L=MX1+J-1 7 X(J)=SX(L) DO 9 I=1,MY DO 8 J=1,MX K=MXY1+(J-1)*MY+I-1 8 FX(J)=A02(K) CALL SPLIN(X,FX,1,MX) DO 9 J=1,MX K=MXY1+(J-1)*MY+I-1 9 A22(K)=FX(J) C RETURN END C C********************************************************* C SUBROUTINE DISC(Y,DEP) C C DETERMINATION OF DEPTH OF 3D INTERFACES AND ITS DERIVATIVES C FOR BICUBIC POLYNOMIAL APPROXIMATION C SAVE B,DX,DY,llay,ibb,iu,il,ju,jl DIMENSION Y(18),DEP(6),B(16,2),DX(2),DY(2) C COMMON /AUXI/ IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOUT, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori COMMON /AUXX/ MMX(20),MMY(20),MMXY(20) COMMON /INTRF/ Z(1000),SX(350),SY(350),NX(20),NY(20),BRD(6),NINT, 1 XINTA COMMON/ZCOEF/ A02(1000),A20(1000),A22(1000) C IBB=0 LB=2 IF(INTR.EQ.LAY)LB=1 MX=NX(INTR) MY=NY(INTR) DO 1 I=2,MX K=MMX(INTR)+I-1 IF(Y(1).LT.SX(K))GO TO 2 1 CONTINUE 2 I1=K DO 3 I=2,MY K=MMY(INTR)+I-1 IF(Y(2).LT.SY(K))GO TO 4 3 CONTINUE 4 J1=K IF(MAUX.EQ.0) GOTO 8 IF(LB.EQ.2) GOTO 5 IF(I1.EQ.IU.AND.J1.EQ.JU.AND.LLAY.EQ.LAY) GOTO 10 GOTO 8 5 IF(I1.EQ.IL.AND.J1.EQ.JL.AND.LLAY.EQ.LAY) GOTO 10 IL=I1 JL=J1 GOTO 9 8 IU=I1 JU=J1 9 LLAY=LAY I=I1-MMX(INTR) J=J1-MMY(INTR) DX(LB)=SX(I1-1) DY(LB)=SY(J1-1) MM=MMXY(INTR)-1 K=MM+(I-1)*MY+J B20=A20(K) B02=A02(K) B22=A22(K) B00=Z(K) K=MM+I*MY+J C20=A20(K) C02=A02(K) C22=A22(K) C00=Z(K) K=MM+(I-1)*MY+J+1 D20=A20(K) D02=A02(K) D22=A22(K) D00=Z(K) K=MM+I*MY+J+1 E20=A20(K) E02=A02(K) E22=A22(K) E00=Z(K) HX=SX(I1)-DX(LB) HY=SY(J1)-DY(LB) XA=3.*HX YA=3.*HY XB=HX/3. YB=HY/3. D32=(E22-D22)/XA D30=(E20-D20)/XA B30=(C20-B20)/XA B32=(C22-B22)/XA D12=(E02-D02)/HX-XB*(E22+2.*D22) D10=(E00-D00)/HX-XB*(E20+2.*D20) B10=(C00-B00)/HX-XB*(C20+2.*B20) B12=(C02-B02)/HX-XB*(C22+2.*B22) B03=(D02-B02)/YA B13=(D12-B12)/YA B23=(D22-B22)/YA B33=(D32-B32)/YA B01=(D00-B00)/HY-YB*(D02+2.*B02) B11=(D10-B10)/HY-YB*(D12+2.*B12) B21=(D20-B20)/HY-YB*(D22+2.*B22) B31=(D30-B30)/HY-YB*(D32+2.*B32) MAUX=1 B(1,LB)=B00 B(2,LB)=B01 B(3,LB)=B02 B(4,LB)=B03 B(5,LB)=B10 B(6,LB)=B11 B(7,LB)=B12 B(8,LB)=B13 B(9,LB)=B20 B(10,LB)=B21 B(11,LB)=B22 B(12,LB)=B23 B(13,LB)=B30 B(14,LB)=B31 B(15,LB)=B32 B(16,LB)=B33 IBB=1 10 AX=Y(1)-DX(LB) AZ=Y(2)-DY(LB) IF(IBB.EQ.1) GOTO 11 B00=B(1,LB) B01=B(2,LB) B02=B(3,LB) B03=B(4,LB) B10=B(5,LB) B11=B(6,LB) B12=B(7,LB) B13=B(8,LB) B20=B(9,LB) B21=B(10,LB) B22=B(11,LB) B23=B(12,LB) B30=B(13,LB) B31=B(14,LB) B32=B(15,LB) B33=B(16,LB) 11 AUX1=((B33*AZ+B32)*AZ+B31)*AZ+B30 AUX2=((B23*AZ+B22)*AZ+B21)*AZ+B20 AUX3=((B13*AZ+B12)*AZ+B11)*AZ+B10 AUX4=((B03*AZ+B02)*AZ+B01)*AZ+B00 DEP(1)=((AUX1*AX+AUX2)*AX+AUX3)*AX+AUX4 IF(NDER.EQ.0) RETURN DEP(2)=(3.*AUX1*AX+2.*AUX2)*AX+AUX3 IF(NDER.EQ.1)GO TO 7 DEP(4)=6.*AUX1*AX+2.*AUX2 7 AUX1=(3.*B33*AZ+2.*B32)*AZ+B31 AUX2=(3.*B23*AZ+2.*B22)*AZ+B21 AUX3=(3.*B13*AZ+2.*B12)*AZ+B11 AUX4=(3.*B03*AZ+2.*B02)*AZ+B01 DEP(3)=((AUX1*AX+AUX2)*AX+AUX3)*AX+AUX4 IF(NDER.EQ.1) RETURN DEP(5)=(3.*AUX1*AX+2.*AUX2)*AX+AUX3 AUX1=3.*B33*AZ+B32 AUX2=3.*B23*AZ+B22 AUX3=3.*B13*AZ+B12 AUX4=3.*B03*AZ+B02 DEP(6)=2.*(((AUX1*AX+AUX2)*AX+AUX3)*AX+AUX4) RETURN END C C********************************************************* C SUBROUTINE SPLIN(X,FX,NMIN,NMAX) C DIMENSION A(200),B(200),H(200),F(200),X(200),FX(200) C C CUBIC SPLINE INTERPOLATION WITH ZERO CURVATURES AT C THE END POINTS C IF((NMAX-NMIN).EQ.1)GO TO 4 NMIN1=NMIN+1 NMAX1=NMAX-1 H(NMIN1)=X(NMIN1)-X(NMIN) D2=(FX(NMIN1)-FX(NMIN))/H(NMIN1) DO 1 I=NMIN1,NMAX1 H(I+1)=X(I+1)-X(I) D1=D2 D2=(FX(I+1)-FX(I))/H(I+1) B(I)=H(I)+H(I+1) FX(I)=3.*(D2-D1)/B(I) A(I)=H(I)/B(I) 1 B(I)=H(I+1)/B(I) 4 FX(NMIN)=0. FX(NMAX)=0. IF((NMAX-NMIN).EQ.1)RETURN H(NMIN)=0. F(NMIN)=0. DO 2 I=NMIN1,NMAX1 XPOM=2.+A(I)*H(I-1) H(I)=-B(I)/XPOM 2 F(I)=(FX(I)-A(I)*F(I-1))/XPOM DO 3 I=NMIN,NMAX1 J=NMAX1-(I-NMIN) 3 FX(J)=H(J)*FX(J+1)+F(J) RETURN END C C ********************************************************* C SUBROUTINE FACETS(N1,N2,NSRF) INTEGER LU,N1,N2,NSRF C C Subroutine FACETS writes the index file listing the vertices of each C tetragon covering the structural interface. The vertices are assumed C to be stored in a separate file, with inner loop over N1 points along C the first horizontal axis, middle loop over N2 points along the second C horizontal axis and outer loop over the surfaces. The vertices are C indexed by positive integers according to their order in the vertex C file. C C Input: C LU... Logical unit number connected to the output file to be C written by this subroutine. C N1... Number of points along the first horizontal axis. C N2... Number of points along the second horizontal axis. C NSRF... Number of interfaces. C The input parameters are not altered. C C No output. C C Output index file with the tetragons: C For each tetragon, a line containing I1,I2,I3,I4,/ C I1,I2,I3,I4... Indices of the vertices of the tetragon. C The vertices are indexed by positive integers according to C their order in the respective vertex file. C /... List of vertices is terminated by a slash. C C Date: 1999, October 4 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: CHARACTER*9 FORMAT INTEGER I1,I2,ISRF COMMON/VRML/LUBRD,LUGRD,LU,LURAY C IF(LU.EQ.0)RETURN C Setting output format: FORMAT='(4(I0,A))' I1=INT(ALOG10(FLOAT(N1*N2*NSRF)+0.5))+1 FORMAT(5:5)=CHAR(ICHAR('0')+I1) C C Writing the file: DO 33 ISRF=0,N1*N2*(NSRF-1),N1*N2 DO 32 I2=ISRF,ISRF+N1*(N2-2),N1 DO 31 I1=I2+1,I2+N1-1 WRITE(LU,FORMAT) I1,' ',I1+1,' ',I1+1+N1,' ',I1+N1,' /' 31 CONTINUE 32 CONTINUE 33 CONTINUE C RETURN END C C ********************************************************* C SUBROUTINE NEWB(VS,B11,B12,B22) REAL VS,B11,B12,B22 C C Subroutine to slightly modify matrix (7) of Psencik (1998): C Quasi-shear waves in the zero-order approximation of the C quasi-isotropic approach. Preliminary results. In: Seismic Waves C in Complex 3-D Structures, Report 7, pp.225-266. C The modification makes the coupling equations dependent on the C reference isotropic model just in terms of the reference ray path. C C Input: C VS... S-wave velocity in the reference isotropic model. C B11,B12,B22... Components of symmetric 2*2 matrix B, C B(M,N) = DeltaA_ijkl E(M)_i E(3)_j E(3)_k E(N)_l C C Output: C B11,B12,B22... Components of symmetric 2*2 matrix C 2*{E-[E+B*VS**(-2)]**(-1/2)}*VS**2, C where E is the unit matrix and B is the input matrix. C C Date: 1999, March 13 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C REAL VS2,A11,A12,A22,AD,AT,AS,BD,BT,BS C C B=E+B*VS**(-2): VS2=VS**2 B11=B11/VS2+1. B12=B12/VS2 B22=B22/VS2+1. C C A=SQRT(B): BD=B11*B22-B12*B12 BT=B11+B22 BS=B11-B22 AD=SQRT(BD) AT=SQRT(BT+AD+AD) AS=BS/AT A11=0.5*(AT+AS) A12=B12/AT A22=0.5*(AT-AS) C C B=2*[A**(-1)-E]*VS**2: B11=2.*(1.-A22/AD)*VS2 B12=2.*( A12/AD)*VS2 B22=2.*(1.-A11/AD)*VS2 C RETURN END