module FourierTransform
	implicit none
	!private calcfft
	private calcfft2
	
	interface
		recursive function calcfft(x, sgn) result (c)
			use kind_defs
			implicit none
			complex(dbl), intent(in) :: x(:)
			integer, intent(in) :: sgn
			complex(dbl) :: c(size(x))
		end function
		
		function calcfft2(x, sgn) result (c)
			use kind_defs
			implicit none
			complex(dbl), intent(in) :: x(:,:)
			integer, intent(in) :: sgn
			complex(dbl) :: c(size(x,1),size(x,2))
		end function
		
		logical function ispower2(N)
			integer, intent(in) :: N
		end function
		
		integer function nextpow2(m)
			integer, intent(in) :: m
		end function
	end interface
	
	interface fft
		module procedure fftc
		module procedure fftd
	end interface
	interface ifft
		module procedure ifftc
		module procedure ifftd
	end interface
	interface fft2
		module procedure fft2c
		module procedure fft2d
	end interface
	interface ifft2
		module procedure ifft2c
		module procedure ifft2d
	end interface
contains
	! Fast Fourier Transform:
	function fftc(x)
		use kind_defs
		implicit none
		complex(dbl), intent(in) :: x(:)
		complex(dbl) :: fftc(size(x))
		
		fftc = calcfft(x, -1)
	end function
	function fftd(x)
		use kind_defs
		implicit none
		real(dbl), intent(in) :: x(:)
		complex(dbl) :: fftd(size(x))
	
		fftd = calcfft(cmplx(x,kind=dbl), -1)
	end function
	
	! Inverse FFT:
	function ifftc(x)
		use kind_defs
		implicit none
		complex(dbl), intent(in) :: x(:)
		complex(dbl) :: ifftc(size(x))
	
		ifftc = calcfft(x, 1)
	end function
	function ifftd(x)
		use kind_defs
		implicit none
		real(dbl), intent(in) :: x(:)
		complex(dbl) :: ifftd(size(x))
	
		ifftd = calcfft(cmplx(x,kind=dbl), 1)
	end function
	
	! 2D-FFT:
	function fft2c(x)
		use kind_defs
		implicit none
		complex(dbl), intent(in) :: x(:,:)
		complex(dbl) :: fft2c(size(x,1),size(x,2))
		
		fft2c = calcfft2(x, -1)
	end function
	function fft2d(x)
		use kind_defs
		implicit none
		real(dbl), intent(in) :: x(:,:)
		complex(dbl) :: fft2d(size(x,1),size(x,2))
		
		fft2d = calcfft2(cmplx(x,kind=dbl), -1)
	end function
	
	! Inverse 2D-FFT:
	function ifft2c(x)
		use kind_defs
		implicit none
		complex(dbl), intent(in) :: x(:,:)
		complex(dbl) :: ifft2c(size(x,1),size(x,2))
		
		ifft2c = calcfft2(x, 1)
	end function
	function ifft2d(x)
		use kind_defs
		implicit none
		real(dbl), intent(in) :: x(:,:)
		complex(dbl) :: ifft2d(size(x,1),size(x,2))
		
		ifft2d = calcfft2(cmplx(x,kind=dbl), 1)
	end function
end module

logical function ispower2(N)
	implicit none
	integer, intent(in) :: N
	
	ispower2 = (IAND(N, N-1) == 0)
end function

integer function nextpow2(m)
	implicit none
	integer, intent(in) :: m
	
	nextpow2 = nint( log10(real(m))/log10(2.0) )
end function

recursive function calcfft(x, sgn) result (c)
	use kind_defs
	implicit none
	complex(dbl), intent(in) :: x(:)
	integer, intent(in) :: sgn
	complex(dbl) :: c(size(x))
	integer :: M, N, k, j
	complex(dbl), dimension(size(x)/2) :: odd, even, co, ce
	complex(dbl) :: w
	
	N = size(x)
	w = exp(sgn*2.0_dbl*pi*complexi/real(N))
	
	if (mod(N,2) == 0) then
		! Fast Fourier Transformation:
		! Split the data in to odd and even:
		M = N/2
		do k = 1,M
			odd(k) = x(2*k)
			even(k) = x(2*k-1)
		enddo
		
		! Calculate the FFT of the even and odd part:
		co = calcfft(odd, sgn)
		ce = calcfft(even, sgn)
		
		! Calculate the coefficients:
		do k = 1,M
			c(k) = ce(k) + co(k) * w**(k-1)
		enddo
		do k = M+1,N
			c(k) = ce(k-M) + co(k-M) * w**(k-1)
		enddo
		c = c/sqrt(2.0_dbl)
	else
		! Slow Fourier Transformation:
		do k = 1,N
			c(k) = (0.0_dbl,0.0_dbl)
			do j = 1,N
				c(k) = c(k) + x(j) * w**((j-1)*(k-1))
			enddo
		enddo
		if (sgn == 1) c = c/(real(N))
	endif
end function

function calcfft2(x, sgn) result (c)
	use kind_defs
	use FourierTransform, only : calcfft
	implicit none
	complex(dbl), intent(in) :: x(:,:)
	integer, intent(in) :: sgn
	complex(dbl) :: c(size(x,1),size(x,2))
	integer :: N1, N2, i

	! Get matrix dimensions:
	N1 = size(x,1)
	N2 = size(x,2)
	
	! Transform by rows:
	do i = 1,N1
		c(i,:) = calcfft(x(i,:), sgn)
	enddo
	
	! Transform by cols:
	do i = 1,N2
		c(:,i) = calcfft(c(:,i), sgn)
	enddo
end function