module fft5d_mod
!! 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 fftwflags
  use sep3d_struct_mod
  use utilities_mod
  implicit none
  integer*8,private,save :: planF,planI
  integer, private, save :: nlast(5)
  complex, dimension(:,:,:),     pointer ,private:: sft
 

contains

logical function fft_init(struct,data,first)
type(sep3d) :: struct
complex :: data(:,:,:,:,:)
logical :: first
integer :: iz,ihx,ihy,i
integer, external :: sep_thread_num
real    :: arg,scale,pi


if(any(nlast /= struct%n)) then
 nlast=struct%n
 if(.not. first) then
   call sfftw_destroy_plan(planF)
   call sfftw_destroy_plan(planI)
   deallocate(sft)
 end if
 allocate(sft(struct%n(1),struct%n(2),struct%n(3)))
 call sfftw_plan_dft(planF,5,struct%n,data(1,1,1,1,1),data(1,1,1,1,1),&
    ffTW_FORWARD,fftW_ESTIMATE)
 call sfftw_plan_dft(planI,5,struct%n,data(1,1,1,1,1),data(1,1,1,1,1),&
    ffTW_BACKWARD,fftW_ESTIMATE)
  if(struct%o(1)/=0. .or. struct%o(2)/=0. .or. struct%o(3)/=0 ) then
    pi = acos(-1.)
    scale=1.
    do i=1,5
      scale=scale/struct%n(i)
    end do
    scale=sqrt(scale)
     do ihx=1,struct%n(3)
       do ihy=1,struct%n(2)
         do iz=1,struct%n(1)
              arg = 2*pi &
              * ( (iz -1.)/ struct%n(1) * struct%o(1) / struct%d(1) &
              +   (ihx-1.)/struct%n(2) * struct%o(2)/struct%d(2) &
              +   (ihy-1.)/struct%n(3) * struct%o(3)/struct%d(3) )
              
              sft(iz,ihx,ihy) = &
              cmplx(cos(arg),sin(arg))*scale
           end do
        end do
     end do
 	else 
    scale=1.
    do i=1,5
      scale=scale/struct%n(i)
    end do
    sft=cmplx(sqrt(scale),0.)
  end if
end if
  fft_init=.true.
end function

logical function fft_forward(struct,data)
type(sep3d) :: struct
complex     :: data(:,:,:,:,:)
integer     :: i4,i5 
call sfftw_execute(planF)
do i5=1,size(data,5)
  do i4=1,size(data,4)
    data(:,:,:,i4,i5)=data(:,:,:,i4,i5)*conjg(sft)
  end do
end do
fft_forward=.true.
end function


logical function fft_backward(struct,data)
type(sep3d) :: struct
complex     :: data(:,:,:,:,:)
integer     :: i4,i5 
do i5=1,size(data,5)
  do i4=1,size(data,4)
    data(:,:,:,i4,i5)=data(:,:,:,i4,i5)*sft
  end do
end do
call sfftw_execute(planI)
fft_backward=.true.
end function



end module
