C 
C     P R O G R A M   P O L A R P L O T
C     *********************************
C
C     PROGRAM POLARPLOT IS DESIGNED FOR THE PLOTTING OF POLAR
C     DIAGRAMS FROM THE FILES GENERATED BY PROGRAMS SYNTAN OR
C     BPLOT
C
C     ************************************************************
C
      CHARACTER*80 TEXT,TXT,PSTEXT,FILEIN,FILEOU,FILE1,FILE2
      DIMENSION SEISA(3001),SEISB(3001),IEP(100)
C
C
C
C**************************************************
C
      LIN=5
      LOU=6
      LU4A=1
      LU4B=2
      FILEIN='polar.dat'
      FILEOU='polar.out'
      FILE1='lu4a.dat'
      FILE2='lu4b.dat'
      WRITE(*,'(2A)') ' SPECIFY NAMES OF INPUT AND OUTPUT FILES',
     1' LIN, LOU, LU4A, LU4B: '
      READ(*,*) FILEIN,FILEOU,FILE1,FILE2
      IF(FILE1.EQ.' ') GO TO 99
      IF(FILE2.EQ.' ') GO TO 99
      OPEN(LIN,FILE=FILEIN,FORM='FORMATTED',STATUS='OLD')
      OPEN(LOU,FILE=FILEOU,FORM='FORMATTED')
      OPEN(LU4A,FILE=FILE1,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU4B,FILE=FILE2,FORM='FORMATTED',STATUS='OLD')
C
C**************************************************
C
      IRUN=0
      PSTEXT=' '
      IPRINT=0
      XSHIFT=3.
      YSHIFT=6.
      READ(LIN,*)IPRINT,XSHIFT,YSHIFT,PSTEXT
      WRITE(LOU,107)IPRINT,XSHIFT,YSHIFT,PSTEXT
      IF(IPRINT.LT.0)THEN
        CALL PLOTN(PSTEXT,0)
        IPRINT=-IPRINT
      END IF
      CALL PLOTS(LDUM1,LDUM2,7)
      SHF=XSHIFT
      CALL PLOT(SHF,YSHIFT,-3)
C
    2 MCONT=0
      MEPIC=0
      NTICX=1
      NTEXT=0
      NVER=0 
      READ(LIN,*)MCONT,MEPIC,NTICX,NTEXT,NVER
      WRITE(LOU,101)MCONT,MEPIC,NTICX,NTEXT,NVER
      IF(MCONT.EQ.0)GO TO 99
      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
      SC=1.
      TSTART=0.
      TFIN=10.
      AMP=0.
      B1=1.      
      TXT='POLAR'
      READ(LIN,*)XLEN,DTICX,SC,TSTART,TFIN,AMP,B1
      WRITE(LOU,102)XLEN,DTICX,SC,TSTART,TFIN,AMP,B1
      READ(LIN,*)TXT
      WRITE(LOU,100)TXT
      REWIND LU4A
      READ(LU4A,100)TEXT
      WRITE(LOU,100)TEXT
      READ(LU4A,105)MDISTA,MRED,MCOMPA,ILOC,VRED,RSTEP,XSOUR,YSOUR,DT
      WRITE(LOU,105)MDISTA,MRED,MCOMPA,ILOC,VRED,RSTEP,XSOUR,YSOUR,DT
      READ(LU4A,104)XMXA,AMAXIM
      WRITE(LOU,104)XMXA,AMAXIM
      REWIND LU4B
      READ(LU4B,100)TEXT
      WRITE(LOU,100)TEXT
      READ(LU4B,105)MDISTB,MRED,MCOMPB,ILOC,VRED,RSTEP,XSOUR,YSOUR,DT
      WRITE(LOU,105)MDISTB,MRED,MCOMPB,ILOC,VRED,RSTEP,XSOUR,YSOUR,DT
      READ(LU4B,104)XMXB,BMAXIM
      WRITE(LOU,104)XMXB,BMAXIM
      IF(MDISTA.NE.MDISTB)WRITE(LOU,108)
      IF(MDISTA.NE.MDISTB)GO TO 99
      SMAXIM=AMAXIM
      IF(SMAXIM.LT.BMAXIM)SMAXIM=BMAXIM
