fourier3d.f90

Ground state of the time-independent Gross-Pitaevskii equation

!***
! The following code is for use with the IMSL library
!***

module fft_data
  integer, parameter :: dp=kind(1.d0)
  real(dp) :: fact
  real(dp), dimension(:,:), allocatable :: wsave
end module fft_data


subroutine fft_init(ng)
  !
  ! Initialize the Fourier transform
  !
  use fft_data
  implicit none

  integer, dimension(0:3), intent(in) :: ng

  integer :: m
  
  allocate(wsave(2*ng(0)+15,3))

  do m=1,3
     call dffti(ng(m),wsave(1,m))
  end do

  fact = 1.d0/real(ng(1)*ng(2)*ng(3),dp)

  return
end subroutine fft_init


subroutine fourier3D(ng,fin,fout,direction)
  !
  ! 3D Fourier transform
  !
  use fft_data
  implicit none

  integer, dimension(0:3), intent(in) :: ng
  real(dp), dimension(ng(1),ng(2),ng(3)), intent(in) :: fin
  real(dp), dimension(ng(1),ng(2),ng(3)), intent(out) :: fout
  integer, intent(in) :: direction

  integer :: i,j,k
  real(dp), dimension(ng(0)) :: f

  if(direction >= 0) then
     ! Forward transform
     fout = fin

     ! transform in x
     do k = 1,ng(3)
        do j = 1,ng(2)
           call dfftf(ng(1),fout(1,j,k),wsave(1,1))
        end do
     end do

     ! transform in y
     do k = 1,ng(3)
        do i = 1,ng(1)
           f(1:ng(2)) = fout(i,:,k)
           call dfftf(ng(2),f,wsave(1,2))
           fout(i,:,k) = f(1:ng(2))
        end do
     end do

     ! transform in z
     do j = 1,ng(2)
        do i = 1,ng(1)
           f(1:ng(3)) = fout(i,j,:)
           call dfftf(ng(3),f,wsave(1,3))
           fout(i,j,:) = f(1:ng(3))
        end do
     end do

  else
     ! Backward transform
     fout = fin

     ! transform in x
     do k = 1,ng(3)
        do j = 1,ng(2)
           call dfftb(ng(1),fout(1,j,k),wsave(1,1))
        end do
     end do

     ! transform in y
     do k = 1,ng(3)
        do i = 1,ng(1)
           f(1:ng(2)) = fout(i,:,k)
           call dfftb(ng(2),f,wsave(1,2))
           fout(i,:,k) = f(1:ng(2))
        end do
     end do

     ! transform in z
     do j = 1,ng(2)
        do i = 1,ng(1)
           f(1:ng(3)) = fout(i,j,:)
           call dfftb(ng(3),f,wsave(1,3))
           fout(i,j,:) = f(1:ng(3))
        end do
     end do
     
     fout = fout*fact

  end if


  return
end subroutine fourier3D