C
C     PROGRAM FRESAN
C     **************
C
      CHARACTER*80 MTEXT,ITEXT,FILEIN,FILEOU,FILE1,FILE2
      COMPLEX IMAG,CTIME,CAUX,cax,cay,caz,cax1,cay1,caz1,csour,aaux
      DIMENSION ISEL(100),iep(100),III(100),AMFR(1000),AMFI(1000),
     1aaux(3),p(3),pol(3),pol1(3)
      COMMON/PLOT1/RESP(30000),IFREQ(100),AMPL(100),DST(100),MCOMP
      COMMON/SOUR/ROS
      COMMON/SH/MSH,NSH
C
C
C**************************************************
C
      LIN=5
      LOU=6
      LU2=1
      LU8=2
      FILEIN='fresan.dat'
      FILEOU='fresan.out'
      FILE1='lu2.out'
      FILE2='lu8.out'
      WRITE(*,'(2A)') ' Specify names of input and output files',
     1' LIN, LOU, LU2, LU8: '
      READ(*,*) FILEIN,FILEOU,FILE1,FILE2
      IF(FILE1.EQ.' ') GO TO 99
      IF(FILE2.EQ.' ') LU8=0
      OPEN(LIN,FILE=FILEIN,FORM='FORMATTED',STATUS='OLD')
      OPEN(LOU,FILE=FILEOU,FORM='FORMATTED')
      OPEN(LU2,FILE=FILE1,FORM='FORMATTED',STATUS='OLD')
      IF(LU8.NE.0)OPEN(LU8,FILE=FILE2,FORM='FORMATTED')
C
C**************************************************
C
      CALL PLOTS(LDUM1,LDUM2,7)
C
      READ(lin,104)ITEXT
      READ(lin,101)IPRINT,JPRINT
      WRITE(lou,105)ITEXT
      WRITE(lou,101)IPRINT,JPRINT
      REWIND LU2
      IF(LU8.NE.0)REWIND LU8
C
      NEPIC=0
      NSPR=30000
      IMAG=(0.,1.)
      PI=4.*ATAN(1.)
C
C     READ INPUT PARAMETERS
C
      READ(lin,101)MCOMP,MABS,MSOUR,MSELec,mepic,mfr,mplot
      WRITE(lou,101)MCOMP,MABS,MSOUR,MSELec,mepic,mfr,mplot
      READ(lin,100)RFR,GM,qred,freqm,arot
      WRITE(lou,100)RFR,GM,qred,freqm,arot
      csrt=cos(arot)
      snrt=sin(arot)
      IF(ABS(rfr).LT.0.00001)rfr=1.
      IF(ABS(qred).LT.0.00001)qred=1.
C
C     FREQUENCY filtering DATA
C
      READ(lin,107)NFREQ,FL,FR,FD
      WRITE(lou,107)NFREQ,FL,FR,FD
      IF(NFREQ.GT.0)FD=1./FD/FLOAT(NFREQ)
      NF1=INT(FL/FD-0.5)
      NF=INT(FR/FD+0.5)-NF1
      FL=FD*FLOAT(NF1)+FD
      FR=FL+FD*FLOAT(NF)-FD
      FDA=6.2831853*FD
C
c     selection of receivers and elementary waves
c
      if(mepic.eq.0)go to 1
      read(lin,101)nepic,(iep(i),i=1,nepic)
      write(lou,101)nepic,(iep(i),i=1,nepic)
    1 if(mselec.eq.0)go to 2
      read(lin,101)nselec,(isel(i),i=1,nselec)
      write(lou,101)nselec,(isel(i),i=1,nselec)
    2 continue
c
C     READING FROM LU2,GENERAL QUANTITIES
C
      READ(LU2,104)MTEXT
      read(lu2,106)ndst,ksh,iloc
      READ(LU2,100)XSOUR,ysour,ZSOUR,TSOUR,rstep,ROS
      if(ndst.ne.1)READ(LU2,100)(dst(i),i=1,ndst)
      if(ndst.eq.1)read(LU2,100)xrec,yrec,xprf,xatan
      WRITE(lou,104)MTEXT
      write(lou,101)ndst,ksh,iloc
      WRITE(lou,100)XSOUR,ysour,ZSOUR,TSOUR,rstep,ROS
      if(ndst.ne.1)write(lou,100)(dst(i),i=1,ndst)
      if(ndst.eq.1)write(lou,100)xrec,yrec,xprf,xatan
      if(ndst.eq.1)dst(1)=xprf
      nrdst=ndst-nepic