C
C     LOOP FOR THE RECEIVER POSITIONS
C
      DO 10 I=1,MDISTA
      READ(LU4A,110)XX,SMAXIA,TMINA,NPTSA
      READ(LU4A,109)(SEISA(M),M=1,NPTSA)
      READ(LU4B,110)XX,SMAXIB,TMINB,NPTSB
      READ(LU4B,109)(SEISB(M),M=1,NPTSB)
      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
    5 SAUX=SMAXIA/999.
      if(nver.eq.1.and.mcompb.eq.0)saux=-saux
      DO 22 M=1,NPTSA
   22 SEISA(M)=SEISA(M)*SAUX
      SAUX=SMAXIB/999.
      DO 23 M=1,NPTSB
   23 SEISB(M)=SEISB(M)*SAUX
C
C     PLOT OF FRAME
      XMER=.5*XLEN
      DDX=XMER
      IF(IRUN.NE.0)SHF=SHF+2.*XMER+XSHIFT
      IF(IRUN.NE.0)CALL PLOT(2.*XMER+XSHIFT,0.,-3)
      IRUN=1
      IF(NTEXT.NE.0)CALL BORDER(XLEN,DTICX,XLEN,DTICX,SC,TXT,nver,-1.,
     11.,-1.,1.,NTICX,NTICX,0,0)
      IF(NTEXT.EQ.0)CALL BORDER(XLEN,DTICX,XLEN,DTICX,SC,TEXT,nver,
     1-1.,1.,-1.,1.,NTICX,NTICX,0,0)
      ELM=.45*SC
      T=.5*(XLEN-6.*ELM)
      IF(MCOMPA.EQ.1)
     1CALL SYMBOL(T,-1.6*SC,ELM,'RADIAL',0.,6)
      T=T-ELM
      IF(MCOMPA.EQ.0)
     1CALL SYMBOL(T,-1.6*SC,ELM,'VERTICAL',0.,8)
      T=T-ELM
      IF(MCOMPA.EQ.2)
     1CALL SYMBOL(T,-1.6*SC,ELM,'TRANSVERSE',0.,10)
      U=-(1.6+.4)*SC
      IF(MCOMPB.EQ.2)
     1CALL SYMBOL(U,T,ELM,'TRANSVERSE',90.,10)
      T=T+ELM
      IF(MCOMPB.EQ.0)
     1CALL SYMBOL(U,T,ELM,'VERTICAL',90.,8)
      T=T+ELM
      IF(MCOMPB.EQ.1)
     1CALL SYMBOL(U,T,ELM,'RADIAL',90.,6)
      CALL PLOT(0.,0.,3)
      if(iloc.ne.1.and.ntext.ge.0)
     1CALL SYMBOL(ELM,XLEN+SC,ELM,'X= ',0.,3)
      if(iloc.eq.1.and.ntext.ge.0)
     1CALL SYMBOL(ELM,XLEN+SC,ELM,'Z= ',0.,3)
      if(ntext.ge.0)then
        CALL NUMBER(4*ELM,XLEN+SC,ELM,XX,0.,2)
        CALL SYMBOL(10.*ELM,XLEN+SC,ELM,'KM, ',0.,4)
        CALL SYMBOL(14.*ELM,XLEN+SC,.5*ELM,'T1,T2= ',0.,7)
        CALL NUMBER(17.5*ELM,XLEN+SC,.5*ELM,TSTART,0.,2)
        CALL NUMBER(20.*ELM,XLEN+SC,.5*ELM,TFIN,0.,2)
      end if
C
C
      SMAXI=SMAXIA
      IF(SMAXI.LT.SMAXIB)SMAXI=SMAXIB
      TMIN=TMINA
      IF(TMIN.GT.TMINB)TMIN=TMINB
      TMAXA=TMINA+(NPTSA-1)*DT
      TMAXB=TMINB+(NPTSB-1)*DT
      TMAX=TMAXA
      IF(TMAX.LT.TMAXB)TMAX=TMAXB
      if(iprint.ge.2)WRITE(LOU,102)TMIN,TMAX,SMAXI,SMAXIM
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
      IF(AMP.LT.(-0.00001))FACTOR=B1*DDX/SMAXIM
      IF(AMP.GT.0.00001)FACTOR=B1
      GO TO 8
    7 FACTOR=0.
    8 CONTINUE
      SFMAX=FACTOR*SMAXI
      SF1=.003*SFMAX
      IF(IPRINT.ge.1)WRITE(LOU,103)XX,SMAXI,FACTOR,SFMAX
