C
C	  Calcul d'une transformee de Fourier rapide FFT sur un
C     tableau de complexes. Inspire de AL Garcia (Numerical
C     Methods for Physics).
	  SUBROUTINE FFT(A, N)
      IMPLICIT NONE

C     Description de l'interface de la subroutine
      INTEGER N    ! Nombre de points du vecteur complexe d'entree
                   ! ATTENTION : ce nombre doit etre une puissance de 2
      COMPLEX A(N) ! Vecteur complexe d'entree
C
C     En entree le tableau de complexes donnee contient les donnees
C     En sortie le tableau de complexes donnee contient les transformees

C     Declaration des constantes et parametres
	  REAL PI
      PARAMETER (PI = 3.141592654)

C     Declaration des variables
      INTEGER half_N, Nm1,i,j,k, m,ki,ki1,ip
      COMPLEX T, U, W
      REAL angle

C     Initisalisation des variables
      half_N = N/2
      Nm1 = N-1
	  m = INT(LOG(REAL(N))/LOG(2.0) + 0.5)  ! N = 2**m

C     Bit-scramble du vecteur d'entree
      j = 1
      DO i=1,Nm1
         IF (i.LT.j) THEN    ! on swappe A(i) et A(j)
           T = A(j)
           A(j) = A(i)
           A(i) = T
         ENDIF
         k = half_N
         DO WHILE (k.LT.j)
            j = j-k
            k = k/2
         ENDDO
         j = j+k
      ENDDO

C     Calcul de la transformee (Butterfly operation)
	  DO k=1,m
         ki = 2**k
         ki1 = ki/2
         U = (1.0, 0.0)
         angle = PI/REAL(ki1)
         W = CMPLX(COS(angle),-SIN(angle))
         DO j=1,ki1
            DO i=j,N, ki
              ip = i+ki1
              T = A(ip)*U
			  A(ip) = A(i)-T
              A(i) = A(i)+T
            ENDDO
            U = U*W
         ENDDO
      ENDDO

	  RETURN
	  END