C
C     WRITING ON LU8
C
      WRITE(LU8,104)MTEXT
      WRITE(LU8,104)ITEXT
      WRITE(LU8,152)XSOUR,ysour,ZSOUR,TSOUR,RSTEP,FL,FD
      WRITE(LU8,106)NrDST,NF,NF1,MCOMP,iloc
C
C     check of the dimensions of the resp field
C     INITIAL SPECIFICATION OF the resp field
C
      NSUM=2*NDST*NF
      IF(NSUM.GT.NSPR)WRITE(lou,153)
      IF(NSUM.GT.NSPR)GO TO 99
      DO 4 I=1,NSUM
    4 RESP(I)=0.
C
C     INITIAL PARAMETERS FOR RADIATION PATTERNS
C     AND FOR FREQUENCY DEPENDENT FACTORS
C
      if(msour.ne.0)CALL SOURCE(lin,lou,0,0,MSOUR,p,pol,AMSOUR,PHSOUR)
      IF(MFR.NE.0)CALL FDEP(0,G,FL,FD,NF,AMFR,AMFI)
C
C
C     S T A R T   O F   S U M M A T I O N
C     O F   F R E Q U E N C Y  R E S P O N S E
C     ****************************************
C
C     LOOP over termination points of all ELEMENTARY WAVES
C     at all receivers
c
    6 continue
      READ(LU2,103,END=84)NC,idist,time,cax,cay,caz,tast
      if(nc.lt.0)READ(LU2,109)cax1,cay1,caz1
      read(lu2,110)(p(i),i=1,3)
      read(lu2,110)(pol(i),i=1,3)
      if(nc.lt.0)read(lu2,110)(pol1(i),i=1,3)
      caux=cax*csrt+cay*snrt
      cay=-cax*snrt+cay*csrt
      cax=caux
      caux=cax1*csrt+cay1*snrt
      cay1=-cax1*snrt+cay1*csrt
      cax1=caux
      nc1=iabs(nc)
      if(mselec.eq.0)go to 8
C
C     SELECTION OF WAVES
C
      DO 7 j=1,NSELEC
      IF(ISEL(j).eq.nc1)GO TO 6
    7 CONTINUE
    8 if(mepic.eq.0)go to 10
c
c     selection of receivers
c
      do 9 j=1,nepic
      if(iep(j).eq.idist)go to 6
    9 continue
   10 continue
c
C
      if(msour.ne.0)then
        call source(lin,lou,1,1,msour,p,pol,amdir,phdir)
        csour=amdir*cexp(imag*phdir)
      end if
      if(msour.eq.0)csour=(1.,0.)
      aaux(1)=cax*csour
      aaux(2)=cay*csour
      aaux(3)=caz*csour
      if(nc.lt.0)then
        if(msour.ne.0)then
          call source(lin,lou,1,1,msour,p,pol1,amdir,phdir)
          csour=amdir*cexp(imag*phdir)
        end if
        if(msour.eq.0)csour=(1.,0.)
        aaux(1)=aaux(1)+cax1*csour
        aaux(2)=aaux(2)+cay1*csour
        aaux(3)=aaux(3)+caz1*csour
      end if
      if(mcomp.eq.0)k=3
      if(mcomp.eq.1)k=1
      if(mcomp.eq.2)k=2
      rw=real(aaux(k))
      cw=aimag(aaux(k))
      ar=sqrt(rw*rw+cw*cw)
      if(abs(rw).lt..000001.and.abs(cw).lt..000001)then
        ph=0.
      else
        ph=atan2(cw,rw)
      end if
c
      if(iprint.eq.1)write(lou,102)NC,idist,time,ar,ph,amdir,tast
C
C     absorption effects
C
      IF(MABS.EQ.0)GO TO 67
      AUX=Qred*tast
      abs1=.5*aux
      IF(MABS.le.2)GO TO 67
      abs2=1.+alog(freqm/rfr)
      ph=ph+2.*aux*freqm
   67 CONTINUE
C
C     FREQUENCY DEPENDENT EFFECTS
C
      KFR=0
      IF(MABS.LE.1.AND.MFR.EQ.0)GO TO 74
      KFR=1
      DO 73 JF=1,NF
      AMFR(JF)=1.
   73 AMFI(JF)=0.
