SUBROUTINE FFT(g,p)
*     fast Fourier transform of complex input data set g
*     transform returned in g
      IMPLICIT NONE
      REAL*8 Wp(2), angle, factor(2), pi, temp(2), g(2,*)
      PARAMETER(pi = 3.141592653589793D0)
      INTEGER L, N, N2, del_i, i, ip, j, jmax, k, p, Itmp1_
      N = 2**p                              ! number of data points
      N2 = N/2
*     rearrange input data according to bit reversal
      j = 1
      DO i = 1, N-1
*        g(i,1) is real part of f((i-1)*del_t)
*        g(i,2) is imaginary part of f((i-1)*del_t)
*        set g(i,2) = 0 if data real
         IF(i .LT. j) THEN
            temp(1) = g(1,j)
            temp(2) = g(2,j)
            g(1,j) = g(1,i)
            g(2,j) = g(2,i)
            g(1,i) = temp(1)
            g(2,i) = temp(2)
         END IF
         k = N2
         DO WHILE(k .LT. j)
            j = j-k
            k = k/2
         END DO
         j = j + k
      END DO
*     begin Danielson-Lanczos algorithm
      jmax = 1
      DO L = 1, p
         del_i = 2*jmax
         Wp(1) = 1                          !  Wp initialized at W^0
         Wp(2) = 0
         angle = pi/jmax
         factor(1) = cos(angle)             ! ratio of new to old W^p    
         factor(2) = -sin(angle)
         DO j = 1, jmax
            DO i = j, N, del_i
*              calculate transforms of length 2^L
               ip = i + jmax
               temp(1) = g(1,ip)*Wp(1) - g(2,ip)*Wp(2)
               temp(2) = g(1,ip)*Wp(2) + g(2,ip)*Wp(1)
               g(1,ip) = g(1,i) - temp(1)
               g(2,ip) = g(2,i) - temp(2)
               g(1,i) = g(1,i) + temp(1)
               g(2,i) = g(2,i) + temp(2)
            END DO
*           find new W^p
            temp(1) = Wp(1)*factor(1) - Wp(2)*factor(2)
            temp(2) = Wp(1)*factor(2) + Wp(2)*factor(1)
            DO Itmp1_ = 1, 2
               Wp(Itmp1_) = temp(Itmp1_)
            END DO
         END DO
         jmax = del_i
      END DO
      END