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: '
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.
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
WRITE(LOU,101)ICONT,NDST,ILOC
IF(ICONT.EQ.0)WRITE(LOU,107)LU1
IF(ICONT.EQ.0)GO TO 99
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'
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
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.
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
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
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
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
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
IF(NRAY.NE.0) WRITE(LOU,104)(NR(I),I = 1,NRAY)
I = 1
K = 1
IPL = 0
IF(KPRINT.GE.2)WRITE(LOU,106)N,IND
IF(N.EQ.0)GOTO 24
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
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
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.
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
WRITE(LOU,104)NEXP
IF(NEXP.EQ.0)GOTO 14
NS=NEXP
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
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
WRITE(LOU,104)NEXP
IF(NEXP.EQ.0)GOTO 18
NS=NEXP
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 'calcops.for'
C     calcops.for
C
C=======================================================================
C                                                                 ```