      SUBROUTINE fft(a,m,inv)
! Paul L. DeVries, Department of Physics, Miami University

! This subroutine performs the Fast Fourier Transform by
! the method of Cooley and Tukey

!        start:     1965
!  last modified:   1993
!

!  Parameter description by H. Sormann 2000
!  PARAMETERS:   a     vector of complex data which are to be
!                      Fourier-transformed
!                      this vector is VARIABLY dimensioned (see below)
!                      On return, a  contains the FT data.
!
!                m     2**m  is the number of data points
!
!                inv equal 1        Fourier transform
!                inv not equal 1    INVERSE Fourier transform
!
      IMPLICIT NONE

      INTEGER                               :: m,inv,n,nd2,i,j,k,l,le,le1,ip
      COMPLEX(KIND=KIND(0.d0))              :: u,w,t
      COMPLEX(KIND=KIND(0.d0)),DIMENSION(*) :: a    ! * = variable dimension
      DOUBLE PRECISION             :: ang,pi
!      PARAMETER (pi=3.141592653589793d0)
      pi=4.d0*atan(1.d0)

      n=2**m
      nd2=n/2
      j=1
      DO i=1,n-1
        IF(i .lt. j)THEN
          t=a(j)
          a(j)=a(i)
          a(i)=t
        ENDIF
        k=nd2
100     IF(k .lt. j)THEN
          j=j-k
          k=k/2
          goto 100
        ENDIF 
        j=j+k
      END DO
      le=1
      DO l=1,m
        le1=le
        le=le+le
        u=(1.d0,0.d0)
        ang=pi/dble(le1)
        w=cmplx(cos(ang), -sin(ang),KIND(0.d0))
        IF(inv .eq. 1)w=conjg(w)
        DO j=1,le1
          DO i=j,n,le
            ip=i+le1
            t=a(ip)*u
            a(ip)=a(i)-t
            a(i)=a(i)+t
          END DO
          u=u*w
        END DO
      END DO
      IF(inv .ne. 1)THEN 
        DO i=1,n
          a(i)=a(i)/dble(n)
        END DO
      ENDIF
      END SUBROUTINE fft
