program FT1d_
  !! 1-D FFT
  !! Paul Sava (paul@sep.stanford.edu)
  !! 12/6/03

  !! complex input/output
  !! input  is  a1-a2-a3-a4-a5
  !! output is ka1-a2-a3-a4-a5

  use sep
  use util_
  use fftwflags

  integer, parameter :: dbl=selected_real_kind(10,50)
  integer(kind=dbl)  :: plan

  logical :: adj, verb
  integer :: stat
  real    :: scale, pi, arg, sign

  type(myaxis):: a1,a2,a3,a4,a5
  integer     :: i1,i2,i3,i4,i5

  complex, dimension(:,:) :: data
  complex, dimension(:)   :: tmpi,tmpo

  logical :: shift
  complex, dimension(:)   :: sft

  logical :: center
  complex, dimension(:)   :: cen

  integer :: ntot      !! total n2...n5
  integer :: ninc,iinc !! in core
  integer :: nblk,iblk !! no of blocks

  call sep_init()

  call from_param("adj"   ,adj,   .false.)
  call from_param("verb"  ,verb,  .false.)
  call from_param("shift" ,shift, .true. )
  call from_param("center",center,.false.)

  call hinaxis(a1,1);call houaxis(a1,1)
  call hinaxis(a2,2);call houaxis(a2,2)
  call hinaxis(a3,3);call houaxis(a3,3)
  call hinaxis(a4,4);call houaxis(a4,4)
  call hinaxis(a5,5);call houaxis(a5,5)
  if(verb) then
     call printaxis(a1)
     call printaxis(a2)
     call printaxis(a3)
     call printaxis(a4)
     call printaxis(a5)
  end if

  !! how many traces to read in core
  call from_param("ninc",ninc,1)
  ntot = a2%n*a3%n*a4%n*a5%n
  if(verb) then
     write(0,*) "ntot=",ntot
     write(0,*) "ninc=",ninc
  end if

  if(ntot/ninc * ninc == ntot) then
     nblk = ntot/ninc
  else
     if(verb) write(0,*) "WARNING: ninc does not divide ntot; reseting..."
     do while(ninc>=1 .and. ntot/ninc * ninc /= ntot)
        ninc=ninc-1
     end do
     nblk = ntot/ninc
     if(verb) write(0,*) "ninc=",ninc
  end if
  if(verb) write(0,*) "nblk=",nblk

  call sep_close()

  !----------------------------------------------------------------
  allocate(data(a1%n,ninc), stat=stat)
  if(stat /=0) call seperr("cant allocate data")
  allocate( tmpi(a1%n), stat=stat)
  if(stat /=0) call seperr("cant allocate tmpi")
  allocate( tmpo(a1%n), stat=stat)
  if(stat /=0) call seperr("cant allocate tmpo")
  !----------------------------------------------------------------
  if(adj) then
     call fftw_f77_create_plan(plan,a1%n,FFTW_BACKWARD,FFTW_ESTIMATE)
     sign = +1.
  else
     call fftw_f77_create_plan(plan,a1%n,FFTW_FORWARD, FFTW_ESTIMATE)
     sign = -1.
  end if
  !----------------------------------------------------------------
  scale = 1./sqrt(1.*a1%n)
  !----------------------------------------------------------------
  allocate(sft(a1%n), stat=stat)
  if(stat /=0) call seperr("cant allocate sft")
  sft=cmplx(1.,0.)

  if(shift) then
     if(a1%o/=0.) then
        pi = acos(-1.)
        do i1=1,a1%n
           arg = 2*pi &
           * (i1-1.)/a1%n * a1%o/a1%d
           sft(i1) = &
           cmplx(cos(arg),(sign)*sin(arg))
        end do

     end if
  end if
  !----------------------------------------------------------------
  allocate(cen(a1%n), stat=stat)
  if(stat /=0) call seperr("cant allocate cen")
  cen=cmplx(1.,0.)

  if(center) then
     do i1=1,a1%n,2
        cen(i1)=-1
     end do
     do i1=2,a1%n,2
        cen(i1)=+1
     end do
  end if
  !----------------------------------------------------------------
  if(adj) then
     do iblk=1,nblk
        if(verb) write(0,*) iblk,'/',nblk
        call sep_read(data)

        do iinc=1,ninc
           tmpi = data(:,iinc) * scale
           if(shift) tmpi = tmpi * sft !! origin shift
           call fftw_f77_one(plan,tmpi,tmpo)
           data(:,iinc) = tmpo * cen
        end do

        call sep_write(data)
     end do
  else
     do iblk=1,nblk
        if(verb) write(0,*) iblk,'/',nblk
        call sep_read(data)

        do iinc=1,ninc
           tmpi = data(:,iinc) * cen
           call fftw_f77_one(plan,tmpi,tmpo)
           if(shift) tmpo = tmpo * sft !! origin shift
           data(:,iinc) = tmpo * scale
        end do

        call sep_write(data)
     end do
  end if
  !----------------------------------------------------------------
  call exit(0)
end program FT1d_