C
C
      K=0
      IA=0
      IB=0
      XNEW=0.
      YNEW=0.
      TST=TSTART
      TEND=TFIN
      IF(TST.LT.TMIN)TST=TMIN
      IF(TEND.GT.TMAX)TEND=TMAX
      IF(TST.LT.TMINA)XNEW=0.
      IF(TST.LT.TMINA)GO TO 14
      IA=(TST-TMINA)/DT+1
      T=TMINA+DT*FLOAT(IA-1)
      TM=T
      AMPL=SEISA(IA)+(SEISA(IA+1)-SEISA(IA))*(TST-T)/DT
      XNEW=FACTOR*AMPL
      IF(ABS(XNEW).GT.XMER)GO TO 15
   14 IF(TST.LT.TMINB)YNEW=0.
      IF(TST.LT.TMINB)GO TO 12
      IB=(TST-TMINB)/DT+1
      T=TMINB+DT*FLOAT(IB-1)
      TM=T
      BMPL=SEISB(IB)+(SEISB(IB+1)-SEISB(IB))*(TST-T)/DT
      YNEW=FACTOR*BMPL
      IF(ABS(YNEW).GT.XMER)GO TO 15
   12 CONTINUE
      XNEW=XNEW+XMER
      YNEW=YNEW+XMER
      if(iprint.ge.2)WRITE(LOU,102)XNEW,YNEW,AMPL,BMPL
      CALL PLOT(XNEW,YNEW,3)
   15 CONTINUE
      IF(IA.NE.0)IA=IA+1
      IF(IB.NE.0)IB=IB+1
      IF(ABS(T-TMINA).LT..0001)IA=1
      IF(ABS(T-TMINB).LT..0001)IB=1
      XNEW=0.
      IF(IA.GT.0.AND.IA.LE.NPTSA)XNEW=FACTOR*SEISA(IA)
      YNEW=0.
      IF(IB.GT.0.AND.IB.LE.NPTSB)YNEW=FACTOR*SEISB(IB)
      IF(ABS(XNEW).GT.XMER)GO TO 15
      IF(ABS(YNEW).GT.XMER)GO TO 15
      XNEW=XNEW+XMER
      YNEW=YNEW+XMER
      K=K+1
      T=TM+K*DT
      IF(T.GT.TEND)GO TO 13
      if(iprint.ge.2)WRITE(LOU,106)K,IA,IB,T,XNEW,YNEW
      CALL PLOT(XNEW,YNEW,2)
      GO TO 15
   13 XNEW=0.
      IF(T.GT.TMAXA)GO TO 11
      AMPL=SEISA(IA-1)+(SEISA(IA)-SEISA(IA-1))*(T-TEND)/DT
      XNEW=FACTOR*AMPL
      IF(ABS(XNEW).GT.XMER)GO TO 10
   11 YNEW=0.
      IF(T.GT.TMAXB)GO TO 9
      BMPL=SEISB(IB-1)+(SEISB(IB)-SEISB(IB-1))*(T-TEND)/DT
      YNEW=FACTOR*BMPL
      IF(ABS(YNEW).GT.XMER)GO TO 10
    9 XNEW=XNEW+XMER
      YNEW=YNEW+XMER
      CALL PLOT(XNEW,YNEW,2)
   10 CONTINUE
      call plot(-shf,0.,-3)
C
C     END OF THE LOOP FOR RECEIVER POSITIONS
C
      GO TO 2
C
C
  100 FORMAT(A)
  101 FORMAT(16I5)
  102 FORMAT(8F10.5)
  103 FORMAT(2X,4E15.5)
  104 FORMAT(22X,F10.5,9X,E15.9)
  105 FORMAT(4I5,5F10.5)
  106 FORMAT(3I5,4F10.5)
  107 FORMAT(I5,2F10.5,1X,A)
  108 FORMAT(/1X,'DIFFERENT SELECTION OF RANGES ON THE AXES,
     1 COMPUTATION TERMINATED'//)
  109 FORMAT(20F4.0)
  110 FORMAT(F10.5,E15.8,F10.5,I5)
   99 CALL PLOT(0.,0.,999)
C
C
      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