C
C     FREQUENCY DEPENDENT EFFECTS: SUBROUTINE FDEP
C
      IF(MFR.NE.0)CALL FDEP(1,G,FL,FD,NF,AMFR,AMFI)
C
C     FREQUENCY DEPENDENT EFFECTS:CAUSAL ABSORPTION
C
      IF(MABS.LE.1.or.mabs.eq.3)GO TO 74
      DO 75 JF=1,NF
      FRJ=FD*FLOAT(JF+NF1)
      abs4=2.*aux*frj*alog(frj/rfr)
      caux=(0.,0.)
      IF(MABS.EQ.2)CAUX=CEXP(-IMAG*abs4)
      AUXR=REAL(CAUX)
      AUXI=AIMAG(CAUX)
      AUX1=AMFR(JF)*AUXR-AMFI(JF)*AUXI
      AUX2=AMFR(JF)*AUXI+AMFI(JF)*AUXR
      AMFR(JF)=AUX1
      AMFI(JF)=AUX2
   75 CONTINUE
C
c     FFR computation
c
   74 continue
      sph=sin(ph)
      cph=cos(ph)
      are=ar*cph
      aim=ar*sph
c
c     print of real and imaginary parts of the complex amplitude
c
      if(iprint.eq.2)write(lou,154)are,aim
      ctime=time
      if(mabs.eq.0)go to 83
      ctime=ctime+abs1*imag
      if(mabs.eq.3)ctime=ctime-abs2*aux/pi
   83 tr=real(ctime)
      ti=aimag(ctime)
      TID=TI*FDA
      TRD=TR*FDA
      TIE=TID*FLOAT(NF1)
      TRE=TRD*FLOAT(NF1)
      AUX=EXP(-TIE)
      ARC=AUX*COS(TRE)
      ARS=AUX*SIN(TRE)
      AA1=ARC*ARE-ARS*AIM
      AA2=ARC*AIM+ARS*ARE
      AUX=exp(-tid)
      AR=AUX*COS(TRD)
      AI=AUX*SIN(TRD)
c
c     ar, ai - constant factors for recurent response computation
c     aa1, aa2 - real and imaginary parts of the amplitude
c     corresponding to the lowest considered frequency
c
      IF(IPRINT.EQ.3)WRITE(lou,100)TR,TI,AR,AI,AA1,AA2
C
C     LOOP FOR FREQUENCIES
C
      N=2*(Idist-1)*NF-1
      DO 82 JF=1,nF
      A11=AR*AA1-AI*AA2
      AA2=AR*AA2+AI*AA1
      AA1=A11
      aux=aa1
      auxx=aa2
      IF(KFR.EQ.0)GO TO 81
      AUX=AA1*AMFR(JF)-AA2*AMFI(JF)
      Auxx=AA1*AMFI(JF)+AA2*AMFR(JF)
   81 CONTINUE
      N=N+2
      RESP(N)=RESP(N)+Aux
      RESP(N+1)=RESP(N+1)+Auxx
   82 CONTINUE
C
C     END OF LOOP FOR FREQUENCIES
C
      go to 6
c
c     E N D   O F   C O M P U T A T I O N   O F
C     F R E Q U E N C Y   R E S P O N S E
C     *****************************************
C
C     printing and storing THE FREQUENCY RESPONSE
C
   84 CONTINUE
      DO 95 IR=1,NDST
      K=2*IR*NF
      J=K-2*NF+1
      A=0.1E-20
      DO 91 I=J,K
   91 A=AMAX1(A,ABS(RESP(I)))
      WRITE(LU8,160)DST(IR),A
      if(iprint.eq.1)WRITE(lou,108)DST(IR),A
      A=A/99999.1
   92 DO 93 I=1,12
      III(I)=RESP(J)/A
      IF(J.GE.K)GO TO 94
   93 J=J+1
      WRITE(LU8,161)(III(I),I=1,12)
      IF(IPRINT.EQ.2)WRITE(lou,161)(III(I),I=1,12)
      GO TO 92
   94 WRITE(LU8,161)(III(K),K=1,I)
      IF(IPRINT.EQ.2)WRITE(lou,161)(III(K),K=1,I)
   95 CONTINUE
      REWIND LU8
      REWIND LU2
