c   MFINIT
c
c   computes the radial power spectrum and the residual
c   power spectrum of a standard gridded file.  The
c   input grid must be complete (no dvals), and its
c   column dimension must be less than 2049.
c
c   This is the first program in the matched filtering
c   sequence.  It is followed by mfdesign, and mffilter.
c
c   input:  a specified standard grid
c   output: mfilt.6 fourier transform parameters
c           mfilt.8 fourier transform
c           mfilt.9 log averaged radial power spectrum
c           mfilt.40 log residual power spectrum
c           mfilt.50 extended input grid
c
c   link with SFFTMG and SFOURT
c
c   Jeff Phillips June 1992
c
      dimension a1(2,2048,16),work(4096)
      common/mane/id(14),pgm(2),nz,yo,xo,dx,dy,iw,kr,ny,ir
      data kr/11/,ir/5/,iw/6/
      call prep(a1,work,work(2049))
      print *,'Extended input grid written to mfilt.50'
c  read header record of extended input file.
  180 open(kr,status='old',form='unformatted',file='mfilt.50')
c  note that x & y have been switched from normal (usgs) grid specficati
c  ...convention used is x-north y-east and z-down
      read(kr) id,pgm,n2,n1,nz,yo,dy,xo,dx
      pgm(1)='mfin'
      pgm(2)='it  '
      call get_dims(n1,n2,nx,ny,nri,id2,nxa)
      write(iw,275) nri
  275 format(' blocking factor for rows =',i3)
      write(14) nri,n1,n2,id2
      close(14)
c  input parameters for fft are set: call main program.
      call mfftfil(a1,work,nri,n1,n2,id2)
      stop
      end
      subroutine prep(g,g1,g2)
c
c  Open grid, tilt, and cosine taper to reduce FFT wraparound.
c  The new nrow, ncol should be chosen from the list of rich-in-
c  -factors-of-2 numbers that SFFTMG likes.
c
c  Lin Cordell Jan 1990.
c  Revised June 1991.
c  Modified by Jeff Phillips Feb 1992.
c
	character*50 ifile,ofile
      dimension g(1),g1(1),g2(1),id(14),pgm(2)
      index(i,j,nc) = ((j - 1) * nc) + i
        pi=3.1415927
      print*, 'Enter input grid name (NO DVALS !):'
	print 101
101	format(' *',$)
	read(5,1)ifile
1	format(a50)
      open(10,file=ifile,status='old',form='unformatted')
        read(10)id,pgm,nc,nr,nz,xo,dx,yo,dy
        call fftdims(nc,nr)
	print*, 'Enter (new) ncol and nrow:'
	print 101
      read*, ncol,nrow
c  Check to make sure SFFTMG will accept these.
c  Added 19 June 91.
      kr=11
      iw=6
      ir=5
c      nadd=0
c  note that x & y have been switched from normal (usgs) grid specficati
c  ...convention used is x-north y-east and z-down
      n1=nrow
      n2=ncol
      call get_dims(n1,n2,nx,ny,nri,id2,nxa)
      itest=n1-nx
      if(itest.ne.0) then
        print*,'Select another nrow; SFFTMG cannot handle this one.'
        stop
      endif
      if(nc.gt.ncol.or.nr.gt.nrow)then
        print*,'New ncol, nrow are too small.'
        stop
      endif
c  Tilt grid.
      call edge(alpha,beta,gamma,nc,nr,10)
      close(10)
c        print*, ' alpha = ',alpha,', beta = ',beta,' ,gamma = ',gamma
	open(10,file=ifile,status='old',form='unformatted')
	read(10)id,pgm,nc,nr,nz,xo,dx,yo,dy
c  Output extended file.
        ofile='mfilt.50'
      open(11,file=ofile,status='unknown',form='unformatted')
	pgm(1)='prep'
	pgm(2)='    '
	write(11)id,pgm,ncol,nrow,nz,xo,dx,yo,dy
      open(14,file='mfilt.6',status='unknown',form='unformatted')
        write(14)id,pgm,nc,nr,nz,xo,dx,yo,dy
        write(14)alpha,beta,gamma,ncol,nrow
c  First row (saved).
        j2=1
        call rowio(g,yr,nz,nc,10,1)
        y=float(j2)
      do 211 j1=1,nc
        x=float(j1)
        g(j1)=g(j1)-alpha*x-beta*y-gamma
