c  FFT_ARAS.FOR
c
c  To calculate area-weighted radial amplitude spectrum = sum of|F| in dw
c  band, quadrants I and II only.
c
c  Note that |F| is scaled, ie divided by n1*n2.
c
c  Program adapted from F_COF by Lin Cordell, 8 Mar 89.
c
c  In complex frequency domain, reads the fft coefficient file "fftfil.cof"
c  (obtained  from FFTFIL using the icoef option).  This program follows
c  the blocking of 'NEWFFTFIL'.
c
c  All user-supplied changes occur in subroutine MODIFY, at several locations
c  indicated by "SUPPLY".  The main program just accommodates FFTFIL's
c  restructuring of the grid. 
c
c  See comments in subroutine MODIFY regarding switch of coordinate
c  axes from x=east, y=north; to x=north, y=east.
c
c  Lin Cordell, 9 Dec '87
c  Revised for PC 21 Aug '90.
c  Repairs 1 Feb 91.
c
      character title*56,pgm*8
      character*50 ingrd
      dimension g(2,2048),sum(1000),number(1000)
      data lnri/16/
      ierror=0
      write(6,10)
10    format(1x,'Enter grid from which Fourier coeff.s created:'/1x'*'$)
      read(5,100)ingrd
100   format(a50)
      open(10,file=ingrd,status='old',form='unformatted')
      read(10)title,pgm,n2,n1,nz,yo,dy,xo,dx
c  Note: n2,n1 and dx,dy rotate sense of rows and cols. Now x=north,
c   y=east.
      close(10)
      write(6,7)
7     format(' Enter no. of rows & cols added to grid(nadd):'/1x'*'$)
      read(5,*)nadd
      icmd=10
      l=lnri+1
      lnri21=lnri/2+1
c  set no. of rows for fft: need m=l*2**k, m.gt. or .eq. nx+2*nadd.
c  m=no. of rows, l=no. from 9-16, k=interger.
      n1=n1+2*nadd
  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 n2=n2+2*nadd
      write(6,275)nri
  275 format(' blocking rows =',i3)
      if(n1.gt.2048.or.n2.gt.2048) go to 320
      nxa=nxa+2*nadd
      id2=n1
      if(n2.gt.n1)id2=n2
      id2=2*id2
c  Blocking parameters for fft are set; call subroutine.
      n22=2*n2
      open(14,access='direct',status='old',
     1 form='unformatted',file='fftfil.cof',
     2 recl=n22*4)
      call modify(g,sum,number,n2,n1,n22,dy,dx)
      go to 99
320   write(6,321)
321   format(1x,'Extended rows/columns exceeds 2048; run aborted.')
99    continue
      close(14)
      close(12)
      stop
      end
      subroutine modify(g,a,n,n2,n1,n22,dy,dx)
c  Note:  In contrast with standard-grid notation, the x,dx,v1
c  direction is north; y,dy,v2 direction is east, and angles are
c  measured clockwise from the +x (north) direction.
      dimension g(n22),a(1),n(1)
      character*50 ofile
      ncp1=n2+1
      cn=1./(n2*dy)
      rn=1./(n1*dx)
      scale=float(n1)*float(n2)
      nrnq=float(n1)/2.+1.0000001
      ncnq=float(n2)/2.+1.0000001
      pi=3.1415927e0
c  If necessary, SUPPLY parameter read here.
      write(6,10)
10    format(1x,'Enter delta-w and w-max:'/' *'$)
      read(5,*)dw,wmax
      dwh=0.5*dw
      rdx=1./dx
      rdy=1./dy
      rmin=rdx
      if(rdy.lt.rdx)rmin=rdy
      if(2.*wmax.ge.rmin)go to 900
      nrowmax=1+wmax/rn
      ncolmax=1+wmax/cn
      write(6,14)
14    format(' Enter output file name:'/' *'$)
      read(5,15)ofile
15    format(a50)
      open(15,file=ofile,status='new',form='formatted')
      write(15,50)
50    format(6x,' w',11x,'sum |F|',4x,'n',5x,'mean |F|',4x,
     1'log mean |F|')
      do 11 kw=1,1000
      a(kw)=0.0
11    n(kw)=0
c
c  Read fftfil.coef file and modify.
c   First row, DC term treated separately.
      jr=1
      read(14,rec=jr)g
c   DC term.
      a(1)=a(1)+abs(g(1))
      n(1)=n(1)+1
c
c   Rest of row 1.
      do 1 i=2,n2
      ic=2*i
      ii=i-1
      if(i.gt.ncnq) ii=-(ncp1-i)
      v2=float(ii)*cn
c    (If necessary, here consider effect of sign of v2.)
      w=abs(v2)
      Gr=g(ic-1)
      Gi=g(ic)
      if(w.gt.wmax)go to 1
      kw=1+w/dw
      a(kw)=a(kw)+sqrt(Gr*Gr+Gi*Gi)
      n(kw)=n(kw)+1
1     continue
c
c  Read rest of fftfil.cof file, row-wise, and modify.
      do 3 jr=2,nrowmax
      jj=jr-1
      v1=float(jj)*rn
      v1sq=v1*v1
      read(14,rec=jr)g
      do 2 i=1,n2
      ic=2*i
      ii=i-1
      if(i.gt.ncnq) ii=-(ncp1-i)
      v2=float(ii)*cn
      v2sq=v2*v2
      w=sqrt(v1sq+v2sq)
c  (If necessary, consider effect of sgn of v1, v2.)
      Gr=g(ic-1)
      Gi=g(ic)
      if(w.gt.wmax)go to 2
      kw=1+w/dw
       a(kw)=a(kw)+sqrt(Gr*Gr+Gi*Gi)
      n(kw)=n(kw)+1
2     continue
3     continue
c
      c=1./(n1*n2)
      kwmax=1+wmax/dw
      do 12 kw=1,kwmax
      rkw=kw
      w=rkw*dw-dwh
      aras=a(kw)*c
      rasm=-1
      rasmlog=-1
      if(n(kw).ne.0)rasm=aras/n(kw)
      if(rasm.gt.0) rasmlog=log(rasm)
      write(15,13)w,aras,n(kw),rasm,rasmlog
13    format(1x,2e13.5,i5,2e13.5)
12    continue
      go to 999
900   write(6,901)
901   format(' Cant handle it.')
999   close(15)
      close(14)
      return
      end