C
C     PRINTING OF FREQUENCY DEPENDENT AMPLITUDE CURVES
C
      IF(JPRINT.EQ.0)GO TO 350
      READ(lin,101)MTAB,(IFREQ(I),I=1,MTAB)
      WRITE(lou,101)MTAB,(IFREQ(I),I=1,MTAB)
      DO 310 I=1,MTAB
      IFF=IFREQ(I)
      AMAX=0.
      IF(IFF.GT.NF)GO TO 310
      FREQ=FD*FLOAT(NF1+IFF)
      DO 305 IR=1,NDST
      A=RESP(2*(IR-1)*NF+2*IFF-1)
      B=RESP(2*(IR-1)*NF+2*IFF)
      AMPL(IR)=SQRT(A*A+B*B)
      IF(AMPL(IR).GT.AMAX)AMAX=AMPL(IR)
  305 CONTINUE
      B=AMAX/999.1
      DO 306 IJ=1,NDST
      III(IJ)=AMPL(IJ)/B
  306 CONTINUE
      WRITE(lou,170)FREQ,AMAX
      WRITE(lou,171)(III(K),K=1,NDST)
  310 CONTINUE
  350 CONTINUE
C
C     PLOTTING OF FREQUENCY DEPENDENT AMPLITUDE CURVES
C
      IF(MPLOT.EQ.0)GO TO 400
      CALL ZPLOT(NF1,NF,FD,NDST,MPLOT,lin,lou)
  400 CONTINUE
C
C
  100 FORMAT(8F10.5)
  101 FORMAT(16i5)
  102 FORMAT(2I3,F10.5,2e12.6,3f10.6)
  103 FORMAT(2I3,F12.7,6e12.6,f10.6)
  104 FORMAT(A)
  105 FORMAT(1x,A)
  106 FORMAT(26I3)
  107 FORMAT(I5,6F10.5)
  108 format(f10.3,e12.5)
  109 format(6e12.6)
  110 format(3f10.5)
  152 FORMAT(5F10.5,2E15.7)
  153 FORMAT(1X,'DIMENSION OF RESP EXCEEDED')
  154 FORMAT(3E16.8)
  160 FORMAT(F10.3,E12.5)
  161 FORMAT(12I6)
  170 FORMAT(/,1X,'FREQUENCY=',F10.5,2X,'(HZ),  MAXIMUM AMPLITUDE=',
     11E15.6)
  171 FORMAT(20I4)
C
C
   99 CONTINUE
      STOP
      END
C
C
      SUBROUTINE FDEP(N,G,FL,DF,NF,AMFR,AMFI)
      DIMENSION AMFR(1000),AMFI(1000)
      IF(N.EQ.0)RETURN
      DO 1 I=1,NF
      AMFR(I)=1.
      AMFI(I)=0.
    1 CONTINUE
      RETURN
      END
C
C
      SUBROUTINE ZPLOT(NF1,NF,FD,NDST,MPLOT,lin,lou)
      character*80 mtext
      COMMON/PLOT1/RESP(30000),IFREQ(100),AMPL(100),DST(100),MCOMP
C
      irun=0
      shift=5.
C
    2 READ(lin,104)mTEXT
      WRITE(lou,105)mTEXT
      READ(lin,100)XMIN,XMAX,XLEN,DTICX,SC
      WRITE(lou,100)XMIN,XMAX,XLEN,DTICX,SC
      IF(ABS(XLEN).LT..00001)GO TO 99
      IF(ABS(SC).LT..0001)SC=1.
      READ(lin,100)AMIN,AMAX,ALEN,DTICA
      WRITE(lou,100)AMIN,AMAX,ALEN,DTICA
      READ(lin,101)NLOG,NTICX,NTICA,NDX,NDA
      WRITE(lou,101)NLOG,NTICX,NTICA,NDX,NDA
      XMER=XLEN/(XMAX-XMIN)
      AMER=ALEN/(AMAX-AMIN)
      call plot(5.,5.,-3)
C
C     LOOP FOR PLOTS
C
    1 READ(lin,101) MTAB,(IFREQ(I),I=1,MTAB)
      WRITE(lou,101)MTAB,(IFREQ(I),I=1,MTAB)
      IF(MTAB.EQ.0)GO TO 2
      IF(IRUN.EQ.1)CALL PLOT(XLEN+SHIFT,0.,-3)
      IRUN=1
      CALL BORDER(XLEN,DTICX,ALEN,DTICA,SC,mTEXT,0,XMIN,XMAX,
     1AMIN,AMAX,NTICX,NTICA,NDX,NDA)