211   g1(j1)=g(j1)
c  Cosine taper right side
        s=0.5*(g(nc)+g(1))
	r=0.5*(g(nc)-g(1))
	lim=ncol-nc
	rlamb=float(lim+1)
	arg=pi/rlamb
	do 201 k=1,lim
	rk=float(k)
	ge=s+r*cos(arg*rk)
	g(nc+k)=ge
201     g1(nc+k)=ge
        call rowio(g,yr,nz,ncol,11,2)
c  Rows down to nr
        do 204 j2=2,nr
        y=float(j2)
	call rowio(g,yr,nz,nc,10,1)
	do 202 j1=1,nc
        x=float(j1)
202     g(j1)=g(j1)-alpha*x-beta*y-gamma
c  Cosine taper right side.
        s=0.5*(g(nc)+g(1))
	r=0.5*(g(nc)-g(1))
	do 203 k=1,lim
	rk=float(k)
203     g(nc+k)=s+r*cos(arg*rk)
        call rowio(g,yr,nz,ncol,11,2)
204	continue
c  Last row
c  Extend right side
c  Rows nr+1 to nrow.
        do 205 j1=1,ncol
205     g2(j1)=g(j1)
        lim=nrow-nr
	rlamb=float(lim+1)
	arg=pi/rlamb
	do 207 k=1,lim
	yr=float(nr+k-1)*dx+yo
	rk=float(k)
	do 206 j1=1,ncol
	s=0.5*(g2(j1)+g1(j1))
	r=0.5*(g2(j1)-g1(j1))
206     g(j1)=s+r*cos(arg*rk)
        call rowio(g,yr,nz,ncol,11,2)
207	continue
	go to 999
900     print*, "Can't handle it."
	close(10)
999	close(11)
        return
	end
      subroutine rowio(a,y,nz,ncol,ld,key)
      dimension a(nz,ncol)
      go to (1,2),key
1     read(ld) y,a
      go to 90
2     write(ld)y,a
90    return
      end
      subroutine edge(a,b,c,nc,nr,iunit)
c
c     returns the coefficients a,b,c of a planar surface that best fits the
c     edges of an (nc by nr) standard grid:
c                         edge(i,j) = a*i + b*j + c
c     The standard grid must be opened as unit iunit and the header must be
c     read prior to calling this subroutine.
c
c     Jeff Phillips 1992
c
      dimension g(2048)
      nz=1
      n=0
      c=0.
      si2=0.
      sj2=0.
      sih=0.
      sjh=0.
      avi=float(nc+1)/2.
      avj=float(nr+1)/2.
      call rowio(g,y,nz,nc,iunit,1)
      aj=1.0-avj
      do 10 i=1,nc
      ai=float(i)-avi
      n=n+1
      c=c+g(i)
      si2=si2+ai*ai
      sj2=sj2+aj*aj
      sih=sih+ai*g(i)
      sjh=sjh+aj*g(i)
   10 continue
      ai1=float(1)-avi
      ainc=float(nc)-avi
      do 20 j=2,nr-1
      aj=float(j)-avj
      call rowio(g,y,nz,nc,iunit,1)
      n=n+2
      c=c+g(1)+g(nc)
      si2=si2+ai1*ai1+ainc*ainc
      sj2=sj2+2*aj*aj
      sih=sih+ai1*g(1)+ainc*g(nc)
      sjh=sjh+aj*(g(1)+g(nc))
   20 continue
      call rowio(g,y,nz,nc,iunit,1)
      aj=float(nr)-avj
      do 30 i=1,nc
      ai=float(i)-avi
      n=n+1
      c=c+g(i)
      si2=si2+ai*ai
      sj2=sj2+aj*aj
      sih=sih+ai*g(i)
      sjh=sjh+aj*g(i)
   30 continue
      c=c/float(n)
      d=si2*sj2
      a=(sj2*sih)/d
      b=(si2*sjh)/d
      c=c-avi*a-avj*b
      return
      end
      subroutine fftdims(ncol, nrow)
c
c     suggests new dimensions for fourier transform
c
c     Jeff Phillips 1992
c
      dimension nlist(54), ngood(54)
      data nlist / 9, 10, 11, 12, 13, 16
     &, 18, 20, 22, 24, 26, 32
     &, 36, 40, 44, 48, 52, 56, 64
     &, 72, 80, 88, 96, 104, 112, 128
     &, 144, 160, 176, 192, 208, 224, 256
     &, 288, 320, 352, 384, 416, 448, 512
     &, 576, 640, 704, 768, 832, 896, 1024
     &, 1152, 1280, 1408, 1536, 1664, 1792, 2048 /
