c   MFFILTER
c
c   performs matched filtering
c
c   This is the third and final program in the matched
c   filtering sequence.  It is preceeded by, and uses
c   output from, the programs mfinit and mfdesign.
c
c   input:  mfilt.6 fourier transform parameters
c           mfilt.7 filter parameters (from mfdesign)
c           mfilt.8 fourier transform of original grid
c           mfilt.40 residual spectrum
c                       (from mfinit)
c   output: mfilt.31-39 bandpass filtered grids
c           mfilt.41-49 azimuthally filtered grids
c           mfilt.51-59 downward continued grids
c
c   Jeff Phillips September 1992
c
      common /params/ zp(9), bp(9), ilam(9), phi(361), wj(181)
      common /mane/ id(14),pgm(2),nz,yo,xo,dx,dy,iw,kr,nc,nr,ir
      dimension a1(2,2048,16),work(4096)
      character coname*50,suffix*2,drive*1
      character ida*56, pgma*8, ans*1, ask_yn*1
c
      open(10,file='mfilt.6',form='unformatted',status='old',err=990)
      read(unit=10, err=990) ida, pgma, nc, nr, nz, yo, dy, xo, dx
      write(id,'(a56)') ida
      write(pgm,'(a8)') pgma
      write(unit=6, fmt=200) ida, pgma, nc, xo, dx, nr, yo, dy
  200 format(/,6h id = ,a56,7h pgm = ,a8,/,6h nc = ,i5,5x,5hxo = 
     &,g10.3,5x,5hdx = ,g10.3,/,6h nr = ,i5,5x,5hyo = ,g10.3,5x,5hdy = 
     &,g10.3,/)
      read(unit=10) alpha, beta, gamma, nx, ny
      read(10) nri,n1,n2,id2
      close(10)
c
      nnx = (nx / 2) + 1
      nny = (ny / 2) + 1
      kmax = max0(nnx,nny)
      open(7,file='mfilt.7',form='formatted',status='old',err=991)
      do 32 i = 1, 9
      read(unit=7, fmt=*, err=33, end=33) zp(i), bp(i), ilam(i)
   32 continue
   33 close(unit=7) 
      kk = i - 1
      print'(a\)',' Enter (RAM) drive letter for temporary files: '
      read(5,199) drive
  199 format(a)
      write(unit=6, fmt=201) kk
c
  201 format(35h Filtering can be performed within ,i2,6h bands)
      idev = 30
      idev2 = 40
      idev3 = 50
      l = kk + 1
      n22=2*nx
      open(10,file='mfilt.8',form='unformatted',access='direct',
     1 status='old',recl=n22*4)
      open(12,file=drive//':\mfilt.tmp',form='unformatted',access=
     1'direct',status='unknown',recl=n22*4)
      open(15,file=drive//':\mfilt.tm2',form='unformatted',access=
     1'direct',status='unknown',recl=n22*4)

    1 idev = idev + 1
      idev2 = idev2 + 1
      idev3 = idev3 + 1
      l = l - 1
      if (l .lt. 1) goto 99
      write(unit=6, fmt=202) l, zp(l), bp(l), ilam(l)
  202 format(/,' Parameters for band ',i2,1h:,/,'       depth = ',f8.3,
     &'  amplitude = ',f8.3,'  layer code = ',i3)
      ans = ask_yn(' Do you want to filter in this band? : ')
      if (ans .eq. 'n') goto 26
