c  PROFFILT.FOR
c
c  This program takes 2D (profile) FFT coefficients in PROFFT format and
c  applies various filters:
c
c    1) lowpass
c    2) highpass
c    3) swing-the-tail (experimental)
c    4) upward (downward) continuation
c    5) threshold filter (experimental)
c    6) strip filter
c    7) horizontal derivative
c    8) vertical derivative
c    9) pseudo gravity
c
c  Output is a modified coefficient file.
c  Tien Grauch, modified by Lin Cordell
      dimension line(20),slope(4),b(4),r(4)
      character*50 infile,ofile
1     write(*,'(a,$)')
     1 ' Enter input coef file (from profft) or print stop  '
      read 809,infile
809   format(a50)
101   format(' *',$)
      if(infile.eq.'stop') go to 300
      open(10,file=infile,form='formatted',status='old',err=295)
      write(*,'(a,$)')' Enter output filtered coef file name  '
      read 809,ofile
      open(11,file=ofile,form='formatted',status='unknown')
      read(10,800) ncol,nrow,nz,x0,dx,y0,dy
800   format(7x,i4,2(1x,i1),4e16.5)
      write(11,800) ncol,nrow,nz,x0,dx,y0,dy
      ncnq=float(ncol)*.5e0+1.0000001
      do 5 i=1,2
      read(10,801) line
801   format(20a4)
      write(11,801) line
5      continue
      pi=3.1415927e0