c
      write(unit=6, fmt=901) ncol, nrow
  901 format(' ncol = ',i5,' nrow = ',i5)
      n = 0
      do 10 i = 1, 54
      if (nlist(i) .lt. ncol) goto 10
      if (nlist(i) .gt. ncol+ncol/2) goto 20
      n = n + 1
      ngood(n) = nlist(i)
   10 continue
   20 write(unit=6, fmt=902)
  902 format(' Suggested new values:')
      write(unit=6, fmt=101) (ngood(i),i = 1, n)
  101 format('  ncol = ',14i5)
      n = 0
      do 30 i = 1, 54
      if (nlist(i) .lt. nrow) goto 30
      if (nlist(i) .gt. nrow+nrow/2) goto 40
      n = n + 1
      ngood(n) = nlist(i)
   30 continue
   40   write(unit=6, fmt=102) (ngood(i),i = 1, n)
  102 format('  nrow = ',14i5)
      return 
      end
c***********************************************************************
c   subroutine 'mfftfil' checks for flagged values in grid, transforms
c   data, computes log average radial power spectrum and log residual
c   spectrum.
c***********************************************************************
          subroutine mfftfil(a1,work[huge],nri,n1,n2,id2)
          dimension a1(2,n2,nri),work(id2)
          dimension g(1025),j1(1025)
          character drive*1, tempfile*80
          common/mane/id(14),pgm(2),nz,yo,xo,dx,dy,iw,kr,ny,ir
          data iunita/10/,iunitb/12/
          data pi2/6.28318531/,iunitc/15/
c
      print*, 'Enter (RAM) drive letter for temporary files:'
        print 101
101     format(' *',$)
        read(5,1)drive
1       format(a)

          n22=2*n2
   20     do 30 i=1,2
          do 30 j=1,n2
   30     a1(2,j,i)=0.0
         tempfile=drive//':\slave1.tmp'
         open(10,access='direct',status='unknown',form='unformatted',
     1   file=tempfile,recl=n22*4)
c   read & check input data.
          do 450 iwr=1,n1
          read(kr)dum,(a1(1,i,1),i=1,n2)
          do 320 ii=1,n2
          if(a1(1,ii,1).ge.1.0e37) then
            write(iw,340)jr,ii
  340 format(' #flagged value row=',i4,3x,'col=',i4,' may be more?')
            close(11)
            close(10)
            return
          endif
  320     continue
          call swrda(iunita,iwr,a1,n22)
  450 continue
          ih=iw
          rn1=1./(float(n1)*dx)
          rn2=1./(float(n2)*dy)
  930     print *,'Computing fourier transform'
  980     isign=-1
          tempfile=drive//':\slave2.tmp'
          open(12,access='direct',status='unknown',form='unformatted',
     1    file=tempfile,recl=n22*4)
          call sfftmg(iunita,iunitb,n2,n1,a1,nri,isign,work)
          open(15,access='direct',status='unknown',form='unformatted',
     1    file='mfilt.8',recl=n22*4)
c
          do 990 jr=1,n1
          call srdda(iunita,jr,a1,n22)
  990     call swrda(iunitc,jr,a1,n22)
      print *,'Fourier transform written to mfilt.8'
 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'
          do 31 i=1,kmax
          j1(i)=0
          g(i)=0.
   31     continue
