c  FFT : Fast Fourier Transform
c  Author: Richard Jones, richard.t.jones@uconn
c  Verion: 1.0, August 1, 2005
c
c  Purpose: FFT reads a sampled function from a 1d HBOOK histogram, computes
c	the FFT and stores it in another histogram.  The input function is
c	assumed to be sampled at the midpoints of the histogram bins.  The
c	FFT is computed at frequencies n*(2pi/L) where L is the range of
c	the input histogram and n is an integer in the interval [0,N/2]
c	for N the total number of bins in the input histogram.
c
c  Notes:
c    1. The output histogram is defined on the frequency interval from
c	n = -(N-1)/2 to n = N/2.  The negative frequency data points hold
c	the imaginary part of the FFT and the non-negative points contain
c       the real part.  If N is even then there is one less data point
c	below zero than above zero because the FFT is pure real at n=N/2.
c 	The number of bins in the output FFT histogram is the same as the
c	input histogram.


      subroutine fft(id1,id2)
      integer id1,id2
      character*80 title1,title2
      integer nbinx,nbiny,nwt,loc
      real xlo,xhi,ylo,yhi
      integer MAXBINS
      parameter (MAXBINS=99999)
      real fun(MAXBINS)
      real x,dx,f,df,ftr,fti
      real fhi,fmax
      integer nf
      integer i,n
      real PI
      parameter (PI=3.1415926)
      logical hexist
      external hexist
      if (.not.hexist(id1)) then
         print *, 'fft error: input histogram id',id1,' not found'
         return
      endif
      call hgive(id1,title1,nbinx,xlo,xhi,nbiny,ylo,yhi,nwt,loc)
      if (nbinx.gt.MAXBINS) then
         print *, 'fft error: current maximum is',MAXBINS,' data points'
         return
      endif
      dx=(xhi-xlo)/nbinx
      df=2*PI/(xhi-xlo)
      fmax=PI/dx
      nf=fmax/df*1.00001
      fhi=df*(nf+0.5)
 1000 format('FFT of ',a80)
      write(title2,1000) title1
      call hbook1(id2,title2,2*nf+1,-fhi,fhi,0.)
      call hunpak(id1,fun,'hist',0)
      do n=0,nf
         ftr=0
         fti=0
         f=n*df
         do i=1,nbinx
            x=xlo+dx*(i-0.5)
            ftr=ftr+fun(i)*cos(f*x)
            fti=fti+fun(i)*sin(f*x)
         enddo
         call hfill(id2,f,0.,ftr/nbinx)
         call hfill(id2,-f,0.,fti/nbinx)
      enddo
      end

      subroutine ffti(id1,id2,xlo)
      integer id1,id2
      character*80 title1,title2
      integer nbinx,nbiny,nwt,loc
      real xlo,xhi,flo,fhi,ylo,yhi
      integer MAXBINS
      parameter (MAXBINS=99999)
      real fft(MAXBINS)
      real x,dx,f,df,fun
      integer nf
      integer i,n
      real fmax
      real PI
      parameter (PI=3.1415926)
      complex weight
      external weight
      complex amp
      complex i_
      parameter (i_=(0.,1.))
      logical hexist
      external hexist
      if (.not.hexist(id1)) then
         print *, 'ffti error: input histogram id',id1,' not found'
         return
      endif
      call hgive(id1,title1,nbinx,flo,fhi,nbiny,ylo,yhi,nwt,loc)
      if (nbinx.gt.MAXBINS) then
         print *, 'ffti error: current max is',MAXBINS,' data points'
         return
      endif
      df=(fhi-flo)/nbinx
      fmax=fhi-df/2
      nf=nbinx/2
      dx=PI/fmax
      xhi=xlo+dx*nbinx
 1000 format('Inverse FFT of ',a80)
      write(title2,1000) title1
      call hbook1(id2,title2,nbinx,xlo,xhi,0.)
      call hunpak(id1,fft,'hist',0)
      do i=1,nbinx
         x=xlo+dx*(i-0.5)
         fun=0
         fun=fun+weight(f)*fft(nf+1)
         do n=1,nf
            f=n*df
            amp=2*weight(f)*(fft(nf+n+1)+i_*fft(nf-n+1))
            fun=fun+real(amp)*cos(f*x)
     +            +aimag(amp)*sin(f*x)
         enddo
         call hfill(id2,x,0.,fun)
      enddo
      end

      complex function weight(f)
      real f
      complex i_
      parameter (i_=(0.,1.))
** set the weight to unity to obtain the unfiltered function
      weight=-((f+i_*0.0129878074)*(f+i_*0.0763085186))
*     weight=-((f+i_*0.0107076885)*(f+i_*0.0025207386))
     +            *exp(-(f/0.2)**2)
*     weight=1.
      end