PROGRAM FFTdemo
      IMPLICIT NONE
      CHARACTER*80 Stmp1_
      REAL*8 g, fn, gr, gi
      REAL pi, GWaspect
      PARAMETER(pi = 3.141592)
      INTEGER IR_, N, p
      EXTERNAL fn, gr, gi
      COMMON /ARRY/g(2,1024), N
      CALL GWopen(IR_, 0)
      CALL GWindow(IR_, -4.0*GWaspect(2), -2.0, 4.0*GWaspect(2), 2.0)
      CALL GWsetpen(IR_, 0, -1, -1, -1)
      CALL GWline(IR_, -pi*1.1, 0.0, pi*1.1, 0.0)
      CALL GWline(IR_, 0.0, -1.5, 0.0, 1.5)
      CALL plotf(fn, 0)
      CALL GWinput(IR_, 'p = ', Stmp1_)
      IF(IR_ .NE. 0) THEN
         READ(Stmp1_,*) p
         N = 2**p
         CALL coefficients(fn)
         CALL FFT(g, p)
         CALL plotf(gr, 13)
         CALL plotf(gi, 16)
      END IF
      CALL GWquit(IR_)
      END
C
      SUBROUTINE coefficients(f)
*     transform returned in g
      IMPLICIT NONE
      REAL*8 f, g, pi, t, dt, pi2
      PARAMETER(pi = 3.141592653589793D0, pi2 = 2*pi)
      INTEGER N, k
      COMMON /ARRY/g(2,1024), N
      t = 0
      dt = pi2/N
      DO k = 1, N
         g(1, k) = f(t)
         g(2, k) = 0
         t = t + dt
      END DO
      END
C
      SUBROUTINE plotf(f, c)
      IMPLICIT NONE
      INTEGER IR_, N, c, nplot
      REAL*8 f, g, pi, t, dt
      PARAMETER(pi = 3.141592653589793D0)
      COMMON /ARRY/g(2,1024), N
      DATA nplot/100/
      dt = pi/nplot
      t = -pi
      CALL GWsetpen(IR_, c, -1, -1, -1)
      CALL GWmove2(IR_, REAL(t), REAL(f(t)))
      DO WHILE(t .LT. pi)
         t = t + dt
         CALL GWline2(IR_, REAL(t), REAL(f(t)))
      END DO
      CALL GWline2(IR_, REAL(pi), REAL(f(pi)))
      END
C
      FUNCTION fn(xx)
      IMPLICIT NONE
      REAL*8 fn, x, xx, pi, pi2, eps
      PARAMETER(pi = 3.141592653589793D0, pi2 = 2*pi, eps = 1.0D-13)
      x = xx
      IF(x .LT. 0.0D0) x = x + pi2
      IF(ABS(x) .LT. eps) THEN
         fn = 0.0D0
      ELSE IF(x .LT. pi - eps) THEN
         fn = 1.0D0
      ELSE IF(x .GT. pi + eps) THEN
         fn = -1.0D0
      ELSE
         fn = 0.0D0
      END IF
      END
C
      FUNCTION gr(xx)
      IMPLICIT NONE
      REAL*8 g, gr, x, xx, kx, pi, pi2
      PARAMETER(pi = 3.141592653589793D0, pi2 = 2*pi)
      INTEGER N, k
      COMMON /ARRY/g(2,1024), N
      x = xx
      IF(x .LT. 0.0D0) x = x + pi2
      gr = 0
      kx = -N/2*x
      IF(kx .GT. 0.0D0) THEN
         kx = mod(kx, pi2)
      ELSE IF(kx .LT. 0.0D0) THEN
         kx = -mod(-kx, pi2) + pi2
      ENDIF
      DO k = -N/2, N/2-1
         IF(kx .GT. pi2) kx = kx - pi2
         gr = gr + g(1, mod(k+N,N)+1)*cos(kx)
     +           - g(2, mod(k+N,N)+1)*sin(kx)
         kx = kx + x
      END DO
      gr = gr/N
      END
C
      FUNCTION gi(xx)
      IMPLICIT NONE
      REAL*8 g, gi, x, xx, kx, pi, pi2
      PARAMETER(pi = 3.141592653589793D0, pi2 = 2*pi)
      INTEGER N, k
      COMMON /ARRY/g(2,1024), N
      x = xx
      IF(x .LT. 0.0D0) x = x + pi2
      gi = 0
      kx = -N/2*x
      IF(kx .GT. 0.0D0) THEN
         kx = mod(kx, pi2)
      ELSE IF(kx .LT. 0.0D0) THEN
         kx = -mod(-kx, pi2) + pi2
      ENDIF
      DO k = -N/2, N/2-1
         IF(kx .GT. pi2) kx = kx - pi2
         gi = gi + g(1, mod(k+N,N)+1)*sin(kx)
     +           + g(2, mod(k+N,N)+1)*cos(kx)
         kx = kx + x
      END DO
      gi = gi/N
      END