c--------------------------------------------------------
c $Id$
c Description: testfft
c Author:  Murray S. Daw <daw@clemson.edu>
c Created: 22 Jan 03
c--------------------------------------------------------
c  this is to test (& demonstrate) a simple fft program
c  obtained from netlib (www.netlib.org)
c  it reads in a (real) array from standard input,
c      does the FT, writes it to standard output
c
c  note that, to use the canned routine, i have
c  inserted "real*8" and "complex*16"
c  when importing canned routines, it is important to
c  verify that the precision of the arguments is compatible
c  between the caller and callee!
c
c  also notice that the length of the array must be a power
c  of 2 for this FFT routine to work      
      
      program testfft

      implicit real*8 (a-h,o-z)   ! see note above
      parameter (nmax=6000)
      dimension x(nmax)
      complex*16 cx(nmax)   ! fft routine wants a complex array
      dimension power(nmax)  ! power spectrum

c     set up the data
      n = 0
 10   read(5,*,end=1000)xin
      n = n + 1
      if ( n .le. nmax ) x(n) = xin
      goto 10
 1000 continue
      if ( n .gt. nmax ) then
         print *,' whoops! '
         print *,' n = ',n,' should be <= nmax = ',nmax
         print *,' you need to increase the allocation '
         stop
      endif
c      print *,' array has ',n,' elements '
c     n must be a power of 2 for this routine to work
      m = 0
      nmod = n

      do while ( nmod .gt. 1 )
         modn = mod(nmod,2)
         if ( modn .ne. 0 ) then
            print *,' whoops! '
            print *,' n = ',n,' must be a power of 2! '
            stop
         else
            nmod = nmod/2
            m = m + 1
         endif
      enddo
      if ( n .ne. 2**m ) then
         print *,' darn it ! '
         stop
      endif

      

c     do the fft
c     the fft routine expects a COMPLEX array, so we
c     copy the real array to a complex one
      do i=1,n
         cx(i) = cmplx(x(i),0.)
      enddo
      k = 1  ! this requests a forward transform
      call fft(cx,n,k)
c     now the power spectrum is complex
c     take the abs & copy back to real array
      do i=1,n
         power(i) = abs(cx(i))
      enddo

c     now write
      do i=1,n
         write(6,9000) power(i)
      enddo
 9000 format(1x,g20.10)

      stop
      end

c     
c  here i paste in the fft from netlib (with a minor change noted)
c      
c=================================================================
c  It runs and computes correctly although it is not in any 
c  sense sophisticated. However, it is free, and useful if 
c  you are in a hurry and don't want to think things out for yourself.
c   
c  Arthur Wouk (wouk@brl-vgr)
c  
c  *******************
c  Fast Fourier Transform
c  *******************
c  
	SUBROUTINE FFT (X, N, K)
C	FFT COMPUTES THE (FAST) FOURIER TRANSFORM OF THE VECTOR X
C	(A COMPLEX ARRAY OF DIMENSION N). SOURCE: Ferziger; Numerical
C       methods for engineering applications.
C
C	X = DATA TO BE TRANSFORMED; ON RETURN IT CONTAINS THE TRANSFORM.
C	N = SIZE OF VECTOR. MUST BE A POWER OF 2 (<32769).
C	K = 1 FOR FORWARD TRANSFORM.
C	K = -1 FOR INVERSE TRANSFORM.
C
	IMPLICIT INTEGER (A-Z)
	INTEGER SBY2,S

c     msd: i have changed these to real*8 and complex*16
c	REAL GAIN, PI2, ANG, RE, IM
c	COMPLEX X(N), XTEMP, T, U(16), V, W
	REAL*8 GAIN, PI2, ANG, RE, IM
	COMPLEX*16 X(N), XTEMP, T, U(16), V, W


        
	LOGICAL NEW
	DATA	PI2,GAIN,NO,KO /6.283185307, 1., 0, 0/
C
C	TEST FIRST CALL?
C
	NEW = ( NO .NE. N)
	IF ( .NOT. NEW) GO TO 2
C
C	IF FIRST CALL COMPUTE LOG2 (N).
C
	L2N = 0
	NO = 1
    1	L2N = L2N + 1
	NO = NO + NO
	IF (NO .LT. N) GO TO 1
	GAIN = 1./N
	ANG = PI2*GAIN
	RE = COS (ANG)
	IM = SIN (ANG)
C
C	COMPUTE COMPLEX EXPONENTIALS IF NOT FIRST CALL
C
    2	IF (.NOT. NEW .AND. K*KO .GE. 1) GO TO 4
	U(1) = CMPLX (RE, -SIGN(IM, FLOAT(K)))
	DO 3 I = 2,L2N
  	   U(I) = U(I-1)*U(I-1)
    3	CONTINUE
	K0 = K
C
C	MAIN LOOP
C
    4	SBY2 = N
	DO 7 STAGE = 1,L2N
  	   V = U(STAGE)
  	   W = (1., 0.)
  	   S = SBY2
  	   SBY2 = S/2
  	   DO 6 L = 1,SBY2
		DO 5 I = 1,N,S
		   P = I + L- 1
		   Q = P + SBY2
		   T = X(P) + X(Q)
		   X(Q) = ( X(P) - X(Q))*W
		   X(P) =T
    5 		CONTINUE
		W = W*V
    6	   CONTINUE
    7	CONTINUE
C
C	REORDER THE ELEMENTS BY BIT REVERSAL
C
	DO 9 I = 1,N
	   INDEX = I-1
	   JNDEX = 0
	   DO 8 J = 1,L2N
		JNDEX = JNDEX+JNDEX
		ITEMP = INDEX/2
		IF (ITEMP+ITEMP .NE. INDEX) JNDEX = JNDEX + 1
		INDEX = ITEMP
    8	   CONTINUE
	   J = JNDEX + 1
	   IF (J .LT. I) GO TO 9
	   XTEMP = X(J)
	   X(J) = X(I)
	   X(I) = XTEMP
    9	CONTINUE
C
C	FORWARD TRANSFORM DONE
C
	IF (K .GT. 0) RETURN
C
C	INVERSE TRANSFORM
C
	DO 10 I = 1,N
	   X(I) = X(I)*GAIN
   10  	CONTINUE
	RETURN
	END