module fft
  implicit none
  real, parameter :: pi=3.141592654

  contains

subroutine fourier (a, m, n, sig)
!  Fourier transform of a complex variable.
!  The original data is in the vector a(0:n-1)
!  The length of the vector, n, must be 
!  a power of two, n=2**m
!  The original vector a(0..n-1) is overwritten
!  byt the transform.
!  If sig < 0, calculate the inverse transform

  implicit none
  integer, intent(in) :: m, n, sig
  complex, dimension(0:n), intent(inout) :: a
  integer i,j,k,l, half_n, le,le1,ip
  real  angle 
  complex x,u,w,t 

  half_n = n/2

  ! if inverse transform, conjugate a
  if (sig < 0) a = conjg(a)

  ! rearrange the vector so that a(s)=a(j), where s
  ! is the binary representation of j in reverse order
  j=1 
  do i=1, n-1
    if (i < j)  then 
      x = a(j-1) ; a(j-1) = a(i-1) ;  a(i-1) = x ; 
    end if
    k = half_n 
    do while (k < j)  
       j = j-k ;  k = k/2 
    end do
    j = j+k
  end do

  ! transform pieces of length le=2,4,8,...
  le=2
  do l=1, m
    le1 = le/2 
    u = cmplx(1,0)
    angle = pi/le1
    w = cmplx(cos(angle),sin(angle))
    do j=1,le1
      do i=j,n,le
        ip=i+le1 
        x = a(ip-1) 
	t = x * u
	a(ip-1) = a(i-1)-t
	a(i-1) = a(i-1)+t
      end do
      x = u * w
      u = x
    end do
    le = 2*le
  end do

  ! if inverse transform, conjugate a bck
  ! and normalize
  if (sig < 0) a=conjg(a)/n

end subroutine

end module