10    print 810
810   format(' 1 - lowpass',15x,'2 - highpass',/,' 3 - swing-the-tail',
     1 8x,'4 - up (or down)ward continue',/,' 5 - threshhold filter',
     1/,' 6 - strip filter',/' 7 - horizontal derivative'/' 8 - vertical
     1 derivative',/' 9 - pseudo gravity', //,' Enter option: ',$)
      read*,ifilt
      if(ifilt.lt.1.or.ifilt.gt.9) go to 10
      go to (11,11,50,75,100,200,400,401,409),ifilt
c
c -- LOWPASS & HIPASS OPTIONS---
c
11    print*
      write(*,'(a,$)')
     1 ' Enter two frequency cut-offs f1,f2 (enter 0,0 for help)  '
      read*, f1,f2
      if(f2.lt.f1) go to 11
      wncnq=.5e0/dx
      if(f1.ge.0.0.and.f1.le.wncnq.and.f2.gt.0.0.and.f2.le.wncnq) go to
     & 12
      print 806
806   format(/' f1<f2 and are of the form k/(ncol*dx) where k is a ',
     1 'wavenumber in the frequency',/,' domain.  Program filters along 
     2 a linear taper between f1 and f2 so that')
      if(ifilt.eq.1) print 807
      if(ifilt.eq.2) print 808
807   format(' all wavelengths (peak-to-peak) less than 1/f2 will be ',
     1 'filtered out.')
808   format(' all wavelengths (peak-to-peak) greater than 1/f2 will be'
     1 ,' filtered out')
	print 875, dx,wncnq
875	format(' f1 and f2 should be between 0 and Nyquist.'/
     1 ' For dx=',g14.5,', Nyquist frequency=',g14.5)
      go to 11
c
c  Read coef file and modify for low and hi pass
c
12    delta=f2-f1
      if(delta.eq.0.0) delta=1.0
      delta=1.e0/delta
c --lowpass parameters--
      xadd=1.0
      filt1=0.0
      filt2=1.0
      if(ifilt.eq.1) go to 15
c --highpass parameters--
      delta=-1.e0*delta
      xadd=0.0
      filt1=1.0
      filt2=0.0
15    do 30 i=1,ncnq
      read(10,803) k,w,fr,fi,amp
803   format(1x,i3,4e16.5)
      filt=filt2
      if(w.gt.f2) filt=filt1
      if(w.le.f2.and.w.ge.f1) filt=xadd-delta*(w-f1)
      fr=fr*filt
      fi=fi*filt
      amp=amp*filt
      write(11,803) k,w,fr,fi,amp
30    continue
      go to 300
c
c -----SWING-THE-TAIL OPTION----
c
c  note that the following values asked for are read directly
c off the amplitude spectrum graph--do not take logs even if the
c plot was semi-log
50    write(6,111)
111   format(1x'Enter intercept of signal and of noise (directly ',
     1 'from graph):'/1x'*'$)
      read*,bs,bn
4     write(6,2)
2     format(' Enter freq. and amp. of pt. where noise line crosses ',
     1 'signal line:'/1x'*'$)
      read*,wc,y
      if(wc.lt.0.0) go to 4
c
c  if cutoff >= nyquist, no filtering
      if(wc.ge.0.5/dx) go to 299
c
c  Read coef file and modify - swing-the-tail
c
      as=log(y/bs)
      as=as/wc
      an=log(y/bn)
      an=an/wc
      c1=as-an
      do 60 i=1,ncnq
      read(10,803) k,w,fr,fi,amp
      s=1.0
      if(w.gt.wc) s=exp(c1*(w-wc))
      fr=fr*s
      fi=fi*s
      amp=amp*s
      write(11,803) k,w,fr,fi,amp
60    continue
      go to 300
c
c ------UPWARD AND DOWNWARD CONTINUATION------
c
75    print*,'Enter depth to continue down (<0 for upward cont.)'
      read*,d
      const=2.0e0*pi*d
      if(d.gt.0.0) print 811
      if(d.lt.0.0) print 812
811   format(/'  ...downward continuing...'/)
812   format(/'  ...upward continuing...'/)
c --read coefs and modify--
      do 80 i=1,ncnq
       read(10,803) k,w,fr,fi,amp
      filt=exp(const*w)
      fr=fr*filt
      fi=fi*filt
      amp=amp*filt
      write(11,803) k,w,fr,fi,amp
80    continue
      go to 300
c
c -----THRESHHOLD FILTER-----
c
c  note that the following values asked for are read directly
c off the amplitude spectrum graph--do not take logs even if the
c plot was semi-log
100   call signal(slope,b,iseg,r)
c  Modify coefs.  DC left unchanged
      read(10,803) k,w,fr,fi,amp
      write(11,803) k,w,fr,fi,amp
      if(iseg.eq.1) go to 150
c  modify for more than one segment defining signal
      do 120 i=1,ncnq-1
      read(10,803) k,w,fr,fi,amp
c  check which segment applicable
      nseg=iseg
      do 105 k=1,iseg-1
      if(w.ge.r(iseg-k)) go to 110
      nseg=nseg-1
105   continue
110   sig=b(nseg)*exp(slope(nseg)*w)
c  change amp if spectrum is above signal
      s=1.0
      if(amp.gt.sig) s=sig/amp
      fr=fr*s
      fi=fi*s
      amp=amp*s
      write(11,803) k,w,fr,fi,amp
120   continue
      go to 300
c  for signal having only one segment
150   do 180 i=1,ncnq-1
      read(10,803) k,w,fr,fi,amp
      sig=b(1)*exp(slope(1)*w)
c  change amp if above signal
      s=1.0
      if(amp.gt.sig) s=sig/amp
      fr=fr*s
      fi=fi*s
      amp=amp*s
      write(11,803) k,w,fr,fi,amp
180   continue
      go to 300
200   continue
c ------STRIP FILTER------
c
201   print*,' Enter half and three-quarter cut frequencies (cy/km),'
      print*,' vh, vtq; (vh < vtq < nyquist):'
      read*,vh,vtq
      if(vtq.gt.(1./(2.*dx)).or.vh.ge.vtq)then
	print*,'Cant handle it.'
	go to 201
      endif
      z=1.098612/(2.*pi*(vtq-vh))
      rs=exp(2.*pi*vh*z)
      print*,'R =',rs,' , z =',z
      const=2.0e0*pi*z
c --read coefs and modify--
      do 202 i=1,ncnq
      read(10,803) k,w,fr,fi,amp
      filt=rs/(rs+exp(const*w))
      fr=fr*filt
      fi=fi*filt
      amp=amp*filt
      write(11,803) k,w,fr,fi,amp
202   continue
      go to 300
400   continue
c ------HORIZONTAL DERIVATIVE------
c
c --read coefs and modify--
      cnst=2.*pi
      do 402 i=1,ncnq
      read(10,803) k,w,fr,fi,amp
      tr=0.
      ti=cnst*w
      frn=-fi*ti
      fin=fr*ti
      fr=frn
      fi=fin
      amp=sqrt(fr*fr+fi*fi)
      write(11,803) k,w,fr,fi,amp
402   continue
      go to 300
401   continue
c ------VERTICAL DERIVATIVE------
c
c --read coefs and modify--
      cnst=2.*pi
      do 403 i=1,ncnq
      read(10,803) k,w,fr,fi,amp
      filt=cnst*w
      fr=fr*filt
      fi=fi*filt
      amp=amp*filt
      write(11,803) k,w,fr,fi,amp
403   continue
      go to 300
409   continue
c ------PSEUDO GRAVITY------
c
      print*,'Enter density and susceptibility:'
      print 101
      read*,rho,susp
      print*,'Enter inclination, declination, and field strength:'
      print 101
      read*,rincl,decl,field
      print*,'Enter profile azimuth (+ east, - west, from north):'
      print 101
      read*,azimuth
      rad=pi/180.
      t1=cos((azimuth-decl)*rad)*cos(rincl*rad)
      t3=sin(rincl*rad)
      c1=t3*t3-t1*t1
      c2=2.*t1*t3
      c1sq=c1*c1
      c2sq=c2*c2
      c3=c1sq+c2sq
      if(c3.eq.0.) then
        print*, 'Cant handle it.'
        stop
      endif
      alpha=(6.67*rho)/(2.*pi*susp*field)
c --read coefs and modify--
c -- DC term set equal to zero arbitrarily.
      read(10,803) k,w,fr,fi,amp
      fr=0.
      fi=0.
      amp=0.
      write(11,803) k,w,fr,fi,amp
      do 404 i=2,ncnq
      read(10,803) k,w,fr,fi,amp
      den=w*w*c3
      tr=(alpha*(abs(w)*c1))/den
      ti=(-alpha*w*c2)/den
      frn=fr*tr-fi*ti
      fin=fr*ti+tr*fi
      fr=frn
      fi=fin
      amp=sqrt(fr*fr+fi*fi)
      write(11,803) k,w,fr,fi,amp
404   continue
      go to 300
295    print*,'Input file not found or not ascii.  Try again'
       go to 1
299   print*,'you''ve asked for no filtering'
300   close(10)
      close(11)
      stop
      end
      subroutine signal(slope,b,iseg,r)
      dimension slope(4),b(4),r(4)
c
1     print*,'How many line segments define signal of spectrum? (max
     1 4)'
      read*,iseg
      if(iseg.le.4.and.iseg.gt.0) go to 5
      print*,'min=1,max=4'
      go to 1
5     print*,'Enter y-intercept of line segment 1'
      read*,y0
      r1=0.
      y1=y0
c  ask for line segment endpoints and calc equation of lines, i.e.,
c         lny=lnb+slope*r
c
      do 10 i=1,iseg
      print 801,i
801   format(' Enter freq. and amp. value of endpoint of segment ',i1)
      read*,r(i),y2
      if(r(i).gt.r1) go to 8
      print*,'Define segments in ascending freq. order'
      go to 5
8     slope(i)=log(y2/y1)/(r(i)-r1)
      b(i)=y2*exp(-slope(i)*r(i))
      y1=y2
      r1=r(i)
10    continue
      return
      end