C
C     LOOP FOR FREQUENCIES
C
      DO 10 I=1,MTAB
      IFF=IFREQ(I)
      AMAXIM=0.
      IF(IFF.GT.NF)GO TO 10
      FREQ=FD*FLOAT(NF1+IFF)
      DO 5 IR=1,NDST
      A=RESP(2*(IR-1)*NF+2*IFF-1)
      B=RESP(2*(IR-1)*NF+2*IFF)
      AMPL(IR)=SQRT(A*A+B*B)
      IF(AMPL(IR).GT.AMAXIM)AMAXIM=AMPL(IR)
      IF(NLOG.EQ.1)AMPL(IR)=ALOG10(AMPL(IR))
    5 CONTINUE
      IF(NLOG.EQ.1)AMAXIM=ALOG10(AMAXIM)
      IEXP=3
      DO 8 IR=1,NDST
      XNEW=(DST(IR)-XMIN)*XMER
      IF(XNEW.LT.0..OR.XNEW.GT.XLEN)GO TO 8
      A=AMPL(IR)
      IF(MPLOT.EQ.2)A=A/AMAXIM
      YNEW=(A-AMIN)*AMER
      IF(IR.EQ.1)GO TO 11
      IF(YNEW.LT.0..OR.YNEW.GT.ALEN)GO TO 6
      IF(YOLD.GE.0..AND.YOLD.LE.ALEN)GO TO 7
      IF(YOLD.LT.0.)YB=0.
      IF(YOLD.GT.ALEN)YB=ALEN
      IEXP=2
    4 AUX1=(DST(IR)-DST(IR-1))*XMER
      AUX2=YNEW-YOLD
      XOLD=XNEW-(YNEW-YB)*AUX1/AUX2
      CALL PLOT(XOLD,YB,IEXP)
      IF(IEXP.EQ.3)GO TO 12
      IEXP=3
      GO TO 7
    6 IF(YOLD.GT.ALEN)GO TO 12
      IF(YNEW.LT.0.)YB=0.
      IF(YNEW.GT.ALEN)YB=ALEN
      IEXP=3
      GO TO 4
   11 IF(YNEW.LT.0..OR.YNEW.GT.ALEN)GO TO 12
    7 CALL PLOT(XNEW,YNEW,IEXP)
   12 YOLD=YNEW
      IEXP=2
    8 CONTINUE
      U=(4.5+1.2*FLOAT(I))*SC
      UU=(4.5+1.8*FLOAT(I))*SC
      CALL NUMBER(U,ALEN+.2*SC,.3*SC,FREQ,0.,1)
      IF(MPLOT.EQ.2)CALL NUMBER(UU,ALEN+.7*SC,.3*SC,AMAXIM,0.,5)
   10 CONTINUE
      CALL SYMBOL(.2*SC,ALEN+.2*SC,.3*SC,'FREQUENCIES(HZ)=',0.,16)
      ELM=.45*SC
      T=.5*(XLEN-6.3*SC)
      CALL SYMBOL(T,-1.6*SC,ELM,'DISTANCE IN KM',0.,14)
      T=.5*(ALEN-7.65*SC)
      U=-(1.6+.4*NDA)*SC
      CALL SYMBOL(U,T,ELM,'A M P L I T U D E',90.,17)
      IF(MCOMP.EQ.0)
     1CALL SYMBOL(.2*SC,ALEN+1.5*SC,ELM,'VERTICAL',0.,8)
      IF(MCOMP.EQ.1)
     1CALL SYMBOL(.2*SC,ALEN+1.5*SC,ELM,'RADIAL',0.,6)
      IF(MCOMP.EQ.2)
     1CALL SYMBOL(.2*SC,ALEN+1.5*SC,ELM,'TRANSVERSE',0.,10)
      IF(MPLOT.EQ.2)
     1CALL SYMBOL(.2*SC,ALEN+.7*SC,.3*SC,'REDUCTION FACTOR=',0.,17)
      GO TO 1
C
C     THE END OF THE LOOP FOR FREQUENCIES
C
  100 FORMAT(8F10.5)
  101 FORMAT(16I5)
  104 FORMAT(A)
  105 FORMAT(1X,A)
   99 CALL PLOT(0.,0.,999)
      RETURN
      END
C
C=======================================================================
C
      INCLUDE 'source.for'
C     source.for
      INCLUDE 'border.for'
C     border.for
      INCLUDE 'calcops.for'
C     calcops.for
C
C=======================================================================
C