module dimfft_
  !! image perturbation FFTs

  use utilities
  use fftwflags
  implicit none

  integer, parameter,                 private :: dbl=selected_real_kind(10,50)
  integer(kind=dbl),                  private :: planF, planI
  real,                               private :: signF, signI
  complex, dimension(:,:,:), pointer, private :: sftF, sftI
  real,                               private :: scale

  type(paxis) :: laz,lahx,lahy,lamx,lamy

contains
  !----------------------------------------------------------------
  subroutine dimfft_init(az,ahx,ahy,amx,amy)
    type(paxis) :: az,ahx,ahy,amx,amy
    integer     :: iz,ihx,ihy,imx,imy
    integer :: stat
    integer, dimension(5) :: ni
    real    :: pi, arg, sign

    laz=az
    lahx=ahx
    lahy=ahy
    lamx=amx
    lamy=amy

    ni = (/az%n,ahx%n,ahy%n,amx%n,amy%n/)
    pi = acos(-1.)

    scale = 1./sqrt(1.*product(ni))

    !----------------------------------------------------------------
    call fftwnd_f77_create_plan(planF,5,ni,FFTW_FORWARD &
    ,    FFTW_ESTIMATE + FFTW_THREADSAFE + FFTW_IN_PLACE)
    signF = +1.

    allocate(sftF(az%n,ahx%n,ahy%n), stat=stat)
    if(stat /=0) call seperr("cant allocate sftF")
    sftF=cmplx(1.,0.)

!!$    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 )
!!$
!!$             sftF(iz,ihx,ihy) = &
!!$             cmplx(cos(arg),(sign)*sin(arg))
!!$
!!$          end do
!!$       end do
!!$    end do
    !----------------------------------------------------------------
    call fftwnd_f77_create_plan(planI,5,ni,FFTW_BACKWARD &
    ,    FFTW_ESTIMATE + FFTW_THREADSAFE + FFTW_IN_PLACE)
    signI = -1.

    allocate(sftI(az%n,ahx%n,ahy%n), stat=stat)
    if(stat /=0) call seperr("cant allocate sftI")
    sftI=cmplx(1.,0.)

!!$    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 )
!!$
!!$             sftI(iz,ihx,ihy) = &
!!$             cmplx(cos(arg),(sign)*sin(arg))
!!$
!!$          end do
!!$       end do
!!$    end do
    !----------------------------------------------------------------    
  end subroutine dimfft_init
  !----------------------------------------------------------------  
  subroutine dimfft(cdata)
    complex, dimension(:,:,:,:,:), pointer :: cdata
    integer :: imx,imy

    call fftwnd_f77_one(planF,cdata,0)
    cdata = cdata * scale

!    do imx=1,lamx%n
!       do imy=1,lamy%n
!          cdata(:,:,:,imx,imy) = &
!          cdata(:,:,:,imx,imy) * sftF
!       end do
!    end do

  end subroutine dimfft
  !----------------------------------------------------------------
  subroutine dimift(cdata)
    complex, dimension(:,:,:,:,:), pointer :: cdata
    integer :: imx,imy

!    do imx=1,lamx%n
!       do imy=1,lamy%n
!          cdata(:,:,:,imx,imy) = &
!          cdata(:,:,:,imx,imy) * sftI
!       end do
!    end do

    call fftwnd_f77_one(planI,cdata,0)
    cdata = cdata * scale

  end subroutine dimift
  !----------------------------------------------------------------
end module dimfft_











