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)
   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=1,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 'calcops.for'
C     calcops.for
C
C=======================================================================
C