a2.for 0100666 0000765 0000765 00000065313 11023551462 011452 0 ustar bulant bulant C
SUBROUTINE COEF(STU,UC,ITRANS) C C ROUTINE FOR THE COMPUTATION OF REFLECTION/TRANSMISSION COEFFICIENTS C FOR GENERAL ISOTROPIC/ANISOTROPIC MEDIA C C INPUT PARAMETERS: C STU(1-3)... CARTESIAN COMPONENTS OF THE DISPLACEMENT VECTOR C STU(4-6)... CARTESIAN COMPONETS OF THE TRACTION VECTOR CALCULATED C IN THIS ROUTINE C ITRANS... ITRANS=O: REFLECTION C ITRANS=1: TRANSMISSION C ITRANS=999: SURFACE CONVERSION C C OUTPUT PARAMETERS: C UC(1-3)... RAY-CENTERED COMPONENTS OF THE DISPLACEMENT VECTOR C COMPLEX A(6,6),C(6,6),COA(3,3),COC(3,3),AA(6),TAU(3),P(3),U(3), 1CP1,CP2,CP3,DETA,DETC,UC(3),STU(6),CAB DIMENSION Y(18),UN(3) 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 INTEGER CODE COMMON /COD/CODE(50,2),KREF,KC,ITYPE COMPLEX PS COMMON /RAY/ AY(28,2000),DS(20,50),KINT(50),HHH(3,3),TMAX, 1 PS(3,7,50),IS(8,50),N,IREF,IND,IND1 COMMON /ZERO/ RNULL C INDI=0 IF(ITRANS.EQ.999) THEN ITRANS=0 INDI=1 END IF I1=1 I2=3 LAY1=IS(1,IREF) LAY=IS(7,IREF) DO 910 K=1,6 IF(K.GT.3) GOTO 900 UN(K)=DS(K,IREF) 900 Y(K)=AY(K+1,N) 910 CONTINUE IF(ICOEF.GT.0) 1WRITE(LOUT,'(A,/,3(3F14.7,/),3I5)')'COEF: Y(1)-Y(6),UN,IREF,LAY, 2ITYPE',(Y(K),K=1,6),UN,IREF,LAY,ITYPE CALL PARDIS (Y,0) C C DISPLACEMENT AND STRESS VECTOR OF INCIDENT WAVE C BSTU=SQRT(REAL(STU(1)*CONJG(STU(1))+STU(2)*CONJG(STU(2)) 1 +STU(3)*CONJG(STU(3)))) BST=STU(1)*HHH(ITYPE,1)+STU(2)*HHH(ITYPE,2) 1 +STU(3)*HHH(ITYPE,3) BSTU=SIGN(1.,BST)*BSTU DO 950 K=1,3 P(K)=PS(K,7,IREF) STU(K)=STU(K)/BSTU U(K)=STU(K) 950 CONTINUE CALL STRESS (P,UN,U,TAU) DO 960 K=1,3 STU(K+3)=TAU(K) 960 CONTINUE IF(ICOEF.GT.0) 1WRITE(LOUT,'(A,/,6(2F14.7,/))') 'COEF: INCIDENT U,TAU',U,TAU 1000 IF(IANI(LAY).EQ.0) GOTO 3000 C C COMPUTATION OF STRESS-DISPLACEMENT-VECTORS FOR GENERATED WAVES C C ANISOTROPIC CASE C 1010 IF(I1.GT.3)CALL PARDIS(Y,0) DO 1020 K=I1,I2 DO 1030 L=1,3 P(L)=PS(L,K,IREF) 1030 CONTINUE C C EVALUATION OF DISPLACEMENT VECTOR (GENERALLY COMPLEX-VALUED) C CALL DISPL(P,U,K) CALL STRESS(P,UN,U,TAU) IF(ICOEF.GT.0)WRITE(LOUT,'(A,/,9(2F14.7,/))')'COEF: P,U,TAU' 1,P,U,TAU DO 1060 L=1,3 A(K,L)=U(L) A(K,L+3)=TAU(L) 1060 CONTINUE 1020 CONTINUE IF(INTR.EQ.1.AND.INDI.EQ.1) GOTO 5000 LAY=IS(8,IREF) IF(LAY.EQ.0) THEN LAY=IS(7,IREF) GO TO 5000 END IF IF(I2.EQ.6) GOTO 4000 I1=4 I2=6 IF(IANI(LAY).NE.0) GOTO 1010 C C ISOTROPIC CASE C 3000 IF(I1.GT.3) CALL PARDIS(Y,0) DO 3010 L=1,3 P(L)=PS(L,I1,IREF) 3010 CONTINUE CP1=P(1) CP2=P(2) IVERT=0 IF(CABS(CP1).LT..0000001.AND.CABS(CP2).LT..0000001)IVERT=1 IF(IVERT.EQ.1)THEN U(1)=1. U(2)=0. U(3)=0. GO TO 3020 END IF CP3=-(CP1*CP1+CP2*CP2) CP1=CP1*P(3) CP2=CP2*P(3) CAB=SQRT(REAL(CP1*CONJG(CP1)+CP2*CONJG(CP2)+CP3*CONJG(CP3))) U(1)=CP1/CAB U(2)=CP2/CAB U(3)=CP3/CAB 3020 CONTINUE CALL STRESS(P,UN,U,TAU) IF(ICOEF.GT.0)WRITE(LOUT,'(A,/,9(2F14.7/))')'COEF: P,U,TAU' 1,P,U,TAU DO 3030 L=1,3 A(I1,L)=U(L) A(I1,L+3)=TAU(L) 3030 CONTINUE IF(IVERT.EQ.1)THEN U(1)=0. U(2)=1. IF(REAL(P(3)).LT.0.)U(2)=-1. U(3)=0. GO TO 3040 END IF CP1=-P(2) CP2=P(1) CP3=CMPLX(0.,0.) CAB=SQRT(REAL(CP1*CONJG(CP1)+CP2*CONJG(CP2)+CP3*CONJG(CP3))) U(1)=CP1/CAB U(2)=CP2/CAB U(3)=CP3/CAB 3040 CONTINUE CALL STRESS (P,UN,U,TAU) IF(ICOEF.GT.0)WRITE(LOUT,'(A,/,9(2F14.7/))')'COEF: P,U,TAU' 1,P,U,TAU DO 3050 L=1,3 A(I1+1,L)=U(L) A(I1+1,L+3)=TAU(L) 3050 CONTINUE DO 3210 L=1,3 P(L)=PS(L,I2,IREF) 3210 CONTINUE CP1=P(1) CP2=P(2) CP3=P(3) CAB=SQRT(REAL(CP1*CONJG(CP1)+CP2*CONJG(CP2)+CP3*CONJG(CP3))) U(1)=CP1/CAB U(2)=CP2/CAB U(3)=CP3/CAB CALL STRESS (P,UN,U,TAU) IF(ICOEF.GT.0)WRITE(LOUT,'(A,/,9(2F14.7/))')'COEF: P,U,TAU' 1,P,U,TAU DO 3220 L=1,3 A(I2,L)=U(L) A(I2,L+3)=TAU(L) 3220 CONTINUE IF(INTR.EQ.1.AND.INDI.EQ.1) GOTO 5000 LAY=IS(8,IREF) IF(LAY.EQ.0) THEN LAY=IS(7,IREF) GO TO 5000 END IF IF(I2.EQ.6)GOTO 4000 I1=4 I2=6 GOTO 1000 C C REVERSE OF SIGNS FOR REFLECTIONS C 4000 DO 4010 L=1,3 DO 4010 K=1,6 A(L,K)=-A(L,K) 4010 CONTINUE IF(ICOEF.EQ.0) GOTO 4012 WRITE(LOUT,'(//,A)') 'COEF: MATRIX A' DO 4011 K=1,6 WRITE(LOUT,'(6(12F10.5))') (A(L,K),L=1,6) 4011 CONTINUE C C SOLUTION OF SYSTEM OF LINEAR EQUATIONS (KRAMER'S METHOD) C 4012 CALL DET6(A,DETA) DO 4020 K=1,3 UC(K)=CMPLX(0.,0.) 4020 CONTINUE IF(CABS(DETA).LT.0.0000001) THEN WRITE(LOUT,'(A)') 1 'COEF: MATRIX FOR R/T COEF. SINGULAR: AMPLITUDE NOT EVALUATED' IND=11 RETURN END IF IF(INDI.EQ.1) THEN L1=1 L2=3 GO TO 4041 END IF IF(IANI(LAY1).NE.0.OR.ITYPE.EQ.3)THEN L1=ITYPE IF(ITRANS.NE.0)L1=L1+3 L2=L1 END IF IF(IANI(LAY1).EQ.0.AND.ITYPE.NE.3)THEN L1=1 IF(ITRANS.NE.0)L1=L1+3 L2=L1+1 END IF 4041 CONTINUE IF(ICOEF.NE.0) WRITE(LOUT,'(A,2I5)')'COEF: L1,L2',L1,L2 DO 4044 L=L1,L2 DO 4043 J=1,6 DO 4042 K=1,6 C(J,K)=A(J,K) IF(L.EQ.J)C(J,K)=STU(K) 4042 CONTINUE 4043 CONTINUE CALL DET6(C,DETC) AA(L)=DETC/DETA IF(ICOEF.GT.0)WRITE(LOUT,'(A,2F12.6)') 'COEF: AA=',AA(L) 4044 CONTINUE IF(INDI.EQ.1)GO TO 5045 C DO 4100 L=L1,L2 LL=L IF(ITRANS.NE.0)LL=L-3 UC(LL)=AA(L)*BSTU 4100 CONTINUE RETURN C C INTERACTION WITH THE FREE SURFACE C 5000 DO 5010 L=1,3 DO 5010 K=1,3 A(L,K)=-A(L,K) COA(L,K)=-A(L,K+3) 5010 CONTINUE IF(ICOEF.EQ.0)GO TO 5012 WRITE(LOUT,'(//,A)')'COEF: MATRIX COA' WRITE(LOUT,'(3(6F10.5/))')((COA(L,K),L=1,3),K=1,3) C C SOLUTION OF SYSTEM OF 3 LINEAR EQUATIONS (CRAMER'S METHOD) C 5012 CALL DET3(COA,DETA) IF(CABS(DETA).LT.0.0000001)THEN WRITE(LOUT,'(A)') 1 'COEF: MATRIX FOR R/T COEF. SINGULAR: AMPLITUDE NOT EVALUATED' IND=11 RETURN END IF DO 5044 L=1,3 DO 5043 J=1,3 DO 5042 K=1,3 COC(J,K)=COA(J,K) IF(L.EQ.J)COC(J,K)=STU(K+3) 5042 CONTINUE 5043 CONTINUE CALL DET3(COC,DETC) AA(L)=DETC/DETA IF(ICOEF.GT.0) 1WRITE(LOUT,'(A,2F12.6)')'COEF: AA=',AA(L) 5044 CONTINUE IF(INDI.NE.1)GO TO 5200 C C CONVERSION AT AN INTERFACE C 5045 CONTINUE IF(MREG.EQ.2) THEN C C CALCULATION OF CONVERSION FOR PRESSURE C CP1=STU(1) CP2=STU(2) CP3=STU(3) ARE=REAL(CP1) AIM=AIMAG(CP1) APHI=ATAN2(AIM,ARE) ARE=SQRT(REAL(CP1*CONJG(CP1)+CP2*CONJG(CP2)+CP3*CONJG(CP3))) UC(1)=ARE*(1.+AA(3))*CMPLX(COS(APHI),SIN(APHI)) UC(2)=(0.,0.) UC(3)=(0.,0.) IF(ICOEF.GT.0)WRITE(LOUT,'(A,4F10.5)') 'COEF: UC(1),ARE,APHI', 1 UC(1),ARE,APHI GO TO 5120 END IF C DO 5100 K=1,3 U(K)=CMPLX(0.,0.) DO 5050 J=1,3 U(K)=-AA(J)*A(J,K)+U(K) 5050 CONTINUE UC(K)=U(K)+STU(K) 5100 CONTINUE 5120 CONTINUE IF(ICOEF.GT.0)WRITE(LOUT,'(A,/,3(2F12.5/))') 'COEF: 1 NORMALIZED CONVERSION-VECTOR',UC DO 5150 K=1,3 UC(K)=UC(K)*BSTU 5150 CONTINUE IF(ICOEF.GT.0)WRITE(LOUT,'(A,/,3(2F12.5/))') 'COEF: 1 CONVERSION-VECTOR',UC RETURN C C REFLECTION FROM THE FREE SURFACE C 5200 L1=1 IF(ITYPE.EQ.3)L1=3 L2=2 IF(ITYPE.EQ.3)L2=3 DO 5250 K=1,3 UC(K)=CMPLX(0.,0.) 5250 CONTINUE DO 5300 L=L1,L2 UC(L)=AA(L)*BSTU 5300 CONTINUE RETURN END C C ********************************************************* C SUBROUTINE DET6(A,D6) COMPLEX A(6,6),B(5,5),D6,D D6=(0.,0.) DO 10 J=1,6 IF(A(1,J).NE.(0.,0.)) THEN CALL SET(A,B,1,J,6) CALL DET5(B,D) C=(-1)**(J+1) D6=D6+CMPLX(C,0.)*D*A(1,J) ENDIF 10 CONTINUE RETURN END C C ********************************************************* C SUBROUTINE DET5(A,D5) COMPLEX A(5,5),B(4,4),D5,D D5=(0.,0.) DO 10 J=1,5 IF(A(1,J).NE.(0.,0.)) THEN CALL SET(A,B,1,J,5) CALL DET4(B,D) C=(-1)**(J+1) D5=D5+CMPLX(C,0.)*D*A(1,J) ENDIF 10 CONTINUE RETURN END C C *********************************************************** C SUBROUTINE DET4(A,D4) COMPLEX A(4,4),B(3,3),D4,D D4=(0.,0.) DO 10 J=1,4 IF(A(1,J).NE.(0.,0.)) THEN CALL SET(A,B,1,J,4) CALL DET3(B,D) C=(-1)**(J+1) D4=D4+CMPLX(C,0.)*D*A(1,J) ENDIF 10 CONTINUE RETURN END C C ********************************************************** C SUBROUTINE DET3(A,D3) COMPLEX A(3,3),D3 D3=A(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(3,2)) 1+ A(1,2)*(A(2,3)*A(3,1)-A(2,1)*A(3,3)) 2+ A(1,3)*(A(2,1)*A(3,2)-A(2,2)*A(3,1)) RETURN END C C ************************************************************ C SUBROUTINE SET(A,B,I,J,K) COMPLEX A(K,K),B(K-1,K-1) IFL=0 JFL=0 DO 20 I1=1,K IF(I1.EQ.I) THEN IFL=1 GOTO 20 ENDIF DO 10 J1=1,K IF(J1.EQ.J) THEN JFL=1 GOTO 10 ENDIF IF(IFL.EQ.0.AND.JFL.EQ.0) THEN B(I1,J1)=A(I1,J1) ELSE IF(IFL.EQ.0.AND.JFL.EQ.1) THEN B(I1,J1-1)=A(I1,J1) ELSE IF(IFL.EQ.1.AND.JFL.EQ.0) THEN B(I1-1,J1)=A(I1,J1) ELSE B(I1-1,J1-1)=A(I1,J1) ENDIF 10 CONTINUE JFL=0 20 CONTINUE 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 DISPL(P,U,I1) C C ROUTINE FOR THE COMPUTATION OF DISPLACEMENT VECTOR OF GENERATED C WAVES AT AN INTERFACE. IT WORKS EVEN FOR A COMPLEX REFRACTION VECTOR. C ROUTINES CHRM AND POLAR WORK ONLY FOR REAL REFRACTION VECTORS C C INPUT PARAMETERS: C P(1-3)... SLOWNESS VECTOR C I1... TYPE OF THE WAVE WITH THE SLOWNESS VECTOR P C I1=1... REFLECTED S WAVE (SV COMPONENT IN ISOTROPIC CASE) C I1=2... REFLECTED S WAVE (SH COMPONENT IN ISOTROPIC CASE) C I1=3... REFLECTED P WAVE C I1=4... TRANSMITTED S WAVE (SV COMPONENT IN ISOTROPIC CASE) C I1=5... TRANSMITTED S WAVE (SH COMPONENT IN ISOTROPIC CASE) C I1=6... TRANSMITTED P WAVE C C OUTPUT PARAMETERS: C U(1-3)... CORRESPONDING DISPLACEMENT VECTOR C complex uxn(3,3) COMPLEX U,U1,U2,U3,P,P1,P2,P3,P2P3,P1P2,P1P3,P1P1,P2P2,P3P3,C11, 1 C12,C13,C22,C23,C33,C11N,C22N,C33N,C23SQ,C13SQ,C12SQ,CD11,CD12, 2 CD13,CD22,CD23,CD33,CZ01,CZ02,CZ03,CDTR DIMENSION U(3),P(3),xx(3) COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33, 1 A34,A35,A36,A44,A45,A46,A55,A56,A66, 1 DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23, 1 DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44, 1 DXA45,DXA46,DXA55,DXA56,DXA66, 1 DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23, 1 DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44, 1 DYA45,DYA46,DYA55,DYA56,DYA66, 1 DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23, 1 DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44, 1 DZA45,DZA46,DZA55,DZA56,DZA66, 1 A2546,A1266,A1355,A1456,A3645,A2344 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 COMPLEX PS COMMON /RAY/ AY(28,2000),DS(20,50),KINT(50),HHH(3,3),tmax, 1 PS(3,7,50),IS(8,50),N,IREF,IND,IND1 C C C THIS PART CORRESPONDS TO CHRM ,BUT FOR COMPLEX REFRACTION VECTORS C 900 P1=P(1) P2=P(2) P3=P(3) 910 continue IF(ICOEF.GT.0)WRITE(LOUT,'(A,3(/2F12.6))')'DISPL: START P', 1p1,p2,p3 P2P3=P2*P3 P1P2=P1*P2 P1P3=P1*P3 P1P1=P1*P1 P2P2=P2*P2 P3P3=P3*P3 C11=P1P1*A11+P2P2*A66+P3P3*A55 1+2.*P2P3*A56+2.*P1P3*A15+2.*P1P2*A16 C22=P1P1*A66+P2P2*A22+P3P3*A44 1+2.*P2P3*A24+2.*P1P3*A46+2.*P1P2*A26 C33=P1P1*A55+P2P2*A44+P3P3*A33 1+2.*P2P3*A34+2.*P1P3*A35+2.*P1P2*A45 C23=P1P1*A56+P2P2*A24+P3P3*A34 1 +P2P3*A2344+P1P3*A3645+P1P2*A2546 C13=P1P1*A15+P2P2*A46+P3P3*A35 1 +P2P3*A3645+P1P3*A1355+P1P2*A1456 C12=P1P1*A16+P2P2*A26+P3P3*A45 1 +P2P3*A2546+P1P3*A1456+P1P2*A1266 C11N=C11-1. C22N=C22-1. C33N=C33-1. C23SQ=C23*C23 C13SQ=C13*C13 C12SQ=C12*C12 CD11=C22N*C33N-C23SQ CD22=C11N*C33N-C13SQ CD33=C11N*C22N-C12SQ CD12=C13*C23-C12*C33N CD13=C12*C23-C13*C22N CD23=C12*C13-C23*C11N CDTR=CD11+CD22+CD33 IF(ICOEF.GT.0) 1WRITE(LOUT,'(A,4(/4F12.6))')'DISPL: CDIJ' 1,CD11,CD12,CD13,CD22,CD23,CD33,CDTR c if(abs(cdtr).lt..000001)write(lout,'(A)')'DISPL: SHEAR WAVE c 1 SINGULARITY IN CALCULATION OF R/T COEFFICIENT' C C THIS PART CORRESPONDS TO ROUTINE POLAR BUT FOR COMPLEX POLARIZATION C VECTORS C U1=CD11 U2=CD12 U3=CD13 uxn(1,1)=u1 uxn(1,2)=u2 uxn(1,3)=u3 XN=real(U1*CONJG(U1)+U2*CONJG(U2)+U3*CONJG(U3)) xx(1)=xn IF(ICOEF.GT.0)WRITE(LOUT,'(A,F14.7)')'DISPL: XN',XN U1=CD12 U2=CD22 U3=CD23 uxn(2,1)=u1 uxn(2,2)=u2 uxn(2,3)=u3 XN=real(U1*CONJG(U1)+U2*CONJG(U2)+U3*CONJG(U3)) xx(2)=xn IF(ICOEF.GT.0)WRITE(LOUT,'(A,F14.7)')'DISPL: XN',XN U1=CD13 U2=CD23 U3=CD33 uxn(3,1)=u1 uxn(3,2)=u2 uxn(3,3)=u3 XN=real(U1*CONJG(U1)+U2*CONJG(U2)+U3*CONJG(U3)) xx(3)=xn IF(ICOEF.GT.0)WRITE(LOUT,'(A,F14.7)')'DISPL: XN',XN xn=0. do 920 l=1,3 if(xn.ge.xx(l))go to 920 xn=xx(l) nx=l 920 continue xn=1./sqrt(xn) u1=uxn(nx,1)*xn u2=uxn(nx,2)*xn u3=uxn(nx,3)*xn U(1)=U1 U(2)=U2 U(3)=U3 IF(ICOEF.GT.0)WRITE(LOUT,'(A,3(/2F12.6))')'DISPL: U',U C C CHECK OF PRECISION C CZ01=C11N*u1+C12*u2+C13*u3 CZ02=C12*u1+C22N*u2+C23*u3 CZ03=C13*u1+C23*u2+C33N*u3 IF(ICOEF.GT.0) 1WRITE(LOUT,'(A,3(/2F14.7))')'DISPL: PRECISSION OF DISPL' 1,CZ01,CZ02,CZ03 RETURN END C C ********************************************************** C SUBROUTINE DMAT C C EVALUATES ELEMENTS OF THE MATRIX DIJ C INTEGER CODE COMMON /COD/ CODE(50,2),KREF,KC,ITYPE COMMON/DJK/D11,D12,D13,D22,D23,D33,DTR COMMON/GAM/G11,G12,G13,G22,G23,G33 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 COMPLEX PS COMMON /RAY/ AY(28,2000),DS(20,50),KINT(50),HHH(3,3),tmax, 1 PS(3,7,50),IS(8,50),N,IREF,IND,IND1 C G1=G11-1. G2=G22-1. G3=G33-1. D11=G2*G3-G23*G23 D22=G1*G3-G13*G13 D33=G1*G2-G12*G12 D12=G13*G23-G3*G12 D13=G12*G23-G2*G13 D23=G12*G13-G1*G23 DTR=D11+D22+D33 IF(ABS(DTR).LT.0.0001)THEN write(lout,'(A)')'DMAT: SHEAR WAVE SINGULARITY' IND=10 END IF RETURN END C C ********************************************************** C SUBROUTINE DDMAT(DG,g) C C EVALUATES DERIVATIVES OF ELEMENTS OF THE MATRIX DIJ C DIMENSION DG(3,3) COMMON/DDJK/DD11,DD12,DD13,DD22,DD23,DD33,DDTR COMMON/GAM/G11,G12,G13,G22,G23,G33 C DG11=DG(1,1) DG22=DG(2,2) DG33=DG(3,3) DG12=DG(1,2) DG13=DG(1,3) DG23=DG(2,3) G1=G11-1. G2=G22-1. G3=G33-1. DD11=G3*DG22+G2*DG33-2.*G23*DG23+g*(2.-g22-g33) DD22=G3*DG11+G1*DG33-2.*G13*DG13+g*(2.-g33-g11) DD33=G2*DG11+G1*DG22-2.*G12*DG12+g*(2.-g11-g22) DD12=G23*DG13+G13*DG23-DG12*G3-G12*DG33+g*g12 DD13=G23*DG12+G12*DG23-DG13*G2-G13*DG22+g*g13 DD23=G13*DG12+G12*DG13-DG23*G1-G23*DG11+g*g23 DDTR=DD11+DD22+DD33 RETURN END C C ********************************************************** C FUNCTION TR(G) C C TRACE OF PRODUCT OF MATRICES G AND D C DIMENSION G(3,3) COMMON/DJK/ D11,D12,D13,D22,D23,D33,DTR C TR=G(1,1)*D11+G(2,2)*D22+G(3,3)*D33+ 1 2.*(G(1,2)*D12+G(1,3)*D13+G(2,3)*D23) RETURN END C C ********************************************************** C FUNCTION TRD(G) C C TRACE OF PRODUCT OF MATRICES G AND DD C DIMENSION G(3,3) COMMON/DDJK/ DD11,DD12,DD13,DD22,DD23,DD33,DDTR C TRD=G(1,1)*DD11+G(2,2)*DD22+G(3,3)*DD33+ 1 2.*(G(1,2)*DD12+G(1,3)*DD13+G(2,3)*DD23) RETURN END C C ********************************************************** C SUBROUTINE FCT(X,Y,DERY) C C COMPUTATION OF THE RIGHT HAND SIDES OF THE RAY TRACING EQUATIONS C save g,gpp,gpx,gxx,gx,gp,gx1,gx2,gx3,gp1,gp2,gp3 DIMENSION DERY(18),Y(18),AUX(2),AUXX(3,2),VX(3) DIMENSION G(3,3),GPP(3,3),GPX(3,3),GXX(3,3),GX(3),GP(3) DIMENSION GX1(3,3),GX2(3,3),GX3(3,3),GP1(3,3),GP2(3,3),GP3(3,3) COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33, 1 A34,A35,A36,A44,A45,A46,A55,A56,A66, 1 DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23, 1 DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44, 1 DXA45,DXA46,DXA55,DXA56,DXA66, 1 DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23, 1 DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44, 1 DYA45,DYA46,DYA55,DYA56,DYA66, 1 DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23, 1 DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44, 1 DZA45,DZA46,DZA55,DZA56,DZA66, 1 A2546,A1266,A1355,A1456,A3645,A2344 COMMON /APROX1/ E(21),EX(21),EY(21),EZ(21),EXX(21),EXY(21), 1 EXZ(21),EYY(21),EYZ(21),EZZ(21) 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 INTEGER CODE COMMON /COD/ CODE(50,2),KREF,KC,ITYPE COMPLEX PS COMMON /RAY/ AY(28,2000),DS(20,50),KINT(50),HHH(3,3),tmax, 1 PS(3,7,50),IS(8,50),N,IREF,IND,IND1 COMMON /DJK/ D11,D12,D13,D22,D23,D33,DTR COMMON /DDJK/ DD11,DD12,DD13,DD22,DD23,DD33,DDTR COMMON /ZERO/ RNULL C IF(IANI(LAY).NE.0) GOTO 20 CALL PARDIS(Y,0) IF(ITYPE.LT.3)ITT=16 IF(ITYPE.EQ.3)ITT=1 C C COMPUTATION OF RIGHT HAND SIDES OF RAY TRACING EQUATIONS C IN AN ISOTROPIC LAYER C C P-WAVE FOR ITT=1, S-WAVE FOR ITT=16 C V=E(ITT) VX(1)=EX(ITT) VX(2)=EY(ITT) VX(3)=EZ(ITT) DO 1 I=1,3 DERY(I)=V*Y(3+I) DERY(3+I)=-.5*VX(I)/V 1 CONTINUE IF(NDER.EQ.1)RETURN VXX=EXX(ITT) VXY=EXY(ITT) VXZ=EXZ(ITT) VYY=EYY(ITT) VYZ=EYZ(ITT) VZZ=EZZ(ITT) DO 3 J=1,2 AUX(J)=0. DO 2 I=1,3 II=3+I+3*J AUX(J)=AUX(J)+VX(I)*Y(II) 2 CONTINUE JJ=4+3*J AUXX(1,J)=VXX*Y(JJ)+VXY*Y(JJ+1)+VXZ*Y(JJ+2) AUXX(2,J)=VXY*Y(JJ)+VYY*Y(JJ+1)+VYZ*Y(JJ+2) AUXX(3,J)=VXZ*Y(JJ)+VYZ*Y(JJ+1)+VZZ*Y(JJ+2) 3 CONTINUE DO 4 I=1,3 DERY(I+6)=AUX(1)*Y(I+3)+V*Y(I+12) DERY(I+9)=AUX(2)*Y(I+3)+V*Y(I+15) DERY(I+12)=.5*(VX(I)*AUX(1)-V*AUXX(I,1))/V/V DERY(I+15)=.5*(VX(I)*AUX(2)-V*AUXX(I,2))/V/V 4 CONTINUE RETURN C C C COMPUTATION OF RIGHT HAND SIDES OF RAY TRACING EQUATIONS C IN AN ANISOTROPIC LAYER C C DETERMINATION OF PARAMETERS OF THE MEDIUM C 20 CALL PARDIS(Y,0) CALL CHRM2(Y,G,1) CALL DMAT IF(IND.EQ.10)RETURN CALL PCHRM(Y,GP1,1,1) CALL PCHRM(Y,GP2,2,1) CALL PCHRM(Y,GP3,3,1) CALL CHRM2(Y,GX1,2) CALL CHRM2(Y,GX2,3) CALL CHRM2(Y,GX3,4) GP(1)=TR(GP1)/DTR GP(2)=TR(GP2)/DTR GP(3)=TR(GP3)/DTR GX(1)=TR(GX1)/DTR GX(2)=TR(GX2)/DTR GX(3)=TR(GX3)/DTR DO 5 I=1,3 DERY(I)=.5*GP(I) DERY(I+3)=-.5*GX(I) 5 CONTINUE IF(NDER.EQ.1)RETURN CALL CHRM2(Y,G,1) CALL DDMAT(GP1,gp(1)) GPP(1,1)=TRD(GP1)-GP(1)*DDTR GPP(1,2)=TRD(GP2)-GP(2)*DDTR GPP(1,3)=TRD(GP3)-GP(3)*DDTR GPX(1,1)=TRD(GX1)-GX(1)*DDTR GPX(1,2)=TRD(GX2)-GX(2)*DDTR GPX(1,3)=TRD(GX3)-GX(3)*DDTR CALL DDMAT(GP2,gp(2)) GPP(2,1)=TRD(GP1)-GP(1)*DDTR GPP(2,2)=TRD(GP2)-GP(2)*DDTR GPP(2,3)=TRD(GP3)-GP(3)*DDTR GPX(2,1)=TRD(GX1)-GX(1)*DDTR GPX(2,2)=TRD(GX2)-GX(2)*DDTR GPX(2,3)=TRD(GX3)-GX(3)*DDTR CALL DDMAT(GP3,gp(3)) GPP(3,1)=TRD(GP1)-GP(1)*DDTR GPP(3,2)=TRD(GP2)-GP(2)*DDTR GPP(3,3)=TRD(GP3)-GP(3)*DDTR GPX(3,1)=TRD(GX1)-GX(1)*DDTR GPX(3,2)=TRD(GX2)-GX(2)*DDTR GPX(3,3)=TRD(GX3)-GX(3)*DDTR CALL DDMAT(GX1,gx(1)) GXX(1,1)=TRD(GX1)-GX(1)*DDTR GXX(1,2)=TRD(GX2)-GX(2)*DDTR GXX(1,3)=TRD(GX3)-GX(3)*DDTR CALL DDMAT(GX2,gx(2)) GXX(2,1)=TRD(GX1)-GX(1)*DDTR GXX(2,2)=TRD(GX2)-GX(2)*DDTR GXX(2,3)=TRD(GX3)-GX(3)*DDTR CALL DDMAT(GX3,gx(3)) GXX(3,1)=TRD(GX1)-GX(1)*DDTR GXX(3,2)=TRD(GX2)-GX(2)*DDTR GXX(3,3)=TRD(GX3)-GX(3)*DDTR C DO 11 L=1,3 DO 11 M=L,3 CALL PPCHRM(G,L,M,1) AUX1=TR(G) GPP(M,L)=(GPP(M,L)+AUX1)/DTR IF(L.NE.M)GPP(L,M)=(GPP(L,M)+AUX1)/DTR 11 CONTINUE DO 12 L=1,3 CALL PCHRM(Y,G,L,2) GPX(L,1)=(GPX(L,1)+TR(G))/DTR CALL PCHRM(Y,G,L,3) GPX(L,2)=(GPX(L,2)+TR(G))/DTR CALL PCHRM(Y,G,L,4) GPX(L,3)=(GPX(L,3)+TR(G))/DTR 12 CONTINUE CALL CHRM2(Y,G,5) GXX(1,1)=(GXX(1,1)+TR(G))/DTR CALL CHRM2(Y,G,8) GXX(2,2)=(GXX(2,2)+TR(G))/DTR CALL CHRM2(Y,G,10) GXX(3,3)=(GXX(3,3)+TR(G))/DTR CALL CHRM2(Y,G,6) AUX1=(GXX(1,2)+TR(G))/DTR GXX(1,2)=AUX1 GXX(2,1)=AUX1 CALL CHRM2(Y,G,7) AUX1=(GXX(1,3)+TR(G))/DTR GXX(1,3)=AUX1 GXX(3,1)=AUX1 CALL CHRM2(Y,G,9) AUX1=(GXX(2,3)+TR(G))/DTR GXX(2,3)=AUX1 GXX(3,2)=AUX1 CALL CHRM2(Y,G,1) C DO 13 I=1,3 DERY(I+6)=0. DERY(I+9)=0. DERY(I+12)=0. DERY(I+15)=0. DO 13 K=1,3 DERY(I+6)=DERY(I+6)+.5*(GPX(I,K)*Y(K+6)+GPP(I,K)*Y(K+12)) DERY(I+9)=DERY(I+9)+.5*(GPX(I,K)*Y(K+9)+GPP(I,K)*Y(K+15)) DERY(I+12)=DERY(I+12)-.5*(GXX(I,K)*Y(K+6)+GPX(K,I)*Y(K+12)) DERY(I+15)=DERY(I+15)-.5*(GXX(I,K)*Y(K+9)+GPX(K,I)*Y(K+15)) 13 CONTINUE RETURN END Ca3.for 0100666 0000765 0000765 00000067511 11023551526 011456 0 ustar bulant bulant C
SUBROUTINE INIT(PN,V) C C ROUTINE FOR THE CALCULATION OF THE PHASE VELOCITIES FOR THE C SPECIFIED UNIT NORMAL TO THE PHASE FRONT C SOLUTION OF THE EIGENVALUE PROBLEM C DIMENSION C(3,3),PN(3),V(3) DOUBLE PRECISION A1,A2,A3,R,S,T,Q,P,TH COMMON /AUXI/ IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,MORI C CALL CHRM1(C,PN,PN) C C COMPUTATION OF THE ROOTS OF THE CHARACTERISTIC EQUATION WITH C CARDANS FORMULA C C COMPUTATION OF THE COEFFICIENTS OF THE CHARACTERISTIC EQUATION C X**3+R*X**2+S*X+T=0 C C11=C(1,1) C12=C(1,2) C13=C(1,3) C22=C(2,2) C23=C(2,3) C33=C(3,3) A1=C22*C33-C23*C23 A2=C11*C33-C13*C13 A3=C11*C22-C12*C12 R=-(C11+C22+C33) S=A1+A2+A3 A2=C12*C33-C13*C23 A3=C12*C23-C13*C22 T=-C11*A1+C12*A2-C13*A3 C C SOLUTION OF CUBIC EQUATION FOLLOWING THE ALGORITHM C FROM NUMERICAL RECIPES C Q=(R*R-3.*S)/9. P=(2.*R*R*R-9.*R*S+27.*T)/54. A1=Q*Q*Q-P*P PI=3.14159265 IF(A1.LE.0.00001)THEN TH=PI END IF IF(A1.GT.0.00001)THEN TH=Q*Q*Q TH=P/DSQRT(TH) IF(TH.LE.(-1.))TH=-1. IF(TH.GE.(1.))TH=1. TH=DACOS(TH) END IF A1=-2.*DSQRT(Q) A2=-R/3. A3=TH/3. V(1)=A1*DCOS(A3)+A2 A3=(TH+2.*PI)/3. V(2)=A1*DCOS(A3)+A2 A3=(TH+4.*PI)/3. V(3)=A1*DCOS(A3)+A2 X=V(1) AUX1=X**3+R*X**2+S*X+T X=V(2) AUX2=X**3+R*X**2+S*X+T X=V(3) AUX3=X**3+R*X**2+S*X+T IF(IPRINT.GT.2)WRITE(LOU,'(A/1X,3E15.6)')'INIT: PRECISION OF SOLUT *IONS OF CUBIC EQUATION ', AUX1,AUX2,AUX3 RETURN END C C ********************************************************* C SUBROUTINE INV(A,AINV,DET) C C ROUTINE INV DETERMINESS AN ADJOINT MATRIX OF A REAL 3*3 MATRIX C DIMENSION A(3,3),AINV(3,3) 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 /ZERO/ RNULL C DET=A(1,1)*(A(2,2)*A(3,3)-A(2,3)*A(3,2))-A(2,1)*(A(1,2)*A(3,3) 1-A(3,2)*A(1,3))+A(3,1)*(A(1,2)*A(2,3)-A(2,2)*A(1,3)) INULL=0 IF(ABS(DET).GT.RNULL) GOTO 100 DO 10 K=1,3 DO 10 J=1,3 IF(ABS(A(J,K)).GT.RNULL) INULL=1 AINV(J,K)=0.0 10 CONTINUE IF(INULL.EQ.0) THEN IF(IPRINT.GT.2)WRITE(LOUT,1000) A ELSE IF(IPRINT.GT.2)WRITE(LOUT,1010) END IF 100 AINV(1,1)=(A(2,2)*A(3,3)-A(3,2)*A(2,3)) AINV(2,1)=-(A(2,1)*A(3,3)-A(2,3)*A(3,1)) AINV(3,1)=(A(2,1)*A(3,2)-A(2,2)*A(3,1)) AINV(1,2)=-(A(1,2)*A(3,3)-A(1,3)*A(3,2)) AINV(2,2)=(A(1,1)*A(3,3)-A(1,3)*A(3,1)) AINV(3,2)=-(A(1,1)*A(3,2)-A(1,2)*A(3,1)) AINV(1,3)=(A(1,2)*A(2,3)-A(1,3)*A(2,2)) AINV(2,3)=-(A(1,1)*A(2,3)-A(1,3)*A(2,1)) AINV(3,3)=(A(1,1)*A(2,2)-A(1,2)*A(2,1)) IF(IPRINT.LE.2)RETURN WRITE(LOUT,1030)A WRITE(LOUT,1031)AINV WRITE(LOUT,1032)DET DO 30 K=1,3 E1=0.0 E2=0.0 E3=0.0 DO 20 J=1,3 E1=E1+A(J,K)*AINV(1,J) E2=E2+A(J,K)*AINV(2,J) E3=E3+A(J,K)*AINV(3,J) 20 CONTINUE WRITE(LOUT,1033)E1,E2,E3 30 CONTINUE C 1000 FORMAT(1x,' INV: ALL ELEMENTS OF MATRIX A ZERO',/,3(3F12.6,/)) 1010 FORMAT(1x,' INV :DETERMINANT OF MATRIX A ZERO') 1030 FORMAT(1X,' INV: A'/3(3F16.8/)) 1031 FORMAT(1X,' INV: AINV'/3(3F16.8/)) 1032 FORMAT(1X,' INV: DET ',F16.8) 1033 FORMAT(1X,' INV: ',3F16.8) RETURN END C C ********************************************************* C SUBROUTINE OUT(X,Y,DERY,IHLF,NDIM,PRMT) C DIMENSION Y(18),DERY(18),UN(3),POLD(3),PNEW(3),YA(3),YB(3),DEP(6), *PRMT(5),DN(3,3),DNG(3,2),XK(3),XKG(3),YK(3),YKG(3),TM(2), *XDAUX(3),YDAUX(3) COMMON /AUXI/ IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,MORI INTEGER CODE COMMON /COD/ CODE(50,2),KREF,KC,ITYPE COMMON /DENS/ RHO(20) COMMON /INTRF/ Z(1000),SX(350),SY(350),NX(20),NY(20),BRD(6),NINT, 1 XINTA COMPLEX PS COMMON /RAY/ AY(28,2000),DS(20,50),KINT(50),HHH(3,3),tmax, 1 PS(3,7,50),IS(8,50),N,IREF,IND,IND1 COMMON /RAY2/ DRY(3,2000) COMMON /ZERO/ RNULL COMMON /VSP/XVSP,YVSP,XNRM,YNRM,ICOD,IVSP COMMON /DYN/XDYN(3,3),ydyn(3,3) common/appr/ xold,xnew,yold(18),dold(18),ynew(18),dnew(18) COMMON /KM/KMAH,SPROLD,TSTOLD C N=N+1 NTR=20 IF(IND.EQ.10)GO TO 25 IF(N.GE.2000) GOTO 152 IF(X.GE.TMAX)GO TO 155 C C TEST OF SATISFACTION OF THE EIKONAL EQUATION AND OF C THE DERIVATIVES OF THE EIKONAL EQUATION C IF(IPREC.EQ.1.AND.NDER.GT.1)THEN TEST0=DERY(1)*Y(4)+DERY(2)*Y(5)+DERY(3)*Y(6) IF(ABS(TEST0).GT.RNULL)THEN DO 80 I=4,6 Y(I)=Y(I)/TEST0 80 CONTINUE END IF V2=1./(Y(4)*Y(4)+Y(5)*Y(5)+Y(6)*Y(6)) TEST3=Y(7)*Y(4)+Y(8)*Y(5)+Y(9)*Y(6) IF(ABS(TEST3).GT.RNULL)THEN DO 81 I=4,6 Y(I+3)=Y(I+3)-V2*TEST3*Y(I) 81 CONTINUE END IF ALP1=DERY(4)*Y(7)+DERY(5)*Y(8)+DERY(6)*Y(9) ALP2=DERY(1)*Y(13)+DERY(2)*Y(14)+DERY(3)*Y(15) IF(ABS(ALP1-ALP2).GT.RNULL)THEN IF(ABS(ALP1).GT..000001)THEN ALP=ALP2/ALP1 DO 82 I=7,9 Y(I)=Y(I)*ALP 82 CONTINUE ELSE DO 86 I=4,6 DERY(I+3)=DERY(I+3)-V2*ALP2*Y(I) 86 CONTINUE END IF END IF TEST4=Y(10)*Y(4)+Y(11)*Y(5)+Y(12)*Y(6) IF(ABS(TEST4).GT.RNULL)THEN DO 83 I=4,6 Y(I+6)=Y(I+6)-V2*TEST4*Y(I) 83 CONTINUE END IF ALP1=DERY(4)*Y(10)+DERY(5)*Y(11)+DERY(6)*Y(12) ALP2=DERY(1)*Y(16)+DERY(2)*Y(17)+DERY(3)*Y(18) IF(ABS(ALP1-ALP2).GT.RNULL)THEN IF(ABS(ALP1).GT..000001)THEN ALP=ALP2/ALP1 DO 84 I=10,12 Y(I)=Y(I)*ALP 84 CONTINUE ELSE DO 87 I=4,6 DERY(I+6)=DERY(I+6)-V2*ALP2*Y(I) 87 CONTINUE END IF END IF END IF IF(IPREC.EQ.2.AND.NDER.GT.1)THEN TEST0=DERY(1)*Y(4)+DERY(2)*Y(5)+DERY(3)*Y(6) TEST1=DERY(1)*Y(13)+DERY(2)*Y(14)+DERY(3)*Y(15) 1 -DERY(4)*Y(7)-DERY(5)*Y(8)-DERY(6)*Y(9) TEST2=DERY(1)*Y(16)+DERY(2)*Y(17)+DERY(3)*Y(18) 1 -DERY(4)*Y(10)-DERY(5)*Y(11)-DERY(6)*Y(12) TEST3=Y(7)*Y(4)+Y(8)*Y(5)+Y(9)*Y(6) TEST4=Y(10)*Y(4)+Y(11)*Y(5)+Y(12)*Y(6) IF(ABS(TEST0-1.).GT.RNULL)WRITE(LOU,100)TEST0 IF(ABS(TEST1).GT.RNULL.OR.ABS(TEST2).GT.RNULL) 1 WRITE(LOU,101)TEST1,TEST2 IF(ABS(TEST3).GT.RNULL.OR.ABS(TEST4).GT.RNULL) 1 WRITE(LOU,102)TEST3,TEST4 100 FORMAT(1X,'P_I*V_I=',E15.5) 101 FORMAT(1X,'DER.EIKON.EQ.=',2E15.5) 102 FORMAT(1X,'P_I*Q_{IJ}=',2E15.5) END IF C C CHECK THE POSITION OF THE RAY WITH RESPECT TO BOUNDARIES OF THE MODEL C NTR=21 IF(N.EQ.1)GO TO 13 IF(IREF.GT.1.AND.(N-1).EQ.KINT(IREF-1)) GOTO 13 NTR=22 IF(Y(1).LT.BRD(1).OR.Y(1).GT.BRD(2).OR.Y(2).LT.BRD(3).OR.Y(2). 1GT.BRD(4)) GOTO 20 C C CHECK WETHER THE RAY CROSSED THE BOREHOLE C IF(ICOD.EQ.0)GO TO 4 IF(ICOD.GT.IREF)GO TO 4 FXYZA=(AY(2,N-1)-XVSP)*XNRM+(AY(3,N-1)-YVSP)*YNRM FXYZB=(Y(1)-XVSP)*XNRM+(Y(2)-YVSP)*YNRM IF(FXYZA*FXYZB.GT.0.)GO TO 4 C C THE RAY CROSSED THE VERTICAL PLANE PERPENDICULAR TO THE LINE C SOURCE-BOREHOLE AND CONTAINING THE BOREHOLE C DETERMINATION OF THE POINT OF INTERSECTION C XNEW=X KDIM=6 IF(NDER.GT.1)KDIM=18 DO 41 I=1,KDIM YNEW(I)=Y(I) DNEW(I)=DERY(I) 41 CONTINUE XA=XOLD XB=XNEW C C DETERMINATION OF THE POINT OF INTERSECTION WITH THE VERTICAL PLANE C PERPENDICULAR TO THE LINE SOURCE-BOREHOLE C CALL ROOT(XA,FXYZA,XB,FXYZB,X,PRMT,-1) IF(IAC.GE.10)GO TO 153 CALL APPROX(X,Y,DERY,KDIM) IND=43 4 CONTINUE C C CHECK WETHER THE RAY CROSSED AN INTERFACE C NTR=23 INTR=LAY+1 CALL DISC(Y,DEP) ZPL=DEP(1) INTR=LAY CALL DISC(Y,DEP) ZPU=DEP(1) NTR=24 IF(ZPL.LE.ZPU) THEN C C TWO NEIGHBOURHOOD INTERFACES CROSS EACH OTHER C WRITE(LOU,'(A)') ' MAUX, LAY, IND,IND1,X,Y, Z,ZU,ZL' WRITE(LOU,'(4I5,5F12.5,/)')MAUX,LAY,IND,IND1,Y(1),Y(2),Y(3),ZPU, 1ZPL GOTO 150 END IF NTR=25 IS(3,IREF)=LAY IF(Y(3).LT.ZPU.OR.Y(3).GT.ZPL) GOTO 30 NTR=27 C C THE RAY DID NOT CROSS AN INTERFACE C C IF(IND.EQ.1.OR.IND.EQ.2)GO TO 148 IF(IND.EQ.30.OR.IND.EQ.31)GO TO 148 13 CONTINUE AY(1,N)=X XOLD=X DO 15 I=1,6 AY(I+1,N)=Y(I) YOLD(I)=Y(I) DOLD(I)=DERY(I) 15 CONTINUE DRY(1,N)=DERY(1) DRY(2,N)=DERY(2) DRY(3,N)=DERY(3) IF(NDER.GT.1)THEN C C DETERMINATION OF KMAH INDEX C IF(N.GT.1)THEN SPR1=Y(8)*Y(12)-Y(9)*Y(11) SPR2=Y(9)*Y(10)-Y(7)*Y(12) SPR3=Y(7)*Y(11)-Y(8)*Y(10) SPR=SPR1*DERY(1)+SPR2*DERY(2)+SPR3*DERY(3) IF(SPR*SPROLD.LT.0..AND.(X-TSTOLD).GT.0.)KMAH=KMAH+1 SPROLD=SPR TSTOLD=X END IF DO 14 I=1,12 YOLD(i+6)=y(i+6) DOLD(I+6)=DERY(I+6) 14 CONTINUE END IF NTR=99 IF(IND.EQ.10)GO TO 25 IF(IND.EQ.43)GO TO 25 C C WRITING IN AY FIELD C CALL PARDIS (Y,1) RETURN C C THE RAY CROSSED ONE OF THE VERTICAL BOUNDARIES OF THE MODEL C 20 IF(Y(1).LT.BRD(1)) IND=1 IF(Y(1).GT.BRD(2)) IND=2 IF(Y(2).LT.BRD(3)) IND=3 IF(Y(2).GT.BRD(4)) IND=4 IF(IND.EQ.3.OR.IND.EQ.4) GOTO 23 AUX=(BRD(IND)-AY(2,N-1))/(Y(1)-AY(2,N-1)) TR=X X=AY(1,N-1)+AUX*(X-AY(1,N-1)) K=6 IF(NDER.GT.1)K=18 DO 21 I=1,K 21 Y(I)=Y(I)+DERY(I)*(X-TR) Y(1)=BRD(IND) AY(1,N)=X DO 3 I=1,6 3 AY(I+1,N)=Y(I) GO TO 4 23 AUX=(BRD(IND)-AY(3,N-1))/(Y(2)-AY(3,N-1)) TR=X X=AY(1,N-1)+AUX*(X-AY(1,N-1)) K=6 IF(NDER.GT.1)K=18 DO 16 I=1,K 16 Y(I)=Y(I)+DERY(I)*(X-TR) Y(2)=BRD(IND) AY(1,N)=X DO 17 I=1,6 17 AY(I+1,N)=Y(I) IF(IND.EQ.3) IND=30 IF(IND.EQ.4) IND=31 GO TO 4 C C THE RAY CROSSED AN INTERFACE C 30 CONTINUE IF(IREF.LE.1.OR.KC.GT.0) GOTO 39 IF(N-KINT(IREF-1).EQ.2) THEN IND=9 IND1=LAY GO TO 25 END IF C C WHICH INTERFACE WAS CROSSED C 39 INTR=LAY IF(Y(3).GE.ZPL) INTR=LAY+1 XNEW=X KDIM=6 IF(NDER.GT.1)KDIM=18 DO 40 I=1,KDIM YNEW(I)=Y(I) DNEW(I)=DERY(I) 40 CONTINUE DO 35 I=1,3 YA(I)=YOLD(I) YB(I)=YNEW(I) 35 CONTINUE CALL DISC(YA,DEP) FXYZA=DEP(1)-YA(3) CALL DISC(YB,DEP) FXYZB=DEP(1)-YB(3) C C DETERMINATION OF THE POINT OF INTERSECTION WITH THE INTERFACE C XA=XOLD XB=XNEW CALL ROOT(XA,FXYZA,XB,FXYZB,XINT,PRMT,1) IF(IAC.GE.10)GO TO 153 NTR=10 C C THE QUANTITIES OF RAY TRACING AND DYNAMIC RAY TRACING C AT THE POINT OF INCIDENCE WILL BE DETERMINED C IND1=INTR CALL APPROX(XINT,Y,DERY,KDIM) AY(1,N)=XINT DO 63 I=1,6 AY(I+1,N)=y(I) 63 CONTINUE CALL PARDIS(Y,1) CALL FCT(X,Y,DERY) DO 64 I=1,3 DRY(I,N)=DERY(I) 64 CONTINUE IF(NDER.GT.1)THEN IF(N.GT.1)THEN SPR1=Y(8)*Y(12)-Y(9)*Y(11) SPR2=Y(9)*Y(10)-Y(7)*Y(12) SPR3=Y(7)*Y(11)-Y(8)*Y(10) SPR=SPR1*DERY(1)+SPR2*DERY(2)+SPR3*DERY(3) IF(SPR*SPROLD.LT.0..AND.(X-TSTOLD).GT.0.)KMAH=KMAH+1 SPROLD=SPR TSTOLD=X END IF DO 62 I=1,3 XDYN(I,1)=Y(I+6) XDYN(I,2)=Y(I+9) XDYN(I,3)=DERY(I) YDYN(I,1)=Y(I+12) YDYN(I,2)=Y(I+15) YDYN(I,3)=DERY(I+3) XK(I)=DERY(I) YK(I)=DERY(I+3) 62 CONTINUE END IF C C DETERMINATION OF THE UNIT NORMAL TO THE INTERFACE AT THE POINT OF C INCIDENCE C CALL DISC (Y,DEP) DFX=DEP(2) DFY=DEP(3) ROO=1+DFX*DFX+DFY*DFY UN3=ROO**(-0.5) UN(3)=-UN3 UN(2)=UN3*DFY UN(1)=UN3*DFX C C SCALAR PRODUCT OF GROUP VELOCITY WITH UNIT NORMAL VECTOR C PN=DERY(1)*UN(1)+DERY(2)*UN(2)+DERY(3)*UN(3) C C ORIENTATION OF THE UNIT NORMAL SO THAT IT POINTS TO THE MEDIUM IN C WHICH THE ENERGY OF INCIDENT WAVE PROPAGATES C IF(PN.LT.0) GOTO 120 UN(1)=-UN(1) UN(2)=-UN(2) UN(3)=-UN(3) UN3=-UN3 120 CONTINUE IF(NDER.GT.1)THEN XIINC=Y(4)*UN(1)+Y(5)*UN(2)+Y(6)*UN(3) UN33=UN3*UN3*UN3 DFXX=DEP(4) DFXY=DEP(5) DFYY=DEP(6) TM(1)=0. TM(2)=0. DO 126 J=1,2 AUX=0. DO 125 I=1,3 TM(J)=TM(J)+UN(I)*XDYN(I,J) AUX=AUX+UN(I)*XK(I) 125 CONTINUE TM(J)=-TM(J)/AUX 126 CONTINUE C C DETERMINATION OF DERIVATIVES OF THE UNIT NORMAL TO INTERFACE C AUX1=1.+DFX*DFX AUX2=1.+DFY*DFY AUX3=DFX*DFY DN(1,1)=(DFXX*AUX2-AUX3*DFXY)*UN33 DN(1,2)=(DFXY*AUX2-AUX3*DFYY)*UN33 DN(1,3)=0. DN(2,1)=(DFXY*AUX1-AUX3*DFXX)*UN33 DN(2,2)=(DFYY*AUX1-AUX3*DFXY)*UN33 DN(2,3)=0. DN(3,1)=(DFX*DFXX+DFY*DFXY)*UN33 DN(3,2)=(DFX*DFXY+DFY*DFYY)*UN33 DN(3,3)=0. DO 129 L=1,3 DO 128 M=1,2 DNG(L,M)=0. DO 127 K=1,3 DNG(L,M)=DNG(L,M)+DN(L,K)*(XDYN(K,M)+XK(K)*TM(M)) 127 CONTINUE 128 CONTINUE 129 CONTINUE END IF NTR=20 IF(KRE.EQ.1) GOTO 24 NTR=30 IF(KRE.EQ.0) GOTO 130 C C MULTIPLE REFLECTED WAVE C IF((IREF+1).GT.KRE.AND.INTR.EQ.INT1) GOTO 151 C C TERMINATION OF THE RAY AT AN INNER INTERFACE OR AT THE FREE SURFACE C OR AT THE BOTTOM OF THE MODEL C IF((IREF+1).GT.KRE) IRR=IREF NTR=35 IF((IREF+1).GT.KRE) GOTO 148 C C NCD : INDICATOR WHETHER A REFLECTION OR TRANSMISSION TAKES PLACE C WITH RESPECT TO THE CODE OF THE WAVE C NCD=CODE(IREF+1,1)-CODE(IREF,1) C C NCD1 : INDICATOR FOR THE TYPE OF THE WAVE IN THE NEW LAYER C NCD1=CODE(IREF+1,2)-CODE(IREF,2) NTR=40 IF(KC.GT.0.AND.IREF.EQ.1.AND.INTR.EQ.LAY) GOTO 151 NTR=50 IF(KC.LT.0.AND.IREF.EQ.1.AND.INTR.NE.LAY) GOTO 151 NTR=60 IF(IREF.EQ.1) GOTO 170 NTR=70 C IF(INTR.EQ.INT1.AND.NCD.NE.0.OR.NCD1.NE.0) GOTO 151 IF(INTR.EQ.INT1.AND.NCD.NE.0)GOTO 151 NTR=75 IF(INTR.EQ.INT1.AND.NCD1.NE.0)GOTO 151 C C REFLECTION OR TRANSMISSION OF THE RAY C NTR=80 IF(INTR.NE.INT1.OR.NCD.NE.0.OR.NCD1.NE.0) GOTO 170 C C COMPOUND ELEMENT OF THE RAY C IREFR=1 KINT(IREF)=0 IREF=IREF+1 NCD=CODE(IREF+1,1)-CODE(IREF,1) NCD1=CODE(IREF+1,2)-CODE(IREF,2) C C TERMINATION AT AN INNER INTERFACE C IRR=IREF NTR=90 IF((IREF+1).GT.KRE) GOTO 24 NTR=100 GOTO 170 C C REFRACTED WAVE C C ORDINARY TERMINATION AT UPPER BOUNDARY C 130 NTR=110 IF(INTR.EQ.LAY.AND.LAY.EQ.1) GOTO 148 NCD=1 NCD1=0 C C SPECIFICATION OF THE LAYER OF THE GENERATED WAVE C 170 INT1=INTR IF(IRHO.NE.0)DS(8,IREF)=RHO(LAY) IRR=IREF IREF=IREF+1 NTR=120 IF(NCD.EQ.0) GOTO 200 NTR=130 IF(INTR.EQ.LAY) GOTO 190 C C TRANSMISSION AT THE LOWER INTERFACE C NTR=140 IF(NCD.LT.0) GOTO 151 C C ORDINARY TERMINATION AT LOWER BOUNDARY C NTR=150 IF(KRE.LE.1.AND.INTR.EQ.NINT) GOTO 148 NTR=160 IF(INTR.EQ.NINT) GOTO 151 LAY=LAY+1 NTR=170 GOTO 200 C C TRANSMISSION AT THE UPPER INTERFACE C 190 NTR=180 IF(NCD.GT.0.AND.KRE.GT.1) GOTO 151 NTR=190 IF(KRE.LE.1.AND.LAY.EQ.1) GOTO 24 NTR=200 IF(LAY.EQ.1) GOTO 151 LAY=LAY-1 C C TRANSFORMATION OF THE SLOWNESS VECTOR C 200 IF(INTR.EQ.NINT.AND.MDIM.NE.0) GOTO 154 DO 210 I=1,3 POLD(I)=Y(3+I) 210 CONTINUE ITRANS=0 C C CHECK WHETHER A REFLECTION OR TRANSMISSION TAKES PLACE C IF(KC.EQ.0) GOTO 217 IF(CODE(IREF,1).EQ.CODE(IREF-1,1)) GOTO 218 217 ITRANS=1 218 CALL TRANSL (Y,POLD,PNEW,UN,ITRANS,1) IF(IND.EQ.9.OR.IND.EQ.10) GOTO 25 DO 220 I=1,3 Y(I+3)=PNEW(I) 220 CONTINUE IF(NDER.GT.1)THEN CALL FCT(X,Y,DERY) DO 225 I=1,3 XKG(I)=DERY(I) YKG(I)=DERY(I+3) 225 CONTINUE XI=PNEW(1)*UN(1)+PNEW(2)*UN(2)+PNEW(3)*UN(3) DXTN=XKG(1)*UN(1)+XKG(2)*UN(2)+XKG(3)*UN(3) DO 228 M=1,2 AUX=XKG(1)*DNG(1,M)+XKG(2)*DNG(2,M)+XKG(3)*DNG(3,M) DO 224 I=1,3 XDAUX(I)=XDYN(I,M) YDAUX(I)=YDYN(I,M) 224 CONTINUE DO 227 I=1,3 AUXX=(DNG(I,M)-UN(I)*AUX/DXTN)*(XI-XIINC) DO 226 K=1,3 AUXX=AUXX-UN(I)*(XKG(K)*(YDaux(K)+YK(K)*TM(M))- 1 YKG(K)*(XDaux(K)+XK(K)*TM(M)))/DXTN 226 CONTINUE XDYN(I,M)=XDYN(I,M)+(XK(I)-XKG(I))*TM(M) YDYN(I,M)=YDYN(I,M)+(YK(I)-YKG(I))*TM(M)+AUXX 227 CONTINUE 228 CONTINUE DO 229 I=1,3 Y(I+6)=XDYN(I,1) Y(I+9)=XDYN(I,2) Y(I+12)=YDYN(I,1) Y(I+15)=YDYN(I,2) 229 CONTINUE END IF PRMT(5)=2. RETURN C C ORDINARY TERMINATION OF THE RAY C 148 IF(IND.EQ.1.OR.IND.EQ.2.OR.IND.EQ.30.OR.IND.EQ.31) GOTO 25 24 IND=INTR+100 IF((MREG.EQ.0.OR.MREG.GE.1).AND.MDIM.NE.0) THEN POLD(1)=Y(4) POLD(2)=Y(5) POLD(3)=Y(6) IS(3,IREF)=LAY IREF=IREF+1 CALL TRANSL(Y,POLD,PNEW,UN,0,1) IREF=IREF-1 END IF IF(KRE.LE.1) IRR=IREF IF(INTR.EQ.1) IND=3 IF(INTR.EQ.NINT) IND=4 KINT(IRR)=N IF(IND.NE.3) GOTO 25 C C COMPUTATION OF REFLECTED WAVES AT THE EARTH SURFACE FOR COEFFICIENTS C OF CONVERSION C GOTO 25 150 IND=20 GOTO 25 151 IND=8 GOTO 25 152 IND=13 GOTO 25 153 IND=39 GOTO 25 154 IND=15 GOTO 25 155 IND=12 25 PRMT(5)=3.0 IF(IND.NE.3.AND.IND.NE.43) RETURN ITYPE=CODE(IREF,2) CALL PARDIS(Y,2) RETURN END C C ********************************************************* C SUBROUTINE POLAR (N1,N2,NN,I) C C ROUTINE COMPUTING POLARIZATION VECTORS C C IF(IPOL.EQ.0) POLARIZATION VECTORS ARE NOT PRINTED C IF(IPOL.EQ.1) POLARIZATION VECTORS ARE PRINTED ONLY C FOR THE POINTS OF INCIDENCE OF A RAY AT INTERFACES C IF(IPOL.EQ.2) POLARIZATION VECTORS ARE PRINTED FOR ALL C COMPUTED POINTS ALONG THE RAY C N1,N2 - SUCCESSIVE NUMBERS OF THE FIRST AND THE LAST POINT C OF AN ELEMENT OF A RAY (IN CASE IPOL.EQ.0, N1=N2 AND THIS C NUMBER CORRESPONDS TO THE POINT OF INCIDENCE OF A RAY AT C AN INTERFACE; NEGATIVE N2 - INDICATION OF COMPOSED ELEMENT C OF THE RAY C NN - TOTAL NUMBER OF POINTS ALONG THE RAY C I - NUMBER OF THE ELEMENT OF THE RAY C C CALLED FROM ROUTINES: AMPL C ROUTINES CALLED: CHRM C DIMENSION PXN(3,3),Y(6),B1(3),B2(3),XX(3),E(21),E1(3),E2(3),T(3) COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33, 1 A34,A35,A36,A44,A45,A46,A55,A56,A66, 1 DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23, 1 DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44, 1 DXA45,DXA46,DXA55,DXA56,DXA66, 1 DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23, 1 DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44, 1 DYA45,DYA46,DYA55,DYA56,DYA66, 1 DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23, 1 DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44, 1 DZA45,DZA46,DZA55,DZA56,DZA66, 1 A2546,A1266,A1355,A1456,A3645,A2344 COMMON /AUXI/ IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,MORI INTEGER CODE COMMON /COD/ CODE(50,2),KREF,KC,ITYPE COMMON /DJK/ D11,D12,D13,D22,D23,D33,DTR COMMON /GAM/ C11,C12,C13,C22,C23,C33 COMPLEX PS COMMON /RAY/ AY(28,2000),DS(20,50),KINT(50),HHH(3,3),TMAX, 1 PS(3,7,50),IS(8,50),N,IREF,IND,IND1 COMMON /ZERO/ RNULL C EQUIVALENCE(E(1),A11),(E(2),A12),(E(3),A13),(E(4),A14),(E(5),A15) 1 ,(E(6),A16),(E(7),A22),(E(8),A23),(E(9),A24) 2 ,(E(10),A25),(E(11),A26),(E(12),A33),(E(13),A34),(E(14),A35) 2 ,(E(15),A36),(E(16),A44),(E(17),A45),(E(18),A46),(E(19),A55) 2 ,(E(20),A56),(E(21),A66) C IF(N2.LT.0)THEN LAY=IS(3,I-1) N2=-N2 ELSE LAY=IS(3,I) END IF IF(IPOL.GT.0) 1WRITE(LOU,'(A,6I5)')' LAY,I,N1,N2,ITYPE,IANI',LAY,I,N1,N2,ITYPE, 2IANI(LAY) IF(IANI(LAY).NE.0) GO TO 2000 C C ISOTROPIC CASE C IF(NN.EQ.N2) GO TO 1005 IF(ITYPE.NE.3) GO TO 1020 C C P-WAVE POLARIZATION VECTOR C 1005 DO 1010 J=N1,N2 V=1./SQRT(AY(5,J)*AY(5,J)+AY(6,J)*AY(6,J)+AY(7,J)*AY(7,J)) T(1)=V*AY(5,J) T(2)=V*AY(6,J) T(3)=V*AY(7,J) IF((IPOL.EQ.1.AND.J.EQ.N2).OR.IPOL.EQ.2) 1WRITE(LOU,'(A,/,6F10.5)') ' X1,X2,X3,TX,TY,TZ', 2AY(2,J),AY(3,J),AY(4,J),T 1010 CONTINUE DO 1015 K=1,3 HHH(1,K)=0. HHH(2,K)=0. HHH(3,K)=T(K) 1015 CONTINUE IF(N2.EQ.NN) GO TO 1020 RETURN C C S-WAVE POLARIZATION VECTORS C 1020 CONTINUE DO 1100 J=N1,N2 V=SQRT(AY(12,J)) V2=2.*V M1=5 M2=6 M3=7 IF(ABS(V2).GT..000001)THEN DVS1=AY(13,J)/V2 DVS2=AY(14,J)/V2 ELSE DVS1=0. DVS2=0. END IF 1030 P1=AY(M1,J) P2=AY(M2,J) P3=AY(M3,J) XL=SQRT(P1*P1+P2*P2) IF(XL.LT.RNULL)THEN IF(M1.EQ.4)THEN M1=6 M2=7 M3=5 IF(ABS(V2).GT..000001)THEN DVS1=AY(14,J)/V2 DVS2=AY(15,J)/V2 ELSE DVS1=0. DVS2=0. END IF GO TO 1030 END IF IF(M1.EQ.5)THEN M1=7 M2=5 M3=6 IF(ABS(V2).GT..000001)THEN DVS1=AY(15,J)/V2 DVS2=AY(13,J)/V2 ELSE DVS1=0. DVS2=0. END IF GO TO 1030 END IF END IF C C COMPUTATION OF ANGLE XI SPECIFYING POLARIZATION VECTORS C IF(IPOL.EQ.2)WRITE(LOU,'(A,/,6F12.6)') 1' P1,P2,P3,XL,D1,D2',P1,P2,P3,XL,DVS1,DVS2 IF(J.EQ.N1)THEN X1=0. TIME=0. XO=(P3*(DVS1*P2-DVS2*P1))/(XL*XL) XI=0. GO TO 1040 END IF TIME=AY(1,J)-AY(1,J-1) X1=(P3*(DVS1*P2-DVS2*P1))/(XL*XL) XI=.5*(X1+XO)*TIME+XI 1040 IF(IPOL.EQ.2)WRITE(LOU,'(A,/,6F12.6)') 1' XI0,XI1,TSTEP,XISTEP',XO,X1,TIME,XI XO=X1 P1=AY(5,J) P2=AY(6,J) P3=AY(7,J) XL=SQRT(P1*P1+P2*P2) IF(XL.LT..0001)THEN B1(1)=1. B1(2)=0. B1(3)=0. B2(1)=0. B2(2)=1. IF(P3.LT.0.)B2(2)=-1. B2(3)=0. ELSE U=SQRT(P1*P1+P2*P2+P3*P3) P13=P1*P3/XL P23=P2*P3/XL B1(1)=P13/U B1(2)=P23/U B1(3)=-xl/U B2(1)=-P2/XL B2(2)=P1/XL B2(3)=0. END IF SN=SIN(XI) CS=COS(XI) DO 1050 K=1,3 E1(K)=CS*B1(K)-SN*B2(K) E2(K)=SN*B1(K)+CS*B2(K) 1050 CONTINUE AY(17,J)=E1(1) AY(18,J)=E1(2) AY(19,J)=E1(3) AY(20,J)=E2(1) AY(21,J)=E2(2) AY(22,J)=E2(3) IF((IPOL.EQ.1.AND.J.EQ.N2).OR.IPOL.EQ.2) 1WRITE(LOU,'(A,/,4F10.5/6F10.5)') 2' TIME,X1,X2,X3/ E1X,E1Y,E1Z,E2X,E2Y,E2Z',(AY(K,J),K=1,4),E1,E2 1100 CONTINUE DO 1110 K=1,3 HHH(1,K)=E1(K) HHH(2,K)=E2(K) IF(N2.EQ.NN) GO TO 1110 HHH(3,K)=0. 1110 CONTINUE RETURN C C ANISOTROPIC CASE C C C COMPUTATION OF POLARIZATION VECTORS C 2000 DO 2300 J=N1,N2 DO 2110 K=1,21 E(K)=AY(K+7,J) 2110 CONTINUE DO 2120 K=1,6 Y(K)=AY(K+1,J) 2120 CONTINUE A2546=A25+A46 A1266=A12+A66 A1355=A13+A55 A1456=A14+A56 A3645=A36+A45 A2344=A23+A44 CALL CHRM(Y) 2125 IF(ABS(DTR).LT..000001)GO TO 2135 PXN(1,1)=D11 PXN(1,2)=D12 PXN(1,3)=D13 XX(1)=D11*D11+D12*D12+D13*D13 PXN(2,1)=D12 PXN(2,2)=D22 PXN(2,3)=D23 XX(2)=D12*D12+D22*D22+D23*D23 PXN(3,1)=D13 PXN(3,2)=D23 PXN(3,3)=D33 XX(3)=D13*D13+D23*D23+D33*D33 XN=0. DO 2130 L=1,3 IF(XN.GE.XX(L))GO TO 2130 XN=XX(L) NX=L 2130 CONTINUE XN=SQRT(XN) IF(XN.GT.RNULL)GO TO 2140 2135 CONTINUE WRITE(LOU,'(A,A,5I5)') 1' XN.LT.RNULL IN POLAR - SHEAR-WAVE SINGULARITY ???', 2' LAY,N,N1,N2',LAY,J,J,N1,N2 IND=10 RETURN 2140 P1=PXN(NX,1)/XN P2=PXN(NX,2)/XN P3=PXN(NX,3)/XN IF(J.NE.N1)THEN XN=P1OLD*P1+P2OLD*P2+P3OLD*P3 IF(XN.LT.0.)THEN P1=-P1 P2=-P2 P3=-P3 ENDIF ENDIF P1OLD=P1 P2OLD=P2 P3OLD=P3 C C CHECK OF PRECISION C Z01=(C11-1.)*P1+C12*P2+C13*P3 Z02=C12*P1+(C22-1.)*P2+C23*P3 Z03=C13*P1+C23*P2+(C33-1.)*P3 Z01=ABS(Z01) Z02=ABS(Z02) Z03=ABS(Z03) IF(Z01.GT.0.01.OR.Z02.GT.0.01.OR.Z03.GT.0.01) THEN WRITE(LOU,'(A,4I5)') ' CHRISTOFFEL EQUATION IS SATISFIED 1 WITH PRECISION LESS THAN 0.01: LAY,ITYPE,RAY ELEMENT,NPOINT ', 2 LAY,ITYPE,I,J END IF IF((IPOL.EQ.1.AND.J.EQ.N2).OR.IPOL.EQ.2) 1WRITE(LOU,'(A,/,6F10.5)') 2' X1,X2,X3,GX,GY,GZ',AY(2,J),AY(3,J),AY(4,J),P1,P2,P3 IF(J.EQ.N2)THEN HHH(ITYPE,1)=P1 HHH(ITYPE,2)=P2 HHH(ITYPE,3)=P3 END IF 2300 CONTINUE RETURN END C C ********************************************************* C SUBROUTINE POLRT(XCOF,COF,M,ZERO,IER) DIMENSION XCOF(7),COF(7) COMPLEX ZERO(6) DOUBLE PRECISION XO,YO,X,Y,XPR,YPR,UX,UY,V,YT,XT,U,XT2,YT2,SUMSQ, / DX,DY,TEMP,ALPHA,XCOF,COF,DABS C IFIT=0 N=M IER=0 IF(XCOF(N+1)) 10,25,10 10 IF(N) 15,15,32 C 15 IER=1 20 RETURN C 25 IER=4 GO TO 20 C 30 IER=2 GO TO 20 32 IF(N-36) 35,35,30 35 NX=N NXX=N+1 N2=1 KJ1=N+1 DO 40 L=1,KJ1 MT=KJ1-L+1 40 COF(MT)=XCOF(L) C 45 XO=-.00500101D0 YO=0.01000101D0 C IN=0 50 X=XO C XO=-10.D0*YO YO=-10.D0*X C X=XO Y=YO IN=IN+1 GO TO 59 55 IFIT=1 XPR=X YPR=Y C 59 ICT=0 60 UX=0.D0 UY=0.D0 V =0.D0 YT=0.D0 XT=1.D0 U=COF(N+1) IF(U) 65,130,65 65 DO 70 I=1,N L=N-I+1 TEMP=COF(L) XT2=X*XT-Y*YT YT2=X*YT+Y*XT U=U+TEMP*XT2 V=V+TEMP*YT2 FI=I UX=UX+FI*XT*TEMP UY=UY-FI*YT*TEMP XT=XT2 70 YT=YT2 SUMSQ=UX*UX+UY*UY IF(SUMSQ) 75,110,75 75 DX=(V*UY-U*UX)/SUMSQ X=X+DX DY=-(U*UY+V*UX)/SUMSQ Y=Y+DY 78 IF(DABS(DY)+DABS(DX)-1.0D-12) 100,80,80 C 80 ICT=ICT+1 IF(ICT-500) 60,85,85 85 IF(IFIT) 100,90,100 90 IF(IN-5) 50,95,50 C 95 IER=3 GO TO 20 100 DO 105 L=1,NXX MT=KJ1-L+1 TEMP=XCOF(MT) XCOF(MT)=COF(L) 105 COF(L)=TEMP ITEMP=N N=NX NX=ITEMP IF(IFIT) 120,55,120 110 IF(IFIT) 115,50,115 115 X=XPR Y=YPR 120 IFIT=0 IF(DABS(Y)-1.0D-10*DABS(X))135,125,125 125 ALPHA=X+X SUMSQ=X*X+Y*Y N=N-2 GO TO 140 130 X=0.D0 NX=NX-1 NXX=NXX-1 135 Y=0.D0 SUMSQ=0.0 ALPHA=X N=N-1 140 COF(2)=COF(2)+ALPHA*COF(1) 145 DO 150 L=2,N 150 COF(L+1)=COF(L+1)+ALPHA*COF(L)-SUMSQ*COF(L-1) 155 ZERO(N2)=CMPLX(X,Y) N2=N2+1 IF(SUMSQ) 160,165,160 160 Y=-Y SUMSQ=0.D0 GO TO 155 165 IF(N) 20,20,45 END Ca4.for 0100666 0000765 0000765 00000122320 11023551576 011452 0 ustar bulant bulant C
SUBROUTINE PROFIL(XSOUR,YSOUR,ZSOUR,TSOUR,PSI0,PAZM,RANG, 1XXX,YYY,ZZZ,TTT,DT,AC,ASTART,ASTEP,AFIN,ITMAX,MOUT,NCODE, 2METHOD,ITPR,indr1) C C 3-D INITIAL VALUE RAY TRACING AND RAY TRACING FROM THE SOURCE C TO A PRESCRIBED PROFILE PASSING THROUGH THE EPICENTER C COMMON /AUXI/ IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori INTEGER CODE COMMON /COD/ CODE(50,2),KREF,KC,ITYPE COMMON /DIST/ XDST(200),NDSTX,XREPS,DST(2),NDST,REPS,LNDST, 1XPRF,YPRF,ILOC COMPLEX PS COMMON /RAY/ AY(28,2000),DS(20,50),KINT(50),HHH(3,3),tmax, 1 PS(3,7,50),IS(8,50),N,IREF,IND,IND1 COMMON /ZERO/ RNULL COMMON/VSP/XVSP,YVSP,XNRM,YNRM,ICOD,IVSP COMMON/DYN/XDYN(3,3),ydyn(3,3) COMMON/KM/KMAH,SPROLD,TSTOLD C iwave=0 itp=code(1,2) RANG=0. XXX=0. ZZZ=0. TTT=0. REPS1=.05*REPS DD=dst(1) xcos=cos(dd) xsin=sin(dd) dd=0. X=0. ITER=0 II=0 LNDST=0 C AA=ASTART-ASTEP INDEX=0 INUM=0 ICLS=0 ISUC=0 INDS=ISOUR C C LOOP IN AZIMUTH, FROM VALUE ASTART TO AFIN, WITH THE STEP C ASTEP C 1 AA=AA+ASTEP PNEW=AA IF(ASTEP.GT.0..AND.AA.GT.AFIN)GO TO 99 IF(ASTEP.LT.0..AND.AA.LT.AFIN)GO TO 99 IND=INDS NDER=1 IF(MDIM.GE.1)NDER=2 SPROLD=0. CALL RAYA(XSOUR,YSOUR,ZSOUR,TSOUR,PSI0,AA,PX,PY,PZ,XX,YY,ZZ,T, 1DT,AC) NDER=1 IF(IND.EQ.14)RETURN x=(yprf-yy)*xcos-(xprf-xx)*xsin IF(NDSTX.EQ.0)GO TO 65 IF(IND.EQ.ITPR)XAX=X IF(IND.EQ.ITPR)PNEW=AA IF(MOUT.EQ.3)WRITE(LOU,100)IND,IND1,X,XX,YY,ZZ,T,AA,PSI0 IF(INUM.NE.0)GO TO 2 AOLD=AA IOLD=IND IOLD1=IND1 XOLD=X TOLD=T INUM=1 GO TO 1 C C PARAMETERS FOR THE PRECEDING RAY: AA=AOLD, X=XOLD, IND=IOLD C PARAMETERS FOR THE NEW RAY: AA=ANEW, X=XNEW, IND=INEW C 2 INEW=IND INEW1=IND1 ANEW=AA XNEW=X TNEW=T IF(INEW.EQ.ITPR.AND.IOLD.EQ.ITPR)GO TO 21 IF(INEW.EQ.ITPR)GO TO 50 IF(IOLD.EQ.ITPR)GO TO 55 IF(INEW.EQ.9.AND.IOLD.NE.9.AND.IOLD.NE.8)GO TO 30 IF(INEW.NE.9.AND.INEW.NE.8.AND.IOLD.EQ.9)GO TO 35 GOTO 3 21 IF(IOLD1.NE.INEW1)then if(inew1.eq.indr1)go to 50 if(iold1.eq.indr1)go to 55 else GO TO 40 end if C C NO ITERATIONS, TAKE A NEW RAY IN THE LOOP C 3 CONTINUE if(isuc.eq.0)ind=0 IF(IOLD.NE.INEW)IND=0 IOLD=INEW IOLD1=INEW1 XOLD=XNEW AOLD=ANEW TOLD=TNEW GO TO 1 C C REGULAR CASE: IOLD=3, INEW=3 C 40 XXNEW=XNEW XXOLD=XOLD AANEW=ANEW AAOLD=AOLD TTNEW=TNEW TTOLD=TOLD IBNEW=0 IBOLD=0 41 IF(XXNEW.GT.XXOLD.AND.DST(2).GT.DST(1))GO TO 46 IF(XXNEW.LT.XXOLD.AND.DST(2).LT.DST(1))GO TO 46 C C REGULAR CASE: XXNEW.LE.XXOLD, ITREND=-1 (REVERSE BRANCH) C DX=XXOLD IF(IBOLD.EQ.1) DX=DX+REPS IF(DD.GE.DX) GO TO 3 DX=XXNEW IF(IBNEW.EQ.1) DX=DX-REPS IF(DD.LT.DX) GOTO 3 II=1 GO TO 43 C C REGULAR CASE: XXNEW.GT.XXOLD, ITREND=1 (REGULAR BRANCH) C 46 continue DX=XXOLD IF(IBOLD.EQ.1) DX=DX-REPS IF(DD.LE.DX) GO TO 3 DX=XXNEW IF(IBNEW.EQ.1) DX=DX+REPS IF(DD.GT.DX) GOTO 3 II=1 43 P1=AAOLD P2=AANEW X1=XXOLD X2=XXNEW T1=TTOLD T2=TTNEW C C I T E R A T I O N S C ITER=0 ISIGN=1 IPR1=0 IPR2=0 ISUC=0 GO TO 61 68 XAX=X PAX=PNEW 61 ITER=ITER+1 IF(ITER.GT.ITMAX)GO TO 80 ISIGN=-ISIGN AAUX=0.5*(P1+P2) IF(METHOD.LE.1.AND.IND.EQ.ITPR.and.iter.gt.1)GO TO 62 GO TO 69 62 if(mori.eq.0)AUX=(XDYN(1,1)*xsin-XDYN(2,1)*xcos)*cos(psi0) if(mori.ne.0)AUX=(XDYN(1,2)*xsin-XDYN(2,2)*xcos)*cos(pnew) IF(ABS(AUX).LT..00001)GO TO 69 AAUX=PNEW+(DD-X)/AUX 69 PNEW=AAUX 71 IND=INDS SPROLD=0. XOLD=0. IF(MDIM.GE.1)NDER=2 CALL RAYA(XSOUR,YSOUR,ZSOUR,TSOUR,PSI0,PNEW,PX,PY,PZ,XX,YY,ZZ,T, 1DT,AC) if(iwave.eq.1)code(1,2)=itp NDER=1 XE=XX-Xprf YE=YY-Yprf RPRF=SIGN(1.,XE*XCOS+YE*XSIN)*SQRT(XE*XE+YE*YE) x=(yprf-yy)*xcos-(xprf-xx)*xsin IF(MOUT.EQ.3)WRITE(LOU,101) 1IND,IND1,ITER,KMAH,DD,X,XX,YY,T,PNEW,PSI0 C C TESTING WHETHER THE RAY OF A QS WAVE DOES NOT TERMINATE C OUTSIDE THE RANGE IN WHICH PREVIOUS RAYS TERMINATED; IF YES, C A RAY OF THE OTHER QS WAVE WITH THE SAME INITIAL PARAMETERS C IS CALCULATED C if((x-x1)*(x-x2).gt.0..and.itype.eq.1)then iter=iter+1 if(iter.gt.itmax)go to 80 code(1,2)=2 iwave=1 go to 71 end if if((x-x1)*(x-x2).gt.0..and.itype.eq.2)then iter=iter+1 if(iter.gt.itmax)go to 80 code(1,2)=1 iwave=1 go to 71 end if IF(ICLS.NE.0)GO TO 70 IF(IND.NE.ITPR)P2=PNEW IF(IND.NE.ITPR)GO TO 61 IF(ABS(X-XAX).LT..000001)GO TO 67 IF(ABS(X-DD).LT.REPS)GO TO 65 IF(X1.LT.X2.AND.DD.GT.X)GO TO 63 IF(X1.GT.X2.AND.DD.LT.X)GO TO 63 IF(ABS(X-X1).LT..000001)GO TO 67 P2=PNEW X2=X T2=T IPR2=1 GO TO 68 63 IF(ABS(X-X2).LT..000001)GO TO 67 P1=PNEW X1=X T1=T IPR1=1 GO TO 68 67 IF(ABS(PNEW-PAX).GT..000001)ITER=ITMAX AX1=X1-DD AX2=X2-DD IF((IPR1*IPR2).EQ.0)ITER=ITMAX X=X1 PNEW=P1 IF(ABS(AX1).GT.ABS(AX2))PNEW=P2 IF(ABS(AX1).GT.ABS(AX2))X=X2 IF(ITER.EQ.ITMAX)GO TO 61 ICLS=1 GO TO 69 70 ICLS=0 GO TO 65 C C SUCCESSFUL ITERATIONS C 65 INDEX=INDEX+1 isuc=1 RANG=rPRF XXX=XX YYY=YY ZZZ=ZZ TTT=T PAZM=PNEW XAX=X IF(MOUT.EQ.3)WRITE(LOU,113)DD,X,XX,YY,ZZ,T,PNEW,PSI0, 1IND,IND1,ITER,II,INDEX GO TO 98 C 80 continue P1=PNEW X1=X T1=T IF(ITER.GT.ITMAX)P1=AAOLD IF(ITER.GT.ITMAX)X1=XXOLD IF(ITER.GT.ITMAX)T1=TTOLD P2=AANEW X2=XXNEW T2=TTNEW GO TO 3 C C E N D O F I T E R A T I O N S C C BOUNDARY RAYS: CASE IOLD.NE.ITPR, INEW=ITPR C OR CASE IOLD=ITPR, INEW=ITPR BUT IOLD1.NE.INEW1 C (IOLD1.NE.INDR1, INEW1=INDR1) C 50 XXNEW=XNEW TTNEW=TNEW AANEW=ANEW IBNEW=0 P1=AOLD P2=ANEW 54 IRES=0 ITER=0 51 PNEW=0.5*(P1+P2) ITER=ITER+1 IND=INDS NDER=1 CALL RAYA(XSOUR,YSOUR,ZSOUR,TSOUR,PSI0,PNEW,PX,PY,PZ,XX,YY,ZZ,T, 1DT,AC) x=(yprf-yy)*xcos-(xprf-xx)*xsin IF(MOUT.EQ.3)WRITE(LOU,103)IND,IND1,ITER,X,XX,YY,T,PNEW,PSI0 IF(IND.EQ.ITPR.AND.IND1.EQ.Indr1)GO TO 52 P1=PNEW if((x-dd)*(xnew-dd).gt.0.)iter=itmax GO TO 53 52 XXOLD=X AAOLD=PNEW TTOLD=T IBOLD=1 if((x-dd)*(xnew-dd).lt.0.)iter=itmax IF(ABS(X-XAX).LT.REPS1)ITER=ITMAX IRES=1 XAX=X T1=T P2=PNEW 53 IF(ITER.LT.ITMAX)GO TO 51 IF(MOUT.EQ.3)WRITE(LOU,107)X,ZZ,XX,YY,T,PNEW,IND,IND1,IRES IF(IRES.EQ.1) GOTO 41 GO TO 3 C C BOUNDARY RAYS: CASE IOLD=ITPR, INEW.NE.ITPR C OR CASE IOLD=ITPR, INEW=ITPR BUT IOLD1.NE.INEW1 C (IOLD1=INDR1, INEW1.NE.INDR1) C 55 XXOLD=XOLD AAOLD=AOLD TTOLD=TOLD IBOLD=0 P1=AOLD P2=ANEW 59 IRES=0 ITER=0 56 PNEW=0.5*(P1+P2) ITER=ITER+1 IND=INDS NDER=1 CALL RAYA(XSOUR,YSOUR,ZSOUR,TSOUR,PSI0,PNEW,PX,PY,PZ,XX,YY,ZZ,T, 1DT,AC) x=(yprf-yy)*xcos-(xprf-xx)*xsin IF(MOUT.EQ.3)WRITE(LOU,103)IND,IND1,ITER,X,XX,YY,T,PNEW,PSI0 IF(IND.EQ.ITPR.AND.IND1.EQ.Indr1)GO TO 57 P2=PNEW if((x-dd)*(xold-dd).gt.0.)iter=itmax GO TO 58 57 XXNEW=X AANEW=PNEW TTNEW=T IBNEW=1 if((x-dd)*(xold-dd).lt.0.)iter=itmax IF(ABS(X-XAX).LT.REPS1.AND.IRES.EQ.1) ITER=ITMAX IRES=1 XAX=X T2=T P1=PNEW 58 IF(ITER.LT.ITMAX)GO TO 56 IF(MOUT.EQ.3)WRITE(LOU,107)X,ZZ,XX,YY,T,PNEW,IND,IND1,IRES IF(IRES.EQ.1)GOTO 41 GO TO 3 C C CRITICAL RAYS. CASE IOLD.NE.9, IOLD.NE.3, INEW=9 C 30 ITER=0 XCR=XNEW P1=AOLD P2=ANEW IRES=0 31 PNEW=0.5*(P1+P2) ITER=ITER+1 IND=INDS NDER=1 CALL RAYA(XSOUR,YSOUR,ZSOUR,TSOUR,PSI0,PNEW,PX,PY,PZ,XX,YY,ZZ,T, 1DT,AC) x=(yprf-yy)*xcos-(xprf-xx)*xsin IF(MOUT.EQ.3)WRITE(LOU,104)IND,IND1,ITER,X,XX,YY,T,PNEW,PSI0 IF(IND.EQ.9)GO TO 32 IF(IND.EQ.ITPR)GO TO 33 P1=PNEW GO TO 34 32 CONTINUE C 32 IF(IND1.NE.INEW1)P1=PNEW C IF(IND1.NE.INEW1) GOTO 34 P2=PNEW IF(ABS(X-XCR).LT.0.01.AND.KC.NE.0.AND.IRES.EQ.1) GOTO 89 XCR=X GOTO 34 89 ITER=ITMAX-1 GO TO 31 33 IF(ABS(X-XAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX IRES=1 XAX=X T2=T P1=PNEW PAP=PNEW 34 IF(ITER.LT.ITMAX)GO TO 31 IF(MOUT.EQ.3)WRITE(LOU,111)X,ZZ,XX,YY,T,PNEW,IND,IND1,IRES IF(IRES.EQ.0) GOTO 3 P2=PAP XXNEW=XAX AANEW=P2 TTNEW=T2 IBNEW=1 P1=AOLD GO TO 54 C C CRITICAL RAYS. CASE IOLD=9, INEW.NE.9, INEW.NE.3. C 35 ITER=0 P1=AOLD P2=ANEW IRES=0 36 PNEW=0.5*(P1+P2) ITER=ITER+1 IND=INDS NDER=1 CALL RAYA(XSOUR,YSOUR,ZSOUR,TSOUR,PSI0,PNEW,PX,PY,PZ,XX,YY,ZZ,T, 1DT,AC) x=(yprf-yy)*xcos-(xprf-xx)*xsin IF(MOUT.EQ.3)WRITE(LOU,104)IND,IND1,ITER,X,XX,YY,T,PNEW,PSI0 IF(IND.EQ.9)GO TO 37 IF(IND.EQ.ITPR)GO TO 38 P2=PNEW GO TO 39 37 IF(IND1.NE.IOLD1)P2=PNEW IF(IND1.NE.IOLD1) GO TO 39 P1=PNEW IF(ABS(X-XCR).LT.0.01.AND.KC.NE.0.AND.IRES.EQ.1) GO TO 94 XCR=X GO TO 39 94 ITER=ITMAX-1 GOTO 36 38 IF(ABS(X-XAX).LT.REPS1.AND.IRES.EQ.1) ITER=ITMAX IRES=1 XAX=X P2=PNEW PAP=PNEW T1=T 39 IF(ITER.LT.ITMAX)GO TO 36 IF(MOUT.EQ.3)WRITE(LOU,111)X,ZZ,XX,YY,T,PNEW,IND,IND1,IRES IF(IRES.EQ.0) GOTO 3 P1=PAP XXOLD=XAX AAOLD=P1 TTOLD=T1 IBOLD=1 P2=ANEW GO TO 59 C C 100 FORMAT('*',2I3,5F10.5,2F15.10) 101 FORMAT(1X,'*','ITERATIONS',5X,4I3,5F10.5,2F15.10) 103 FORMAT(2X,'*','BOUNDARY RAY',5X,3I3,4F10.5,2F15.10) 104 FORMAT(2X,'*','CRITICAL RAY',5X,3I3,4F10.5,2F15.10) 107 FORMAT(10X,'*',5F10.5,F15.10,3I3,5X,'BOUNDARY RAY') 111 FORMAT(10X,'*',5F10.5,F15.10,3I3,5X,'CRITICAL RAY') 113 FORMAT('*',7F10.5,F15.10,5I3) C C 98 CONTINUE LNDST=1 99 CONTINUE RETURN END C C ********************************************************* C SUBROUTINE RAYA (X0,Y0,Z0,T0,FI0,PSI0,PX,PY,PZ,XX,YY,ZZ,T,DT,AC) C C INITIAL-VALUE RAY TRACING BY THE RUNGE-KUTTA METHOD C DIMENSION Y(18),DEP(6),PRM(5),DERY(18),AUX(8,18),DIN(18),VSQ(3) dimension pn(3) COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33, 1 A34,A35,A36,A44,A45,A46,A55,A56,A66, 1 DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23, 1 DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44, 1 DXA45,DXA46,DXA55,DXA56,DXA66, 1 DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23, 1 DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44, 1 DYA45,DYA46,DYA55,DYA56,DYA66, 1 DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23, 1 DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44, 1 DZA45,DZA46,DZA55,DZA56,DZA66, 1 A2546,A1266,A1355,A1456,A3645,A2344 COMMON /AUXI/ IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,mscon,lou, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori INTEGER CODE COMMON /COD/ CODE(50,2),KREF,KC,ITYPE COMMON /INTRF/ Z(1000),SX(350),SY(350),NX(20),NY(20),BRD(6),NINT, 1 XINTA COMPLEX PS COMMON /RAY/ AY(28,2000),DS(20,50),KINT(50),HHH(3,3),tmax, 1 PS(3,7,50),IS(8,50),N,IREF,IND,IND1 COMMON /ZERO/ RNULL common/dyn/xdyn(3,3),ydyn(3,3) common/appr/ xold,xnew,yold(18),dold(18),ynew(18),dnew(18) COMMON/KM/KMAH,SPROLD,TSTOLD EXTERNAL FCT,OUT C IAC=0 kmah=0 Y(1)=X0 Y(2)=Y0 Y(3)=Z0 IREFR=0 KRE=KREF IF(KC.EQ.0) KRE=0 N=0 IREF=1 IF(IND.GT.0) GOTO 6 C C FOR IND=-1 DETERMINATION OF THE NUMBER OF THE LAYER IN WHICH THE C INITIAL POINT OF THE RAY IS SITIUATED C IF(Y(1).LT.BRD(1).OR.Y(1).GT.BRD(2).OR.Y(2).LT.BRD(3).OR.Y(2).GT. 1 BRD(4)) GOTO 4 INTR=1 1 CALL DISC (Y,DEP) ZINT=DEP(1) IF(ABS(Y(3)).GT.0.00001) GOTO 2 ISOUR=1 C C REDEFINITION OF Z-COORDINATES FOR A SOURCE ON THE EARTH SURFACE C Z0=ZINT GOTO 5 2 IF(Y(3).LT.ZINT.AND.INTR.EQ.1) GOTO 4 IF(Y(3).LT.ZINT) GOTO 5 IF(ABS(Y(3)-ZINT).LT.RNULL.AND.INTR.EQ.NINT) GOTO 5 IF(INTR.EQ.NINT) GOTO 4 ISOUR=INTR INTR=INTR+1 GOTO 1 4 WRITE(lou,'(A,/,6F10.5,/,4F10.5)')' Y,BRD',Y,BRD IND=50 RETURN C C DETERMINATION OF INITIAL CONDITIONS FOR THE RUNGE-KUTTA PROCEDURE C 5 IF(IND.GE.0) GOTO 6 IND=ISOUR RETURN 6 LAY=ISOUR is(3,1)=lay INT1=ISOUR IF(ISOUR.NE.CODE(1,1)) IND=14 IF(ISOUR.NE.CODE(1,1)) RETURN ITYPE=CODE(1,2) CALL PARDIS(Y,0) C C DETERMINATION OF INITIAL VALUES FOR RAY TRACING C AND DYNAMIC RAY TRACING C kdim=6 if(nder.gt.1)kdim=18 csp=cos(psi0) snp=sin(psi0) csf=cos(fi0) snf=sin(fi0) if(mori.eq.0)then Y(4)=CSP*CSF Y(5)=SNP*CSF Y(6)=SNF else y(4)=csf*csp y(5)=snp y(6)=snf*csp end if do 8 k=1,3 pn(k)=y(k+3) 8 continue if(nder.gt.1)then do 3 k=7,12 y(k)=0. 3 continue if(mori.eq.0)then y(13)=-snp y(14)=csp y(15)=0. y(16)=-csp*snf y(17)=-snp*snf y(18)=csf else y(13)=-snf y(14)=0. y(15)=csf y(16)=-csf*snp y(17)=csp y(18)=-snf*snp end if end if IF(IANI(ISOUR).ne.0)then C C SOURCE LOCATED IN AN ANISOTROPIC LAYER C CALL INIT(pn,VSQ) IF(IPRINT.GT.2)WRITE(lou,'(a,3F14.6)')' V1,V2,V3=', VSQ VP=AMAX1(VSQ(1),VSQ(2),VSQ(3)) VS1=AMIN1(VSQ(1),VSQ(2),VSQ(3)) VS2=VSQ(1)+VSQ(2)+VSQ(3)-VP-VS1 IF(ITYPE.EQ.3)V=SQRT(VP) IF(ITYPE.EQ.1)V=SQRT(VS1) IF(ITYPE.EQ.2)V=SQRT(VS2) do 7 i=4,6 y(i)=y(i)/v 7 continue if(nder.gt.1)then nder=1 call fct(0.,y,dery) nder=2 vg=sqrt(dery(1)*dery(1)+dery(2)*dery(2)+dery(3)*dery(3)) if(mori.eq.0)then auxp=-dery(1)*snp+dery(2)*csp auxf=-dery(1)*csp*snf-dery(2)*snp*snf+dery(3)*csf y(13)=y(13)-auxp*csp*csf/v y(14)=y(14)-auxp*snp*csf/v y(15)=y(15)-auxp*snf/v y(16)=y(16)-auxf*csp*csf/v y(17)=y(17)-auxf*snp*csf/v y(18)=y(18)-auxf*snf/v else auxp=-dery(1)*snf+dery(3)*csf auxf=-dery(1)*csf*snp+dery(2)*csp-dery(3)*snf*snp y(13)=y(13)-auxp*csf*csp/v y(14)=y(14)-auxp*snp/v y(15)=y(15)-auxp*snf*csp/v y(16)=y(16)-auxf*csf*csp/v y(17)=y(17)-auxf*snp/v y(18)=y(18)-auxf*snf*csp/v end if do 11 i=13,18 y(i)=y(i)/v 11 continue C C DETERMINATION OF THE SOURCE INDEX IN ANISOTROPIC MEDIUM C call fct(0.,y,dery) aaa=y(4)*dery(1)+y(5)*dery(2)+y(6)*dery(3) if(abs(aaa-1.).gt.1.0e-02)then ind=10 return end if el=-(dery(7)*y(13)+dery(8)*y(14)+dery(9)*y(15))/vg em=-(dery(10)*y(13)+dery(11)*y(14)+dery(12)*y(15))/vg en=-(dery(10)*y(16)+dery(11)*y(17)+dery(12)*y(18))/vg ee=y(13)*y(13)+y(14)*y(14)+y(15)*y(15) ff=y(13)*y(16)+y(14)*y(17)+y(15)*y(18) gg=y(16)*y(16)+y(17)*y(17)+y(18)*y(18) egf=(ee*gg-ff*ff)/v/v be=el*gg+en*ee-2.*em*ff ce=el*en-em*em if(egf.gt.0.)then if(ce.lt.0.)kmah=-1 if(ce.gt.0.)then if(be.lt.0.)kmah=0 if(be.gt.0.)kmah=-2 end if end if if(egf.lt.0.)then if(ce.gt.0.)kmah=-1 if(ce.lt.0.)then if(be.gt.0.)kmah=0 if(be.lt.0.)kmah=-2 end if end if end if end if C C SOURCE LOCATED IN AN ISOTROPIC LAYER C IF(IANI(ISOUR).eq.0)then IF(ITYPE.EQ.3)V=SQRT(A11) IF(ITYPE.NE.3)V=SQRT(A44) do 9 i=4,kdim y(i)=y(i)/v 9 continue end if C IND=0 IND1=0 PRM(1)=T0 PRM(2)=TMAX PRM(3)=DT IF(ITYPE.NE.3) PRM(3)=DT*1.7 PRM(4)=AC T=PRM(1) 20 CONTINUE DO 10 I=1,3 auxx=y(4)*y(4)+y(5)*y(5)+y(6)*y(6) auxx=sqrt(auxx) DIN(I)=auxx din(i+3)=prm(3)/auxx 10 CONTINUE do 25 i=7,kdim din(i)=0. 25 continue DO 30 I=1,kdim DERY(I)=DIN(I) 30 CONTINUE C C COMPUTATION OF THE RAY C CALL RKGS(PRM,Y,DERY,kdim,IHLF,FCT,OUT,AUX) IF(IHLF.EQ.11) IND=5 IF(IHLF.EQ.12) IND=6 IF(IHLF.EQ.13) IND=7 IF(IND.GE.5.AND.IND.LE.7) RETURN IF(ABS(PRM(5)).GT.0.0001) GOTO 35 IF(IND.eq.12) GOTO 70 GOTO 60 35 CONTINUE XX=Y(1) YY=Y(2) ZZ=Y(3) T=AY(1,N) IF(ABS(PRM(5)-2.).GT.0.0001) GOTO 80 C C INTEGRATION FROM THE POINT OF REFLECTION/TRANSMISSION TO THE CLOSEST C POINT OF THE RAY CORRESPONDING TO REGULAR TIME STEP C PRM(1)=AY(1,N) I=INT((PRM(1)-T0)/DT) PRM(2)=FLOAT(I+1)*DT+T0 PRM(3)=DT GOTO 20 60 PRM(1)=PRM(2) PRM(2)=TMAX PRM(3)=DT IF(ITYPE.NE.3) PRM(3)=1.7*DT N=N-1 GOTO 20 70 CONTINUE XX=Y(1) YY=Y(2) ZZ=Y(3) 80 continue c if(kmah.ne.0)ind1=ind1+50 IF(IREFR.EQ.1) IND1=-IND1 PX=Y(4) PY=Y(5) PZ=Y(6) if(nder.gt.1)then do 90 i=1,3 xdyn(i,1)=y(i+6) xdyn(i,2)=y(i+9) xdyn(i,3)=dery(i) ydyn(i,1)=y(i+12) ydyn(i,2)=y(i+15) ydyn(i,3)=dery(i+3) 90 continue end if RETURN END C C ********************************************************* C SUBROUTINE RECEIV(XSOUR,YSOUR,ZSOUR,TSOUR,DT,AC,ITMAX,ASTART, 1ASTEP,AFIN,BMIN,BSTEP,BMAX,MOUT,LU1,LU2,METHOD,ITPR,NCD) C C TWO-POINT RAY TRACING C COMPLEX AMPX1,AMPX2,AMPY1,AMPY2,AMPZ1,AMPZ2,APX,APY,APZ DIMENSION JC(50,2),YDD(2),DEP(6) DIMENSION TIME(500),XCOOR(500),ZCOOR(500),INDI(500),AMPX1(500), 1AMPY1(500),AMPZ1(500),AMPX2(500),AMPY2(500),AMPZ2(500), 2p(500,3),pol(500,3),pol1(500,3),APX(2),APY(2),APZ(2) COMMON /AUXI/ IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori INTEGER CODE COMMON /COD/ CODE(50,2),KREF,KC,ITYPE COMMON /DIST/ DST(200),NDST,REPS,PROF(2),NDSTP,PREPS,LNDST, 1XPRF,YPRF,ILOC COMPLEX PS,CKMAH COMMON /RAY/ AY(28,2000),DS(20,50),KINT(50),HHH(3,3),tmax, 1 PS(3,7,50),IS(8,50),N,IREF,IND,IND1 COMMON /RAY2/ DRY(3,2000) COMMON/VSP/XVSP,YVSP,XNRM,YNRM,ICOD,IVSP COMMON/DYN/XDYN(3,3),ydyn(3,3) COMMON/KM/KMAH,SPROLD,TSTOLD COMMON/VRML/LUBRD,LUGRD,LUIND,LURAY C ITER=0 indr1=0 II=0 DD=0. C AA=ASTART-ASTEP BMIN1=BMIN REPS1=.05*REPS INDEX=0 INUM=0 IK1=0 ICR=0 INDS=ISOUR C C LOOP IN DECLINATION, FROM VALUE ASTART TO AFIN, WITH THE STEP C ASTEP C 1 AA=AA+ASTEP PNEW=AA IF(ASTEP.GT.0..AND.AA.GT.AFIN)GO TO 99 IF(ASTEP.LT.0..AND.AA.LT.AFIN)GO TO 99 4 IND=INDS CALL PROFIL(XSOUR,YSOUR,ZSOUR,TSOUR,AA,PAZM,r,XX,YY,Zz,T, 1DT,AC,BMIN,BSTEP,BMAX,ITMAX,MOUT,NCD,METHOD,ITPR,indr1) IF(IND.EQ.14)RETURN IF(ITPR.EQ.43)r=Zz IF(NDST.EQ.0.AND.LNDST.EQ.0)THEN BMIN=BMIN1 GO TO 1 END IF IF(NDST.EQ.0)GO TO 65 IF(IND.EQ.ITPR)rAX=r IF(IND.EQ.ITPR)PNEW=AA IF(MOUT.GE.2)WRITE(LOU,100)IND,IND1,KMAH,r,T,AA IF(INUM.NE.0)GO TO 2 dOLD=AA aold=pazm IOLD=IND IOLD1=IND1 rOLD=r xold=xx yold=yy zold=zz TOLD=T INUM=1 GO TO 1 C C PARAMETERS FOR THE PRECEDING RAY: C DOLD (AA), AOLD (PAZM), ROLD (R), IOLD (IND) C PARAMETERS FOR THE NEW RAY: C DNEW (AA), ANEW (PAZM), RNEW (R), INEW (IND) C 2 INEW=IND INEW1=IND1 dNEW=AA anew=pazm rNEW=r xnew=xx ynew=yy znew=zz TNEW=T IF(INEW.EQ.ITPR.AND.IOLD.EQ.ITPR)GO TO 21 IF(INEW.EQ.ITPR)GO TO 50 IF(IOLD.EQ.ITPR)GO TO 55 IF(INEW.EQ.9.AND.IOLD.NE.9.AND.IOLD.NE.8)GO TO 30 IF(INEW.NE.9.AND.INEW.NE.8.AND.IOLD.EQ.9)GO TO 35 c indr1=ind1 GO TO 3 21 IF(IOLD1.NE.INEW1)IK1=2 IF(IK1.EQ.2)GO TO 55 c indr1=ind1 GO TO 40 C C NO ITERATIONS, TAKE A NEW RAY IN THE LOOP C 3 CONTINUE IOLD=INEW IOLD1=INEW1 rOLD=rNEW xold=xnew yold=ynew zold=znew dOLD=dNEW aold=anew TOLD=TNEW GO TO 1 C C REGULAR CASE: IOLD=3, INEW=3 C 40 rrNEW=rNEW xxnew=xnew yynew=ynew zznew=znew rrOLD=rOLD xxold=xold yyold=yold zzold=zold ddNEW=dNEW aanew=anew ddOLD=dOLD aaold=aold TTNEW=TNEW TTOLD=TOLD IBNEW=0 IBOLD=0 41 IF(rrNEW.GT.rrOLD.AND.DST(2).GT.DST(1))GO TO 46 IF(rrNEW.LT.rrOLD.AND.DST(2).LT.DST(1))GO TO 46 C C REGULAR CASE: RRNEW.LE.RROLD, ITREND=-1 (REVERSE BRANCH) C ITREND=-1 DO 42 I=1,NDST NDD=NDST-I+1 DD=DST(NDD) dr=rrOLD IF(IBOLD.EQ.1)dr=dr+REPS IF(DD.GE.dr)GO TO 42 dr=rrNEW IF(IBNEW.EQ.1)dr=dr-REPS IF(DD.LT.dr.AND.IK1.NE.0)GO TO 6 IF(DD.LT.dr)GO TO 3 II=NDD GO TO 43 42 CONTINUE IF(IK1.NE.0)GO TO 6 GO TO 3 C C REGULAR CASE: RRNEW.GT.RROLD, ITREND=1 (REGULAR BRANCH) C 46 ITREND=1 DO 44 I=1,NDST DD=DST(I) dr=rrOLD IF(IBOLD.EQ.1)dr=dr-REPS IF(DD.LE.dr)GO TO 44 dr=rrNEW IF(IBNEW.EQ.1)dr=dr+REPS IF(DD.GT.dr.AND.IK1.NE.0)GO TO 6 IF(DD.GT.dr)GO TO 3 II=I GO TO 43 44 CONTINUE IF(IK1.NE.0)GO TO 6 GO TO 3 43 d1=ddOLD a1=aaold d2=ddNEW a2=aanew x1=xxold y1=yyold z1=zzold x2=xxnew y2=yynew z2=zznew T1=TTOLD T2=TTNEW C C I T E R A T I O N S C 60 continue C C TRANSFORMATION OF COORDINATES OF A RECEIVER FROM CYLINDRICAL C TO CARTESIAN COORDINATES C IF(ITPR.NE.43)THEN XD=Xprf+DD*COS(PROF(1)) YD=Yprf+DD*SIN(PROF(1)) YDD(1)=XD YDD(2)=YD INTR=1 IF(ITPR.GT.100)INTR=ITPR-100 CALL DISC(YDD,DEP) ZD=DEP(1) END IF IF(ITPR.EQ.43)THEN XD=XVSP YD=YVSP ZD=DD END IF DELX=XD-X1 DELY=YD-Y1 DELZ=ZD-Z1 dr1=sqrt(delx*delx+dely*dely+delz*delz) DELX=XD-X2 DELY=YD-Y2 DELZ=ZD-Z2 dr2=sqrt(delx*delx+dely*dely+delz*delz) c ITER=0 GO TO 61 68 rAX=r 61 ITER=ITER+1 IF(ITER.GT.ITMAX)GO TO 80 C C PREPARATION FOR ITERATIVE SOLUTION OF TWO-POINT RAY TRACING C if(method.eq.2.or.itpr.ne.ind)go to 69 C C PARAXIAL RAY APPROXIMATION C AUX1=XDYN(2,1)*XDYN(3,2)-XDYN(3,1)*XDYN(2,2) AUX2=XDYN(3,1)*XDYN(1,2)-XDYN(1,1)*XDYN(3,2) AUX3=XDYN(1,1)*XDYN(2,2)-XDYN(2,1)*XDYN(1,2) DET=AUX1*XDYN(1,3)+AUX2*XDYN(2,3)+AUX3*XDYN(3,3) IF(ABS(DET).LT..0000001)GO TO 69 AUX1=DELY*XDYN(3,2)-DELZ*XDYN(2,2) AUX2=DELZ*XDYN(1,2)-DELX*XDYN(3,2) AUX3=DELX*XDYN(2,2)-DELY*XDYN(1,2) IF(mori.eq.0)CSF=COS(PNEW) if(mori.ne.0)csf=cos(pazm) IF(ABS(csf).LT..0000001)GO TO 69 aux11=(AUX1*XDYN(1,3)+AUX2*XDYN(2,3)+AUX3*XDYN(3,3))/DET/csf AUX1=DELZ*XDYN(2,1)-DELY*XDYN(3,1) AUX2=DELX*XDYN(3,1)-DELZ*XDYN(1,1) AUX3=DELY*XDYN(1,1)-DELX*XDYN(2,1) aux22=(AUX1*XDYN(1,3)+AUX2*XDYN(2,3)+AUX3*XDYN(3,3))/DET if(mori.eq.0)pazm=pazm+aux11 if(mori.ne.0)pazm=pazm+aux22 if(mori.eq.0)pnew=pnew+aux22 if(mori.ne.0)pnew=pnew+aux11 go to 72 C C HALVING OF INTERVAL C 69 pnew=0.5*(d1+d2) pazm=0.5*(a1+a2) C C INITIAL ANGLES FOR A NEW RAY WERE DETERMINED C 72 ind=inds rOLD=0. SPROLD=0. IF(MDIM.GE.1)NDER=2 CALL RAYA(XSOUR,YSOUR,ZSOUR,TSOUR,PNEW,PAZM,PX,PY,PZ,XX,YY,Zz, 1T,DT,AC) NDER=1 XE=XX-Xprf YE=YY-Yprf r=SQRT(XE*XE+YE*YE) delx=xd-xx dely=yd-yy delz=zd-zz drs=sqrt(delx*delx+dely*dely+delz*delz) IF(ITPR.EQ.43)r=Zz IF(MOUT.GE.2)WRITE(LOU,101)IND,IND1,ITER,KMAH,DD,r,T,PNEW,PAZM if(mout.eq.4)write(lou,120)xd,yd,zd,xx,yy,zz 120 format(1x,'(x,y,z) receiver',3F15.8/1x,'(xx,yy,zz) ray',3F15.8) IF(IND.NE.ITPR)then d2=PNEW a2=pazm x2=xx y2=yy z2=zz dr2=drs GO TO 61 end if IF(drs.LT.REPS)GO TO 65 IF(dr2.LT.dr1)GO TO 63 d2=PNEW a2=pazm x2=xx y2=yy z2=zz dr2=drs T2=T GO TO 68 63 continue d1=PNEW a1=pazm x1=xx y1=yy z1=zz dr1=drs T1=T GO TO 68 C C SUCCESSFUL ITERATIONS C 65 INDEX=INDEX+1 IF(NDST.EQ.0.AND.IND.NE.3)GO TO 900 IF(MDIM.NE.0)CALL AMPL(APX,APY,APZ,UU) 900 CONTINUE IF(LURAY.NE.0)GO TO 800 IF(LU1.EQ.0.OR.NDST.EQ.0)GO TO 800 TIME(INDEX)=T XCOOR(INDEX)=r ZCOOR(INDEX)=Zz INDI(INDEX)=II 800 CONTINUE rAX=r IF(MOUT.GE.1)WRITE(LOU,113)DD,r,XX,YY,Zz,T,PNEW,PAZM, 1IND,IND1,ITER,II,INDEX IF(LU1.NE.0.AND.NDST.NE.0)WRITE(LU1,105)N,II IF(LU1.NE.0.AND.NDST.NE.0) 1WRITE(LU1,108)(AY(1,I),AY(2,I),AY(3,I),AY(4,I),AY(5,I),AY(6,I), 2AY(7,I),AY(8,I),AY(12,I),AY(16,I),AY(17,I),AY(18,I),AY(19,I), 3AY(20,I),AY(21,I),AY(22,I),I=1,N) IF(MDIM.EQ.0)GO TO 80 IF(IND.NE.ITPR)GO TO 80 SPR=1. CKMAH=(1.,0.) IF(MDIM.EQ.2)THEN IF(KMAH.NE.0)THEN PH=-1.57079632*KMAH CSKMAH=COS(PH) SNKMAH=SIN(PH) CKMAH=CMPLX(CSKMAH,SNKMAH) END IF SPR1=XDYN(2,1)*XDYN(3,2)-XDYN(3,1)*XDYN(2,2) SPR2=XDYN(3,1)*XDYN(1,2)-XDYN(1,1)*XDYN(3,2) SPR3=XDYN(1,1)*XDYN(2,2)-XDYN(2,1)*XDYN(1,2) SPR=XDYN(1,3)*SPR1+XDYN(2,3)*SPR2+XDYN(3,3)*SPR3 VV=ay(5,n)*ay(5,n)+ay(6,n)*ay(6,n)+ay(7,n)*ay(7,n) SPR=SPR*SQRT(VV) SPR=SQRT(ABS(SPR)) IF(MOUT.GE.2)WRITE(LOU,110)XDYN IF(MOUT.GE.2)WRITE(LOU,112)yDYN END IF DO 802 I=1,2 APX(I)=APX(I)*UU*CKMAH/SPR APY(I)=APY(I)*UU*CKMAH/SPR APZ(I)=APZ(I)*UU*CKMAH/SPR 802 CONTINUE IF(MOUT.GE.1) 1WRITE(LOU,'(2X,F8.5,6(2X,E11.5)/10X,6(2X,E11.5),F10.5,I5)') 2UU,(APX(I),APY(I),APZ(I),I=1,2),SPR,KMAH TAUX=T TAST=0. NCC=code(1,2) ncod=ncd IF(iani(isour).eq.0.and.ncc.ne.3)NCOD=-NCD call polar(1,1,1,1) IF(LU2.NE.0.AND.NDST.NE.0)then WRITE(LU2,116)ncod,II,T,APX(1),APY(1),APZ(1),TAST if(ncc.eq.1.and.ncod.lt.0)WRITE(LU2,115)APX(2),APY(2),APZ(2) WRITE(LU2,114)ay(5,1),ay(6,1),ay(7,1) if(ncc.eq.1)WRITE(LU2,114)(hhh(1,i),i=1,3) if(ncc.eq.1.and.ncod.lt.0)WRITE(LU2,114)(hhh(2,i),i=1,3) if(ncc.eq.2)WRITE(LU2,114)(hhh(2,i),i=1,3) if(ncc.eq.3)WRITE(LU2,114)(hhh(3,i),i=1,3) c write(LU2,114)spr end if IF(LURAY.NE.0)GO TO 801 IF(LU1.EQ.0.OR.NDST.EQ.0)GO TO 801 AMPX1(INDEX)=APX(1) AMPY1(INDEX)=APY(1) AMPZ1(INDEX)=APZ(1) if(ncc.eq.1.and.ncod.lt.0)then AMPX2(INDEX)=APX(2) AMPY2(INDEX)=APY(2) AMPZ2(INDEX)=APZ(2) end if p(index,1)=ay(5,1) p(index,2)=ay(6,1) p(index,3)=ay(7,1) if(ncc.eq.1)then pol(index,1)=hhh(1,1) pol(index,2)=hhh(1,2) pol(index,3)=hhh(1,3) end if if(ncc.eq.1.and.ncod.lt.0)then pol1(index,1)=hhh(2,1) pol1(index,2)=hhh(2,2) pol1(index,3)=hhh(2,3) end if if(ncc.eq.2)then pol(index,1)=hhh(2,1) pol(index,2)=hhh(2,2) pol(index,3)=hhh(2,3) end if if(ncc.eq.3)then pol(index,1)=hhh(3,1) pol(index,2)=hhh(3,2) pol(index,3)=hhh(3,3) end if 801 CONTINUE C C GENERATE FILE FOR PLOTTING RAYS C IF(LURAY.NE.0)THEN WRITE(LURAY,119)INDEX DO 803 J=1,N WRITE(LURAY,122)AY(2,J),AY(3,J),AY(4,J) 803 CONTINUE WRITE(LURAY,121) END IF C C 80 IF(NDST.EQ.0.AND.LNDST.EQ.1)THEN BMIN=BMIN+BSTEP GO TO 4 END IF d1=PNEW a1=pazm dr1=drs x1=xx y1=yy z1=zz T1=TAUX d2=ddNEW a2=aanew x2=xxnew y2=yynew z2=zznew T2=TTNEW IF(ITREND.EQ.(-1))GO TO 85 II=II+1 IF(II.GT.NDST.AND.IK1.NE.0)GO TO 6 IF(method.eq.1)then aa=pnew inum=0 end if IF(II.GT.NDST)GO TO 3 DD=DST(II) if(method.eq.1)go to 60 dr=rrNEW IF(IBNEW.EQ.1)dr=dr+REPS IF(DD.GT.dr.AND.IK1.NE.0)GO TO 6 IF(DD.GT.dr)GO TO 3 GO TO 60 85 II=II-1 IF(II.LT.1.AND.IK1.NE.0)GO TO 6 IF(method.eq.1)then aa=pnew inum=0 end if IF(II.LT.1)GO TO 3 DD=DST(II) if(method.eq.1)go to 60 dr=rrNEW IF(IBNEW.EQ.1)dr=dr-REPS IF(DD.LT.dr.AND.IK1.NE.0)GO TO 6 IF(DD.LT.dr)GO TO 3 GO TO 60 C C 6 CONTINUE IF(IK1.EQ.1)GO TO 7 IK1=1 d1=dNEW a1=anew d2=ddNEW a2=aanew IOLD1=INEW1 indr1=iold1 GO TO 59 7 IK1=0 rrOLD=rrNEW xxold=xxnew yyold=yynew zzold=zznew ddOLD=ddNEW aaold=aanew TTOLD=TTNEW IBOLD=IBNEW rrNEW=rNEW xxnew=xnew yynew=ynew zznew=znew ddNEW=dNEW aanew=anew TTNEW=TNEW IBNEW=0 GO TO 41 C C E N D O F I T E R A T I O N S C C C BOUNDARY RAYS. CASE IOLD.NE.3, INEW=3 C 50 rrNEW=rNEW xxnew=xnew yynew=ynew zznew=znew TTNEW=TNEW ddNEW=dNEW aanew=anew IBNEW=0 d1=dOLD d2=dNEW a1=dOLD a2=dNEW 54 IRES=0 ITER=0 51 PNEW=0.5*(d1+d2) ITER=ITER+1 IND=INDS CALL PROFIL(XSOUR,YSOUR,ZSOUR,TSOUR,PNEW,PAZM,r,XX,YY,Zz,T, 1DT,AC,BMIN,BSTEP,BMAX,ITMAX,MOUT,NCD,METHOD,ITPR,indr1) IF(ITPR.EQ.43)r=Zz IF(MOUT.GE.2)WRITE(LOU,103)IND,IND1,ITER,r,T,PNEW IF(IND.EQ.ITPR.AND.LNDST.EQ.1)GO TO 52 d1=PNEW a1=pazm GO TO 53 52 rrOLD=r xxold=xx yyold=yy zzold=zz ddOLD=PNEW aaold=pazm TTOLD=T IBOLD=1 IF(ABS(r-rAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX IRES=1 rAX=r ZAX=Zz IAX=IND IAX1=IND1 T1=T d2=PNEW a2=pazm 53 IF(ITER.LT.ITMAX)GO TO 51 IF(MOUT.GE.1.AND.IRES.EQ.1) 1WRITE(LOU,107)rAX,ZAX,T1,d2,IAX,IAX1,IRES IF(MOUT.GE.1.AND.IRES.EQ.0) 1WRITE(LOU,107)r,Zz,T,PNEW,IND,IND1,IRES IF(IRES.EQ.1)GO TO 41 GO TO 3 C C BOUNDARY RAYS. CASE IOLD=3, INEW.NE.3 C OR CASE IOLD=3, INEW=3 BUT IOLD1.NE.INEW1 C 55 rrOLD=rOLD xxold=xold yyold=yold zzold=zold ddOLD=dOLD aaold=aold TTOLD=TOLD IBOLD=0 d1=dOLD d2=dNEW a1=aOLD a2=aNEW 59 IRES=0 ITER=0 56 PNEW=0.5*(d1+d2) ITER=ITER+1 IND=INDS CALL PROFIL(XSOUR,YSOUR,ZSOUR,TSOUR,PNEW,PAZM,r,XX,YY,Zz,T, 1DT,AC,BMIN,BSTEP,BMAX,ITMAX,MOUT,NCD,METHOD,ITPR,indr1) IF(ITPR.EQ.43)r=Zz IF(MOUT.GE.2)WRITE(LOU,103)IND,IND1,ITER,r,T,PNEW IF(IND.EQ.ITPR.AND.IBOLD.EQ.1.AND.LNDST.EQ.1)GO TO 57 IF(IND.EQ.ITPR.AND.IND1.EQ.IOLD1.AND.LNDST.EQ.1)GO TO 57 d2=PNEW a2=pazm GO TO 58 57 rrNEW=r xxnew=xx yynew=yy zznew=zz ddNEW=PNEW aanew=pazm TTNEW=T IBNEW=1 IF(ABS(r-rAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX IRES=1 rAX=r ZAX=Zz IAX=IND IAX1=IND1 T2=T d1=PNEW a1=pazm 58 IF(ITER.LT.ITMAX)GO TO 56 IF(MOUT.GE.1.AND.IRES.EQ.1) 1WRITE(LOU,107)rAX,ZAX,T2,d1,IAX,IAX1,IRES IF(MOUT.GE.1.AND.IRES.EQ.0) 1WRITE(LOU,107)r,Zz,T,PNEW,IND,IND1,IRES IF(IRES.EQ.1.AND.IK1.EQ.1)GO TO 7 IF(IRES.EQ.1)GO TO 41 GO TO 3 C C CRITICAL RAYS. CASE IOLD.NE.9, IOLD.NE.3, INEW=9 C 30 ITER=0 d1=dOLD d2=dNEW IRES=0 31 PNEW=0.5*(d1+d2) ITER=ITER+1 IND=INDS CALL PROFIL(XSOUR,YSOUR,ZSOUR,TSOUR,PNEW,PAZM,r,XX,YY,Zz,T, 1DT,AC,BMIN,BSTEP,BMAX,ITMAX,MOUT,NCD,METHOD,ITPR,indr1) IF(ITPR.EQ.43)r=Zz IF(MOUT.GE.2)WRITE(LOU,104)IND,IND1,ITER,r,T,PNEW IF(IND.EQ.9)GO TO 32 IF(IND.EQ.ITPR)GO TO 33 d1=PNEW GO TO 34 32 IF(IND1.NE.INEW1)d1=PNEW IF(IND1.NE.INEW1)GO TO 34 d2=PNEW IF(ITER.EQ.ITMAX.AND.KC.NE.0.AND.IRES.EQ.1)GO TO 89 GO TO 34 89 DO 86 K=1,KREF DO 86 L=1,2 86 JC(K,L)=CODE(K,L) IRF3=IREF+3 DO 87 K=IRF3,KREF DO 87 L=1,2 87 CODE(K-2,L)=JC(K,L) KREF1=KREF KREF=KREF-2 ICR=1 ITER=ITMAX-1 GO TO 31 33 IF(ABS(r-rAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX IRES=1 rAX=r rrnew=r xxnew=xx yynew=yy zznew=zz ZAX=Zz IAX=IND IAX1=IND1 T2=T d1=PNEW PAP=PNEW 34 IF(ITER.LT.ITMAX)GO TO 31 IF(MOUT.GE.1.AND.IRES.EQ.1) 1WRITE(LOU,111)rAX,ZAX,T2,PAP,IAX,IAX1,IRES IF(MOUT.GE.1.AND.IRES.EQ.0) 1WRITE(LOU,111)r,Zz,T,PNEW,IND,IND1,IRES IF(IRES.EQ.0)GO TO 3 d2=PAP rrNEW=rAX ddNEW=d2 aanew=anew TTNEW=T2 IBNEW=1 d1=dOLD IF(ICR.EQ.0)GO TO 54 ICR=0 KREF=KREF1 DO 88 K=1,KREF DO 88 L=1,2 88 CODE(K,L)=JC(K,L) GO TO 54 C C CRITICAL RAYS. CASE IOLD=9, INEW.NE.9, INEW.NE.3. C 35 ITER=0 d1=dOLD d2=dNEW IRES=0 36 PNEW=0.5*(d1+d2) ITER=ITER+1 IND=INDS CALL PROFIL(XSOUR,YSOUR,ZSOUR,TSOUR,PNEW,PAZM,r,XX,YY,Zz,T, 1DT,AC,BMIN,BSTEP,BMAX,ITMAX,MOUT,NCD,METHOD,ITPR,indr1) IF(ITPR.EQ.43)r=Zz IF(MOUT.GE.2)WRITE(LOU,104)IND,IND1,ITER,r,T,PNEW IF(IND.EQ.9)GO TO 37 IF(IND.EQ.ITPR)GO TO 38 d2=PNEW GO TO 39 37 IF(IND1.NE.INEW1)d2=PNEW IF(IND1.NE.INEW1)GO TO 39 d1=PNEW IF(ITER.EQ.ITMAX.AND.KC.NE.0.AND.IRES.EQ.1)GO TO 94 GO TO 39 94 DO 91 K=1,KREF DO 91 L=1,2 91 JC(K,L)=CODE(K,L) IRF3=IREF+3 DO 92 K=IRF3,KREF DO 92 L=1,2 92 CODE(K-2,L)=JC(K,L) KREF1=KREF KREF=KREF-2 ICR=1 ITER=ITMAX-1 GO TO 36 38 IF(ABS(r-rAX).LT.REPS1.AND.IRES.EQ.1)ITER=ITMAX IRES=1 rAX=r rrold=r xxold=xx yyold=yy zzold=zz ZAX=Zz IAX=IND IAX1=IND1 d2=PNEW PAP=PNEW T1=T 39 IF(ITER.LT.ITMAX)GO TO 36 IF(MOUT.GE.1.AND.IRES.EQ.1) 1WRITE(LOU,111)rAX,ZAX,T1,PAP,IAX,IAX1,IRES IF(MOUT.GE.1.AND.IRES.EQ.0)WRITE(LOU,111)r,Zz,T,PNEW,IND,IND1,IRES IF(IRES.EQ.0)GO TO 3 d1=PAP rrOLD=rAX ddOLD=d1 aaold=aold TTOLD=T1 IBOLD=1 d2=dNEW IF(ICR.EQ.0)GO TO 59 ICR=0 KREF=KREF1 DO 93 K=1,KREF DO 93 L=1,2 93 CODE(K,L)=JC(K,L) GO TO 59 C C 100 FORMAT(3I3,2F10.5,F15.10) 101 FORMAT(1X,'ITERATIONS',5X,4I3,F10.5,4F15.10) 102 FORMAT(5X,6I3,3F10.5,F15.10) 103 FORMAT(2X,'BOUNDARY RAY',5X,3I3,2F10.5,F15.10) 104 FORMAT(2X,'CRITICAL RAY',5X,3I3,2F10.5,F15.10) 105 FORMAT(2I5) 107 FORMAT(10X,3F10.5,F15.10,3I3,5X,'BOUNDARY RAY') 108 FORMAT(16E15.5) 109 FORMAT(I5,9E15.5) 110 FORMAT(1X,'XDYN',3F10.5) 111 FORMAT(10X,3F10.5,F15.10,3I3,5X,'CRITICAL RAY') 112 FORMAT(1X,'YDYN',3F10.5) 113 FORMAT(7F10.5,F15.10,5I3) 114 FORMAT(3F10.5) 115 FORMAT(6E12.6) 116 FORMAT(2I3,F12.7,6E12.6,F10.6) 117 FORMAT(6E15.5) 119 FORMAT(2H'R,I3,1H',1X,'/') 121 FORMAT('/') 122 FORMAT(3(F10.5,1X),'/') C C 99 N=0 NAUX=0 IF(LURAY.NE.0)RETURN IF(LU1.NE.0.AND.NDST.NE.0)WRITE(LU1,100)N,NAUX IF(NCC.EQ.1.AND.NCOD.LT.0)INDEX=-INDEX IF(LU1.NE.0.AND.NDST.NE.0)WRITE(LU1,100)INDEX IF(INDEX.EQ.0)RETURN IF(INDEX.LT.0)INDEX=-INDEX IF(LU1.NE.0.AND.NDST.NE.0)THEN DO 200 I=1,INDEX WRITE(LU1,109)INDI(I),XCOOR(I),ZCOOR(I),TIME(I), 1 AMPX1(I),AMPY1(I),AMPZ1(I) IF(NCC.EQ.1.AND.NCOD.LT.0) 1 WRITE(LU1,117)AMPX2(I),AMPY2(I),AMPZ2(I) WRITE(LU1,108)(P(I,K),K=1,3) WRITE(LU1,108)(POL(I,K),K=1,3) IF(NCC.EQ.1.AND.NCOD.LT.0) 1 WRITE(LU1,108)(POL1(I,K),K=1,3) 200 CONTINUE END IF RETURN END C C ********************************************************* C SUBROUTINE RKGS(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) DIMENSION Y(18),DERY(18),AUX(8,18),A(4),B(4),C(4),PRMT(5) EXTERNAL FCT,OUTP DO 1 I=1,NDIM 1 AUX(8,I)=.0666667*DERY(I) X=PRMT(1) XEND=PRMT(2) H=PRMT(3) PRMT(5)=0. CALL FCT(X,Y,DERY) IF(H*(XEND-X))38,37,2 2 A(1)=.5 A(2)=.2928932 A(3)=1.707107 A(4)=.1666667 B(1)=2. B(2)=1. B(3)=1. B(4)=2. C(1)=.5 C(2)=.2928932 C(3)=1.707107 C(4)=.5 DO 3 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=0. 3 AUX(6,I)=0.0 IREC=0 H=H+H IHLF=-1 ISTEP=0 IEND=0 4 IF((X+H-XEND)*H)7,6,5 5 H=XEND-X 6 IEND=1 7 CALL OUTP(X,Y,DERY,IREC,NDIM,PRMT) IF(PRMT(5))40,8,40 8 ITEST=0 9 ISTEP=ISTEP+1 J=1 10 AJ=A(J) BJ=B(J) CJ=C(J) DO 11 I=1,NDIM R1=H*DERY(I) R2=AJ*(R1-BJ*AUX(6,I)) Y(I)=Y(I)+R2 R2=R2+R2+R2 11 AUX(6,I)=AUX(6,I)+R2-CJ*R1 IF(J-4)12,15,15 12 J=J+1 IF(J-3)13,14,13 13 X=X+.5*H 14 CALL FCT(X,Y,DERY) GO TO 10 15 IF(ITEST)16,16,20 16 DO 17 I=1,NDIM 17 AUX(4,I)=Y(I) ITEST=1 ISTEP=ISTEP+ISTEP-2 18 IHLF=IHLF+1 X=X-H H=.5*H DO 19 I=1,NDIM Y(I)=AUX(1,I) DERY(I)=AUX(2,I) 19 AUX(6,I)=AUX(3,I) GO TO 9 20 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)21,23,21 21 CALL FCT(X,Y,DERY) DO 22 I=1,NDIM AUX(5,I)=Y(I) 22 AUX(7,I)=DERY(I) GO TO 9 23 DELT=0. DO 24 I=1,NDIM 24 DELT=DELT+AUX(8,I)*ABS(AUX(4,I)-Y(I)) IF(DELT-PRMT(4))28,28,25 25 IF(IHLF-10)26,36,36 26 DO 27 I=1,NDIM 27 AUX(4,I)=AUX(5,I) ISTEP=ISTEP+ISTEP-4 X=X-H IEND=0 GO TO 18 28 CALL FCT(X,Y,DERY) DO 29 I=1,NDIM AUX(1,I)=Y(I) AUX(2,I)=DERY(I) AUX(3,I)=AUX(6,I) Y(I)=AUX(5,I) 29 DERY(I)=AUX(7,I) CALL OUTP(X-H,Y,DERY,IHLF,NDIM,PRMT) IF(PRMT(5))40,30,40 30 DO 31 I=1,NDIM Y(I)=AUX(1,I) 31 DERY(I)=AUX(2,I) IREC=IHLF IF(IEND)32,32,39 32 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H IF(IHLF)4,33,33 33 IMOD=ISTEP/2 IF(ISTEP-IMOD-IMOD)4,34,4 34 IF(DELT-.02*PRMT(4))35,35,4 35 IHLF=IHLF-1 ISTEP=ISTEP/2 H=H+H GO TO 4 36 IHLF=11 CALL FCT(X,Y,DERY) GO TO 39 37 IHLF=12 GO TO 39 38 IHLF=13 39 CALL OUTP(X,Y,DERY,IHLF,NDIM,PRMT) 40 RETURN END C C ********************************************************* C SUBROUTINE ROOT(XA,FXYZA,XB,FXYZB,XINT,PRMT,ITP) C C ROUTINE FOR THE ITERATIVE COMPUTATION OF THE POINTS OF C INTERSECTION OF RAYS WITH INTERFACES C DIMENSION YINT(3),YDINT(3),PRMT(5),DEP(6) COMMON /VSP/XVSP,YVSP,XNRM,YNRM,ICOD,IVSP COMMON /AUXI/ IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori C IAC=0 ISIGN=-1 XXINT=0. 1 IAC=IAC+1 ISIGN=-ISIGN AUX=FXYZA-FXYZB IF(ABS(AUX).LT..000001)XINT=.5*(XA+XB) IF(ABS(AUX).GE..000001.AND.ISIGN.GT.0)XINT=(FXYZA*XB-FXYZB*XA)/AUX IF(ABS(AUX).GE..000001.AND.ISIGN.LT.0)THEN IF(ITP.GT.0)AUX=DEP(2)*YDINT(1)+DEP(3)*YDINT(2)-YDINT(3) IF(ITP.LT.0)AUX=XNRM*YDINT(1)+YNRM*YDINT(2) IF(ABS(AUX).LT..000001)XINT=.5*(XA+XB) IF(ABS(AUX).GE..000001)XINT=XXINT-AUX3/AUX END IF IF(XINT.LT.XA.OR.XINT.GT.XB)XINT=.5*(XA+XB) IF(IRT.EQ.2) 1WRITE(LOU,'(A,/,4E15.9,I5,/)') ' XA,XB,XINTOLD,XINT,IAC', 2XA,XB,XXINT,XINT,IAC IF(IAC.GT.1.AND.ABS(XINT-XXINT).LT.PRMT(4))RETURN IF(IAC.GE.10)RETURN CALL APPROX(XINT,YINT,YDINT,3) XXINT=XINT IF(ITP.GT.0)THEN CALL DISC(YINT,DEP) AUX3=DEP(1)-YINT(3) ELSE AUX3=(YINT(1)-XVSP)*XNRM+(YINT(2)-YVSP)*YNRM END IF IF((FXYZA*AUX3).GT.0.)THEN XA=XINT FXYZA=AUX3 GO TO 1 ELSE XB=XINT FXYZB=AUX3 GO TO 1 END IF END Ca5.for 0100666 0000765 0000765 00000053555 11023551626 011464 0 ustar bulant bulant C
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 STRESS(P,XN,U,TAU) C C ROUTINE FOR THE COMPUTATION OF NORMAL STRESSES C IMPLICIT COMPLEX (P,U,T,E) REAL N1,N2,N3 DIMENSION XN(3),P(3),U(3),TAU(3) COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33, 1 A34,A35,A36,A44,A45,A46,A55,A56,A66, 1 DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23, 1 DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44, 1 DXA45,DXA46,DXA55,DXA56,DXA66, 1 DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23, 1 DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44, 1 DYA45,DYA46,DYA55,DYA56,DYA66, 1 DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23, 1 DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44, 1 DZA45,DZA46,DZA55,DZA56,DZA66, 1 A2546,A1266,A1355,A1456,A3645,A2344 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/DENS/RHO(20) C P1=P(1) P2=P(2) P3=P(3) U1=U(1) U2=U(2) U3=U(3) N1=XN(1) N2=XN(2) N3=XN(3) IF(IANI(LAY).EQ.0) GO TO 1000 C C ANISOTROPIC CASE C RO=1.7+0.2*SQRT(A11) IF(IRHO.NE.0) RO=RHO(LAY) C C DISPLACEMENT TENSOR C E11=U1*P1 E12=0.5*(U1*P2+U2*P1) E13=0.5*(U1*P3+U3*P1) E22=U2*P2 E23=0.5*(U2*P3+U3*P2) E33=U3*P3 C C STRESS TENSOR C P11=A11*E11+A12*E22+A13*E33+2.*A14*E23+2.*A15*E13+2.*A16*E12 P22=A12*E11+A22*E22+A23*E33+2.*A24*E23+2.*A25*E13+2.*A26*E12 P33=A13*E11+A23*E22+A33*E33+2.*A34*E23+2.*A35*E13+2.*A36*E12 P23=A14*E11+A24*E22+A34*E33+2.*A44*E23+2.*A45*E13+2.*A46*E12 P13=A15*E11+A25*E22+A35*E33+2.*A45*E23+2.*A55*E13+2.*A56*E12 P12=A16*E11+A26*E22+A36*E33+2.*A46*E23+2.*A56*E13+2.*A66*E12 C C NORMAL STRESS C TAU(1)=P11*N1+P12*N2+P13*N3 TAU(2)=P12*N1+P22*N2+P23*N3 TAU(3)=P13*N1+P23*N2+P33*N3 TAU(1)=TAU(1)*RO TAU(2)=TAU(2)*RO TAU(3)=TAU(3)*RO RETURN C C ISOTROPIC CASE C 1000 A12=A11-2.*A44 RO=1.7+0.2*SQRT(A11) IF(IRHO.NE.0) RO=RHO(LAY) IF(ICOEF.EQ.0) GOTO 100 WRITE(LOUT,'(A,4F12.7)') 'STRESS: A11,A44,A12,RO',A11,A44,A12,RO WRITE(LOUT,'(A,/,3(2F12.7,/))') 'STRESS: PSTR',P WRITE(LOUT,'(A,/,3(2F12.7,/))') 'STRESS: USTR',U 100 TETA=U1*P1+U2*P2+U3*P3 IF(ICOEF.GT.0) 1WRITE(LOUT,'(A,/,5F12.7)') 'STRESS: UN,TETA',XN,TETA P11=A12*TETA+2.*A44*P1*U1 P12=A44*(P1*U2+P2*U1) P13=A44*(P1*U3+P3*U1) P22=A12*TETA+2.*A44*P2*U2 P23=A44*(P2*U3+P3*U2) P33=A12*TETA+2.*A44*P3*U3 IF(ICOEF.GT.0) 1WRITE(LOUT,'(A,/,6(2F12.7,/))')'STRESS: P11,P12,P13,P22,P23,P33' 1,P11,P12,P13,P22,P23,P33 TAU(1)=P11*N1+P12*N2+P13*N3 TAU(2)=P12*N1+P22*N2+P23*N3 TAU(3)=P13*N1+P23*N2+P33*N3 TAU(1)=TAU(1)*RO TAU(2)=TAU(2)*RO TAU(3)=TAU(3)*RO RETURN END C C ********************************************************* C SUBROUTINE SURFPL(LIN,LU3) C DIMENSION XX(500),YY(500),ZZ(500) COMMON /AUXI/ IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,MORI INTEGER CODE COMMON /COD/ CODE(50,2),KREF,KC,ITYPE COMMON /DIST/ XDST(200),NDSTX,XREPS,DST(2),NDST,REPS,LNDST, 1XPRF,YPRF,ILOC COMPLEX PS COMMON /RAY/ AY(28,2000),DS(20,50),KINT(50),HHH(3,3),TMAX, 1 PS(3,7,50),IS(8,50),N,IREF,IND,IND1 COMMON /RAY2/ DRY(3,2000) C MORI=0 NDST=1 REPS=.001 TMAX=.01 IND=1 NPAR=1 LAY=1 X0=1. Y0=1. Z0=1. DDELTA=.05 AZIM=0. READ(LIN,*)NPAR,LAY READ(LIN,*)X0,Y0,Z0,DDELTA,AZIM WRITE(LOU,100)NPAR,LAY WRITE(LOU,102)X0,Y0,Z0,DDELTA,AZIM DST(1)=AZIM DST(2)=AZIM+.1 XPRF=X0 YPRF=Y0 KREF=1 ITYPE=0 CODE(1,1)=LAY ISOUR=LAY 1 CONTINUE INUM=0 DELTA=0. ITYPE=ITYPE+1 CODE(1,2)=ITYPE IF(NPAR.LE.2)NDSTX=0 IF(NPAR.EQ.3)THEN NDSTX=1 NDER=2 MDIM=2 END IF 2 CONTINUE INUM=INUM+1 IF(NPAR.LE.2)THEN CALL PROFIL(X0,Y0,Z0,0.,DELTA,PAZM,RANG, 1 X,Y,Z,T,.1,.0001,AZIM,1.,AZIM+.1,10,3,1,0,12,0) XX(INUM)=AY(5,1) YY(INUM)=AY(6,1) ZZ(INUM)=AY(7,1) END IF IF(NPAR.EQ.3)THEN CALL PROFIL(X0,Y0,Z0,0.,DELTA,PAZM,RANG, 1 X,Y,Z,T,.1,.0001,AZIM-.5,1.,AZIM+.6,10,3,1,0,12,0) IF(IND.NE.0)THEN XX(INUM)=DRY(1,1) YY(INUM)=DRY(2,1) ZZ(INUM)=DRY(3,1) END IF END IF DELTA=DELTA+DDELTA IF(INUM.LT.500.AND.DELTA.LT.6.3)GO TO 2 WRITE(LU3,100)ITYPE,NPAR,INUM DELTA=0. DO 3 I=1,INUM WRITE(LU3,101)DELTA,XX(I),YY(I),ZZ(I) DELTA=DELTA+DDELTA 3 CONTINUE IF(ITYPE.LT.3)GO TO 1 ITYPE=0 WRITE(LU3,100)ITYPE C 100 FORMAT(16I5) 101 FORMAT(4E15.5) 102 FORMAT(8F10.5) RETURN END C C ********************************************************* C SUBROUTINE TRANSL(Y,P,PNEW,UN,ITRANS,IASW) C C ROUTINE FOR THE COMPUTATION OF THE TRANSFORMATION OF THE SLOWNESS C VECTOR AT AN INTERFACE C SAVE A,AI,B,BI,C,CI,CD,CDI,XCOE,COE,AT,BT,DETB SAVE XKO,PN,Y1,Z,NDER1,LAY1,ISG,IMAG DOUBLE PRECISION XCOE,COE,XH1,XH2 DIMENSION A(3,3),AI(3,3),B(3,3),BI(3,3),C(3,3),CI(3,3),CD(3,3), * CDI(3,3),AC(3,3),ACI(3,3),XH1(3),XH2(3), * P(3),PNEW(3),UN(3),Y(18),PN(3),Y1(18),DERY(18),XCOE(7), * COE(7),XSI(3),ISG(6),IR(3),IT(3),ICR(3),ICT(3),XKO(7) COMPLEX Z(6),Z0,CSQ,CSI(3),XPR,IMAG COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33, 1 A34,A35,A36,A44,A45,A46,A55,A56,A66, 1 DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23, 1 DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44, 1 DXA45,DXA46,DXA55,DXA56,DXA66, 1 DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23, 1 DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44, 1 DYA45,DYA46,DYA55,DYA56,DYA66, 1 DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23, 1 DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44, 1 DZA45,DZA46,DZA55,DZA56,DZA66, 1 A2546,A1266,A1355,A1456,A3645,A2344 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 INTEGER CODE COMMON /COD/ CODE(50,2),KREF,KC,ITYPE COMMON /DENS/ RHO(20) COMPLEX PS COMMON /RAY/ AY(28,2000),DS(20,50),KINT(50),HHH(3,3),tmax, 1 PS(3,7,50),IS(8,50),N,IREF,IND,IND1 COMMON /ZERO/ RNULL C C IASW: SWITCH INDICATING FROM WHICH PROGRAM TRANSL IS CALLED C IASW=0 CALLED FROM AMPL C IASW=1 CALLED FROM OUT C ISG: ISG=0 COMPLEX-VALUED ROOT OF CHARACTERISTIC EQ. C ISG=1 REAL-VALUED ROOT OF CHARACTERISTIC EQ. RAY HAS C POSITIVE PROJECTION ON THE NORMAL TO INTERFACE C ISG=-1 REAL-VALUED ROOT OF CHARACTERISTIC EQ. RAY HAS C NEGATIVE PROJECTION ON THE NORMAL TO INTERFACE C MR: NUMBER OF REAL-VALUED ROOTS RELATED TO REFLECTED WAVE C MT: NUMBER OF REAL-VALUED ROOTS RELATED TO TRANSMITTED WAVE C MCR: NUMBER OF COMPLEX-VALUED ROOTS RELATED TO REFLECTED WAVE C MTR: NUMBER OF COMPLEX-VALUED ROOTS RELATED TO TRANSMITTED WAVE C IMAG=(0.,1.) IRF=IREF-1 C C FILLING DS FIELD C IF(IASW.GT.0) THEN DS(1,IRF)=UN(1) DS(2,IRF)=UN(2) DS(3,IRF)=UN(3) IF(IRHO.EQ.0)DS(8,IRF)=0.2*SQRT(A11)+1.7 IF(IRHO.NE.0)DS(8,IRF)=RHO(LAY) KINT(IRF)=N IS(1,IRF)=LAY IS(2,IRF)=ITRANS END IF CALL PARDIS(Y,0) IF(ITRANS.EQ.0)IS(7,IRF)=LAY IF(ITRANS.NE.0)IS(8,IRF)=LAY DO 1 K=1,3 IT(K)=0 IR(K)=0 ICT(K)=0 ICR(K)=0 ISG(K)=0 ISG(K+3)=0 1 CONTINUE KDIM=6 IF(NDER.GT.1)KDIM=18 DO 2 K=1,KDIM 2 Y1(K)=Y(K) PUN=P(1)*UN(1)+P(2)*UN(2)+P(3)*UN(3) PU1=PUN*UN(1) PU2=PUN*UN(2) PU3=PUN*UN(3) PN(1)=P(1)-PU1 PN(2)=P(2)-PU2 PN(3)=P(3)-PU3 IF(IASW.EQ.0) GO TO 9 NDER1=NDER NDER=1 LAY1=LAY LAY=IS(3,IRF) CALL FCT(X,Y1,DERY) IF(IND.EQ.10)RETURN DUN=DERY(1)*UN(1)+DERY(2)*UN(2)+DERY(3)*UN(3) LAY=LAY1 CALL FCT(X,Y1,DERY) IF(IND.EQ.10)RETURN NDER=NDER1 DS(4,IRF)=PN(1) DS(5,IRF)=PN(2) DS(6,IRF)=PN(3) DS(7,IRF)=abs(dun) 9 IF(KC.NE.0) ITYPE=CODE(IREF,2) IF(IREF.GT.KREF) ITYPE=CODE(KREF,2) IF(IANI(LAY).EQ.0) GOTO 230 C C TRANSMISSION OR REFLECTION FOR THE ANISOTROPIC CASE C IF(MTRNS.GT.0) THEN WRITE(LOUT,999)UN,PN,(Y(K),K=1,KDIM),P 999 FORMAT(' UN,PN,Y,P',/,6(E12.5,2X)) WRITE(LOUT,'(A)')' ITYPE,IREF,LAY,IANI(LAY),N,KINT(IREF-1)' WRITE(LOUT,'(6I5)')ITYPE,IREF,LAY,IANI(LAY),N,KINT(IREF-1) WRITE(LOUT,'(A,/,10F10.5,/,11F10.5)')' A(I,J,K,L)=',A11,A12,A13, 1 A14,A15,A16,A22,A23,A24,A25,A26,A33,A34,A35,A36,A44,A45,A46,A55, 2 A56,A66 END IF CALL CHRM1(A,UN,UN) AT=A(1,1)+A(2,2)+A(3,3) CALL INV(A,AI,DET) XCOE(7)=DET CALL CHRM1(C,PN,PN) CALL INV(C,CI,DET) DO 10 J=1,3 DO 10 K=1,3 CD(J,K)=C(J,K) IF(J.EQ.K)CD(J,K)=CD(J,K)-1 10 CONTINUE CALL INV(CD,CDI,DET) XCOE(1)=DET CALL CHRM1(B,PN,UN) CALL INV(B,BI,DET) BT=B(1,1)+B(2,2)+B(3,3) DETB=DET DO 20 J=1,3 DO 20 K=1,3 AC(J,K)=A(J,K)+C(J,K) 20 CONTINUE CALL INV(AC,ACI,DET) DO 22 J=1,3 XH1(J)=0.0 DO 21 M=1,3 21 XH1(J)=AI(M,J)*B(J,M)+XH1(J) 22 CONTINUE XCOE(6)=XH1(1)+XH1(2)+XH1(3) XCOE(6)=2.*XCOE(6) DO 30 J=1,3 XH1(J)=0.0 DO 29 M=1,3 29 XH1(J)=AI(M,J)*CD(J,M)+XH1(J) 30 CONTINUE DO 40 J=1,3 XH2(J)=0.0 DO 39 M=1,3 39 XH2(J)=BI(M,J)*A(J,M)+XH2(J) 40 CONTINUE XCOE(5)=0.0 DO 60 K=1,3 XCOE(5)=XH1(K)+4.*XH2(K)+XCOE(5) 60 CONTINUE DO 70 K=1,3 DO 70 J=1,3 ACI(J,K)=ACI(J,K)-AI(J,K)-CI(J,K) 70 CONTINUE DO 80 J=1,3 XH1(J)=0.0 DO 79 M=1,3 79 XH1(J)=ACI(M,J)*B(J,M)+XH1(J) 80 CONTINUE DO 90 J=1,3 XH2(J)=0.0 DO 89 M=1,3 89 XH2(J)=A(M,J)*B(J,M)+XH2(J) 90 CONTINUE XCOE(4)=0.0 DO 100 J=1,3 XCOE(4)=XH1(J)+XH2(J)+XCOE(4) 100 CONTINUE XCOE(4)=4.*DETB+XCOE(4)-AT*BT XCOE(4)=2.*XCOE(4) DO 110 J=1,3 XH1(J)=0.0 DO 109 M=1,3 109 XH1(J)=CDI(M,J)*A(J,M)+XH1(J) 110 CONTINUE DO 120 J=1,3 XH2(J)=0.0 DO 119 M=1,3 119 XH2(J)=BI(M,J)*CD(J,M)+XH2(J) 120 CONTINUE XCOE(3)=0.0 DO 130 J=1,3 XCOE(3)=XH1(J)+4.*XH2(J)+XCOE(3) 130 CONTINUE DO 140 J=1,3 XH1(J)=0.0 DO 139 M=1,3 139 XH1(J)=CDI(M,J)*B(J,M)+XH1(J) 140 CONTINUE XCOE(2)=XH1(1)+XH1(2)+XH1(3) XCOE(2)=2.*XCOE(2) DO 141 K=1,7 M=8-K 141 XKO(M)=XCOE(K) CALL POLRT(XCOE,COE,6,Z,IER) iww=0 142 IF(MTRNS.GT.0)then WRITE(LOUT,1000)XKO WRITE(LOUT,1010) Z 1000 FORMAT(' COEF. OF POLYNOMIAL',/,4(E12.5,4X)) 1010 FORMAT(' ROOTS OF POLYNOMIAL',/,4(E12.5,4X)) end if IF(IER.EQ.131) THEN WRITE(LOUT,1030) 1030 FORMAT(/' NOT ALL ZEROS FOUND , PROGRAM TERMINATES'//) STOP END IF IF(IER.EQ.130) THEN WRITE(LOUT,1040) 1040 FORMAT(/' DEGREE OF POLYNOMAL IS REDUCED %%%%%%%%%%%%%%%%'//) END IF DO 150 J=1,6 ZI=AIMAG(Z(J)) IF(ABS(ZI).GT.0.)THEN ISG(J)=0 GOTO 165 END IF C C CHECK OF REAL-VALUED SLOWNESS VECTORS OF GENERATED WAVES C Z0=Z(J) XPR=XKO(7)+Z0*(XKO(6)+Z0*(XKO(5)+Z0*(XKO(4)+Z0*(XKO(3)+Z0*(XKO(2)+ *Z0*XKO(1)))))) IF(MTRNS.NE.0.OR.(IWW.EQ.1.AND.CABS(XPR).GT..00001))THEN WRITE(LOUT,1041) XPR,Z0 1041 FORMAT(' PRECISSION OF ROOTS',/,4(E12.5,4X)) END IF DO 160 K=1,3 Y1(K+3)=PN(K)+REAL(Z(J))*UN(K) 160 CONTINUE NDER1=NDER NDER=1 CALL FCT(X,Y1,DERY) IF(IND.EQ.10)RETURN NDER=NDER1 AAA=Y1(4)*DERY(1)+Y1(5)*DERY(2)+Y1(6)*DERY(3) IF(MTRNS.GT.0)WRITE(LOUT,2000)AAA,Y1(4),Y1(5),Y1(6) 2000 FORMAT(1X,'GROUP*SLOWNESS, SLOWNESS'/1X,4E15.5) IF(ABS(AAA-1.).GT.1.0E-02)THEN IND=10 WRITE(LOUT,'(A)')'GROUP*SLOWNESS.NE.1' RETURN END IF SG=UN(1)*DERY(1)+UN(2)*DERY(2)+UN(3)*DERY(3) IF(MTRNS.GT.0)WRITE(LOUT,1042)SG 1042 FORMAT(' DIRECTION OF RAY VECTOR (UN(I)*DERY(I),I=1,3',F10.5) IF(MTRNS.GT.0) WRITE(LOUT,1043)(DERY(i),i=1,3) 1043 FORMAT(' RAY VECTOR (DERY)'/3F10.5) IF(SG.GE.0.)ISG(J)=1 IF(SG.LT.0.)ISG(J)=-1 165 CONTINUE 150 CONTINUE MCT=0 MCR=0 C C CHECK OF COMPLEX-VALUED SLOWNESS VECTORS OF GENERATED WAVES C DO 179 K=1,6 IF(ISG(K).NE.0) GOTO 179 DIR=AIMAG(Z(K))*UN(3) IF(UN(3).GT.0.0) GOTO 175 C C INCIDENT WAVE STRIKES THE INTERFACE FROM ABOVE C IF(DIR.LT.0.0)THEN MCR=MCR+1 ICR(MCR)=K END IF IF(DIR.GT.0.0)THEN MCT=MCT+1 ICT(MCT)=K END IF GOTO 179 C C INCIDENT WAVE STRIKES THE INTERFACE FROM BELOW C 175 IF(DIR.GT.0.0)THEN MCR=MCR+1 ICR(MCR)=K END IF IF(DIR.LT.0.0)THEN MCT=MCT+1 ICT(MCT)=K END IF 179 CONTINUE MT=0 DO 180 K=1,6 IF(ISG(K).EQ.(-1))THEN MT=MT+1 IT(MT)=K END IF 180 CONTINUE MR=0 DO 190 K=1,6 IF(ISG(K).EQ.1)THEN MR=MR+1 IR(MR)=K END IF 190 CONTINUE C IF(MTRNS.EQ.0) GOTO 191 WRITE(LOUT,1021)ISG 1021 FORMAT(' INDICATIONS,ISG(6),MCR ICR(3),MCT ICT(3),MR IR(3) 1MT IT(3) ',/,6I5) WRITE(LOUT,1020)MCR,ICR WRITE(LOUT,1020)MCT,ICT WRITE(LOUT,1020)MR,IR WRITE(LOUT,1020)MT,IT 1020 FORMAT(14I5) 191 IS(5,IRF)=MCR IS(6,IRF)=MCT IF(ITRANS.EQ.0)GO TO 200 C C TRANSMISSION C IF(MT.EQ.3) THEN Z1=REAL(Z(IT(1))) Z2=REAL(Z(IT(2))) Z3=REAL(Z(IT(3))) XSI(3)=AMAX1(Z1,Z2,Z3) XSI(1)=AMIN1(Z1,Z2,Z3) XSI(2)=Z1+Z2+Z3-XSI(3)-XSI(1) END IF IF(MT.EQ.2) THEN Z1=REAL(Z(IT(1))) Z2=REAL(Z(IT(2))) XSI(2)=AMAX1(Z1,Z2) XSI(1)=Z1+Z2-XSI(2) CSI(1)=Z(ICT(1)) END IF IF(MT.EQ.1) THEN XSI(1)=REAL(Z(IT(1))) IF(CABS(Z(ICT(1))).GT.CABS(Z(ICT(2)))) THEN CSI(1)=Z(ICT(2)) CSI(2)=Z(ICT(1)) ELSE CSI(1)=Z(ICT(1)) CSI(2)=Z(ICT(2)) END IF END IF IF(MT.EQ.0) THEN CSI(1)=CMPLX(999.,999.) CSI(3)=CMPLX(0.,0.) DO 189 K=1,3 IF(CABS(Z(ICT(K))).LT.CABS(CSI(1))) CSI(1)=Z(ICT(K)) IF(CABS(Z(ICT(K))).GT.CABS(CSI(3))) CSI(3)=Z(ICT(K)) 189 CONTINUE CSI(2)=Z(ICT(1))+Z(ICT(2))+Z(ICT(3))-CSI(1)-CSI(3) END IF IF(MDIM.LT.1) GOTO 196 C C COMPUTATION OF TRANSMITTED SLOWNESS VECTORS FOR EVALUATING C AMPLITUDES IN ROUTINE AMPL C DO 192 K=1,MT DO 193 L=1,3 PS(L,K+3,IRF)=PN(L)+XSI(K)*UN(L) 193 CONTINUE 192 CONTINUE DO 194 K=1,MCT I=ICT(K) M=MT+K DO 195 L=1,3 PS(L,M+3,IRF)=PN(L)+CSI(K)*UN(L) 195 CONTINUE 194 CONTINUE IF(IASW.EQ.0)RETURN C C OVERCRITICAL INCIDENCE C 196 IF(ITYPE.GT.MT) IND=9 IF(IND.EQ.9) RETURN C GOTO 210 C C REFLECTION C 200 IF(MR.EQ.3) THEN Z1=REAL(Z(IR(1))) Z2=REAL(Z(IR(2))) Z3=REAL(Z(IR(3))) XSI(3)=AMIN1(Z1,Z2,Z3) XSI(1)=AMAX1(Z1,Z2,Z3) XSI(2)=Z1+Z2+Z3-XSI(3)-XSI(1) END IF IF(MR.EQ.2) THEN Z1=REAL(Z(IR(1))) Z2=REAL(Z(IR(2))) XSI(2)=AMIN1(Z1,Z2) XSI(1)=Z1+Z2-XSI(2) CSI(1)=Z(ICR(1)) END IF IF(MR.EQ.1) THEN XSI(1)=REAL(Z(IR(1))) IF(CABS(Z(ICR(1))).GT.CABS(Z(ICR(2)))) THEN CSI(1)=Z(ICR(2)) CSI(2)=Z(ICR(1)) ELSE CSI(1)=Z(ICR(1)) CSI(2)=Z(ICR(2)) END IF END IF IF(MR.EQ.0) THEN CSI(1)=CMPLX(0.0) CSI(3)=CMPLX(999.,999.) DO 211 K=1,3 IF(CABS(Z(ICR(K))).GT.CABS(CSI(1))) CSI(3)=Z(ICR(K)) IF(CABS(Z(ICR(K))).LT.CABS(CSI(3))) CSI(1)=Z(ICR(K)) 211 CONTINUE CSI(2)=Z(ICR(1))+Z(ICR(2))+Z(ICR(3))-CSI(1)-CSI(3) END IF C C COMPUTATION OF REFLECTED SLOWNESS VECTORS FOR EVALUATING C AMPLITUDES IN ROUTINE AMPL C IF(MDIM.LT.1) GOTO 209 DO 201 K=1,MR DO 202 L=1,3 PS(L,K,IRF)=PN(L)+XSI(K)*UN(L) 202 CONTINUE 201 CONTINUE DO 203 K=1,MCR I=ICR(K) M=MR+K DO 204 L=1,3 PS(L,M,IRF)=PN(L)+CSI(K)*UN(L) 204 CONTINUE 203 CONTINUE IF(IASW.EQ.0) RETURN C C OVERCRITICAL INCIDENCE C 209 IF(ITYPE.GT.MR) IND=9 IF(IND.EQ.9) RETURN C 210 IF(MTRNS.GT.0) WRITE(LOUT,1111) XSI,Z1,Z2,Z3,ITYPE 1111 FORMAT(' XSI,Z1,Z2,Z3,ITYPE',/,6E14.6,I10) XSSI=XSI(ITYPE) PU1=XSSI*UN(1) PU2=XSSI*UN(2) PU3=XSSI*UN(3) PNEW(1)=PN(1)+PU1 PNEW(2)=PN(2)+PU2 PNEW(3)=PN(3)+PU3 Y(4)=PNEW(1) Y(5)=PNEW(2) Y(6)=PNEW(3) NDER1=NDER NDER=1 c write(lout,'(i5,3f15.10)')itype,y1(4),y1(5),y1(6) CALL FCT(X,Y,DERY) IF(IND.EQ.10)RETURN NDER=NDER1 DUN=DERY(1)*UN(1)+DERY(2)*UN(2)+DERY(3)*UN(3) RHO2=0.2*SQRT(A11)+1.7 IF(IRHO.NE.0) RHO2=RHO(LAY) DS(10,IRF)=abs(dun) DS(11,IRF)=RHO2 IF(MTRNS.GT.0)WRITE(LOUT,1100)PNEW 1100 FORMAT(' PNEW',/,3F12.8) RETURN C C COMPUTATION OF THE TRANSFORMATION OF THE SLOWNESS VECTOR C IN THE ISOTROPIC CASE C 230 IF(ITYPE.EQ.3)CI2=1./A11 IF(ITYPE.NE.3)CI2=1./A44 PP=0.0 DO 250 K=1,3 PP=PP+P(K)*P(K) 250 CONTINUE PUN2=PUN*PUN ROO=CI2-PP+PUN2 ROT1=1./A11-PP+PUN2 ROT2=1./A44-PP+PUN2 IF(MTRNS.GT.0) WRITE(LOUT,1065) ROO,ROT1,ROT2 1065 FORMAT('ROO,ROT1,ROT2',3F10.5) C C ROT1.LT.ROT2 FOR FOR REALISTIC MATERIALS C IF(ITRANS.EQ.0.AND.MDIM.GE.1) THEN IF(ROT1.GT.0.AND.ROT2.GT.0) THEN MR=3 MCR=0 END IF IF(ROT1.LE.0.AND.ROT2.GT.0) THEN MR=2 MCR=1 END IF IF(ROT1.LE.0.AND.ROT2.LE.0) THEN MR=0 MCR=3 END IF IS(5,IRF)=MCR IF(MTRNS.GT.0) WRITE(LOUT,1066)MR,MCR 1066 FORMAT(' MR,MCR',/,2I5) DO 251 K=1,MR IF(K+MCR.EQ.1) RO=ROT1 IF(K+MCR.GT.1) RO=ROT2 SQU=SQRT(RO) IF(MTRNS.GT.0) WRITE(LOUT,1062) SQU 1062 FORMAT(' SQU',/,2F12.6) IS(5,IRF)=MCR J=MR+1-K DO 252 L=1,3 PS(L,J,IRF)=P(L)-(PUN-SQU)*UN(L) 252 CONTINUE 251 CONTINUE DO 253 K=1,MCR M=MR+K IF(M.EQ.3) CSQ=CSQRT(CMPLX(ROT1,0.0)) IF(M.NE.3) CSQ=CSQRT(CMPLX(ROT2,0.0)) IF(MTRNS.GT.0) WRITE(LOUT,1064)CSQ 1064 FORMAT(' CSQ',2F10.5) DO 254 L=1,3 PS(L,M,IRF)=P(L)-(PUN-CSQ)*UN(L) 254 CONTINUE 253 CONTINUE IF(IASW.EQ.0) RETURN END IF IF(ITRANS.EQ.1.AND.MDIM.GE.1) THEN IF(ROT1.GT.0.AND.ROT2.GT.0) THEN MT=3 MCT=0 END IF IF(ROT1.LE.0.AND.ROT2.GT.0) THEN MT=2 MCT=1 END IF IF(ROT1.LE.0.AND.ROT2.LE.0) THEN MCT=3 MT=0 END IF IS(6,IRF)=MCT IF(MTRNS.GT.0) WRITE(LOUT,1061)MT,MCT 1061 FORMAT(' MT,MCT',/,2I5) DO 255 K=1,MT IF(K+MCT.EQ.1) RO=ROT1 IF(K+MCT.NE.1) RO=ROT2 SQU=SQRT(RO) IF(MTRNS.GT.0) WRITE(LOUT,1062)SQU J=MT+1-K DO 256 L=1,3 PS(L,J+3,IRF)=P(L)-(PUN+SQU)*UN(L) 256 CONTINUE 255 CONTINUE DO 257 K=1,MCT M=MT+K IF(M.EQ.3)CSQ=CSQRT(CMPLX(ROT1,0.0)) IF(M.NE.3)CSQ=CSQRT(CMPLX(ROT2,0.0)) IF(MTRNS.GT.0) WRITE(LOUT,1064)CSQ DO 258 L=1,3 PS(L,M+3,IRF)=P(L)-(PUN+CSQ)*UN(L) 258 CONTINUE 257 CONTINUE IF(IASW.EQ.0) RETURN END IF IF(ROO.LE.0.0) GOTO 300 SQU=SQRT(ROO) IF(ITRANS.EQ.1) GOTO 280 C C REFLECTION C PU1=PUN-SQU DO 260 K=1,3 PNEW(K)=P(K)-PU1*UN(K) 260 CONTINUE GOTO 275 C C TRANSMISSION C 280 PU1=PUN+SQU DO 270 K=1,3 PNEW(K)=P(K)-PU1*UN(K) 270 CONTINUE 275 IF(MTRNS.GT.0) WRITE(LOUT,1063)PNEW 1063 FORMAT(' PNEW',/,3F12.6) PP=0.0 DO 290 K=1,3 PP=PP+PNEW(K)*PNEW(K) 290 CONTINUE RHO2=0.2*SQRT(A11)+1.7 IF(IRHO.NE.0) RHO2=RHO(LAY) PUNG=PNEW(1)*UN(1)+PNEW(2)*UN(2)+PNEW(3)*UN(3) DS(10,IRF)=ABS(PUNG)/PP DS(11,IRF)=RHO2 RETURN 300 IND=9 NTR=102 RETURN END Canray.for 0100666 0000765 0000765 00000106312 12154611062 012254 0 ustar bulant bulant C
C PROGRAM A N R A Y, VERSION 4.73 (PRAHA, JUNE 2013) C C******************************************************************* C C PROGRAM ANRAY IS DESIGNED FOR RAY, TRAVEL TIME AND C AMPLITUDE COMPUTATIONS IN 3D GENERAL ANISOTROPIC AND ISOTROPIC C LATERALLY VARYING LAYERED MEDIA. THE PROGRAM MAKES POSSIBLE C COMPUTATION OF RAYS SPECIFIED BY INITIAL ANGLES AT THE SOURCE, C I.E., INITIAL-VALUE RAY TRACING, OR RAYS STARTING FROM THE C SOURCE AND TERMINATING ON A VERTICAL OR SURFACE PROFILE, I.E. C BOUNDARY-VALUE RAY TRACING. RAY AMPLITUDES CAN BE COMPUTED C ALONG RAYS. C C******************************************************************* C C CHARACTER*80 MTEXT,FILEIN,FILEOU,FILE1,FILE2,FILE3,FILE4,FILE5 CHARACTER*80 FILE6,FILE7 DIMENSION Y(18) 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 /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33, 1 A34,A35,A36,A44,A45,A46,A55,A56,A66, 1 DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23, 1 DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44, 1 DXA45,DXA46,DXA55,DXA56,DXA66, 1 DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23, 1 DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44, 1 DYA45,DYA46,DYA55,DYA56,DYA66, 1 DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23, 1 DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44, 1 DZA45,DZA46,DZA55,DZA56,DZA66, 1 A2546,A1266,A1355,A1456,A3645,A2344 INTEGER CODE COMMON /COD/ CODE(50,2),KREF,KC,ITYPE COMMON /DIST/ DST(200),NDST,REPS,PROF(2),NDSTP,PREPS,LNDST, 1XPRF,YPRF,ILOC COMMON /DENS/ RHO(20) COMPLEX PS COMMON /RAY/ AY(28,2000),DS(20,50),KINT(50),HHH(3,3),TMAX, 1 PS(3,7,50),IS(8,50),N,IREF,IND,IND1 COMMON /INTRF/ Z(1000),SX(350),SY(350),NX(20),NY(20),BRD(6),NINT, 1 XINTA COMMON /ZERO/ RNULL COMMON/VSP/XVSP,YVSP,XNRM,YNRM,ICOD,IVSP COMMON/VRML/LUBRD,LUGRD,LUIND,LURAY C C************************************************** C LIN=5 LOU=6 LU1=1 LU2=2 LU3=3 LUBRD=7 LUGRD=8 LUIND=9 LURAY=10 FILEIN='anray.dat' FILEOU='anray.out' FILE1='lu1.dat' FILE2='lu2.dat' FILE3=' ' FILE4=' ' FILE5=' ' FILE6=' ' FILE7=' ' WRITE(*,'(2A)') ' (ANRAY) SPECIFY NAMES OF INPUT AND OUTPUT', 1' FILES LIN, LOU, LU1, LU2, LU3, LUBRD, LUGRD, LUIND, LURAY: ' READ(*,*) FILEIN,FILEOU,FILE1,FILE2,FILE3,FILE4,FILE5,FILE6,FILE7 IF(FILE1.EQ.' ') LU1=0 IF(FILE2.EQ.' ') LU2=0 IF(FILE3.EQ.' ') LU3=0 IF(FILE4.EQ.' ') LUBRD=0 IF(FILE5.EQ.' ') LUGRD=0 IF(FILE6.EQ.' ') LUIND=0 IF(FILE7.EQ.' ') LURAY=0 LOUT=LOU OPEN(LIN,FILE=FILEIN,FORM='FORMATTED',STATUS='OLD') OPEN(LOU,FILE=FILEOU,FORM='FORMATTED') IF(LU1.NE.0)OPEN(LU1,FILE=FILE1,FORM='FORMATTED') IF(LU2.NE.0)OPEN(LU2,FILE=FILE2,FORM='FORMATTED') IF(LU3.NE.0)OPEN(LU3,FILE=FILE3,FORM='FORMATTED') IF(LUBRD.NE.0)OPEN(LUBRD,FILE=FILE4,FORM='FORMATTED') IF(LUGRD.NE.0)OPEN(LUGRD,FILE=FILE5,FORM='FORMATTED') IF(LUIND.NE.0)OPEN(LUIND,FILE=FILE6,FORM='FORMATTED') IF(LURAY.NE.0)OPEN(LURAY,FILE=FILE7,FORM='FORMATTED') C C************************************************** C WRITE(LOU,777) 777 FORMAT(///,'***********************' 1,//,' PROGRAM A N R A Y ',//, 2'***********************',//) NCODE=1 MTEXT='ANRAY' INUL=4 READ(LIN,*)MTEXT WRITE(LOU,115)MTEXT READ(LIN,*)INULL,ISURF IF(INULL.EQ.0)INULL=INUL RNULL=10.**(-INULL) WRITE(LOU,106)INULL,ISURF C C C SPECIFICATION OF THE MODEL C CALL MODEL(MTEXT,LIN) C C GENERATE FILE FOR PLOTTING VARIOUS CHARACTERISTIC SURFACES C IF(LU3.NE.0)CALL SURFPL(LIN,LU3) C C GENERATE FILE FOR VRML PLOTTING BOUNDARIES OF THE MODEL C IF(LUBRD.NE.0)CALL BOX(BRD) C C GENERATE FILE FOR PLOTTING RAYS C IF(LURAY.NE.0)WRITE(LURAY,113) IF(LURAY.NE.0)WRITE(LURAY,105) C C SPECIFICATION OF SYNTHETIC SEISMOGRAMS C 2 ICONT=1 MEP=0 MOUT=0 MDIM=0 METHOD=0 MREG=0 ITMAX=10 IPOL=0 IPREC=0 IRAYPL=0 IPRINT=0 IAMP=0 MTRNS=0 ICOEF=0 IRT=0 ILOC=0 MCOD=0 MORI=0 READ(LIN,*)ICONT,MEP,MOUT,MDIM,METHOD,MREG,ITMAX, 1IPOL,IPREC,IRAYPL,IPRINT,IAMP,MTRNS,ICOEF,IRT,ILOC,MCOD,MORI WRITE(LOU,102)ICONT,MEP,MOUT,MDIM,METHOD,MREG,ITMAX, 1IPOL,IPREC,IRAYPL,IPRINT,IAMP,MTRNS,ICOEF,IRT,ILOC,MCOD,MORI IF(ICONT.EQ.0)GO TO 99 C C c IF(MEP.NE.0.AND.MDIM.EQ.0)MDIM=1 IVSP=0 IF(ILOC.EQ.0)ITPR=3 IF(ILOC.EQ.1)THEN IVSP=1 ITPR=43 MREG=1 END IF IF(ILOC.GT.1)THEN ITPR=ILOC+100 END IF C IF(MEP.EQ.0)THEN NDST=0 END IF C IF(MEP.EQ.1)THEN NDST=1 READ(LIN,*)XREC,YREC WRITE(LOU,104)XREC,YREC GO TO 4 END IF IF(MEP.LT.0)THEN NDST=-MEP PROF(1)=0. XPRF=0. YPRF=0. READ(LIN,*)PROF(1),(DST(I),I=1,NDST),XPRF,YPRF WRITE(LOU,104)PROF(1),(DST(I),I=1,NDST),XPRF,YPRF IF(NDST.EQ.1)RSTEP=1. IF(NDST.EQ.1)DST(2)=DST(1)+1. IF(NDST.EQ.1)GO TO 4 RSTEP=(DST(NDST)-DST(1))/FLOAT(NDST-1) END IF C IF(MEP.GT.0)THEN NDST=MEP READ(LIN,*)PROF(1),RMIN,RSTEP,XPRF,YPRF WRITE(LOU,104)PROF(1),RMIN,RSTEP,XPRF,YPRF DO 13 I=1,MEP 13 DST(I)=RMIN+(I-1)*RSTEP IF(NDST.EQ.1)DST(2)=RMIN+RSTEP END IF PROF(2)=PROF(1)+1. NDSTP=1 C IF(IVSP.EQ.1.AND.NDST.NE.0)THEN READ(LIN,*)XVSP,YVSP WRITE(LOU,104)XVSP,YVSP END IF C 4 TSOUR=0. DT=1. AC=0.0001 REPS=0.05 PREPS=0.05 READ(LIN,*)XSOUR,YSOUR,ZSOUR,TSOUR,DT,AC,REPS,PREPS WRITE(LOU,104)XSOUR,YSOUR,ZSOUR,TSOUR,DT,AC,REPS,PREPS C IF(ABS(XPRF).LT..000001.AND.ABS(YPRF).LT..000001)THEN XPRF=XSOUR YPRF=YSOUR END IF IF(MEP.EQ.1)THEN XE=XREC-XPRF YE=YREC-YPRF RPRF=SQRT(XE*XE+YE*YE) XATAN=ATAN2(YE,XE) PROF(1)=XATAN RMIN=RPRF DST(1)=RMIN WRITE(LOU,104)RPRF,XATAN RSTEP=100. DST(2)=DST(1)+100. PROF(2)=PROF(1)+1. NDSTP=1 END IF C IF(IVSP.EQ.1.AND.NDST.NE.0)THEN XNRM=XVSP-XSOUR YNRM=YVSP-YSOUR AUX=SQRT(XNRM*XNRM+YNRM*YNRM) XNRM=XNRM/AUX YNRM=YNRM/AUX PROF(1)=ATAN2(YNRM,XNRM) PROF(2)=PROF(1)+1. XPRF=XSOUR YPRF=YSOUR END IF IF(MCOD.EQ.0)THEN READ(LIN,*)AMIN,ASTEP,AMAX WRITE(LOU,104)AMIN,ASTEP,AMAX READ(LIN,*)BMIN,BSTEP,BMAX IF(ABS(BSTEP).LT..000001)THEN BMIN=PROF(1)-.3 BMAX=PROF(1)+.4 BSTEP=.6 END IF WRITE(LOU,104)BMIN,BSTEP,BMAX END IF IF((MREG.EQ.0.OR.MREG.EQ.2).AND.MDIM.NE.0) WRITE(LOU,'(/,A,/)') 1 ' COEFFICIENTS OF CONVERSION ARE APPLIED' IF((MREG.NE.0.AND.MREG.NE.2).AND.MDIM.NE.0) WRITE(LOU,'(/,A,/)') 1 ' COEFFICIENTS OF CONVERSION ARE *** NOT *** APPLIED' TMAX=10000. IND=-1 NDER=1 CALL RAYA(XSOUR,YSOUR,ZSOUR,TSOUR,AMIN1,BMIN,PX,PY,PZ,XX,YY,ZZ,T, 1DT,AC) Y(1)=XSOUR Y(2)=YSOUR Y(3)=ZSOUR IF(IND.EQ.50)WRITE(LOU,111)IND IF(IND.EQ.50)GO TO 99 LAY=IND ISOUR=IND ITYPE=3 CALL PARDIS(Y,0) VP=SQRT(A11) IF(IRHO.EQ.0)RO=1.7+.2*VP IF(IRHO.EQ.1)RO=RHO(IND) C C GENERATE FILE LU2 FOR SYNTHETIC SEISMOGRAM COMPUTATIONS C IF(LU2.NE.0.AND.NDST.NE.0)THEN WRITE(LU2,115)MTEXT KSH=2 WRITE(LU2,100)NDST,KSH,ILOC WRITE(LU2,104)XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,RO IF(MEP.NE.1)WRITE(LU2,104)(DST(I),I=1,NDST) IF(MEP.EQ.1)WRITE(LU2,104)XREC,YREC,RPRF,XATAN END IF C C LOOP FOR ELEMENTARY WAVES C 20 READ(LIN,*)KC,KREF,((CODE(I,K),K=1,2),I=1,KREF) WRITE(LOU,100)KC,KREF,((CODE(I,K),K=1,2),I=1,KREF) IF(KREF.EQ.0)GOTO 2 IF(MOUT.NE.0)WRITE(LOU,107) WRITE(LOU,103)NCODE,KC,KREF,((CODE(I,K),K=1,2),I=1,KREF) C IF(MCOD.NE.0)THEN READ(LIN,*)AMIN,ASTEP,AMAX WRITE(LOU,104)AMIN,ASTEP,AMAX READ(LIN,*)BMIN,BSTEP,BMAX IF(ABS(BSTEP).LT..000001)THEN BMIN=PROF(1)-.3 BMAX=PROF(1)+.4 BSTEP=.6 END IF WRITE(LOU,104)BMIN,BSTEP,BMAX END IF C C GENERATE FILE LU1 FOR PLOTTING OF RAY DIAGRAMS, C TIME-DISTANCE AND AMPLITUDE-DISTANCE CURVES C IF(LU1.EQ.0.OR.NDST.EQ.0)GO TO 21 WRITE(LU1,100)ICONT,NDST,ILOC WRITE(LU1,104)RO NPN=2 APN=0. WRITE(LU1,100)NPN,NPN,NPN WRITE(LU1,101)APN,APN,APN,APN,APN WRITE(LU1,101)APN,APN,APN,APN,APN WRITE(LU1,104)Xprf,Yprf,0.0,PROF(1) WRITE(LU1,104)(DST(I),I=1,NDST) 21 CONTINUE C C C SEARCH FOR THE NUMBER OF THE ELEMENT OF THE RAY, STARTING FROM C WHICH THE WAVE DOES UNDERTAKE NEITHER REFLECTION NOR CONVERSION C ICOD=0 IF(IVSP.EQ.0)GO TO 35 DO 34 I=1,KREF ICOD=KREF-I+1 IF(ICOD.EQ.1) GO TO 34 IC1=CODE(ICOD,1) IC2=CODE(ICOD-1,1) IF((IC1-IC2).EQ.0)GO TO 35 IC1=CODE(ICOD,2) IC2=CODE(ICOD-1,2) IF((IC1-IC2).NE.0)GO TO 35 34 CONTINUE 35 CONTINUE IF(MOUT.NE.0)WRITE(LOU,108) C C CALL RECEIV(XSOUR,YSOUR,ZSOUR,TSOUR,DT,AC,ITMAX,AMIN,ASTEP, 1AMAX,BMIN,BSTEP,BMAX,MOUT,LU1,LU2,METHOD,ITPR,NCODE) IF(IND.EQ.14) WRITE(LOU,111) IND NCODE=NCODE+1 GOTO 20 C C END OF LOOP FOR ELEMENTARY WAVES C C 100 FORMAT(26I3) 101 FORMAT(5E15.5) 102 FORMAT(1H0,////,2X,26I3) 103 FORMAT(4X,I4,9X,100I3) 104 FORMAT(8F10.5) 105 FORMAT('/') 106 FORMAT(17I5) 107 FORMAT(//2X,'INT.CODE',5X,'E X T E R N A L C O D E') 108 FORMAT(//) 111 FORMAT(/2X,'IND=',I5,/) 113 FORMAT(6H'RAYS') 115 FORMAT(A) C 99 CONTINUE IF(LURAY.NE.0)WRITE(LURAY,105) IF(LU1.NE.0.AND.NDST.NE.0)WRITE(LU1,100)ICONT,ICONT IF(LU1.NE.0)REWIND LU1 IF(LU2.NE.0)REWIND LU2 C STOP END C C ********************************************************* C SUBROUTINE AMPL (AMPX,AMPY,AMPZ,UU) C C ROUTINE FOR COMPUTING COMPLEX VECTORIAL RAY AMPLITUDES C C OUTPUT PARAMETERS C AMPX(2),AMPY(2),AMPZ(2) - X,Y AND Z COMPONENTS OF COMPLEX C VECTORIAL RAY AMPLITUDES IN THE MODEL COORDINATES. FOR P WAVE C IN ANY MEDIUM AND FOR S WAVES IN AN ANISOTROPIC MEDIUM, I=1. C FOR S WAVE GENERATED IN AN ISOTROPIC MEDIUM, I=1,2. I=1 AND 2 C CORRESPOND TO S WAVES SPECIFIED AT THE SOURCE BY VECTORS E1 C E2. VECTORS E1 AND E2 TOGETHER WITH UNIT VECTOR TANGENT TO C THE RAY FORM A BASIS OF RAY CENTRED COORDINATE SYSTEM. C UU - PRODUCT OF RATIOS OF DENSITIES AND COSINES OF INCIDENCE C AND OF REFLECTION/TRANSMISSION AT POINTS WHERE THE RAY CROSSES c INTERFACES. C C CALLED FROM: RECEIV C ROUTINES CALLED: POLAR,TRANSL,COEF C DIMENSION Y(18),UN(3),POLD(3),PNEW(3) COMPLEX AMPX(2),AMPY(2),AMPZ(2),CR(3),UC(3),STU(6),C1,C2,C3 COMMON /AUXI/ IANI(20),INTR,INT1,IPREC,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori COMMON /DIST/ DST(200),NDST,REPS,PROF(2),NDSTP,PREPS,LNDST, 1XPRF,YPRF,ILOC INTEGER CODE COMMON /COD/ CODE(50,2),KREF,KC,ITYPE COMMON /DENS/ RHO(20) COMPLEX PS COMMON /RAY/ AY(28,2000),DS(20,50),KINT(50),HHH(3,3),tmax, 1 PS(3,7,50),IS(8,50),N,IREF,IND,IND1 COMMON /RAY2/ DRY(3,2000) C KSS=1 ISHEAR=0 ITYPE=CODE(1,2) IF(IANI(ISOUR).EQ.0.AND.ITYPE.NE.3)THEN ISHEAR=1 ITYPE=1 END IF ITP=ITYPE DO 1 I=1,2 AMPX(I)=CMPLX(0.,0.) AMPY(I)=CMPLX(0.,0.) AMPZ(I)=CMPLX(0.,0.) 1 CONTINUE C 3000 NN=N IDD=0 N2=0 N1=1 IRE=IREF AV=1. C C SPECIFICATION OF DISPLACEMENT VECTOR AT SOURCE C IN RAY CENTERED COORDINATES C DO 5 I=1,3 CR(I)=(0.,0.) 5 CONTINUE CR(ITP)=(1.,0.) IREF1=IREF-1 IF(IRE.GT.1)INAUM=CODE(IRE-1,1)-CODE(IRE,1) IF(MREG.GE.1.AND.IRE.GT.1.AND.ILOC.GT.1.AND.INAUM.GE.0)THEN IREF1=IREF1+1 CODE(IREF1+1,2)=3 END IF IF(IREF1.EQ.0) GOTO 100 C C LOOP OVER INTERFACES C DO 10 I=1,IREF1 IREF=I IF(KC.NE.0) ITYPE=CODE(IREF,2) N=KINT(IREF) IF(N.EQ.0) THEN IDD=1 GO TO 10 ELSE N1=N2+1 N2=N IF(IDD.NE.0) N2=-N2 IDD=0 C C COMPUTATION OF POLARIZATION VECTORS C CONSIDERED POLARIZATION VECTOR(S) ARE STORED IN CORRESPONDING C COLUMNS OF THE MATRIX HHH. OTHER COLUMNS ARE ZERO. C CALL POLAR(N1,N2,NN,IREF) END IF DO 20 K=1,6 Y(K)=AY(K+1,N) 20 CONTINUE IF(IAMP.GT.0)WRITE(LOU,'(a,2i5,6f10.5)')' AMPL:I,N,Y',I,N, 1(Y(L),L=1,6) DO 30 K=1,3 POLD(K)=Y(K+3) PS(K,7,IREF)=Y(K+3) 30 CONTINUE DO 40 K=1,3 UN(K)=DS(K,IREF) 40 CONTINUE LAY=IS(1,IREF) ITRANS=IS(2,IREF) ITR1=ITRANS IF(UN(3).GT.0.0) GOTO 50 C C RAY STRIKING THE INTERFACE FROM ABOVE C IF(ITRANS.EQ.0) THEN LAY=LAY+1 ITRANS=1 GOTO 70 END IF IF(ITRANS.GT.0) THEN LAY=LAY-1 ITRANS=0 GOTO 70 END IF C C RAY STRIKING THE INTERFACE FROM BELOW C 50 IF(ITRANS.EQ.0) THEN LAY=LAY-1 ITRANS=1 GOTO 70 END IF IF(ITRANS.GT.0) THEN LAY=LAY+1 ITRANS=0 GOTO 70 END IF C C SLOWNESS VECTORS ON THE SIDE OF THE INTERFACE WHERE GENERATED C WAVE PROPAGATES WERE DETERMINED DURING THE CALL OF TRANSL IN THE C ROUTINE OUT. HERE REMAINING SLOWNESS VECTORS ON THE OTHER SIDE C OF THE INTERFACE ARE DETERMINED C C REDEFINITION OF IREF FOR CALL OF ROUTINE TRANSL C 70 IF(LAY.EQ.0) THEN DO 71 K=4,6 DO 71 L=1,3 PS(L,K,IREF)=CMPLX(0.,0.) 71 CONTINUE GO TO 75 END IF IREF=IREF+1 CALL TRANSL(Y,POLD,PNEW,UN,ITRANS,0) IF(IND.EQ.10)RETURN IREF=IREF-1 75 IF(IAMP.NE.0)THEN WRITE(LOU,'(A)')' REFLECTED/TRANSMITTED SLOWNESS VECTORS' WRITE(LOU,'(6F12.6)')((PS(L,K,IREF),L=1,3),K=1,6) END IF AV1=(DS(11,IREF)*DS(10,IREF))/(DS(8,IREF)*DS(7,IREF)) AV=AV*AV1 IF(IAMP.GT.0) THEN WRITE(LOU,'(A)') 'ROI,ROG,UNVGI,UNVGG,AV1,AV' WRITE(LOU,'(6F10.5)') DS(8,IREF), 1 DS(11,IREF),DS(7,IREF),DS(10,IREF),AV1,AV WRITE(LOU,'(A,/,6F12.5,/,3(3F12.5/))') ' CR,HHH', 2 CR,((HHH(J,K),J=1,3),K=1,3) END IF C C COMPUTATION OF AMPLITUDE COEFFICIENTS OF REFLECTED/TRANSMITTED WAVES C C C COMPUTATION OF CARTESIAN COMPONENTS OF INCIDENT DISPLACEMENT VECTOR C DO 87 K=1,3 STU(K)=CMPLX(0.,0.) DO 87 J=1,3 STU(K)=HHH(J,K)*CR(J)+STU(K) 87 CONTINUE IF(IAMP.GT.0)WRITE(LOU,'(A,6F10.5)') ' STU',(STU(K),K=1,3) IF(KC.NE.0)ITYPE=CODE(IREF+1,2) IF(MREG.GE.1.AND.IRE.GT.1.AND.I.EQ.IREF1.AND.ILOC.GT.1.AND. 1(CODE(IRE,1).LE.CODE(IRE-1,1)))ITR1=1 CALL COEF(STU,CR,ITR1) IF(IND.EQ.11)RETURN BCR=SQRT(REAL(CR(1)*CONJG(CR(1))+CR(2)*CONJG(CR(2)) 1 +CR(3)*CONJG(CR(3)))) IF(BCR.LT.1.E-10) THEN DO 88 K=1,3 UC(K)=(0.,0.) 88 CONTINUE GOTO 130 END IF 10 CONTINUE C C END OF LOOP OVER INTERFACES C C TERMINATION POINT C 100 CONTINUE IF(IRE.GT.1)INAUM=CODE(IRE-1,1)-CODE(IRE,1) IF((MREG.GE.1.AND.IRE.GT.1).AND.ILOC.GT.1.AND.INAUM.GE.0)THEN DO 200 K=1,3 Y(K+3)=REAL(PS(K,6,IREF1)) 200 CONTINUE V=1./SQRT(Y(4)*Y(4)+Y(5)*Y(5)+Y(6)*Y(6)) DO 201 K=1,3 HHH(1,K)=0. HHH(2,K)=0. HHH(3,K)=V*Y(K+3) 201 CONTINUE ELSE N1=N2+1 N2=NN IF(KC.NE.0)ITYPE=CODE(IRE,2) IF(KC.NE.0)IS(7,IRE)=CODE(IRE,1) CALL POLAR(N1,N2,NN,IRE) END IF C C COMPUTATION OF CARTESIAN COMPONENTS OF INCIDENT DISPLACEMENT VECTOR C DO 107 K=1,3 STU(K)=CMPLX(0.,0.) DO 107 J=1,3 STU(K)=HHH(J,K)*CR(J)+STU(K) 107 CONTINUE IF(IAMP.GT.0)WRITE(LOU,'(A,6F10.5)') ' STU',(STU(K),K=1,3) C IF(IRE.GT.1)INAUM=CODE(IRE-1,1)-CODE(IRE,1) IF(MREG.EQ.1.OR.MREG.EQ.3.OR. 1(MREG.EQ.2.AND.IRE.GT.1.AND.INAUM.GE.0))THEN UC(1)=STU(1) UC(2)=STU(2) UC(3)=STU(3) IF(MREG.GT.1) THEN C C CALCULATION OF PRESSURE AT THE TERMINATION POINT C C1=UC(1) C2=UC(2) C3=UC(3) ARE=REAL(C1) IF(ARE.LT.0.)ARE=-ARE AIM=AIMAG(C1) APHI=ATAN2(AIM,ARE) ARE=SQRT(REAL(C1*CONJG(C1)+C2*CONJG(C2)+C3*CONJG(C3))) UC(1)=ARE*CMPLX(COS(APHI),SIN(APHI)) UC(2)=(0.,0.) UC(3)=(0.,0.) IF(IAMP.GT.0)WRITE(LOU,'(A,4F10.5)') ' UC(1),ARE,APHI', 1 UC(1),ARE,APHI END IF GOTO 110 END IF DO 105 K=1,6 Y(K)=AY(K+1,NN) IF(K.LE.3)GO TO 105 PS(K-3,7,IRE)=Y(K) POLD(K-3)=Y(K) UN(K-3)=DS(K-3,IRE) 105 CONTINUE N=NN IF(MREG.EQ.0.OR.MREG.EQ.2) THEN IREF=IREF+1 IF(INTR.EQ.LAY)LAY=LAY-1 IF(INTR.NE.LAY)LAY=LAY+1 CALL TRANSL(Y,POLD,PNEW,UN,1,0) END IF IREF=IRE IF(IAMP.GT.0)THEN WRITE(LOU,'(A)') 1 ' REFLECTED SLOWNESS VECTORS AT TERMINATION POINT' WRITE(LOU,'(6F12.6)')((PS(L,K,IRE),L=1,3),K=1,3) END IF C C COMPUTATION OF CONVERSION COEFFICIENTS C KTR=999 CALL COEF(STU,UC,KTR) IF(IND.EQ.11)RETURN 110 CONTINUE DO 115 K=1,3 Y(K)=AY(K+4,NN) 115 CONTINUE VPEND=1./SQRT(Y(1)*Y(1)+Y(2)*Y(2)+Y(3)*Y(3)) IF(IRE.GT.1)INAUM=CODE(IRE-1,1)-CODE(IRE,1) IF((MREG.GE.1.AND.IRE.GT.1).AND.ILOC.GT.1.AND. 1INAUM.GE.0)VPEND=V DO 120 K=1,3 Y(K)=AY(K+4,1) 120 CONTINUE VP0=1./SQRT(Y(1)*Y(1)+Y(2)*Y(2)+Y(3)*Y(3)) RHO0=0.2*SQRT(AY(8,1))+1.7 IF(IRHO.NE.0) RHO0=RHO(ISOUR) RHEND=0.2*SQRT(AY(8,NN))+1.7 IF(IRHO.NE.0) RHEND=RHO(LAY) AV=AV*VP0*RHO0 AV=AV/(VPEND*RHEND) UU=SQRT(ABS(AV)) IF(IAMP.GT.0) 1WRITE(LOU,'(A,4F12.6)')'VP0,RH0,VPEND,RHEND',VP0,RHO0,VPEND,RHEND 130 CONTINUE N=NN IREF=IRE AMPX(KSS)=UC(1) AMPY(KSS)=UC(2) AMPZ(KSS)=UC(3) IF(MREG.GT.1)AMPX(KSS)=AMPX(KSS)*VPEND*RHEND IF(ISHEAR.NE.0.AND.KSS.NE.2) THEN KSS=2 ITP=2 GOTO 3000 END IF RETURN END C C ********************************************************* C SUBROUTINE APPROX(X,Y,YD,KDIM) C C THE ROUTINE PERFORMS THIRD-ORDER INTERPOLATION BETWEEN POINTS C YOLD AND YNEW PARAMETERIZED BY AN INDEPENDENT VARIABLE X. C DOLD, DNEW ARE THE FIRST DERIVATIVES OF Y WITH RESPECT C TO X AT THE POINTS YOLD AND YNEW. C DIMENSION Y(18),YD(18) COMMON/APPR/ XOLD,XNEW,YOLD(18),DOLD(18),YNEW(18),DNEW(18) C A=(X-XNEW)/(XNEW-XOLD) AUX=A+1. A1=(2.*A+3.)*A*A A2=1.-A1 B1=AUX*A*(X-XNEW) B2=AUX*A*(X-XOLD) AD1=6.*A*AUX/(XNEW-XOLD) AD2=-AD1 BD1=A*(3.*A+2.) BD2=AUX*(3.*A+1.) DO 1 I=1,KDIM Y(I)=A1*YOLD(I)+A2*YNEW(I)+B1*DOLD(I)+B2*DNEW(I) YD(I)=AD1*YOLD(I)+AD2*YNEW(I)+BD1*DOLD(I)+BD2*DNEW(I) 1 CONTINUE 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 ********************************************************* C SUBROUTINE CHRM(Y) C C ROUTINE FOR THE COMPUTATION OF THE ELEMENTS OF THE CHRISTOFFEL C MATRIX FOR AN ARBITRARY ANISOTROPIC MEDIUM C DIMENSION Y(18) COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33, 1 A34,A35,A36,A44,A45,A46,A55,A56,A66, 1 DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23, 1 DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44, 1 DXA45,DXA46,DXA55,DXA56,DXA66, 1 DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23, 1 DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44, 1 DYA45,DYA46,DYA55,DYA56,DYA66, 1 DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23, 1 DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44, 1 DZA45,DZA46,DZA55,DZA56,DZA66, 1 A2546,A1266,A1355,A1456,A3645,A2344 COMPLEX PS COMMON /RAY/ AY(28,2000),DS(20,50),KINT(50),HHH(3,3),tmax, 1 PS(3,7,50),IS(8,50),N,IREF,IND,IND1 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 INTEGER CODE COMMON /COD/ CODE(50,2),KREF,KC,ITYPE COMMON /DJK/ D11,D12,D13,D22,D23,D33,DTR COMMON /GAM/ C11,C12,C13,C22,C23,C33 C P1=Y(4) P2=Y(5) P3=Y(6) P2P3=P2*P3 P1P2=P1*P2 P1P3=P1*P3 P1P1=P1*P1 P2P2=P2*P2 P3P3=P3*P3 C11=P1P1*A11+P2P2*A66+P3P3*A55 1+2.*(P2P3*A56+P1P3*A15+P1P2*A16) C22=P1P1*A66+P2P2*A22+P3P3*A44 1+2.*(P2P3*A24+P1P3*A46+P1P2*A26) C33=P1P1*A55+P2P2*A44+P3P3*A33 1+2.*(P2P3*A34+P1P3*A35+P1P2*A45) C23=P1P1*A56+P2P2*A24+P3P3*A34 1 +P2P3*A2344+P1P3*A3645+P1P2*A2546 C13=P1P1*A15+P2P2*A46+P3P3*A35 1 +P2P3*A3645+P1P3*A1355+P1P2*A1456 C12=P1P1*A16+P2P2*A26+P3P3*A45 1 +P2P3*A2546+P1P3*A1456+P1P2*A1266 C11N=C11-1. C22N=C22-1. C33N=C33-1. C23SQ=C23*C23 C13SQ=C13*C13 C12SQ=C12*C12 D11=C22N*C33N-C23SQ D22=C11N*C33N-C13SQ D33=C11N*C22N-C12SQ D12=C13*C23-C12*C33N D13=C12*C23-C13*C22N D23=C12*C13-C23*C11N DTR=D11+D22+D33 IF(ABS(DTR).LT.0.0000001)THEN WRITE(LOUT,'(A)')'CHRM: SHEAR WAVE SINGULARITY' IND=10 END IF RETURN END C C ********************************************************* C SUBROUTINE CHRM1(C,PN,UN) C C ROUTINE FOR THE COMPUTATION OF THE ELEMENTS OF THE CHRISTOFFEL C MATRIX FOR AN ARBITRARY ANISOTROPIC MEDIUM C DIMENSION C(3,3),PN(3),UN(3) COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33, 1 A34,A35,A36,A44,A45,A46,A55,A56,A66, 1 DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23, 1 DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44, 1 DXA45,DXA46,DXA55,DXA56,DXA66, 1 DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23, 1 DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44, 1 DYA45,DYA46,DYA55,DYA56,DYA66, 1 DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23, 1 DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44, 1 DZA45,DZA46,DZA55,DZA56,DZA66, 1 A2546,A1266,A1355,A1456,A3645,A2344 C P1=PN(1) P2=PN(2) P3=PN(3) U1=UN(1) U2=UN(2) U3=UN(3) P2U3=P2*U3 P3U2=P3*U2 P1U2=P1*U2 P2U1=P2*U1 P1U3=P1*U3 P3U1=P3*U1 P1U1=P1*U1 P2U2=P2*U2 P3U3=P3*U3 C(1,1)=P1U1*A11+P2U2*A66+P3U3*A55 1+(P2U3+P3U2)*A56+(P1U3+P3U1)*A15+(P1U2+P2U1)*A16 C(2,2)=P1U1*A66+P2U2*A22+P3U3*A44 1+(P2U3+P3U2)*A24+(P1U3+P3U1)*A46+(P1U2+P2U1)*A26 C(3,3)=P1U1*A55+P2U2*A44+P3U3*A33 1+(P2U3+P3U2)*A34+(P1U3+P3U1)*A35+(P1U2+P2U1)*A45 C(2,3)=P1U1*A56+P2U2*A24+P3U3*A34 1+0.5*((P2U3+P3U2)*A2344+(P1U3+P3U1)*A3645+(P1U2+P2U1)*A2546) C(1,3)=P1U1*A15+P2U2*A46+P3U3*A35 1+0.5*((P2U3+P3U2)*A3645+(P1U3+P3U1)*A1355+(P1U2+P2U1)*A1456) C(1,2)=P1U1*A16+P2U2*A26+P3U3*A45 1+0.5*((P2U3+P3U2)*A2546+(P1U3+P3U1)*A1456+(P1U2+P2U1)*A1266) C(2,1)=C(1,2) C(3,2)=C(2,3) C(3,1)=C(1,3) RETURN END C C ********************************************************* C SUBROUTINE CHRM2(Y,G,i) C C EVALUATES ELEMENTS OF THE CHRISTOFFEL MATRIX C DIMENSION a(21),Y(18),G(3,3) COMMON/GAM/G11,G12,G13,G22,G23,G33 COMMON /APROX1/ e(21,10) 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 C DO 1 J=1,21 A(J)=E(J,I) 1 CONTINUE P1=Y(4) P2=Y(5) P3=Y(6) P11=P1*P1 P12=P1*P2 P13=P1*P3 P22=P2*P2 P23=P2*P3 P33=P3*P3 G11=A(1)*P11+A(21)*P22+A(19)*P33+ 1 2.*(A(6)*P12+A(5)*P13+A(20)*P23) G22=A(21)*P11+A(7)*P22+A(16)*P33+ 1 2.*(A(11)*P12+A(18)*P13+A(9)*P23) G33=A(19)*P11+A(16)*P22+A(12)*P33+ 1 2.*(A(17)*P12+A(14)*P13+A(13)*P23) G12=A(6)*P11+A(11)*P22+A(17)*P33+ 1 (A(21)+A(2))*P12+(A(20)+A(4))*P13+(A(10)+A(18))*P23 G13=A(5)*P11+A(18)*P22+A(14)*P33+ 1 (A(20)+A(4))*P12+(A(19)+A(3))*P13+(A(17)+A(15))*P23 G23=A(20)*P11+A(9)*P22+A(13)*P33+ 1 (A(10)+A(18))*P12+(A(17)+A(15))*P13+(A(16)+A(8))*P23 G(1,1)=G11 G(1,2)=G12 G(1,3)=G13 G(2,1)=G12 G(2,2)=G22 G(2,3)=G23 G(3,1)=G13 G(3,2)=G23 G(3,3)=G33 RETURN END C C ********************************************************* C SUBROUTINE PCHRM(Y,G,L,I) C C EVALUATES FIRST DERIVATIVES OF ELEMENTS OF CHRISTOFFEL MATRIX C WITH RESPECT TO THE L-TH COMPONENT OF THE SLOWNESS VECTOR C DIMENSION A(21),Y(18),G(3,3) COMMON /APROX1/ E(21,10) C DO 1 J=1,21 A(J)=E(J,I) 1 CONTINUE P1=Y(4) P2=Y(5) P3=Y(6) IF(L.EQ.1)THEN G(1,1)=2.*(A(1)*P1+A(6)*P2+A(5)*P3) G(2,2)=2.*(A(21)*P1+A(11)*P2+A(18)*P3) G(3,3)=2.*(A(19)*P1+A(17)*P2+A(14)*P3) AUX=2.*A(6)*P1+(A(21)+A(2))*P2+(A(20)+A(4))*P3 G(1,2)=AUX G(2,1)=AUX AUX=2.*A(5)*P1+(A(20)+A(4))*P2+(A(19)+A(3))*P3 G(1,3)=AUX G(3,1)=AUX AUX=2.*A(20)*P1+(A(10)+A(18))*P2+(A(17)+A(15))*P3 G(2,3)=AUX G(3,2)=AUX END IF IF(L.EQ.2)THEN G(1,1)=2.*(A(6)*P1+A(21)*P2+A(20)*P3) G(2,2)=2.*(A(11)*P1+A(7)*P2+A(9)*P3) G(3,3)=2.*(A(17)*P1+A(16)*P2+A(13)*P3) AUX=2.*A(11)*P2+(A(21)+A(2))*P1+(A(10)+A(18))*P3 G(1,2)=AUX G(2,1)=AUX AUX=2.*A(18)*P2+(A(20)+A(4))*P1+(A(17)+A(15))*P3 G(1,3)=AUX G(3,1)=AUX AUX=2.*A(9)*P2+(A(10)+A(18))*P1+(A(16)+A(8))*P3 G(2,3)=AUX G(3,2)=AUX END IF IF(L.EQ.3)THEN G(1,1)=2.*(A(5)*P1+A(20)*P2+A(19)*P3) G(2,2)=2.*(A(18)*P1+A(9)*P2+A(16)*P3) G(3,3)=2.*(A(14)*P1+A(13)*P2+A(12)*P3) AUX=2.*A(17)*P3+(A(20)+A(4))*P1+(A(10)+A(18))*P2 G(1,2)=AUX G(2,1)=AUX AUX=2.*A(14)*P3+(A(19)+A(3))*P1+(A(17)+A(15))*P2 G(1,3)=AUX G(3,1)=AUX AUX=2.*A(13)*P3+(A(17)+A(15))*P1+(A(16)+A(8))*P2 G(2,3)=AUX G(3,2)=AUX END IF RETURN END C C ********************************************************* C SUBROUTINE PPCHRM(G,L,M,i) C C EVALUATES SECOND DERIVATIVES OF ELEMENTS OF CHRISTOFFEL MATRIX C WITH RESPECT TO THE L-TH AND M-TH COMPONENTS OF THE SLOWNESS C VECTOR C DIMENSION a(21),G(3,3) COMMON /APROX1/ e(21,10) C do 1 j=1,21 a(j)=e(j,i) 1 continue IF(L.EQ.1.AND.M.EQ.1)THEN G(1,1)=2.*A(1) G(2,2)=2.*A(21) G(3,3)=2.*A(19) AUX=2.*A(6) G(1,2)=AUX G(2,1)=AUX AUX=2.*A(5) G(1,3)=AUX G(3,1)=AUX AUX=2.*A(20) G(2,3)=AUX G(3,2)=AUX END IF IF(L.EQ.2.AND.M.EQ.2)THEN G(1,1)=2.*A(21) G(2,2)=2.*A(7) G(3,3)=2.*A(16) AUX=2.*A(11) G(1,2)=AUX G(2,1)=AUX AUX=2.*A(18) G(1,3)=AUX G(3,1)=AUX AUX=2.*A(9) G(2,3)=AUX G(3,2)=AUX END IF IF(L.EQ.3.AND.M.EQ.3)THEN G(1,1)=2.*A(19) G(2,2)=2.*A(16) G(3,3)=2.*A(12) AUX=2.*A(17) G(1,2)=AUX G(2,1)=AUX AUX=2.*A(14) G(1,3)=AUX G(3,1)=AUX AUX=2.*A(13) G(2,3)=AUX G(3,2)=AUX END IF IF((L.EQ.1.AND.M.EQ.2).OR.(L.EQ.2.AND.M.EQ.1))THEN G(1,1)=2.*A(6) G(2,2)=2.*A(11) G(3,3)=2.*A(17) AUX=A(21)+A(2) G(1,2)=AUX G(2,1)=AUX AUX=A(20)+A(4) G(1,3)=AUX G(3,1)=AUX AUX=A(10)+A(18) G(2,3)=AUX G(3,2)=AUX END IF IF((L.EQ.1.AND.M.EQ.3).OR.(L.EQ.3.AND.M.EQ.1))THEN G(1,1)=2.*A(5) G(2,2)=2.*A(18) G(3,3)=2.*A(14) AUX=A(20)+A(4) G(1,2)=AUX G(2,1)=AUX AUX=A(19)+A(3) G(1,3)=AUX G(3,1)=AUX AUX=A(17)+A(15) G(2,3)=AUX G(3,2)=AUX END IF IF((L.EQ.2.AND.M.EQ.3).OR.(L.EQ.3.AND.M.EQ.2))THEN G(1,1)=2.*A(20) G(2,2)=2.*A(9) G(3,3)=2.*A(13) AUX=A(10)+A(18) G(1,2)=AUX G(2,1)=AUX AUX=A(17)+A(15) G(1,3)=AUX G(3,1)=AUX AUX=A(16)+A(8) G(2,3)=AUX G(3,2)=AUX END IF 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 BOX(BRD) C DIMENSION BRD(6) COMMON/VRML/LUBRD,LUGRD,LUIND,LURAY C WRITE(LUBRD,109) WRITE(LUBRD,105) I=1 WRITE(LUBRD,112)I WRITE(LUBRD,110)BRD(1),BRD(3),BRD(5) WRITE(LUBRD,110)BRD(1),BRD(4),BRD(5) WRITE(LUBRD,110)BRD(1),BRD(4),BRD(6) WRITE(LUBRD,110)BRD(1),BRD(3),BRD(6) WRITE(LUBRD,110)BRD(1),BRD(3),BRD(5) WRITE(LUBRD,105) I=2 WRITE(LUBRD,112)I WRITE(LUBRD,110)BRD(1),BRD(3),BRD(5) WRITE(LUBRD,110)BRD(1),BRD(3),BRD(6) WRITE(LUBRD,110)BRD(2),BRD(3),BRD(6) WRITE(LUBRD,110)BRD(2),BRD(3),BRD(5) WRITE(LUBRD,110)BRD(1),BRD(3),BRD(5) WRITE(LUBRD,105) I=3 WRITE(LUBRD,112)I WRITE(LUBRD,110)BRD(2),BRD(3),BRD(5) WRITE(LUBRD,110)BRD(2),BRD(4),BRD(5) WRITE(LUBRD,110)BRD(2),BRD(4),BRD(6) WRITE(LUBRD,110)BRD(2),BRD(3),BRD(6) WRITE(LUBRD,110)BRD(2),BRD(3),BRD(5) WRITE(LUBRD,105) I=4 WRITE(LUBRD,112)I WRITE(LUBRD,110)BRD(1),BRD(4),BRD(5) WRITE(LUBRD,110)BRD(1),BRD(4),BRD(6) WRITE(LUBRD,110)BRD(2),BRD(4),BRD(6) WRITE(LUBRD,110)BRD(2),BRD(4),BRD(5) WRITE(LUBRD,110)BRD(1),BRD(4),BRD(5) WRITE(LUBRD,105) I=1 WRITE(LUBRD,112)I WRITE(LUBRD,110)BRD(1),BRD(3),BRD(5) WRITE(LUBRD,110)BRD(1),BRD(4),BRD(5) WRITE(LUBRD,110)BRD(2),BRD(4),BRD(5) WRITE(LUBRD,110)BRD(2),BRD(3),BRD(5) WRITE(LUBRD,110)BRD(1),BRD(3),BRD(5) WRITE(LUBRD,105) I=1 WRITE(LUBRD,112)I WRITE(LUBRD,110)BRD(1),BRD(3),BRD(6) WRITE(LUBRD,110)BRD(1),BRD(4),BRD(6) WRITE(LUBRD,110)BRD(2),BRD(4),BRD(6) WRITE(LUBRD,110)BRD(2),BRD(3),BRD(6) WRITE(LUBRD,110)BRD(1),BRD(3),BRD(6) WRITE(LUBRD,105) WRITE(LUBRD,105) C 105 FORMAT('/') 109 FORMAT(25H'BOUNDARIES OF THE MODEL') 110 FORMAT(3(F10.5,1X),'/') 112 FORMAT(6H'BOUND,I1,1H',1X,'/') C RETURN END C C======================================================================= C INCLUDE 'a2.for' C a2.for INCLUDE 'a3.for' C a3.for INCLUDE 'a4.for' C a4.for INCLUDE 'a5.for' C a5.for C C Interpolation method: C Include just one of the following files 'mod*.for': C (a) Isosurface interpolation: INCLUDE 'modis.for' C modis.for C (b) (Bi-)(tri-)cubic B-spline interpolation: * INCLUDE 'modbs.for' C modbs.for C C======================================================================= Canraygse.for 0100666 0000765 0000765 00000012327 12046713172 012762 0 ustar bulant bulant C
C Program ANRAYGSE to read the synthetic seismograms written C in the form of file LU8 of package ANRAY and to write them C in the GSE format. C C Refer to file 'anraygse.htm' for the description C of the input parameters. C C Version: 4.73 C Date: 2012, November 8 C C Coded by: Petr Bulant C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C http://sw3d.cz/staff/bulant.htm C C----------------------------------------------------------------------- C C Subroutines and external functions required: EXTERNAL WGSE1,WGSE2,WGSE3,RSEP1,RSEP3T,ERROR,LENGTH C WGSE1,WGSE2,WGSE3... File 'gse.for'. C RSEP1,RSEP3T... File 'sep.for'. C ERROR... File 'error.for'. C LENGTH... File 'length.for'. C C----------------------------------------------------------------------- C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc INTEGER MT PARAMETER (MT=MRAM/2) INTEGER IS(MT) REAL SEIS(MT) EQUIVALENCE (IS,RAM(1)) EQUIVALENCE (SEIS,RAM(MT+1)) C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER LUSEP,LU8,LUGSE PARAMETER (LUSEP=1,LU8=2,LUGSE=3) CHARACTER*80 FILSEP,FILLU8,FILGSE CHARACTER*6 RECNAM INTEGER I,NT,MCOMP,NDST,ILOC,IREC LOGICAL LREC REAL XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,DT,DF,DST,TO,AREDUC,X1,X2,X3, * TSHIFT CHARACTER*80 MPRINT,IPRINT,STEXT C C....................................................................... C C Reading name of SEP file with input data: WRITE(*,'(A)') '+ANRAYGSE: Enter input filename: ' FILSEP=' ' READ(*,*) FILSEP C C Reading all data from the SEP file into the memory: IF (FILSEP.NE.' ') THEN CALL RSEP1(LUSEP,FILSEP) ELSE C ANRAYGSE-01 CALL ERROR('ANRAYGSE-01: SEP file not given') C Input file in the form of the SEP (Stanford Exploration Project) C parameter or history file must be specified. C There is no default filename. ENDIF C WRITE(*,'(A)') '+ANRAYGSE: Working... ' C C Reading file LU8: CALL RSEP3T('LU8',FILLU8 ,'lu8.out') OPEN(LU8,FILE=FILLU8,STATUS='OLD') C Opening output file GSE: CALL RSEP3T('SS',FILGSE,'ss.gse') OPEN(LUGSE,FILE=FILGSE) C Reading optional time shift: CALL RSEP3R('TSHIFT',TSHIFT,0.) C Reading receiver name generation switch: CALL RSEP3I('IRECNAM',IREC,0) LREC=.FALSE. RECNAM=' ' IF (IREC.EQ.1) THEN LREC=.TRUE. RECNAM='rec ' ENDIF C C Reading and writing the headers of the files: READ(LU8,'(A)') MPRINT READ(LU8,'(A)') IPRINT READ(LU8,'(A)') STEXT READ(LU8,'(5F10.5,2E15.7)') XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,DT,DF READ(LU8,'(16I5)') NDST,NT,MCOMP,ILOC CALL WGSE1(LUGSE,' ') C IF (MCOMP.EQ.0) MCOMP=3 X1=-999. X2=-999. X3=-999. C C Reading and writing seismograms: IREC=0 10 CONTINUE READ(LU8,'(2F10.3,1E12.5,I5)',END=90) DST,TO,AREDUC,NT IREC=IREC+1 IF (NT.GT.MT) THEN C ANRAYGSE-02 CALL ERROR('ANRAYGSE-02: Small arrays IS and SEIS') C The dimension MT of arrays IT and SEIS should be enlarged. ENDIF READ(LU8,'(20I4)') (IS(I),I=1,NT) IF (ILOC.EQ.0) THEN X1=DST X2=DST ELSEIF (ILOC.EQ.1) THEN X3=DST ENDIF TO=TO+TSHIFT IF (LREC) THEN IF (IREC.GE.1000) THEN C ANRAYGSE-03 CALL ERROR('ANRAYGSE-03: Too many receiver names') C This version of the code enables up to 999 receiver names. ENDIF RECNAM(6:6)=CHAR(ICHAR('0')+MOD(IREC,10)) RECNAM(5:5)=CHAR(ICHAR('0')+IREC/10) RECNAM(4:4)=CHAR(ICHAR('0')+IREC/100) ENDIF DO 20, I=1,NT SEIS(I)=(FLOAT(IS(I))/999.1)*AREDUC 20 CONTINUE CALL WGSE2(LUGSE,RECNAM,' ',MCOMP,X1,X2,X3,TO,DT,NT,SEIS) GOTO 10 90 CONTINUE CLOSE(LU8) CALL WGSE3(LUGSE) WRITE(*,'(A)') '+ANRAYGSE: Done. ' STOP END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'sep.for' C sep.for INCLUDE 'gse.for' C gse.for INCLUDE 'length.for' C length.for C C======================================================================= Canraypl.for 0100666 0000765 0000765 00000035343 11404416620 012615 0 ustar bulant bulant C
C PROGRAM A N R A Y P L C ********************* C C THE PROGRAM ANRAYPL IS DESIGNED FOR PLOTTING OF RAY DIAGRAMS, C TRAVEL TIMES AND AMPLITUDES OF SEISMIC BODY WAVES FROM THE C FILE LU1 GENERATED BY PROGRAM ANRAY C C **************************************************************** C CHARACTER*80 TEXT,PSTEXT,FILEIN,FILEOU,FILE1 COMPLEX AX,AY,AZ,AUX,CAUX,CSOUR,IMAG DIMENSION A(30),B(30),C(30),D(30),X1(30),NPNT(20),NR(40) DIMENSION X(500),T(500),XZ(500),AX(2,500),AY(2,500),AZ(2,500), 1Y(500),Z(500),INDI(500),DST(200),DD(12),SV(3),E1(3),E2(3), 2AUX(3),AM(3),PH(3),P(3,500),POL(3,500),POL1(3,500) COMMON/SOUR/ROS C C C************************************************** C LIN=5 LOU=6 LU1=1 FILEIN='anrpl.dat' FILEOU='anrpl.out' WRITE(*,'(2A)') ' (ANRAYPL) SPECIFY NAMES OF INPUT AND OUTPUT', 1' FILES LIN, LOU, LU1: ' READ(*,*) FILEIN,FILEOU,FILE1 IF(FILE1.EQ.' ') GO TO 99 OPEN(LIN,FILE=FILEIN,FORM='FORMATTED',STATUS='OLD') OPEN(LOU,FILE=FILEOU,FORM='FORMATTED') OPEN(LU1,FILE=FILE1,FORM='FORMATTED',STATUS='OLD') C C************************************************** C IMAG=(0.,1.) XB=0. IRUN=0 IPAGE=0 PSTEXT=' ' IPRINT=0 XSHIFT=5. YSHIFT=3. READ(LIN,*)IPRINT,XSHIFT,YSHIFT,PSTEXT WRITE(LOU,116)IPRINT,XSHIFT,YSHIFT,PSTEXT IF(IPRINT.LT.0)THEN CALL PLOTN(PSTEXT,0) KPRINT=-IPRINT END IF CALL PLOTS(LDUM1,LDUM2,7) C REWIND LU1 CALL PLOT(4.,5.,-3) C 200 READ(LU1,101)ICONT,NDST,ILOC WRITE(LOU,101)ICONT,NDST,ILOC IF(ICONT.EQ.0)WRITE(LOU,107)LU1 IF(ICONT.EQ.0)GO TO 99 READ(LU1,102)ROS WRITE(LOU,102)ROS IF(IPAGE.NE.0)THEN CALL PLOT(0.,0.,999) IF(IPRINT.LT.0)CALL PLOTN(PSTEXT,1) CALL PLOTS(LDUM1,LDUM2,7) CALL PLOT(4.,5.,-3) END IF C TEXT='ANRAYPL' READ(LIN,*)TEXT WRITE(LOU,100)TEXT NTICX=0 NTICY=0 NTICH=0 NTICT=0 NTICA=0 IPROF=1 NRAY=0 IBOUND=100 IRED=0 IRS=0 NDX=0 NDY=0 NDH=0 NDT=0 NDA=0 IPAGE=0 READ(LIN,*)NTICX,NTICY,NTICH,NTICT,NTICA,IPROF,NRAY,IBOUND, 1IRED,IRS,NDX,NDY,NDH,NDT,NDA,IPAGE WRITE(LOU,106)NTICX,NTICY,NTICH,NTICT,NTICA,IPROF,NRAY,IBOUND, 1IRED,IRS,NDX,NDY,NDH,NDT,NDA,IPAGE C IF(NTICX.EQ.0)GO TO 99 XB = 0. SC=1. AROT=0. READ(LIN,*)XMIN,XMAX,XLEN,DTICX,SC,AROT WRITE(LOU,102)XMIN,XMAX,XLEN,DTICX,SC,AROT CSRT=COS(AROT) SNRT=SIN(AROT) IF(ABS(XMAX-XMIN).LT..00001)GO TO 32 XMER = XLEN/(XMAX-XMIN) GLH = 1.5 READ(LIN,*)YMIN,YMAX,YLEN,DTICY,HMIN,HMAX,HLEN,DTICH WRITE(LOU,102)YMIN,YMAX,YLEN,DTICY,HMIN,HMAX,HLEN,DTICH IF(NTICY.EQ.0)GO TO 40 YMER = YLEN/(YMAX-YMIN) C C PLOTTING OF BORDER FOR VERTICAL RAY DIAGRAM C IF(IRUN.EQ.1)CALL PLOT(XSHIFT/3.,0.,-3) IRUN=1 CALL BORDER(XLEN,DTICX,YLEN,DTICY,SC,TEXT,1,XMIN,XMAX,YMIN,YMAX, 1NTICX,NTICY,NDX,NDY) TX=.5*(XLEN-6.3*SC) CALL SYMBOL(TX,-1.6*SC,.45*SC,'DISTANCE IN KM',0.,14) TX=.5*(YLEN-4.95*SC) U=-(1.6+.4*NDX)*SC CALL SYMBOL(U,TX,.45*SC,'DEPTH IN KM',90.,11) C 40 IF(NTICH.EQ.0)GO TO 32 HMER = HLEN/(HMAX-HMIN) C C PLOTTING OF BORDER FOR HORIZONTAL RAY DIAGRAM C IF(IRUN.EQ.1.AND.NTICY.NE.0)CALL PLOT(0.,YLEN+YSHIFT,-3) IF(IRUN.EQ.1.AND.NTICY.EQ.0)CALL PLOT(XLEN+XSHIFT,0.,-3) IRUN=1 CALL BORDER(XLEN,DTICX,HLEN,DTICH,SC,TEXT,1,XMIN,XMAX,HMIN,HMAX, 1NTICX,NTICH,NDX,NDH) CALL PLOT(0.,.5*HLEN,3) CALL PLOT(XLEN,.5*HLEN,2) CALL PLOT(0.,0.,3) TX=.5*(XLEN-6.3*SC) CALL SYMBOL(TX,-1.6*SC,.45*SC,'DISTANCE IN KM',0.,14) TX=.5*(HLEN-8.1*SC) U=-(1.6+.4*NDX)*SC CALL SYMBOL(U,TX,.45*SC,'TRANSV.DIST. IN KM',90.,18) IF(NTICY.NE.0)CALL PLOT(0.,-YLEN-YSHIFT,-3) C C PLOTTING OF INTERFACES C 32 READ(LU1,101)NINT,(NPNT(I),I=1,NINT) IF(KPRINT.EQ.3)WRITE(LOU,106)NINT,(NPNT(I),I = 1,NINT) IB=IBOUND IF(IBOUND.LT.0)IB=-IBOUND RD = (XMAX-XMIN)/FLOAT(IB-1) YM = YMAX DO 5 I = 1,NINT N = NPNT(I)-1 READ(LU1,113)(A(J),B(J),C(J),D(J),X1(J),J = 1,N) IF(KPRINT.EQ.3)WRITE(LOU,105)(A(J),B(J),C(J),D(J),X1(J), 1J=1,N) IF(NTICY.EQ.0)GOTO 5 MMM=1 IF(IBOUND.LT.0)YM=YM-.06/YMER IF(IBOUND.LT.0)MMM=3 DO 4 M=1,MMM RB=XMIN IF(IBOUND.LT.0)YM=YM+.03/YMER IPL=0 DO 4 J = 1,IB IF(RB.GT.(XMAX+0.01))GOTO 4 DO 1 K = 1,N IF(K.EQ.N) K1 = K IF(RB.GE.X1(K))GOTO 1 C K1 = K-1 GOTO 2 1 CONTINUE 2 X2 = RB-X1(K1) X33=X1(K1+1) IF(K1.EQ.N)X33=XMAX IF(XMAX.LT.X33)X33=XMAX X3=X33-RB IF(X3.LT.RD)RB=X33 IF(X3.LT.RD)X2=X33-X1(K1) HXX=YM-(((D(K1)*X2+C(K1))*X2+B(K1))*X2+A(K1)-YMIN) IF(HXX.LT.YMIN.OR.HXX.GT.YMAX)IPL=0 IF(HXX.LT.YMIN.OR.HXX.GT.YMAX)GO TO 4 IF(IPL.NE.0)GO TO 3 IPL=1 XSV=(RB-XMIN)*XMER YSV = HXX*YMER CALL PLOT(XSV,YSV,3) GOTO 4 3 XSV = (RB-XMIN)*XMER YSV = HXX*YMER CALL PLOT(XSV,YSV,2) 4 RB = RB+RD IF(IBOUND.LT.0)YM=YM-0.03/YMER 5 CONTINUE 20 READ(LU1,102)X0,Y0,Z0,FI READ(LU1,102)(DST(I),I=1,NDST) IF(KPRINT.EQ.3)WRITE(LOU,102)X0,Y0,Z0,FI IF(KPRINT.EQ.3)WRITE(LOU,102)(DST(I),I=1,NDST) IF(NTICY.EQ.0)GOTO 21 CSF=COS(FI) SNF=SIN(FI) C C PLOTTING OF RAYS C 26 IF(NRAY.NE.0) READ(LIN,*)(NR(I),I = 1,NRAY) IF(NRAY.NE.0) WRITE(LOU,104)(NR(I),I = 1,NRAY) I = 1 K = 1 21 READ(LU1,106)N,IND IPL = 0 IF(KPRINT.GE.2)WRITE(LOU,106)N,IND IF(N.EQ.0)GOTO 24 READ(LU1,103)(T(J),X(J),Y(J),Z(J),(DD(L),L=1,12),J = 1,N) IF(KPRINT.GE.2)WRITE(LOU,102)(T(J),X(J),Y(J),Z(J),J = 1,N) IF(NTICY.EQ.0.AND.NTICH.EQ.0)GO TO 21 IF(NTICY.EQ.0) GO TO 41 IF(NRAY.EQ.0)GOTO 25 IF(NR(K).EQ.I)GOTO 22 25 DO 10 J = 1,N XX=(X(J)-X0)*CSF+(Y(J)-Y0)*SNF YY=Z(J) IF(XX.LT.XMIN.OR.XX.GT.XMAX)GOTO 8 IF(YY.LT.YMIN.OR.YY.GT.YMAX)GOTO 9 IF(IPL.EQ.1)GOTO 6 IF(J.LE.IRS)GO TO 10 IPL = 1 INDEX = 3 GO TO 7 6 INDEX = 2 7 XNEW = (XX-XMIN)*XMER YNEW = (YMAX-YY)*YMER CALL PLOT(XNEW,YNEW,INDEX) GOTO 10 8 IF (IPL.EQ.0) GOTO 10 IF(XX.LT.XMIN) XXX = XMIN IF(XX.GT.XMAX) XXX = XMAX XOL=(X(J-1)-X0)*CSF+(Y(J-1)-Y0)*SNF TG=(YY-Z(J-1))/(XX-XOL) YY = TG*(XXX-XX)+YY XX = XXX IPL=0 GOTO 6 9 IF (IPL.EQ.0) GOTO 10 IF(YY.LT.YMIN) YYY=YMIN IF(YY.GT.YMAX) YYY=YMAX XOL=(X(J-1)-X0)*CSF+(Y(J-1)-Y0)*SNF TG=(XX-XOL)/(YY-Z(J-1)) XX = TG*(YYY-YY)+XX YY = YYY IPL=0 GOTO 6 10 CONTINUE c 41 IF(NTICH.EQ.0)GO TO 23 IF(NTICY.EQ.0)SHF=0. IF(NTICY.NE.0)SHF=YLEN+YSHIFT-HMIN*HMER IPL=0 DO 42 J = 1,N XX=(X(J)-X0)*CSF+(Y(J)-Y0)*SNF HH=-(X(J)-X0)*SNF+(Y(J)-Y0)*CSF IF(XX.LT.XMIN.OR.XX.GT.XMAX)GOTO 45 IF(HH.LT.HMIN.OR.HH.GT.HMAX)GOTO 46 IF(IPL.EQ.1)GOTO 43 IF(J.LE.IRS)GO TO 42 IPL = 1 INDEX = 3 GO TO 44 43 INDEX = 2 44 XNEW = (XX-XMIN)*XMER HNEW = (HMAX-(HH-HMIN))*HMER CALL PLOT(XNEW,HNEW+SHF,INDEX) GOTO 42 45 IF (IPL.EQ.0) GOTO 42 IF(XX.LT.XMIN) XXX = XMIN IF(XX.GT.XMAX) XXX = XMAX XOL=(X(J-1)-X0)*CSF+(Y(J-1)-Y0)*SNF HOL=-(X(J-1)-X0)*SNF+(Y(J-1)-Y0)*CSF TG=(HH-HOL)/(XX-XOL) HH = TG*(XXX-XX)+HH XX = XXX IPL=0 GOTO 43 46 IF (IPL.EQ.0) GOTO 42 IF(HH.LT.HMIN)HHH=HMIN IF(HH.GT.HMAX)HHH=HMAX XOL=(X(J-1)-X0)*CSF+(Y(J-1)-Y0)*SNF HOL=-(X(J-1)-X0)*SNF+(Y(J-1)-Y0)*CSF TG=(XX-XOL)/(HH-HOL) XX = TG*(HHH-HH)+XX HH = HHH IPL=0 GOTO 43 42 CONTINUE GOTO 23 22 IF(K.GE.NRAY)GOTO 23 K = K+1 23 I = I+1 GOTO 21 24 CONTINUE C C PLOTTING OF TIME-DISTANCE CURVE C 11 READ(LU1,101)NS IF(KPRINT.GE.1)WRITE(LOU,101)NS NWTYP=NS IF(NS.LT.0)NS=-NS IF(NS.NE.0)THEN DO 300 I=1,NS READ(LU1,112)INDI(I),X(I),Z(I),T(I),AX(1,I),AY(1,I),AZ(1,I) IF(NWTYP.LT.0)READ(LU1,115)AX(2,I),AY(2,I),AZ(2,I) READ(LU1,111)(P(K,I),K=1,3) READ(LU1,111)(POL(K,I),K=1,3) IF(NWTYP.LT.0)READ(LU1,111)(POL1(K,I),K=1,3) 300 CONTINUE END IF IF(NS.NE.0.AND.KPRINT.GE.1)WRITE(LOU,108)(INDI(I),X(I),Z(I), 1T(I),(AX(J,I),AY(J,I),AZ(J,I),J=1,2),I=1,NS) NSS = NS 35 CONTINUE IF(NTICT.EQ.0)GOTO 15 VRED=6. READ(LIN,*)TMIN,TMAX,TLEN,DTICT,VRED WRITE(LOU,102)TMIN,TMAX,TLEN,DTICT,VRED SHF=0. IF(NTICY.NE.0)SHF=SHF+YLEN+YSHIFT IF(NTICH.NE.0)SHF=SHF+HLEN+YSHIFT IF(IRUN.EQ.1)CALL PLOT(0.,SHF,-3) IRUN=1 IF(NSS.EQ.0)GOTO 30 IEX=0 IF(ILOC.EQ.1)THEN XMER=YLEN/(YMAX-YMIN) DTICX=DTICY XMIN=YMIN XMAX=YMAX XLEN=YLEN NDX=NDY END IF YMER=TLEN/(TMAX-TMIN) NAUX=3 CALL BORDER(XLEN,DTICX,TLEN,DTICT,SC,TEXT,0,XMIN,XMAX,TMIN, 1TMAX,NTICX,NTICT,NDX,NDT) TX=.5*(XLEN-6.3*SC) IF(ILOC.NE.1) 1CALL SYMBOL(TX,-1.6*SC,.45*SC,'DISTANCE IN KM',0.,14) IF(ILOC.EQ.1) 1CALL SYMBOL(TX,-1.6*SC,.45*SC,'DEPTH IN KM',0.,11) TX=.5*(TLEN-8.1*SC) U=-(1.6+.4*NDX)*SC IF(IRED.EQ.0) 1CALL SYMBOL(U,TX,.45*SC,'TRAVEL TIME IN SEC',90.,18) IF(IRED.EQ.0)GO TO 27 CALL SYMBOL(U,TX,.45*SC,'T-D/ ',90.,5) TX=TX+1.8*SC CALL NUMBER(U,TX,.45*SC,VRED,90.,2) TX=TX+2.7*SC CALL SYMBOL(U,TX,.45*SC,'(IN SEC)',90.,8) 27 INDEX=3 12 DO 13 I = 1,NS XNEW=X(I) IF(IEX.EQ.1)XNEW=XZ(I) IF(XNEW.LT.XMIN.OR.XNEW.GT.XMAX)GOTO 13 YNEW=T(I) IF(IRED.NE.0)YNEW=T(I)-ABS(XNEW)/VRED XNEW=(XNEW-XMIN)*XMER IF(YNEW.LT.TMIN.OR.YNEW.GT.TMAX)GOTO 13 YNEW=(YNEW-TMIN)*YMER CALL SYMBOL(XNEW,YNEW,.15,CHAR(INDEX),0.,-1) 13 CONTINUE IF(IEX.EQ.0)GOTO 30 NS=NSS GOTO 14 30 NEXP=0 READ(LIN,*)NEXP WRITE(LOU,104)NEXP IF(NEXP.EQ.0)GOTO 14 NS=NEXP READ(LIN,*)(XZ(I),T(I),I=1,NS) WRITE(LOU,102)(XZ(I),T(I),I=1,NS) IF(NSS.EQ.0)GO TO 15 IEX=1 INDEX=4 GOTO 12 14 CONTINUE CALL PLOT(0.,-SHF,-3) C C PLOTTING OF AMPLITUDE-DISTANCE CURVE C 15 IF(NTICA.EQ.0)GO TO 200 IRUN1=0 ALBACK=0. 19 ICOMP=0 MSOUR=0 READ(LIN,*)AMIN,AMAX,ALEN,DTICA,ICOMP,MSOUR WRITE(LOU,109)AMIN,AMAX,ALEN,DTICA,ICOMP,MSOUR IF(MSOUR.NE.0)CALL SOURCE(LIN,LOU,0,0,MSOUR,SV,E1,AMSOUR,PHSOUR) IF(ALEN.LT..00001)CALL PLOT(XSHIFT+XLEN,-ALBACK,-3) IF(ALEN.LT..00001)GO TO 200 38 IF(NSS.EQ.0)GO TO 36 IF(IRUN.EQ.1.AND.IRUN1.EQ.0)CALL PLOT(XLEN+XSHIFT,0.,-3) IF(IRUN1.EQ.1)CALL PLOT(0.,ALEN+YSHIFT,-3) IF(IRUN1.EQ.1)ALBACK=ALBACK+ALEN+YSHIFT IRUN1=1 IRUN=1 IEX=0 YMER=ALEN/(AMAX-AMIN) NAUX=2 CALL BORDER(XLEN,DTICX,ALEN,DTICA,SC,TEXT,0,XMIN,XMAX,AMIN, 1AMAX,NTICX,NTICA,NDX,NDA) TX=.5*(XLEN-6.3*SC) IF(ILOC.NE.1) 1CALL SYMBOL(TX,-1.6*SC,.45*SC,'DISTANCE IN KM',0.,14) IF(ILOC.EQ.1) 1CALL SYMBOL(TX,-1.6*SC,.45*SC,'DEPTH IN KM',0.,11) TX=.5*(ALEN-7.65*SC) U=-(1.6+.4*NDX)*SC CALL SYMBOL(U,TX,.45*SC,'A M P L I T U D E',90.,17) IF(ICOMP.EQ.0) 1CALL SYMBOL(.45*SC,ALEN+SC,.45*SC,'VERTICAL',0.,8) IF(ICOMP.EQ.1) 1CALL SYMBOL(.45*SC,ALEN+SC,.45*SC,'X-COMPONENT',0.,11) IF(ICOMP.EQ.2) 1CALL SYMBOL(.45*SC,ALEN+SC,.45*SC,'Y-COMPONENT',0.,11) 28 INDEX=3 16 DO 17 I=1,NS IF(MSOUR.NE.0)THEN DO 33 K=1,3 SV(K)=P(K,I) E1(K)=POL(K,I) E2(K)=POL1(K,I) 33 CONTINUE CALL SOURCE(LIN,LOU,1,3,MSOUR,SV,E1,AMSOUR,PHSOUR) CSOUR=AMSOUR*CEXP(IMAG*PHSOUR) END IF IF(MSOUR.EQ.0)CSOUR=(1.,0.) AUX(1)=AX(1,I)*CSOUR AUX(2)=AY(1,I)*CSOUR AUX(3)=AZ(1,I)*CSOUR IF(NWTYP.LT.0)THEN IF(MSOUR.NE.0)THEN CALL SOURCE(LIN,LOU,1,1,MSOUR,SV,E2,AMSOUR,PHSOUR) CSOUR=AMSOUR*CEXP(IMAG*PHSOUR) END IF IF(MSOUR.EQ.0)CSOUR=(1.,0.) AUX(1)=AUX(1)+AX(2,I)*CSOUR AUX(2)=AUX(2)+AY(2,I)*CSOUR AUX(3)=AUX(3)+AZ(2,I)*CSOUR END IF CAUX=AUX(1)*CSRT+AUX(2)*SNRT AUX(2)=-AUX(1)*SNRT+AUX(2)*CSRT AUX(1)=CAUX DO 37 K=1,3 RW=REAL(AUX(K)) CW=AIMAG(AUX(K)) AM(K)=SQRT(RW*RW+CW*CW) IF(ABS(CW).LE..000001.AND.ABS(RW).LE.000001)PH(K)=0. IF(ABS(CW).GT..000001.OR.ABS(RW).GT.000001)PH(K)=ATAN2(CW,RW) 37 CONTINUE IF(KPRINT.EQ.2)WRITE(LOU,114)(AM(K),PH(K),K=1,3) XNEW=X(I) IF(ICOMP.EQ.0)YNEW=AM(3) IF(ICOMP.EQ.1)YNEW=AM(1) IF(ICOMP.EQ.2)YNEW=AM(2) IF(IEX.EQ.0.OR.NEXP.EQ.0)GO TO 31 XNEW=XZ(I) YNEW=T(I) 31 IF(XNEW.LE.XMIN.OR.XNEW.GE.XMAX)GO TO 17 XNEW=(XNEW-XMIN)*XMER IF(ABS(YNEW).LT.1E-10)GO TO 17 YNEW = ALOG10(ABS(YNEW)) IF(YNEW.LT.AMIN.OR.YNEW.GT.AMAX)GOTO 17 YNEW = (YNEW-AMIN)*YMER CALL SYMBOL(XNEW,YNEW,.15,CHAR(INDEX),0.,-1) 17 CONTINUE 36 IF(IEX.EQ.1)GOTO 18 NEXP=0 READ(LIN,*)NEXP WRITE(LOU,104)NEXP IF(NEXP.EQ.0)GOTO 18 NS=NEXP READ(LIN,*)(XZ(I),T(I),I=1,NS) WRITE(LOU,102)(XZ(I),T(I),I=1,NS) IF(NSS.EQ.0)GO TO 18 IEX=1 INDEX=4 GOTO 16 18 IF(IEX.EQ.0)GO TO 29 NS=NSS 29 CONTINUE C C GO TO 19 C C 100 FORMAT(A) 101 FORMAT(26I3) 102 FORMAT(8F10.5) 103 FORMAT(16E15.5) 104 FORMAT(2X,26I3) 105 FORMAT(5E12.6,I5) 106 FORMAT(16I5) 107 FORMAT(2X,'END OF FILE',I3) 108 FORMAT(I5,3F10.5,6E15.9,25x,6e15.9) 109 FORMAT(4F10.5,4I5) 111 FORMAT(3E15.5) 112 FORMAT(I5,9E15.5) 113 FORMAT(5E15.5) 114 FORMAT(2X,3(E15.5,F10.5)) 115 FORMAT(6E15.5) 116 FORMAT(I5,2F10.5,1X,A) C 99 CONTINUE C CALL PLOT(0.,0.,999) STOP END C C======================================================================= C INCLUDE 'source.for' C source.for INCLUDE 'border.for' C border.for INCLUDE 'error.for' C error.for INCLUDE 'length.for' C length.for INCLUDE 'sep.for' C sep.for INCLUDE 'calcops.for' C calcops.for C C======================================================================= Cborder.for 0100666 0000765 0000765 00000004620 11023552062 012414 0 ustar bulant bulant C
SUBROUTINE BORDER(GLX,DX,GLY,DY,SC,TEXT,N,XMIN,XMAX,YMIN,YMAX, 1NX,NY,NDECX,NDECY) C character*80 TEXT C IXX=0 IYY=0 GGG=.15*SC GGH=.225*SC DCX=.2*NDECX DCY=.4*NDECY DNX=ABS(DX)/NX DNY=ABS(DY)/NY XMER=GLX/(XMAX-XMIN) YMER=GLY/(YMAX-YMIN) A=DNX*XMER B=DNY*YMER XST=0. YST=0. IF(DX.GT.0.)GO TO 1 DXX=-DX IX=XMIN/DXX XM=(IX+1)*DXX IXX=(XM-XMIN)/DNX XST=(XM-XMIN)-IXX*DNX XST=XST*XMER GO TO 2 1 XM=XMIN 2 AVST=(XMAX-XM)/DNX+.0001 IX=IXX+AVST+1. CALL PLOT(0.,0.,3) DO 100 I=1,IX I1=I-IXX-1 C=(I-1)*A+XST CALL PLOT(C,0.,2) CALL PLOT(C,GGG,2) IF(I1/NX*NX.NE.I1)GO TO 100 VAL=XM+I1*DNX R=-(.3+DCX)*SC IF(VAL.GE.10.)R=-(.5+DCX)*SC IF(VAL.GE.100.)R=-(.7+DCX)*SC IF(VAL.GE.1000.)R=-(.9+DCX)*SC CALL PLOT(C,GGH,2) CALL NUMBER(C+R,-.7*SC,.4*SC,VAL,0.,NDECX) 100 CALL PLOT(C,0.,3) CALL PLOT(GLX,0.,2) IF(DY.GT.0.)GO TO 3 DY=-DY IY=YMIN/DY YM=(IY+1)*DY IYY=(YM-YMIN)/DNY YST=(YM-YMIN)-IYY*DNY YST=YST*YMER GO TO 4 3 YM=YMIN 4 AVST=(YMAX-YM)/DNY+.0001 IY=IYY+AVST+1. DO 200 I=1,IY IF(N.EQ.1)I1=IY-I-IYY IF(N.NE.1)I1=I-IYY-1 D=(I-1)*B+YST CALL PLOT(GLX,D,2) CALL PLOT(GLX-GGG,D,2) IF(I1/NY*NY.NE.I1)GO TO 200 CALL PLOT(GLX-GGH,D,2) 200 CALL PLOT(GLX,D,3) CALL PLOT(GLX,GLY,2) DO 300 I=1,IX I1=IX-I-IXX E=(IX-I)*A+XST CALL PLOT(E,GLY,2) CALL PLOT(E,GLY-GGG,2) IF(I1/NX*NX.NE.I1)GO TO 300 CALL PLOT(E,GLY-GGH,2) 300 CALL PLOT(E,GLY,3) CALL PLOT(0.,GLY,2) DO 400 I=1,IY IF(N.EQ.1)I1=I-IYY-1 IF(N.NE.1)I1=IY-IYY-I F=(IY-I)*B+YST CALL PLOT(0.,F,2) CALL PLOT(GGG,F,2) IF(I1/NY*NY.NE.I1)GO TO 400 VAL=YM+I1*DNY IF(VAL.LT.10..AND.VAL.GE.0.)R=-(.8+DCY)*SC IF(VAL.LT.0..OR.VAL.GE.10.)R=-(1.2+DCY)*SC IF(VAL.LE.(-10.).OR.VAL.GE.100.)R=-(1.6+DCY)*SC CALL PLOT(GGH,F,2) CALL NUMBER(R,F-.2*SC,.4*SC,VAL,0.,NDECY) 400 CALL PLOT(0.,F,3) CALL PLOT(0.,0.,2) CALL SYMBOL(0.,-3.*SC,.45*SC,TEXT,0.,80) RETURN END Cbplot.for 0100666 0000765 0000765 00000022037 11404426060 012262 0 ustar bulant bulant C
C P R O G R A M B P L O T C C C PROGRAM BPLOT IS DESIGNED FOR THE PLOTTING OF SYNTHETIC C SEISMOGRAMS AND GENERATING A FILE FOR POLARPLOT C C ************************************************************ C CHARACTER*80 TEXT,PSTEXT,MPRINT,IPRINT,TITLE CHARACTER*80 FILEIN,FILEOU,FILE1,FILE2 DIMENSION SMAX(100),SEIS(2048) DIMENSION IS(2048),IEP(100) C C C************************************************** C LIN=5 LOU=6 LU8=1 LU4=2 FILEIN='blot.dat' FILEOU='bplot.out' FILE1='lu8.dat' FILE2='lu4.dat' WRITE(*,'(2A)') ' (BPLOT) SPECIFY NAMES OF INPUT AND OUTPUT', 1' FILES LIN, LOU, LU8, LU4: ' READ(*,*) FILEIN,FILEOU,FILE1,FILE2 IF(FILE1.EQ.' ') GO TO 99 IF(FILE2.EQ.' ') LU4=0 OPEN(LIN,FILE=FILEIN,FORM='FORMATTED',STATUS='OLD') OPEN(LOU,FILE=FILEOU,FORM='FORMATTED') OPEN(LU8,FILE=FILE1,FORM='FORMATTED',STATUS='OLD') IF(LU4.NE.0)OPEN(LU4,FILE=FILE2,FORM='FORMATTED') C C************************************************** C TEXT='BPLOT' PSTEXT=' ' IPR=0 READ(LIN,*)TEXT READ(LIN,*)IPR,PSTEXT WRITE(LOU,108)TEXT WRITE(LOU,104)IPR,PSTEXT IF(IPR.LT.0)THEN CALL PLOTN(PSTEXT,0) IPR=-IPR END IF CALL PLOTS(LDUM1,LDUM2,7) CALL PLOT(4.,6.,-3) C C READING FROM LU8 C IF(LU8.NE.0)REWIND LU8 IF(LU4.NE.0)REWIND LU4 READ(LU8,106)MPRINT READ(LU8,106)IPRINT READ(LU8,106)TITLE READ(LU8,152)XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,DT,DF READ(LU8,101)NDST,NT,MCOMP,ILOC WRITE(LOU,108)MPRINT WRITE(LOU,108)IPRINT WRITE(LOU,108)TITLE WRITE(LOU,152)XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,DT,DF WRITE(LOU,101)NDST,NT,MCOMP C C INPUT DATA C MRED=0 MEPIC=0 NTICX=1 NTICY=1 NDX=0 NDY=0 1 READ(LIN,*)MRED,MEPIC,NTICX,NTICY,NDX,NDY WRITE(LOU,101)MRED,MEPIC,NTICX,NTICY,NDX,NDY IF(MEPIC.EQ.0)GO TO 3 READ(LIN,*)NEPIC,(IEP(I),I=1,NEPIC) WRITE(LOU,101)NEPIC,(IEP(I),I=1,NEPIC) 3 CONTINUE READ(LIN,*)XMIN,XMAX,XLEN,DTICX,YMIN,YMAX,YLEN,DTICY WRITE(LOU,102)XMIN,XMAX,XLEN,DTICX,YMIN,YMAX,YLEN,DTICY NP=IFIX((YMAX-YMIN)/DT) NPMIN=MIN0(NP,NT) VRED=6. AMP=0. B1=1. EPICS=10. EPS=0. SC=1. READ(LIN,*)VRED,AMP,B1,EPICS,EPS,SC WRITE(LOU,102)VRED,AMP,B1,EPICS,EPS,SC SMAXIM=0. XMX=0. C C COMPUTATION OF SMAXIM C DO 40 I=1,NDST READ(LU8,162)XX,TMIN,AREDUC,NUM IF(NUM.EQ.0)GO TO 40 READ(LU8,109)(IS(M),M=1,NUM) IF(XX.LE.XMIN.OR.XX.GE.XMAX)GO TO 40 IF (MEPIC.EQ.0)GO TO 45 DO 46 J=1,NEPIC IF(I.EQ.IEP(J))GO TO 45 46 CONTINUE GO TO 40 45 CONTINUE TSTART=YMIN-TMIN IF(MRED.EQ.1)TSTART=YMIN+ABS(XX-XSOUR)/VRED-TMIN IPOM=IFIX(TSTART/DT)+1 47 IF(IPOM.GT.0)GO TO 41 IPOM=IPOM+NT GO TO 47 41 IF(IPOM.LE.NT)GO TO 42 IPOM=IPOM-NT GO TO 41 42 CONTINUE IIAUX=NT-IPOM DO 43 J=1,NT IF(J.GE.IPOM)JAUX=J-IPOM+1 IF(J.LT.IPOM)JAUX=J+IIAUX+1 43 SEIS(JAUX)=FLOAT(IS(J)) SMAX(I)=0. DO 44 J=1,NPMIN IF(ABS(SEIS(J)).GT.SMAX(I))SMAX(I)=ABS(SEIS(J)) 44 CONTINUE SMAX(I)=SMAX(I)*AREDUC/999.1 IF(SMAX(I).GT.SMAXIM)XMX=XX IF(SMAX(I).GT.SMAXIM)SMAXIM=SMAX(I) 40 CONTINUE WRITE(LOU,103)XMX,SMAXIM WRITE(LOU,102)(SMAX(I),I=1,NDST) REWIND LU8 READ(LU8,106)MPRINT READ(LU8,106)IPRINT READ(LU8,106)TITLE READ(LU8,152)XSOUR,YSOUR,ZSOUR,TSOUR,RSTEP,DT,DF READ(LU8,101)NDST,NT,MCOMP,ILOC C C END OF COMPUTATION OF SMAXIM C 2 CONTINUE IF(LU4.EQ.0)GO TO 12 WRITE(LU4,106)TEXT WRITE(LU4,100)NDST,MRED,MCOMP,ILOC,VRED,RSTEP,XSOUR,YSOUR,DT WRITE(LU4,110)XMX,SMAXIM 12 CONTINUE XMER=XLEN/(XMAX-XMIN) YMER=YLEN/(YMAX-YMIN) DDX=RSTEP*XMER ELM=.45*SC CALL BORDER(XLEN,DTICX,YLEN,DTICY,SC,TEXT,0,XMIN,XMAX, 1YMIN,YMAX,NTICX,NTICY,NDX,NDY) T=.5*(XLEN-6.3*SC) CALL SYMBOL(T,-1.6*SC,ELM,'DISTANCE IN KM',0.,14) T=.5*(YLEN-8.1*SC) U=-(1.6+.4*NDX)*SC IF(MRED.EQ.0) 1CALL SYMBOL(U,T,ELM,'TRAVEL TIME IN SEC',90.,18) IF(MRED.EQ.0)GO TO 9 CALL SYMBOL(U,T,ELM,'T-D/ ',90.,5) T=T+1.8*SC CALL NUMBER(U,T,ELM,VRED,90.,2) T=T+2.7*SC CALL SYMBOL(U,T,ELM,'(IN SEC)',90.,8) 9 CONTINUE IF(MCOMP.EQ.0)CALL SYMBOL(ELM,YLEN+SC,ELM,'VERTICAL',0.,8) IF(MCOMP.EQ.1)CALL SYMBOL(ELM,YLEN+SC,ELM,'X-COMPONENT',0.,11) IF(MCOMP.EQ.2)CALL SYMBOL(ELM,YLEN+SC,ELM,'Y-COMPONENT',0.,11) T=XLEN-7.5*SC CALL NUMBER(T,YLEN+.5*SC,.3*SC,AMP,0.,0) T=T+1.5*SC CALL NUMBER(T,YLEN+.5*SC,.3*SC,B1,0.,2) T=T+1.5*SC CALL NUMBER(T,YLEN+.5*SC,.3*SC,EPS,0.,1) T=T+1.5*SC CALL NUMBER(T,YLEN+.5*SC,.3*SC,SMAXIM,0.,5) C C LOOP FOR THE RECEIVER POSITIONS C 4 DO 10 I=1,NDST READ(LU8,162)XX,TMIN,AREDUC,NUM IF(IPR.GT.1)WRITE(LOU,162)XX,TMIN,AREDUC,NUM IF(NUM.EQ.0)GO TO 10 READ(LU8,109)(IS(M),M=1,NUM) IF(LU4.NE.0)WRITE(LU4,107)XX,AREDUC,TMIN,NUM IF(LU4.NE.0)WRITE(LU4,109)(IS(M),M=1,NUM) IF(XX.LE.(XMIN+0.001).OR.XX.GE.(XMAX-0.001))GO TO 10 IF(MEPIC.EQ.0)GO TO 5 DO 6 J=1,NEPIC IF(I.EQ.IEP(J))GO TO 5 6 CONTINUE GO TO 10 C C SHIFTING SEISMOGRAM INTO A REQUIRED TIME POSITION C 5 TSTART=YMIN-TMIN IF(MRED.EQ.1)TSTART=YMIN+ABS(XX-XSOUR)/VRED-TMIN IPOM=IFIX(TSTART/DT)+1 TL=TMIN+DT*FLOAT(IPOM-1) 37 IF(IPOM.GT.0)GO TO 31 IPOM=IPOM+NT GO TO 37 31 IF(IPOM.LE.NT)GO TO 32 IPOM=IPOM-NT GO TO 31 32 CONTINUE IIAUX=NT-IPOM DO 33 J=1,NT IF(J.GE.IPOM)JAUX=J-IPOM+1 IF(J.LT.IPOM)JAUX=J+IIAUX+1 33 SEIS(JAUX)=AREDUC*FLOAT(IS(J))/999.1 SMAX(I)=0. DO 39 J=1,NPMIN IF(ABS(SEIS(J)).GT.SMAX(I))SMAX(I)=ABS(SEIS(J)) 39 CONTINUE SMAXI=SMAX(I) IF(I.EQ.1)SMAX1=SMAXI C C SYNTHETIC SEISMOGRAM SEIS(J),J=1,NPMIN, CORRESPONDS TO TRAVEL C TIMES (OR REDUCED TRAVEL TIMES) FROM YMIN TO YMAX C C C COMPUTATION OF SCALING FACTORS C IF(SMAXI.LT.0.000001)GO TO 7 IF(ABS(AMP).LT.0.00001)FACTOR=B1*DDX/SMAXI IF(ABS(AMP).LT.0.00001)GO TO 8 7 FACTOR=0. IF(ABS(EPS).GT.0.00001)GO TO 20 IF(AMP.LT.(-0.00001))FACTOR=B1*DDX/SMAXIM IF(AMP.GT.0.00001.AND.AMP.LT.5.)FACTOR=B1 IF(AMP.GT.5.)FACTOR=B1*DDX/SMAX1 GO TO 8 20 IF(AMP.LT.0.00001)FACTOR=B1*DDX*((ABS(XX-XSOUR)/EPICS)**EPS) 1/SMAXIM IF(AMP.GT.0.00001)FACTOR=B1*(ABS(XX-XSOUR)/EPICS)**EPS 8 CONTINUE SFMAX=FACTOR*SMAXI IF(IPR.EQ.1)WRITE(LOU,120)XX,SMAXI,FACTOR,SFMAX IF(IPR.LE.1)GO TO 25 WRITE(LOU,121)XX,TSTART,SMAXI,FACTOR,SFMAX DO 26 J=1,NPMIN 26 IS(J)=IFIX(999.1*SEIS(J)/SMAXI) WRITE(LOU,109)(IS(J),J=1,NPMIN) 25 CONTINUE C C PLOTTING C XM=XMIN X0=(XX-XM)*XMER XNEW=X0 YNEW=0. CALL PLOT(XNEW,YNEW,3) IOLD=1 AMPL=SEIS(1)+(SEIS(2)-SEIS(1))*(TL-YMIN)/DT XNEW=X0-FACTOR*AMPL IF(XNEW.LT.0..OR.XNEW.GT.XLEN)GO TO 15 CALL PLOT(XNEW,YNEW,2) 15 CONTINUE DO 11 J=2,NPMIN INEW=0 U1=FACTOR*SEIS(J) IF(J.EQ.NPMIN)GO TO 50 IF(ABS(U1).LT.0.003*SFMAX)INEW=1 IFUT=0 U11=FACTOR*SEIS(J+1) IF(ABS(U11).LT.0.003*SFMAX)IFUT=1 IF(IOLD.EQ.1.AND.INEW.EQ.1.AND.IFUT.EQ.1)GO TO 11 50 XNEW=X0-U1 YNEW=(TL-YMIN+DT*FLOAT(J-1))*YMER CALL PLOT(XNEW,YNEW,2) 11 IOLD=INEW AMPL=SEIS(NPMIN)+(SEIS(NPMIN+1)-SEIS(NPMIN))*(TL-YMIN)/DT XNEW=X0-FACTOR*AMPL IF(XNEW.LT.0..OR.XNEW.GT.XLEN)GO TO 10 CALL PLOT(XNEW,YLEN,2) 10 CONTINUE C C END OF THE LOOP FOR RECEIVER POSITIONS C 100 FORMAT(4I5,5F10.5) 101 FORMAT(16I5) 102 FORMAT(8F10.5) 103 FORMAT(/,'SMAXIM=',F15.5,' AT XMX=' ,F10.5/) 104 FORMAT(I5,1X,A) 106 FORMAT(A) 107 FORMAT(F10.5,E15.8,F10.5,I5) 108 FORMAT(1X,A) 109 FORMAT(20I4) 110 FORMAT(22X,F10.5,9X,F15.9) 111 FORMAT(20X,F15.9) 112 FORMAT(6F15.9) 120 FORMAT(F10.5,3E15.4) 121 FORMAT(/1X,'SYNTHETIC SEISMOGRAM',2X,'X=',F10.5,2X,'T0=',F10.5, 1/1X,'SMAX=',1E15.6,2X,'FACTOR=',1E15.6,2X,'SFMAX=',F10.5) 152 FORMAT(5F10.5,2E15.7) 154 FORMAT(3E16.8) 162 FORMAT(2F10.3,1E12.5,I5) 99 CALL PLOT(0.,0.,999) STOP END C C======================================================================= C INCLUDE 'border.for' C border.for INCLUDE 'error.for' C error.for INCLUDE 'length.for' C length.for INCLUDE 'sep.for' C sep.for INCLUDE 'calcops.for' C calcops.for C C======================================================================= C