program FFT5d_
  !! 5-D FFT
  !! Paul Sava (paul@sep.stanford.edu)
  !! 12/1/00

  !! complex input/output
  !! input  is  z- hx- hy- mx- my-ro
  !! output is kz-khx-khy-kmx-kmy-ro

  use sep
  use utilities
  use fftwflags

  logical :: adj, verb, reserve
  integer :: stat

  real :: scale, pi, arg, sign
  real :: memory
#ifdef (LINUX)
  integer   :: plan
#else
  integer*8 :: plan
#endif

  type(paxis) :: az,ahx,ahy,amx,amy,ar
  integer     :: iz,ihx,ihy,imx,imy,ir

  !! z-hx-hy-mx-my-nimages
  complex, dimension(:,:,:,:,:), pointer :: data
  complex, dimension(:,:,:),     pointer :: sft
  integer, dimension(5) :: ni
  integer :: n6

  call sep_init()

  call from_param("adj",adj,.false.)
  call from_param("verb",verb,.false.)
  call from_param("reserve",reserve,.false.)

  call hinaxis(az ,1); az%l ="z";  call houaxis(az ,1)
  call hinaxis(ahx,2); ahx%l="hx"; call houaxis(ahx,2)
  call hinaxis(ahy,3); ahy%l="hy"; call houaxis(ahy,3)
  call hinaxis(amx,4); amx%l="mx"; call houaxis(amx,4)
  call hinaxis(amy,5); amy%l="my"; call houaxis(amy,5)
  call hinaxis(ar ,6); ar%l ="ro"; call houaxis(ar ,6)

  call from_history("n6",n6,1)

  memory = 0
  memory = memory + az%n*ahx%n*ahy%n*amx%n*amy%n*8/1024./1024.
  memory = memory + az%n*ahx%n*ahy%n            *8/1024./1024.

  if(verb) then
     write(0,*) "nz =",az%n
     write(0,*) "nhx=",ahx%n
     write(0,*) "nhy=",ahy%n
     write(0,*) "nmx=",amx%n
     write(0,*) "nmy=",amy%n
     write(0,*) "nr =",ar%n,n6

     write(0,*) "memory=",memory,"Mb"
  end if

  call sep_close()

  !----------------------------------------------------------------
  ni = (/az%n,ahx%n,ahy%n,amx%n,amy%n/)
  allocate(data(az%n,ahx%n,ahy%n,amx%n,amy%n), stat=stat)
  if(stat /=0) call seperr("cant allocate data")
  data = 0.
  !----------------------------------------------------------------
  if(reserve) then
     !! reserve enough disk space for the entire output
     if(verb) write(0,*) "reserving disk space"
     do ir=1,ar%n
        if(verb) write(0,*) ir
        call sep_write(data)
     end do
     call sseek('out',0,0)
     if(verb) write(0,*) "OK"
  end if
  !----------------------------------------------------------------
  scale = 1./sqrt(1.*product(ni))

  if(adj) then
     call fftwnd_f77_create_plan(plan,5,ni,FFTW_BACKWARD &
     ,    FFTW_ESTIMATE + FFTW_THREADSAFE + FFTW_IN_PLACE)
     sign = +1.
  else
     call fftwnd_f77_create_plan(plan,5,ni,FFTW_FORWARD  &
     ,    FFTW_ESTIMATE + FFTW_THREADSAFE + FFTW_IN_PLACE)
     sign = -1.
  end if
  !----------------------------------------------------------------
  allocate(sft(az%n,ahx%n,ahy%n), stat=stat)
  if(stat /=0) call seperr("cant allocate sft")
  sft=cmplx(1.,0.)

  if(az%o/=0. .or. ahx%o/=0. .or. ahy%o/=0 ) then

     pi = acos(-1.)
     do iz=1,az%n
        do ihx=1,ahx%n
           do ihy=1,ahy%n

              arg = 2*pi &
              * ( (iz -1.)/ az%n * az%o / az%d &
              +   (ihx-1.)/ahx%n * ahx%o/ahx%d &
              +   (ihy-1.)/ahy%n * ahy%o/ahy%d )
              
              sft(iz,ihx,ihy) = &
              cmplx(cos(arg),(sign)*sin(arg))

           end do
        end do
     end do

  end if
  !----------------------------------------------------------------

  if(adj) then
     do ir=1,ar%n
        if(verb) write(0,*) ir,'/',ar%n
        call sep_read(data)

        do imx=1,amx%n
           do imy=1,amy%n
              data(:,:,:,imx,imy)= data(:,:,:,imx,imy) * sft
           end do
        end do

        call fftwnd_f77_one(plan,data,0)
        data = data * scale

        call sep_write(data)
     end do
  else
     do ir=1,ar%n
        if(verb) write(0,*) ir,'/',ar%n
        call sep_read(data)

        call fftwnd_f77_one(plan,data,0)
        data = data * scale

        do imx=1,amx%n
           do imy=1,amy%n
              data(:,:,:,imx,imy)= data(:,:,:,imx,imy) * sft
           end do
        end do

        call sep_write(data)
     end do
  end if
  !----------------------------------------------------------------

  deallocate(data, stat=stat)
  if(stat /=0) call seperr("cant deallocate data")
  deallocate(sft,  stat=stat)
  if(stat /=0) call seperr("cant deallocate sft")

  call exit(0)
end program FFT5d_







