C
C PROGRAM R A Y P L O T C *********************** C C THE PROGRAM RAYPLOT IS DESIGNED FOR PLOTTING OF RAY DIAGRAMS, C TRAVEL TIMES AND AMPLITUDES OF SEISMIC BODY WAVES FROM THE C FILE LU GENERATED BY PROGRAM SEIS. C C **************************************************************** C CHARACTER*80 TEXT,FILEIN,FILEOU,FILE1 DIMENSION A(30),B(30),C(30),D(30),X1(30),III(30),NPNT(20) DIMENSION X(400),T(400),XZ(400),NR(40),AX(400),AY(400),AZ(400), 1Y(400),TAS(400),ANG(400),PHX(400),PHY(400),PHZ(400),INDI(400) COMMON/SOUR/ROS,VPS,VSS C C *** LIN=5 LOU=6 LU=1 FILEIN='raypl.dat' FILEOU='raypl.out' WRITE(*,'(2A)') ' (RAYPLOT) 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(LU,FILE=FILE1,FORM='FORMATTED',STATUS='OLD') CALL PLOTS(LDUM1,LDUM2,7) C *** C XB = 0. IRUN = 0 C C *** C REWIND LU CALL PLOT(4.,5.,-3) 200 continue NEXT=0 READ(lin,*)NEXT,ISHIFT,IPRINT IF(NEXT.EQ.0)GO TO 99 WRITE(lou,106)NEXT,ISHIFT,IPRINT IF(ISHIFT.EQ.0)ISHIFT=8 SHIFT=FLOAT(ISHIFT) READ(LU,101)ICONT,itpr WRITE(lou,101)ICONT,itpr C C *** C *** C IF(ICONT.EQ.0)WRITE(lou,107)LU IF(ICONT.EQ.0)GO TO 99 READ(lin,*)TEXT WRITE(lou,100)TEXT READ(lin,*)NTICX,NTICY,NTICT,NTICA,INDO,INDT,NRAY,IBOUND, 1IRED,IRS,NDX,NDY WRITE(lou,106)NTICX,NTICY,NTICT,NTICA,INDO,INDT,NRAY,IBOUND, 1IRED,IRS,NDX,NDY C C *** C *** C IF(NTICX.EQ.0)GO TO 99 if(itpr.eq.22)ired=0 IF(INDT.EQ.1.OR.INDT.EQ.2)IRED=0 IF(INDT.EQ.0)INDT=3 XB = 0. IF(IBOUND.EQ.0) IBOUND =100 READ(lin,*)XMIN,XMAX,XLEN,DTICX,SC WRITE(lou,102)XMIN,XMAX,XLEN,DTICX,SC IF(ABS(SC).LT..00001)SC=1. IF(ABS(XMAX-XMIN).LT..00001)GO TO 32 XMER = XLEN/(XMAX-XMIN) GLH = 1.5 IF(NTICY.EQ.0)GO TO 32 READ(lin,*)YMIN,YMAX,YLEN,DTICY WRITE(lou,102)YMIN,YMAX,YLEN,DTICY YMER = YLEN/(YMAX-YMIN) C C PLOTTING OF BORDER FOR RAY DIAGRAM C IF(IRUN.EQ.1)CALL PLOT(XLEN+SHIFT,0.,-3) IRUN=1 72 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 C C PLOTTING OF INTERFACES 32 READ(LU,101)NINT,(NPNT(I),I=1,NINT) IF(IPRINT.EQ.3)WRITE(lou,106)NINT,(NPNT(I),I = 1,NINT) IB=IBOUND IF(IBOUND.LT.0)IB=-IBOUND ib=ib+n RD = (XMAX-XMIN)/FLOAT(IB-1) YM = YMAX DO 5 I = 1,NINT N = NPNT(I)-1 READ(LU,111)(A(J),B(J),C(J),D(J),X1(J),III(J),J=1,N) IF(IPRINT.EQ.3)WRITE(lou,105)(A(J),B(J),C(J),D(J),X1(J),III(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 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 IF(III(K1).EQ.0.OR.III(K1).EQ.(-1))CALL PLOT(XSV,YSV,2) IF(III(K1).EQ.(-2).OR.III(K1).GT.0)CALL PLOT(XSV,YSV,3) RB = RB+RD if((k1+2).le.n.and.rb.gt.x1(k1+2))rb=x1(k1+2) 4 continue IF(IBOUND.LT.0)YM=YM-0.03/YMER 5 CONTINUE 20 READ(LU,102)X0,Y0,ROS,VPS,VSS IF(IPRINT.EQ.3)WRITE(lou,102)X0,Y0,ROS,VPS,VSS IF(NTICY.EQ.0)GOTO 21 C 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(LU,112)N,IND IPL = 0 WRITE(lou,106)N,IND IF(N.EQ.0)GOTO 24 READ(LU,113)(X(J),Y(J),J = 1,N) IF(IPRINT.EQ.3)WRITE(lou,102)(X(J),Y(J),J = 1,N) IF(INDO.NE.0.AND.IND.NE.INDO)GO TO 21 IF(NTICY.EQ.0)GO TO 21 IF(NRAY.EQ.0)GOTO 25 IF(NR(K).EQ.I)GOTO 22 25 DO 10 J = 1,N XX = X(J) YY = Y(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-YMIN))*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 TG=(YY-Y(J-1))/(XX-X(J-1)) 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 TG=(XX-X(J-1))/(YY-Y(J-1)) XX = TG*(YYY-YY)+XX YY = YYY IPL=0 GOTO 6 10 CONTINUE GOTO 23 22 IF(K.GE.NRAY)GOTO 23 K = K+1 23 I = I+1 GOTO 21 24 CONTINUE IF(NTICY.EQ.0)GO TO 11 XNEW = X0-XMIN XNEW = XNEW*XMER YNEW = YMAX-(Y0-YMIN) YNEW = YNEW*YMER IF(XNEW.LT.0..OR.XNEW.GT.XLEN.OR.YNEW.LT.0..OR.YNEW.GT. 1YLEN) GO TO 33 RD=.3*SC CALL PLOT(XNEW+RD,YNEW+RD,3) CALL PLOT(XNEW,YNEW,2) CALL PLOT(XNEW+RD,YNEW-RD,2) CALL PLOT(XNEW-RD,YNEW-RD,3) CALL PLOT(XNEW,YNEW,2) CALL PLOT(XNEW-RD,YNEW+RD,2) CALL PLOT(XNEW,YNEW,3) 33 CONTINUE C C PLOTTING OF TIME-DISTANCE CURVE C 11 READ(LU,112)NS IF(IPRINT.GE.2)WRITE(lou,106)NS NWTYP=NS IF(NS.LT.0)NS=-NS IF(NS.NE.0)READ(LU,114)(INDI(I),X(I),T(I),TAS(I), 1ANG(I),AX(I),AY(I),AZ(I),PHX(I),PHY(I),PHZ(I),I=1,NS) IF(NS.NE.0.AND.IPRINT.GE.2)WRITE(lou,108)(INDI(I),X(I),T(I), 1TAS(I),ANG(I),AX(I),AY(I),AZ(I),PHX(I),PHY(I),PHZ(I),I=1,NS) NSS = NS 35 CONTINUE IF(NTICT.EQ.0)GOTO 15 READ(lin,*)TMIN,TMAX,TLEN,DTICT,VRED WRITE(lou,102)TMIN,TMAX,TLEN,DTICT,VRED IF(IRUN.EQ.1)CALL PLOT(0.,YLEN+SHIFT,-3) IRUN=1 IF(NSS.EQ.0)GOTO 30 IEX=0 if(itpr.ne.22)go to 40 XMER=XLEN/(YMAX-YMIN) DTICX=DTICY XMIN=YMIN XMAX=YMAX NTICX=NTICY 40 CONTINUE YMER = TLEN/(TMAX-TMIN) IF(ABS(VRED).LT..00001) VRED = 6. NAUX=3 IF(INDT.EQ.1.OR.INDT.EQ.2)NAUX=5 C C CALL BORDER(XLEN,DTICX,TLEN,DTICT,SC,TEXT,0,XMIN,XMAX,TMIN, 1TMAX,NTICX,NTICT,NDX,NDY) TX=.5*(XLEN-6.3*SC) if(itpr.ne.22) 1CALL SYMBOL(TX,-1.6*SC,.45*SC,'DISTANCE IN KM',0.,14) if(itpr.eq.22) 1CALL SYMBOL(TX,-1.6*SC,.45*SC,'DEPTH IN KM',0.,11) TX=.5*(YLEN-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 CONTINUE C C INDEX = 9 12 DO 13 I = 1,NS IF(IEX.EQ.0.AND.INDI(I).NE.INDT)GO TO 13 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-X0)/VRED XNEW = (XNEW-XMIN)*XMER IF(YNEW.LT.TMIN.OR.YNEW.GT.TMAX)GOTO 13 YNEW = (YNEW-TMIN)*YMER CALL SYMBOL(XNEW,YNEW,.15*sc,char(INDEX),0.,-index) 13 CONTINUE IF(IEX.EQ.0)GOTO 30 NS = NSS GOTO 14 30 CONTINUE NEXP=0 READ(lin,*)NEXP WRITE(lou,104)NEXP IF(NEXP.EQ.0)GOTO 14 C C 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 = 2 GOTO 12 14 CONTINUE CALL PLOT(0.,-YLEN-SHIFT,-3) C C C PLOTTING OF AMPLITUDE-DISTANCE CURVE C 15 IF(NTICA.EQ.0)GO TO 200 IRUN1=0 ALBACK=0. 19 CONTINUE ALEN=0. READ(lin,*)AMIN,AMAX,ALEN,DTICA,FREQ,KABS,ICOMP,MSOUR IF(ABS(ALEN).LT..00001)CALL PLOT(0.,-ALBACK,-3) IF(ABS(ALEN).LT..00001)GO TO 200 WRITE(lou,109)AMIN,AMAX,ALEN,DTICA,FREQ,KABS,ICOMP,MSOUR CALL SOURCE(lin,lou,0,0,MSOUR,0.,ANGLE,AMSOUR,PHSOUR) IF(KABS.EQ.0)GO TO 38 38 IF(NSS.EQ.0)GO TO 36 IF(INDT.EQ.4)GO TO 36 IF(IRUN.EQ.1.AND.IRUN1.EQ.0)CALL PLOT(XLEN+SHIFT,0.,-3) IF(IRUN1.EQ.1)CALL PLOT(0.,ALEN+SHIFT,-3) IF(IRUN1.EQ.1)ALBACK=ALBACK+ALEN+SHIFT IRUN1=1 IRUN=1 IEX=0 YMER=ALEN/(AMAX-AMIN) NAUX=2 IF(INDT.EQ.1.OR.INDT.EQ.2)NAUX=4 C C CALL BORDER(XLEN,DTICX,ALEN,DTICA,SC,TEXT,0,XMIN,XMAX,AMIN, 1AMAX,NTICX,NTICA,NDX,NDY) TX=.5*(XLEN-6.3*SC) if(itpr.ne.22) 1CALL SYMBOL(TX,-1.6*SC,.45*SC,'DISTANCE IN KM',0.,14) if(itpr.eq.22) 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,'RADIAL',0.,6) IF(ICOMP.EQ.2) 1CALL SYMBOL(.45*SC,ALEN+SC,.45*SC,'TRANSVERSE',0.,10) C C 28 INDEX=9 16 DO 17 I = 1,NS IF(IEX.EQ.0.AND.INDI(I).NE.INDT)GO TO 17 XNEW = X(I) YNEW=AZ(I) IF (ICOMP.EQ.1)YNEW=AX(I) IF (ICOMP.EQ.2)YNEW=AY(I) MWAVE=1 IF(NWTYP.LT.0)MWAVE=2 IF(NWTYP.LT.0.AND.ICOMP.EQ.2)MWAVE=3 C C COMPUTATION OF SOURCE EFFECTS C ***************************** ANGLE=ANG(I) CALL SOURCE(lin,lou,1,MWAVE,MSOUR,0.,ANGLE,AMSOUR,PHSOUR) YNEW=YNEW*AMSOUR C ***************************** C IF(IEX.EQ.0.OR.NEXP.EQ.0)GO TO 41 XNEW=XZ(I) YNEW=T(I) GO TO 31 C C COMPUTATION OF ABSORPTION EFFECTS C ********************************* 41 IF(KABS.GT.0)TAST=TAS(I) IF(KABS.GT.0)YNEW=YNEW*EXP(-3.1415926*FREQ*TAST) C ************************************************** C IF(IPRINT.EQ.1)WRITE(lou,110)X(I),YNEW 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*sc,char(INDEX),0.,-index) 17 CONTINUE 36 IF(IEX.EQ.1)GOTO 18 READ(lin,*)NEXP WRITE(lou,104)NEXP IF(NEXP.EQ.0)GOTO 18 C C 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 IF(INDT.EQ.1.OR.INDT.EQ.2.OR.INDT.EQ.4)GO TO 18 IEX=1 INDEX = 2 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(2X,I2,3(2F7.2,E10.3)) 104 FORMAT(2X,26I3) 105 FORMAT(5E12.6,I5) 106 FORMAT(16I5) 107 FORMAT(2X,'END OF FILE',I3) 108 FORMAT(I5,4F10.5,3E15.9,3F10.5) 109 FORMAT(5F10.5,4I5) 110 FORMAT(F10.5,E15.9) 111 format(5e15.5,i5) 112 format(2i5) 113 format(2e15.5) 114 format(i5,10e15.5) C 99 CONTINUE C C *** C *** C CALL PLOT(0.,0.,999) STOP END C C C======================================================================= C INCLUDE 'source.for' C source.for INCLUDE 'border.for' C border.for INCLUDE 'error.for' C error.for INCLUDE 'sep.for' C sep.for INCLUDE 'length.for' C length.for INCLUDE 'calcops.for' C calcops.for C C======================================================================= C