c===========================================================
c     Test routine for double precision FFT routines.
c===========================================================
      program     tfft
     
      implicit    none

      integer     n
      parameter ( n = 64 )

c-----------------------------------------------------------
c     Minimum length for work array is '2*n+15'
c-----------------------------------------------------------
      real*8      rsave(n),   rfft(n),   w(2*n+15)
      integer     i

c-----------------------------------------------------------
c     Initialize FFT software for specific 'n' (64).
c-----------------------------------------------------------
      call drffti(n,w)

c-----------------------------------------------------------
c     Set up data to be transformed (pure sine wave)
c-----------------------------------------------------------
      do i = 1 , n
         rsave(i) = sin(8.0d0*atan(1.0d0)*(i-1)/n)
         rfft(i)  = rsave(i)
      end do

c-----------------------------------------------------------
c     Compute transform and output original data and FFT
c-----------------------------------------------------------
      call drfftf(n,rfft,w)
      do i = 1 , n
         write(0,*) rsave(i),  rfft(i)
      end do

      write(0,*)

c-----------------------------------------------------------
c     Compute inverse transform of transformed data and 
c     output it (divided by n) and original data.  Note 
c     that 'drfftf' followed by 'drfftb' will generically
c     scale original data by n.
c-----------------------------------------------------------
      call drfftb(n,rfft,w)
      do i = 1 , n
         write(0,*) rsave(i),  rfft(i) / n
      end do

      stop
      end