ccccccccccccccccccccccccc     Program 5.2     cccccccccccccccccccccccccc
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                      c
c Please Note:                                                         c
c                                                                      c
c (1) This computer program is part of the book, "An Introduction to   c
c     Computational Physics," written by Tao Pang and published and    c
c     copyrighted by Cambridge University Press in 1997.               c
c                                                                      c
c (2) No warranties, express or implied, are made for this program.    c
c                                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      SUBROUTINE FFT (AR,AI,N,M)
      implicit none
C
C An example of the fast Fourier transform subroutine with
C N = 2**M. AR and AI are the real and imaginary part of
C data in the input and corresponding Fourier coefficients
C in the output.
C
      REAL*8 AR,AI,PI,A1,A2,Q,U,V
      INTEGER M,N,N2,N1,I,J,K,L,L1,L2

      DIMENSION AR(N),AI(N)
C
      PI = 4.0*ATAN(1.0)
      N2 = N/2
C
      N1 = 2**M
      IF(N1.NE.N) STOP 'Indices do not match'
C
C Rearrange the data to the bit reversed order
C
      L = 1
      DO      150  K = 1, N-1
        IF(K.LT.L) THEN
          A1    = AR(L)
          A2    = AI(L)
          AR(L) = AR(K)
          AR(K) = A1
          AI(L) = AI(K)
          AI(K) = A2
        ENDIF
        J   = N2
        DO      100  WHILE (J.LT.L)
          L = L - J
          J = J/2
  100   END DO
        L = L + J
  150 CONTINUE
C
C Perform additions at all levels with reordered data
C
      L2 = 1
      DO      200  L = 1, M
        Q  =  0.0
        L1 =  L2
        L2 =  2*L1
        DO    190  K = 1, L1
          U   =  COS(Q)
          V   = -SIN(Q)
          Q   =  Q + PI/L1
          DO  180  J = K, N, L2
            I     =  J + L1
            A1    =  AR(I)*U - AI(I)*V
            A2    =  AR(I)*V + AI(I)*U
            AR(I) =  AR(J) - A1
            AR(J) =  AR(J) + A1
            AI(I) =  AI(J) - A2
            AI(J) =  AI(J) + A2
  180     CONTINUE
  190   CONTINUE
  200   CONTINUE
        RETURN
        END