module fft_mod
  implicit none
  integer,       parameter :: dp=selected_real_kind(15,300)
  real(kind=dp), parameter :: pi=3.141592653589793238460_dp
contains
 
  ! In place Cooley-Tukey FFT
  recursive subroutine fft(x)
    complex(kind=dp), dimension(:), intent(inout)  :: x
    complex(kind=dp)                               :: t
    integer                                        :: N
    integer                                        :: i
    complex(kind=dp), dimension(:), allocatable    :: even, odd
 
    N=size(x)
 
    if(N .le. 1) return
 
    allocate(odd((N+1)/2))
    allocate(even(N/2))
 
    ! divide
    odd =x(1:N:2)
    even=x(2:N:2)
 
    ! conquer
    call fft(odd)
    call fft(even)
 
    ! combine
    do i=1,N/2
       t=exp(cmplx(0.0_dp,-2.0_dp*pi*real(i-1,dp)/real(N,dp),kind=dp))*even(i)
       x(i)     = odd(i) + t
       x(i+N/2) = odd(i) - t
    end do
 
    deallocate(odd)
    deallocate(even)
 
  end subroutine fft

! In place Cooley-Tukey inverse FFT
  recursive subroutine ifftwork(x)
    complex(kind=dp), dimension(:), intent(inout)  :: x
    complex(kind=dp)                               :: t
    integer                                        :: N
    integer                                        :: i
    complex(kind=dp), dimension(:), allocatable    :: even, odd
 
    N=size(x)
 
    if(N .le. 1) return
 
    allocate(odd((N+1)/2))
    allocate(even(N/2))
 
    ! divide
    odd =x(1:N:2)
    even=x(2:N:2)
 
    ! conquer
    call ifftwork(odd)
    call ifftwork(even)
 
    ! combine
    do i=1,N/2
       t=exp(cmplx(0.0_dp,2.0_dp*pi*real(i-1,dp)/real(N,dp),kind=dp))*even(i)
       x(i)     = odd(i) + t
       x(i+N/2) = odd(i) - t
    end do
 
    deallocate(odd)
    deallocate(even)
 
  end subroutine ifftwork

  subroutine ifft(x)
    complex(kind=dp), dimension(:), intent(inout)  :: x
    integer                                        :: N
    integer                                        :: i
    call ifftwork(x)
    N=size(x)
    do i=1,N
      x(i)=x(i)/N
    end do
  end subroutine ifft

  function model(c,t)
    real(kind=dp):: model,c(6),t
    model=c(1)*sin(c(2)+c(4)*sin(c(5)+c(6)*t)+c(3)*t)
  end function model
 
end module fft_mod
 
program test
  use fft_mod
  implicit none
  integer,       parameter :: fftsize=1048576
!  complex(kind=dp), dimension(8) :: data = (/1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0/)
  complex(kind=dp), dimension(fftsize):: y
  real(kind=dp), dimension(6):: c=(/1,0,1,1,0,2/)
  real(kind=dp):: t
  integer :: i,t1,t2,count_rate

  do i=1,fftsize
    t=2*pi*(i-1)/fftsize;
    y(i)=model(c,t)
  end do

  call system_clock(t1,count_rate)
  call ifft(y)
  call system_clock(t2)
  write(*,*) 'Elapsed time ',real(t2-t1)/count_rate,' sec'
  do i=1,8
     write(*,'("(", F20.15, ",", F20.15, "i )")') y(i)
  end do
 
end program test