c
      iexc = ask_int('Enter radius (degrees) of angular exclusion zone a
     &bout the cardinal directions: ')
          do 1010 jr=1,nnh1
          jj=jr-1
          x=float(jj)*rn1
          xsq=x*x
          call srdda(iunita,jr,a1,n22)
          do 1090 i=1,n2
          ii=i-1
          if(i.gt.nnh2)ii=-(n2p1-i)
          y=float(ii)*rn2
          xy=sqrt(xsq+y*y)
          k=ifix((xy/dr)+1.5)
          if(xy.eq.0.0) go to 1089
          kp=ifix(atan2(x,y) * 57.29577951 + 0.5)
          if(kp.lt.iexc) go to 1090
          if(kp.gt.(90-iexc).and.kp.lt.(90+iexc)) go to 1090
          if(kp.gt.(180-iexc)) go to 1090
 1089     j1(k)=j1(k)+1
          g(k)=g(k)+a1(1,i,1)**2+a1(2,i,1)**2
 1090     continue
          if(jr.eq.1) go to 1010
 1010 continue
 1110 continue
      open(9,file='mfilt.9',status='unknown',form='formatted')
      dr=dr*pi2
      if(j1(1).eq.0) then
        gzero=0.0
      else
        gzero=alog(g(1)/float(j1(1)))
      endif
c
      do 91 k=1,kmax
      rad=(k-1)*dr
      if(j1(k).eq.0) go to 91
      g(k)=alog(g(k)/float(j1(k)))-gzero
      write(9,*) rad,g(k)
   91 continue
      print *,'Log average radial power spectrum written to mfilt.9'
      open(40,file='mfilt.40',status='unknown',form='unformatted')
      tempfile=drive//':\mfilt.tmp'
      open(41,file=tempfile,access='direct',status='unknown',
     1 form='unformatted',recl=n2*4)
      nz=1
      write(40) id,pgm,n2,n1,nz,-(n2/2)*rn2,rn2,-(n1/2-1)*rn1,rn1
          dr=amin1(rn1,rn2)
c
          do 2010 jr=1,nnh1
          jj=jr-1
          x=float(jj)*rn1
          xsq=x*x
          call srdda(iunita,jr,a1,n22)
          do 2090 i=1,n2
          ii=i-1
          if(i.gt.nnh2)ii=-(n2p1-i)
          y=float(ii)*rn2
          xy=sqrt(xsq+y*y)
          k=ifix((xy/dr)+1.5)
          j=i+nnh2-1
          if(j.gt.n2) j=i-n2/2
          if(k.gt.kmax) then
            work(j)=0.0
          else
            work(j)=alog(a1(1,i,1)**2+a1(2,i,1)**2)-gzero-g(k)
          endif
 2090     continue
          call swrda(41,jr,work,n2)
          if(jr.eq.1) go to 2010
 2010      continue
 2110 continue
      do 32 jr=1,nnh1-2
      call srdda(41,nnh1-jr,work,n2)
      do 34 i=2,nnh2-1
      gzero=work(i)
      work(i)=work(n2-i+2)
   34 work(n2-i+2)=gzero
      call swrbin(40,work,n2)
   32 continue
      do 33 jr=1,nnh1
      call srdda(41,jr,work,n2)
      call swrbin(40,work,n2)
   33 continue
      print *,'Log residual spectrum written to mfilt.40'
      close(40)
      close(41,status='delete')
 1260     close(10,status='delete')
          close(12,status='delete')
          go to 1290
 1270     write(6,1280)
 1280     format(' #entire first row is flagged values-program aborted')
          close(11)
          close(10)
 1290     return
          end

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).
c***********************************************************************
      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***********************************************************************
      subroutine get_dims(n1,n2,nx,ny,nri,id2,nxa)
      data lnri/16/,iw/6/
c
c  nri = blocking factor for rows (l)
c  nxa = number of rows added (n1-nx)
c  id2 = 2*max(n1,n2)
c
      nx=n1
      ny=n2
      l=lnri+1
      lnri21=lnri/2+1
c  set no. of rows for fft: need m=l*2**k, m.gt. or .eq. nx.
c  m=no. of rows, l=no. from 9-16, k=integer.
  190 l=l-1
      if(l.lt.lnri21) go to 260
      mr=n1/l+0.0000001
      k=1
      idiv=2
  200 iquot=mr/idiv+0.0000001
      if(iquot.lt.idiv) go to 210
      k=k+1
      mr=iquot
      go to 200
  210 k=k+1
      m=l*2**k
  220 mtest=l*2**(k-1)
      if(mtest.lt.n1) go to 230
      k=k-1
      m=mtest
      if(k.eq.0) go to 230
      go to 220
  230 lnxa=m-n1
      if(l.ne.lnri) go to 250
      nxa16=lnxa
  240 nri=l
      nxa=lnxa
      go to 190
  250 if(lnxa.ge.nxa) go to 190
      go to 240
  260 n1=n1+nxa
c  check to see if row block size of 16 will be more efficient
      n116=n1-nxa+nxa16
      ntest=0.9*n116
      if(ntest.gt.n1.or.n116.gt.2048) go to 270
      n1=n1-nxa+nxa16
      nxa=nxa16
      nri=16
  270 continue
      if(n1.gt.2048.or.n2.gt.2048) go to 320
      id2=n1
      if(n2.gt.n1)id2=n2
      id2=2*id2
      return
  320 write(iw,330)nx,ny,n1,n2
  330 format(' #no. of extended rows or columns exceeds 2048:'/
     1' input no. of rows and columns='2i4/
     2' no. of rows and columns required for fft =',2i4,/)
      stop
      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
      end
c
      subroutine bell()
      character ding*1
      ding = char(7)
      write(unit=*, fmt=*) ding
      return 
      end







