ccccccccccccccccccccccccc     Program 5.1     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
      PROGRAM DFT_EXAMPLE
C
C Example of the discrete Fourier transform with function
C f(x) = x(1-x) in [0,1]. The inverse transform is also
C performed for comparison.
C
      PARAMETER (N = 128,M=8)
      DIMENSION FR(N),FI(N),GR(N),GI(N)
C
      F0 = 1.0/SQRT(FLOAT(N))
      H  = 1.0/(N-1)
C
      DO      150  I = 1, N
        X     = H*(I-1)
        FR(I) = X*(1.0-X)
        FI(I) = 0.0
  150 CONTINUE
C
      CALL DFT (FR,FI,GR,GI,N)
      DO      200  I = 1, N
        GR(I) = F0*GR(I)
        GI(I) = F0*GI(I)
  200 CONTINUE
C
C     Perform inverse Fourier transform
C
      DO      300  I = 1, N
        GI(I) = -GI(I)
  300 CONTINUE
      CALL DFT (GR,GI,FR,FI,N)
      DO      400  I = 1, N
        FR(I) =  F0*FR(I)
        FI(I) = -F0*FI(I)
  400 CONTINUE
      WRITE (6,999) (H*(I-1),FR(I),I=1,N,M)
      WRITE (6,999)  H*(N-1),FR(N)
      STOP
  999 FORMAT (2F16.8)
      END
C
      SUBROUTINE DFT (FR,FI,GR,GI,N)
C
C Subroutine to perform the discrete Fourier transform with
C FR and FI as the real and imaginary parts of the input and
C GR and GI as the corresponding  output.
C
      DIMENSION FR(N),FI(N),GR(N),GI(N)
      PI = 4.0*ATAN(1.0)
      X  = 2*PI/N
C
      DO      150  I = 1, N
        GR(I) = 0.0
        GI(I) = 0.0
        DO    100  J = 1, N
          Q     = X*(J-1)*(I-1)
          GR(I) = GR(I) + FR(J)*COS(Q) + FI(J)*SIN(Q)
          GI(I) = GI(I) + FI(J)*COS(Q) - FR(J)*SIN(Q)
  100   CONTINUE
  150 CONTINUE
      RETURN
      END