module fkvsp_util
  use fftflags_

  implicit none
  integer, private :: n1,n2
  real,    private :: s1,s2
  complex, dimension(:), pointer, private :: tmp1,tmp2

  integer(kind=dbl), private :: planF1, planF2 !! forward FFT
  integer(kind=dbl), private :: planI1, planI2 !! inverse FFT

contains
  !----------------------------------------------------------------
  subroutine pad2(rdat,cdat)
    real,    dimension(:,:), pointer :: rdat
    complex, dimension(:,:), pointer :: cdat

    cdat(1:size(rdat,1),1:size(rdat,2)) = cmplx(rdat)
  end subroutine pad2
  !----------------------------------------------------------------
  subroutine unpad2(rdat,cdat)
    real,    dimension(:,:), pointer :: rdat
    complex, dimension(:,:), pointer :: cdat

    rdat = real(cdat(1:size(rdat,1),1:size(rdat,2)))
  end subroutine unpad2
  !----------------------------------------------------------------
  integer function pow2(n) result(p)
    integer, intent(in) :: n

    p = 1
    do while( p < n )
       p = p * 2
    end do
  end function pow2
  !----------------------------------------------------------------
  subroutine ftinit(in1,in2)
    integer, intent(in) :: in1,in2

    n1=in1; s1=1./sqrt(1.*n1)
    n2=in2; s2=1./sqrt(1.*n2)

    allocate(tmp1(n1),tmp2(n2))

    call fftwnd_f77_create_plan(planF1,1,n1,FFTW_FORWARD  &
    ,    FFTW_ESTIMATE + FFTW_THREADSAFE + FFTW_IN_PLACE)
    call fftwnd_f77_create_plan(planI1,1,n1,FFTW_BACKWARD &
    ,    FFTW_ESTIMATE + FFTW_THREADSAFE + FFTW_IN_PLACE)

    call fftwnd_f77_create_plan(planF2,1,n2,FFTW_FORWARD  &
    ,    FFTW_ESTIMATE + FFTW_THREADSAFE + FFTW_IN_PLACE)
    call fftwnd_f77_create_plan(planI2,1,n2,FFTW_BACKWARD &
    ,    FFTW_ESTIMATE + FFTW_THREADSAFE + FFTW_IN_PLACE)
  end subroutine ftinit
  !----------------------------------------------------------------
  subroutine fft(cdat)
    complex, dimension(:,:), pointer :: cdat
    integer :: i1,i2

    do i2=1,n2
       do i1=1,n1,2
          cdat(i1,i2) = - cdat(i1,i2)
       end do

       tmp1 = cdat(:,i2)
       call fftwnd_f77_one(planF1,tmp1,0)
       cdat(:,i2) = tmp1*s1
    end do

    do i1=1,n1
       do i2=1,n2,2
          cdat(i1,i2) = - cdat(i1,i2)
       end do

       tmp2 = cdat(i1,:)
       call fftwnd_f77_one(planF2,tmp2,0)
       cdat(i1,:) = tmp2*s2
    end do
  end subroutine fft
  !----------------------------------------------------------------
  subroutine ift(cdat)
    complex, dimension(:,:), pointer :: cdat
    integer :: i1,i2

    do i2=1,n2
       tmp1 = cdat(:,i2)
       call fftwnd_f77_one(planI1,tmp1,0)
       cdat(:,i2) = tmp1*s1

       do i1=1,n1,2
          cdat(i1,i2) = - cdat(i1,i2)
       end do
    end do

    do i1=1,n1
       tmp2 = cdat(i1,:)
       call fftwnd_f77_one(planI2,tmp2,0)
       cdat(i1,:) = tmp2*s2

       do i2=1,n2,2
          cdat(i1,i2) = - cdat(i1,i2)
       end do
    end do
  end subroutine ift
  !----------------------------------------------------------------
end module fkvsp_util
