ccccccccccccccccccccccccc     Program 5.3     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 FFT2D (FR,FI,N1,N2,M1,M2)
C
C Subroutine for the two-dimensional fast Fourier transform
C with N=N1*N2 and N1=2**M1 and N2=2**M2.
C
      DIMENSION FR(N1,N2), FI(N1,N2), GR(N1), GI(N1)
C
C Transformation on the second index
C
      DO        100  I = 1, N1
        CALL FFT (FR(I,1),FI(I,1),N2,M2)
  100 CONTINUE
C
C Transformation on the first index
C
      DO        250  J = 1, N2
        DO      230  I = 1, N1
          GR(I) = FR(I,J)
          GI(I) = FI(I,J)
  230   CONTINUE
        CALL FFT (GR,GI,N1,M1)
        DO      240  I = 1, N1
          FR(I,J) = GR(I)
          FI(I,J) = GI(I)
  240   CONTINUE
  250 CONTINUE
      RETURN
      END
C
      SUBROUTINE FFT (AR,AI,N,M)
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
      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
      END