c apply bandpass filter
      call bpass(a1,work,nri,n1,n2,id2,l,kk,drive)
      ans = ask_yn(' Do you want to output the bandpass filtered grid?:
     &')
c
c  inverse transform
c
      if (ans .eq. 'n') goto 62
      write(unit=6, fmt=204) idev
  204 format(' Writing bandpass filtered grid to mfilt.',i2)
      write(suffix(1:2),'(i2)') idev
      coname='mfilt.'//suffix
          open(11,status='unknown',form='unformatted',file=coname)
          write(11)id,pgm,nc,nr,nz,yo,dy,xo,dx
      call inverse(a1,work,nri,n1,n2,id2,alpha,beta,gamma,drive)
      close(11)
   62 ans = ask_yn(
     &' Do you want to do azimuthal filtering on this bandpass?: ')
c
      if (ans .eq. 'n') goto 25
      iplotr = ask_int('Enter plot device code (5=HP,8=CGA,9=EGA,10=VGA)
     &: ')
      call azimuth(a1,nri,n1,n2,iplotr,drive)
      if(iplotr.eq.-1) go to 26
c
c  inverse transform
c
      write(unit=6, fmt=206) idev2
  206 format(' Writing azimuthally filtered grid for this bandpass to
     1 mfilt.',i2)
      write(suffix(1:2),'(i2)') idev2
      coname='mfilt.'//suffix
          open(11,status='unknown',form='unformatted',file=coname)
          write(11)id,pgm,nc,nr,nz,yo,dy,xo,dx
      call inverse(a1,work,nri,n1,n2,id2,alpha,beta,gamma,drive)
      close(11)
   25 ans = ask_yn(' Do you want to downward continue this band?: ')
      if (ans .eq. 'n') goto 26
      write(unit=6, fmt=207) zp(l)
  207 format(' Maximum source depth in this band is ',f8.3)
      czp = ask_float(' Enter continuation depth: ')
      call downward(a1,nri,n1,n2,czp)
      write(unit=6, fmt=208) idev3
  208 format(' Writing downward continued grid for this bandpass to
     1 mfilt.',i2)
      write(suffix(1:2),'(i2)') idev3
      coname='mfilt.'//suffix
          open(11,status='unknown',form='unformatted',file=coname)
          write(11)id,pgm,nc,nr,nz,yo,dy,xo,dx
      call inverse(a1,work,nri,n1,n2,id2,alpha,beta,gamma,drive)
      close(11)
26    alpha = 0.
      beta = 0.
      gamma = 0.
      goto 1
  990 write(unit=*, fmt=*) 
     &'File not found: run mfinit and mfdesign before mffilter.'
      return 
  991 write(unit=*, fmt=*)
     &'File not found: run mfdesign before mffilter.'
   99 continue
      close(10)
      close(12,status='delete')
      close(15,status='delete')
      end
c
      subroutine bell()
      character ding*1
      ding = char(7)
c      write(unit=*, fmt=*) ding
      print *,ding
      return 
      end
c
      subroutine scold(*)
      call bell
      write(unit=6, fmt=100) 
  100 format(/,38h *** Unacceptable response...Try again)
      return 1
      end
c
      subroutine bpass(a1,work[huge],nri,n1,n2,id2,l,kk,drive)
          dimension a1(2,n2,nri),work(id2)
          character drive*1
      common /weights/ wght(2048)
      common /mane/ id(14),pgm(2),nz,yo,xo,dx,dy,iw,kr,ny,nx,ir
      common /params/ zp(9), bp(9), ilam(9), phi(361), wj(181)
          data iunita/10/,iunitb/12/,dval/0.17014e39/
          data iunitc/15/, iunitx/16/, iunitd/17/
          data pi2/6.2831853/

          n22=2*n2
          rn1=1./(float(n1)*dx)
          rn2=1./(float(n2)*dy)

 1000     n1p1=n1+1
          n2p1=n2+1
          nnh1=(n1/2)+1
          nnh2=(n2/2)+1
          dr=amin1(rn1,rn2)
          kmax=max0(nnh1,nnh2)
          if(kmax.gt.1025) stop 'maximum array dimension > 2048'
          open(iunitx,file=drive//':\mfilt.wgt',status='unknown',
     1         form='unformatted')
c  loop to operate on nnh1 rows of the fourier transform a1
          do 1010 jr=1,nnh1
          jj=jr-1
          x=float(jj)*rn1
          xsq=x*x
c  read a row of the fourier transform a1
          call srdda(iunita,jr,a1,n22)
c  loop to operate on the n2 elements in the row
          do 1090 i=1,n2
          ii=i-1
          if(i.gt.nnh2)ii=-(n2p1-i)
          y=float(ii)*rn2
          xy=sqrt(xsq+y*y)
      wght(i)=1
      if (ilam(l) - 1) 81, 82, 83
   81 bs = bp(l) * xy
      goto 84
   82 bs = bp(l)
      goto 84
   83 if (xy .ne. 0.0) bs = bp(l) / xy
   84 continue
c loop to construct filter weight for this element
      do 53 il = 1, kk
      if (il .eq. l) goto 53
      if (xy .eq. 0.0) then
        if (ilam(l) - 1) 91, 92, 93
   91   if (ilam(il) - 1) 94, 100, 100
   92   if (ilam(il) - 1) 53, 102, 100
   93   if (ilam(il) - 1) 53, 53, 94
   94   br = bp(il) / bp(l)
        goto 104
  100   a1(1,i,1) = 0.0
        a1(2,i,1) = 0.0
        goto 1090
      end if
      if (ilam(il) - 1) 101, 102, 103
  101 br = (bp(il) * xy) / bs
      goto 104
  102 br = bp(il) / bs
      goto 104
  103 br = (bp(il) / xy) / bs
  104 arg = (zp(l) - zp(il)) * xy * pi2
      if ((arg .gt. 88.028) .or. ((alog(br) + arg) .gt. 88.028)) then
      arg = 88.028
      br = amin1(br,1 / br)
      end if
      arg = br * exp(arg)
      if (arg .gt. (dval - wght(i))) goto 53
      wght(i) = wght(i) + arg
   53 continue
c  end of loop to construct filter element
      a1(1,i,1) = a1(1,i,1) / wght(i)
      a1(2,i,1) = a1(2,i,1) / wght(i)
      wght(i)=1./wght(i)
 1090     continue
c  end of loop on on elements of row jr
c  save the bandpass filtered fourier transform on unitb and unitc
          call swrda(iunitb,jr,a1,n22)
          call swrda(iunitc,jr,a1,n22)
          call swrbin(iunitx,wght,n2)
          if(jr.eq.1) go to 1010
          jmr=n1-jr+2
          if(jr.eq.jmr) go to 1110
          a1(2,1,1)=-a1(2,1,1)
          iwr=n2p1
          do 1100 i=2,nnh2
          iwr=iwr-1
          x=a1(2,i,1)
          y=a1(1,i,1)
          a1(2,i,1)=-a1(2,iwr,1)
          a1(1,i,1)=a1(1,iwr,1)
          a1(1,iwr,1)=y
          a1(2,iwr,1)=-x
 1100     continue
          call swrda(iunitb,jmr,a1,n22)
 1010     continue
c  iunitb now holds complex coefficients operated on by filter
c  compute desired anomaly = inv. f.f.t. of a1
 1110     isign=1
          close(iunitx)
          return
c
      entry inverse(a1,work[huge],nri,n1,n2,id2,alpha,beta,gamma,drive)
      isign=1
      open(17,file=drive//':\mfilt.tm1',form='unformatted',access=
     1'direct',status='unknown',recl=n22*4)
      call sfftmg(iunitb,iunitd,n2,n1,a1,nri,isign,work)
          area=1./float(n1*n2)
          dum=0.0
          iwr=0
      n22=2*n2
          do 1210 jr=1,nx
          call srdda(iunitb,jr,a1,n22)
 1170     do 1180 i=1,ny
 1180     a1(1,i,1)=a1(1,i,1)*area+alpha*float(i-1)+
     1              beta*float(jr-1)+gamma
 1190     write(11)dum,(a1(1,i,1),i=1,ny)
 1210     continue
 1220     format(/,1x,'row ',i3)
 1230     format(3x,16f7.1)
 1240     format(1h1,14a4,//' calculated anomaly on level z=',f7.3,//)
          close(17,status='delete')
 1290     return
          end
c
      subroutine azimuth(a1,nri,n1,n2,iplotr,drive)
          dimension a1(2,n2,nri)
          dimension hc(2048)
          character drive*1
      common /weights/ wght(2048)
      dimension iramp(4, 50), phi1(361), isig(2, 50)
      common /mane/ id(14),pgm(2),nz,yo,xo,dx,dy,iw,kr,ny,nx,ir
      common /params/ zp(9), bp(9), ilam(9), phi(361), wj(181)
      dimension xxx(181),yyy(181),ind(181)
      character ans*1, ask_yn*1
          data iunitb/12/
          data iunitc/15/,iunitx/16/
c
      do 63 i = 1, 180
      phi1(i) = 0.0
   63 wj(i) = .0
c
          n22=2*n2
          rn1=1./(float(n1)*dx)
          rn2=1./(float(n2)*dy)
          n2p1=n2+1
          nnh1=(n1/2)+1
          nnh2=(n2/2)+1
          dr=amin1(rn1,rn2)
          kmax=max0(nnh1,nnh2)
c
c compute the phase spectrum within the bandpass
      open(iunitx,file=drive//':\mfilt.wgt',status='old',form=
     1'unformatted')
      open(17,file='mfilt.40',status='old',form='unformatted')
      read(17) id
c  loop to operate on nnh1 rows of the log residual spectrum
          do 64 jr=1,nnh1-2
   64     read(17) dum
          do 3010 jr=1,nnh1
          jj=jr-1
          x=float(jj)*rn1
          xsq=x*x
c  read a row of the log residual spectrum and of the bandpass weights
          call srdbin(17,hc,n2)
          call srdbin(iunitx,wght,n2)
c  loop to operate on the n2 elements in the row
          do 3090 i=1,n2
          ii=i-nnh2
          y=float(ii)*rn2
          xy=sqrt(xsq+y*y)
          ihck = ifix((xy/dr)+1.5)
          if(ihck.gt.kmax) go to 3090
          if(xy.eq.0.0) go to 3090
          kp=ifix(atan2(x,y) * 57.29577951 + 0.5)
          if(kp.eq.0) go to 3090
          j=i+nnh2
          if(j.gt.n2) j=i-nnh2+2
          phi1(kp) = phi1(kp) + hc(i) * wght(j)
          wj(kp) = wj(kp) + wght(j)
 3090     continue
 3010     continue
          close(iunitx)
          close(17)
      do 67 i = 1, 180
      if (wj(i) .le. 0.0) goto 67
      phi1(i) = phi1(i) / wj(i)
   67 continue
   68 n = 0
      cycles = ask_float(' Enter number of cosine cycles (1 or 2): ')
      call cosine(off,amp,pha,cycles,phi1)
      write(unit=6, fmt=210) off, amp, pha
  210 format(21h Cosine term:  off = ,f10.3,8h  amp = ,f10.3,8h  pha =
     &,f10.3)
        ylim1 = off-5.*amp
        ylim2 = off+5.*amp
      ans = ask_yn(' Plot azimuthal spectrum and cosine term?: ')
      if (ans .eq. 'y') then
        ans = ask_yn(' Use range -90 to 90?: ')
        if (ans.eq.'y') then
          call rplot(-90., 90., ylim1, ylim2, 3, iplotr)
          xxx(1) = -90
          xxx(2) = 90.
        else
          call rplot(0., 180., ylim1, ylim2, 3, iplotr)
          xxx(1) = 0.
          xxx(2) = 180.
        endif
        yyy(1) = 0.
        yyy(2) = 0.
        call curv(xxx, yyy, 2, 0)
        do i = 2, 181
        if (ans.eq.'y') then
          xxx(i) = float(i) - 91.
          if(i.lt.92) then
            wj(i) = phi1(89 + i)
          else
            wj(i) = phi1(i - 91)
          endif
        else
          xxx(i) = float(i) - 1.
          wj(i) = phi1(i - 1)
        endif
        yyy(i) = off + (amp * cos((pha + (2. * cycles * xxx(i))) /
     &           57.29577951))
        end do
        if (ans.eq.'y') then
          wj(1) = phi1(90)
          xxx(1) = -90.
        else
          wj(1) = phi1(180)
          xxx(1) = 0.
        endif
        yyy(1) = off + (amp * cos((pha + (2. * cycles * xxx(1))) /
     &           57.29577951))
        call curv(xxx, wj, 181, 300)
        call curv(xxx, yyy, 181, 200)
        call endit
        ans = ask_yn(' Change number of cycles? : ')
        if(ans.eq.'y') go to 68
      end if
c determine ramps for significant azimuths
      n = 0
      sumx = 0.
      sumx2 = 0.
      pmin = 0.
      do 10 i = 1, 180
      if (phi1(i) .eq. 0.0) goto 10
c remove cosine term
      n = n + 1
      xxx(n) = phi1(i) - (amp * cos((pha + (2. * cycles * float(i))) /
     &57.29577951))
      ind(n) = i
      sumx = sumx + xxx(n)
      sumx2 = sumx2 + (xxx(n) * xxx(n))
      if (phi1(i) .lt. pmin) pmin = phi1(i)
   10 continue
      am = sumx / float(n)
      sd = (sumx2 - ((n * am) * am)) / float(n - 1)
      sd = sqrt(sd)
c 67% band
      al = am - sd
      au = am + sd
c 95% band
      all = am - (2. * sd)
      auu = am + (2. * sd)
      write(unit=6, fmt=203) am, sd
  203 format(30h Residual statistics:  mean = ,f10.3,9h  s.d. = ,f10.3)
      ir = 1
      jrr = 1
      isig(1,1) = 0
      isig(2,1) = 0
      iramp(1,1) = -1
c
c if x > au   start a ramp up
c if x > auu  top of ramp up
c if x < auu  start a ramp down
c if x > al   bottom of ramp down
c
      low1 = 0
      open(2,file='temp.out',form='formatted',status='unknown')
      do 30 i = 1, n
      write(2,*) ind(i), xxx(i)
      if (xxx(i) .gt. au) then
        goto (70, 71), ir
   70   iramp(2,jrr) = ind(i)
        iramp(3,jrr) = ind(i)
        iramp(4,jrr) = -1
        ir = 2
        goto 72
   71   iramp(3,jrr) = ind(i)
   72   continue
        if (xxx(i) .gt. auu) then
          isig(2,jrr) = isig(2,jrr) + 1
        else
          isig(1,jrr) = isig(1,jrr) + 1
        end if
        write(2,*) iramp(1,jrr), iramp(2,jrr), iramp(3,jrr),
     &    iramp(4,jrr)
      else if (xxx(i) .lt. al) then
        low2 = ind(i)
        if (low1 .eq. 0) low1 = low2
        goto (75, 76), ir
   75   iramp(1,jrr) = ind(i)
        goto 77
   76   iramp(4,jrr) = ind(i)
        ir = 1
        jrr = jrr + 1
        isig(1,jrr) = 0
        isig(2,jrr) = 0
        iramp(1,jrr) = ind(i)
        iramp(4,jrr) = 0
   77   continue
        write(2,*) iramp(1,jrr), iramp(2,jrr), iramp(3,jrr),
     &    iramp(4,jrr)
      end if
c handle special cases
   30 continue
      if (iramp(4,jrr) .eq. (-1)) then
        if (iramp(1,1) .eq. (-1)) then
          iramp(4,jrr) = 180
          iramp(3,jrr) = 180
          iramp(1,1) = 1
          iramp(2,1) = 1
        else
          iramp(4,jrr) = low1 + 180
        end if
      else
        if (iramp(1,1) .eq. (-1)) iramp(1,1) = low2 - 180
      end if
c
      if (iramp(4,jrr) .eq. 0) jrr = jrr - 1
      ans = ask_yn(' Do you want to remove the cosine term? : ')
      do 35 i = 1, 180
      phi(i) = amp * cos((pha + (2. * cycles * float(i))) /
     &57.29577951)
      if (ans .eq. 'n') phi(i) = 0.0
   35 continue
      if (ans .eq. 'n') amp = 0.0
      ans = ask_yn(
     &' Do you want to remove power within specified azimuth ranges? :
     & ')
      if (ans .eq. 'n') then
        if (amp .ne. 0.0) then
          write(unit=*, fmt=*) 'Removing only cosine term'
          goto 22
        else
          write(unit=*, fmt=*) 'No azimuthal filtering done'
          iplotr=-1
          return
        end if
c
      end if
      write(unit=*, fmt=*) 'Suggested ranges:'
      write(unit=*, fmt=*)
     &' Index    Ramp_on    Top   Ramp_off      67%  95%'
      write(unit=*, fmt=*) '   0        Bypass automatic ranges'
      do 79 i = 1, jrr
   79 write(unit=6, fmt=209) i, (iramp(j,i),j = 1, 4), (isig(j,i),j = 1
     &, 2)
  209 format(i5,4i8,i6,i5)
   23 ir = ask_int('Enter a range index (-1 to exit): ')
      if ((ir .lt. (-1)) .or. (ir .gt. jrr)) call scold(*23)
      if (ir .lt. 0) goto 22
      jrrp=jrr+1
      if(ir.eq.0) ir=jrrp
   24 sd = ask_float(
     &'Enter multiplier (0.25=slight, 0.5=average, 1=severe  filtering):
     & ')
      if (ir .ne. jrrp) goto 21
      print*,'Use azimuths in the range 1 to 180'
  125 iramp(1,jrrp) = ask_int('Enter azimuth for bottom of left ramp: ')
      if ((iramp(1,jrrp) .lt. 1) .or. (iramp(1,jrrp) .gt. 180)) call
     1scold(*125)
  126 iramp(2,jrrp) = ask_int('Enter azimuth for top of left ramp: ')
      if ((iramp(2,jrrp) .lt. 1) .or. (iramp(2,jrrp) .gt. 180)) call
     1scold(*126)
  127 iramp(3,jrrp) = ask_int('Enter azimuth for top of right ramp: ')
      if ((iramp(3,jrrp) .lt. 1) .or. (iramp(3,jrrp) .gt. 180)) call
     1scold(*127)
  128 iramp(4,jrrp) = ask_int('Enter azimuth for bottom of right ramp: '
     1)
      if ((iramp(4,jrrp) .lt. 1) .or. (iramp(4,jrrp) .gt. 180)) call
     1scold(*128)
   21 continue
      denom = iramp(2,ir) - iramp(1,ir)
      if (denom .gt. 1.) then
        denom = denom / sd
        do i = iramp(1,ir) + 1, iramp(2,ir) - 1
        phi(i) = phi(i) + (((phi1(i)-pmin) * (i-iramp(1,ir))) / denom)
        end do
      end if
      do i = iramp(2,ir), iramp(3,ir)
      phi(i) = phi(i) + (sd * (phi1(i) - pmin))
      end do
      denom = iramp(4,ir) - iramp(3,ir)
      if (denom .gt. 1.) then
        denom = denom / sd
        do i = iramp(3,ir) + 1, iramp(4,ir) - 1
        phi(i) = phi(i) + (((phi1(i)-pmin) * (iramp(4,ir)-i)) / denom)
        end do
      end if
      if(ir.ne.jrrp) go to 23
      ans = ask_yn(' Do you want to specify additional ranges? : ')
      if (ans .eq. 'y') goto 24
c write attenuation factors to for092
c
   22 continue
      do 20 i = 1, 180
c phi contains the attenuation factors
   20 phi(i + 180) = phi(i)
c  end of loop on on elements of row jr
c  loop to operate on nnh1 rows of the fourier transform a1
          do 1010 jr=1,nnh1
          jj=jr-1
          x=float(jj)*rn1
          xsq=x*x
c  read a row of the fourier transform a1
          call srdda(iunitc,jr,a1,n22)
c  loop to operate on the n2 elements in the row
          do 1090 i=1,n2
          ii=i-1
          if(i.gt.nnh2)ii=-(n2p1-i)
          y=float(ii)*rn2
          xy=sqrt(xsq+y*y)
          if(xy.eq.0.0) go to 1090
          kp=atan2(x,y) * 57.29577951 + 180.5
          arg=phi(kp)
          if(arg.gt.88.028) arg=88.028
          arg=exp(-arg)
c  end of loop to construct filter element
      a1(1,i,1) = a1(1,i,1) * arg
      a1(2,i,1) = a1(2,i,1) * arg
 1090     continue
c  end of loop on on elements of row jr
c  save the bandpass filtered fourier transform on unitb and unitc
          call swrda(iunitb,jr,a1,n22)
          if(jr.eq.1) go to 1010
          jmr=n1-jr+2
          if(jr.eq.jmr) go to 1110
          a1(2,1,1)=-a1(2,1,1)
          iwr=n2p1
          do 1100 i=2,nnh2
          iwr=iwr-1
          x=a1(2,i,1)
          y=a1(1,i,1)
          a1(2,i,1)=-a1(2,iwr,1)
          a1(1,i,1)=a1(1,iwr,1)
          a1(1,iwr,1)=y
          a1(2,iwr,1)=-x
 1100     continue
          call swrda(iunitb,jmr,a1,n22)
 1010     continue
c  iunitb now holds complex coefficients operated on by filter
c  compute desired anomaly = inv. f.f.t. of a1
 1110     isign=1
      do 85 i = 1, 361
   85 phi(i) = phi1(i) - phi(i)
      ans = ask_yn(
     &' Plot old (dashed) and new (solid) azimuthal spectrum?: ')
      if (ans .eq. 'y') then
        ans = ask_yn(' Use range -90 to 90?: ')
        if (ans.eq.'y') then
          call rplot(-90., 90., ylim1, ylim2, 3, iplotr)
          xxx(1) = -90
          xxx(2) = 90.
        else
          call rplot(0., 180., ylim1, ylim2, 3, iplotr)
          xxx(1) = 0.
          xxx(2) = 180.
        endif
        yyy(1) = 0.
        yyy(2) = 0.
        call curv(xxx, yyy, 2, 0)
        do i = 2, 181
        if (ans.eq.'y') then
          xxx(i) = float(i) - 91.
          if(i.lt.92) then
            wj(i) = phi1(89 + i)
            yyy(i) = phi(89 + i)
          else
            wj(i) = phi1(i - 91)
            yyy(i) = phi(i - 91)
          endif
        else
          xxx(i) = float(i) - 1.
          wj(i) = phi1(i - 1)
          yyy(i) = phi(i - 1)
        endif
        end do
        if (ans.eq.'y') then
          wj(1) = phi1(90)
          yyy(1) = phi(90)
          xxx(1) = -90.
        else
          wj(1) = phi1(180)
          yyy(1) = phi(180)
          xxx(1) = 0.
        endif
        call curv(xxx, wj, 181, 201)
        call curv(xxx, yyy, 181, 300)
        call endit
        ans = ask_yn(' Is this filter acceptable?: ')
        if (ans .eq. 'n') goto 68
      endif
           return
           end
c
c***********************************************************************
      subroutine srdbin(no,dat[huge],n)
c       subroutine 'srdbin' reads standard grids (binary).
c***********************************************************************
      dimension dat(n)
      read(no)dum,dat
      return
      end
c***********************************************************************
      subroutine swrbin(no,dat[huge],n)
c       subroutine 'swrbin' writes standard grids (binary).
      dimension dat(n)
      data dum/0.0/
      write(no)dum,dat
      return
      end
c***********************************************************************
       subroutine srdda(no,ipos,dat[huge],n)
c       subroutine 'srrda' reads keyed sequential files.
c***********************************************************************
      dimension dat(n)
      read(no,rec=ipos)dat
      return
      end
c***********************************************************************
      subroutine swrda(no,ipos,dat[huge],n)
c       subroutine 'swrda' writes keyed sequential files.
c***********************************************************************
      dimension dat(n)
      write(no,rec=ipos)dat
      return
      end
c***********************************************************************
      function ask_yn(request)
      character ask_yn*1
      character request*(*)
    2 write(unit=6, fmt=100) request
  100 format(a,1h $)
      read(unit=5, fmt=102, err=1, end=1) ask_yn
  102 format(a)
      if (ask_yn .eq. 'Y') ask_yn = 'y'
      if (ask_yn .eq. 'N') ask_yn = 'n'
      if ((ask_yn .ne. 'y') .and. (ask_yn .ne. 'n')) call scold(*2)
      return 
    1 call scold(*2)
      end
      function ask_float(request)
      character request*(*)
    2 write(unit=6, fmt=100) request
  100 format(a,$)
      read(unit=5, fmt=*, err=1) ask_float
      return 
    1 write(unit=6, fmt=101) 
  101 format(/,20h **ERROR - try again)
      call bell
      goto 2
c-----------------------------------------------------------------------
c-----
      end
      function ask_int(request)
      character request*(*)
    2 write(unit=6, fmt=100) request
  100 format(1x,a,$)
      read(unit=5, fmt=*, err=1) ask_int
      return 
    1 write(unit=6, fmt=101) 
  101 format(/,20h **ERROR - try again)
      call bell
      goto 2
c-----------------------------------------------------------------------
c---
      end
      subroutine rplot(rmin, rmax, hmin, hmax, ityp, iplotr)
      dimension xp(4), yp(4), dxp(2), dyp(2)
      dimension x(n), y(n)
c      integer str(n)
      call pltset(iplotr, xp(4), yp(4), 1)
      xp(4) = min(xp(4),10.)
      yp(4) = min(yp(4),8.)
      xp(4) = xp(4) - .25
      yp(4) = yp(4) - .25
      dxp(1) = rmin
      dxp(2) = (rmax + .01) - mod(rmax,.01)
      dyp(1) = hmin
      dyp(2) = hmax
      xp(1) = xp(4) - 1.5
      xp(2) = 0.
      xp(3) = 1.25
      yp(1) = yp(4) - 1.25
      yp(2) = 0.
      yp(3) = 1.
      call scale(dxp, dyp, xp, yp, 4, ier)
      if (ier .ne. 0) goto 90
      goto (51, 52, 53), ityp
   51 call yaxis(dyp, dxp, yp, 0.1, 2, .12, '(f4.1)', 4)
      goto 54
   52 call yaxis(dyp, dxp, yp, 1.0, 2, .12, '(f4.0)', 4)
      goto 54
   53 call yaxis(dyp, dxp, yp, 0.5, 2, .12, '(f4.1)', 4)
   54 goto (55, 55, 56), ityp
   55 call xaxis(dxp, dyp, xp, 0.1, 2, .12, '(f4.1)', 4)
      goto 57
   56 call xaxis(dxp, dyp, xp, 5.0, 2, .12, '(f4.0)', 4)
   57 call neatl
      goto (58, 58, 59), ityp
   58 call vchar((dxp(1) + dxp(2)) * .5, dyp(1), 'radians/km', 10, 2, 
     &.12, 0., -.54, -.54)
      goto 60
   59 call vchar((dxp(1) + dxp(2)) * .5, dyp(1), 'azimuth', 7, 2, .12, 
     &0., -.54, -.54)
   60 goto (61, 62, 62), ityp
   61 call vchar(dxp, (dyp(1) + dyp(2)) * .5, 'amplitude', 9, 2, .12, 
     &1.5706, -.48, .54)
      goto 64
   62 call vchar(dxp, (dyp(1) + dyp(2)) * .5, 'log power', 9, 2, .12, 
     &1.5706, -.48, .54)
   64 return 
   90 stop 
      entry curv(x, y, n, m)
      call line(x, y, n, 0, m)
      return 
      entry endit()
      call endpt(ie)
      return 
      end
      subroutine downward(a1,nri,n1,n2,czp)
          dimension a1(2,n2,nri)
      common/mane/id(14),pgm(2),nz,yo,xo,dx,dy,iw,kr,ny,nx,ir
          data iunitb/12/,iunitc/15/
          data pi2/6.2831853/

          n22=2*n2
          rn1=1./(float(n1)*dx)
          rn2=1./(float(n2)*dy)

 1000     n1p1=n1+1
          n2p1=n2+1
          nnh1=(n1/2)+1
          nnh2=(n2/2)+1
          dr=amin1(rn1,rn2)
          kmax=max0(nnh1,nnh2)
          if(kmax.gt.1025) stop 'maximum array dimension > 2048'
c  loop to operate on nnh1 rows of the fourier transform a1
          do 1010 jr=1,nnh1
          jj=jr-1
          x=float(jj)*rn1
          xsq=x*x
c  read a row of the fourier transform a1
          call srdda(iunitc,jr,a1,n22)
c  loop to operate on the n2 elements in the row
          do 1090 i=1,n2
          ii=i-1
          if(i.gt.nnh2)ii=-(n2p1-i)
          y=float(ii)*rn2
          xy=sqrt(xsq+y*y)
      arg = czp * xy * pi2
      if (arg .gt. 88.028) arg = 88.028
      arg=exp(arg)
      a1(1,i,1) = a1(1,i,1) * arg
      a1(2,i,1) = a1(2,i,1) * arg
 1090     continue
c  end of loop on on elements of row jr
c  save the bandpass filtered fourier transform on unitb and unitc
          call swrda(iunitb,jr,a1,n22)
          if(jr.eq.1) go to 1010
          jmr=n1-jr+2
          if(jr.eq.jmr) go to 1110
          a1(2,1,1)=-a1(2,1,1)
          iwr=n2p1
          do 1100 i=2,nnh2
          iwr=iwr-1
          x=a1(2,i,1)
          y=a1(1,i,1)
          a1(2,i,1)=-a1(2,iwr,1)
          a1(1,i,1)=a1(1,iwr,1)
          a1(1,iwr,1)=y
          a1(2,iwr,1)=-x
 1100     continue
          call swrda(iunitb,jmr,a1,n22)
 1010     continue
c  iunitb now holds complex coefficients operated on by filter
c  compute desired anomaly = inv. f.f.t. of a1
 1110     isign=1
          return
          end
c
      subroutine cosine(off,amp,pha,cycles,phi)
c
c     returns the offset, amplitude, and phase of the specified number
c     of cycles of a cosine wave fitting the curve phi
c
c     Jeff Phillips 1992
c
      dimension phi(180)
      n=0
      pha=0.
      si2=0.
      sj2=0.
      sih=0.
      sjh=0.
      do 10 i=1,180
      a1 = cos((2.* cycles * float(i)) / 57.29577951)
      a2 = - sin((2. * cycles * float(i)) / 57.29577951)
      off=off+phi(i)
      si2=si2+a1*a1
      sj2=sj2+a2*a2
      sih=sih+a1*phi(i)
      sjh=sjh+a2*phi(i)
   10 continue
      off=off/180.
      amp=si2*sj2
      a=(sj2*sih)/amp
      b=(si2*sjh)/amp
      amp = sqrt(a**2 + b** 2)
      pha = atan2(b,a) * 57.29577951
      return
      end

