c
c     Second Edition
c     Numerical Analysis:
c     The Mathematics of Scientific Computing
c     D.R. Kincaid & E.W. Cheney
c     Brooks/Cole Publ., 1996
c     ISBN 0-534-33892-5
c     COPYRIGHT (c) 1996
c
c     Section 6.12
c
c     Example of Fast Fourier Transform method
c
c
c     file: fft.f
c
      parameter (m=4,n3=15)
      complex z(0:n3),c(0:n3),d(0:n3),w,f,u,v
c
      print *
      print *,' Fast Fourier Transform example'
      print *,' Section 6.12, Kincaid-Cheney'
      print *
c
      nn=2**m
      m2=nn/2
      pi=4.0*atan(1.0)
      y=2.0*pi/float(nn)
      w=cmplx(0.0 , -y)
      w=cexp(w)
      z(0)=cmplx(1.0,0.0)
      c(0)=f(0.0)
      do 2 k=1,nn-1
         z(k)=w*z(k-1)
         c(k)=f(y*float(k))
 2    continue
      print *,' Initial c-vector'
      do 4 n=0,m-1
         print *,' n = ',n
         print *,c
         print *
         mn=2**(m-n-1)
         n2=2**n
         do 3 k=0,mn-1
            do 3 j=0,n2-1
               u=c(n2*k+j)
               v=z(j*mn)*c(n2*k+j+m2)
               d(2*n2*k+j)=(u+v)*cmplx(0.5,0.0)
               d(2*n2*k+j+n2)=(u-v)*cmplx(0.5,0.0)
 3          continue
            do 4 j=0,nn-1
               c(j)=d(j)
 4          continue
            print *,' Final c-vector'
            print *,c
c
c 5    format (2x,16('(',e13.6,',',e13.6,')',4x))
      stop
      end
c
      function f(x)
      complex f,z,d
      u=cos(x)
      v=sin(x)
      z=cmplx(u,v)
      d=cmplx(16.0,0.0)
      do 2 k=1,15
         y=float(k)
         d=d*z + cmplx(16.0-y,0.0)
 2    continue
      f=d
c
      return
      end
