      PROGRAM FFT
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C  Fast Fourier Transform Program.
C  Written to demonstate the principles of the FFT Method.
C  A data size is defined and a set of Cosine waves of different
C  frequency are defined and generated. 
C  These are added and multiplied by a Gaussian function which has
C  Value 1 at the start and <1 at the end.
C  A FFT is carried out and the result normalised and written to
C  a data file suitable for plotting using  Cricket Graph.
C
C  Henry Rzepa.  September, 1990
C   Modified for SGI Unix Oct. 1993
C
C<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
C       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION XREAL(341024), XIMAG(341024), BROAD(341024), FID(341024),TIT(5)
C  Cricket Graph accepts columns separated by TAB characters.
      CHARACTER  TAB*1, TIT*14
      TAB= CHAR(9)
      TIT(1) = '     Time/Freq'
      TIT(2) = 'Decay Function'
      TIT(3) = '     FID Decay'
      TIT(4) = '     Real Part'
      TIT(5) = 'Imaginary Part'
C	Open files for CricketGraph
      OPEN(3,FILE='cricket.dat',STATUS='UNKNOWN')
 77   WRITE (*,*) 'ENTER No. Points in FFT (2**n) and time delay (CTRL+d to exit)'
      READ(5,*,END=299,ERR=299) N,T0
      IF(N.LT.64 .OR. N.GT.32384) THEN
      GOTO 77
      ENDIF
      	IF (N.EQ.64 .OR. N.EQ. 128 .OR. N.EQ. 256 .OR.
     1 	N.EQ.512 .OR. N.EQ. 1024 .or. N.eq. 32384) THEN
      	DO 2 JK = 1,N
      	XREAL(JK) = 0.0
      	BROAD(JK) = 0.0
      	FID(JK) = 0.0
 2    	XIMAG(JK) = 0.0
 	NP=1
      	ELSE
      	WRITE(*,*) 'N must be a power of 2'
      	GOTO 77
      	ENDIF
 3    Write (*,*) 'Enter frequency and intensity of line (CTRL+d to exit
     1)'
      READ(5,*,END=6,ERR=6) K,A
      	IF (K.GT.N/2) THEN
      	WRITE (*,72) K, N
 72   	FORMAT (' The Nyquist Theorem indicates that'/,I3,
     1' cannot be greater than', I4,'/2')
      	ENDIF
      	IF (K.NE.0) THEN
      	ANP = FLOAT(N)/(2.0*FLOAT(K))
      	ELSE
      	ANP = 1.0
      	ENDIF
      STEP = 3.14159/ANP
      DO 300 KK = 1,N
      PHI = (KK-1-T0)*STEP
      	IF (K.NE.0) THEN
      	FID(KK) = FID(KK) + A*COS(PHI)
      	ELSE
      	FID(KK) = 0.1
      	ENDIF
 300  CONTINUE
      GOTO 3
C Normalise the FID to unity.
  6   AMAX = 0.0
c      DO 301 KK=1,N
c 301  IF (FID(KK).GT. AMAX) AMAX = FID(KK)
c      DO 302 KK=1,N
c 302  FID(KK) = FID(KK)/AMAX
      WRITE (6,*) 'Enter broadening factor (0 for none)'
      BRF = 0.1
      READ(5,*,END=7,ERR=7) BRF
      	IF (BRF.LT. 1.E-7) THEN
      	BRF  = 0.01
      	ENDIF
      	print *, 'value of BRF', BRF
 7    DO 4 JJJ = 1,N
      PART = FLOAT(JJJ-1)/FLOAT(N)
      BROAD(JJJ) = EXP(-BRF*PART)
      FID(JJJ) = FID(JJJ)*BROAD(JJJ)
 4    XREAL(JJJ) = FID(JJJ)
      write (*,*) 'Performing FFT'   
      CALL FASTF(XREAL,XIMAG,N)
      write (*,*) 'Finished FFT'   
C Normalise the Spectrum to unity
      AMAXR = 0.0
      AMAXI = 0.0
      DO 303 KK=1,N
      IF (XIMAG(KK).GT. AMAXI) AMAXI = XIMAG(KK)
 303  IF (XREAL(KK).GT. AMAXR) AMAXR = XREAL(KK)
      DO 304 KK=1,N
      XIMAG(KK) = XIMAG(KK)/AMAXI
 304  XREAL(KK) = XREAL(KK)/AMAXR
      WRITE(3,229) (TIT(J),TAB, J=1,5)
      WRITE(*,229) (TIT(J),TAB, J=1,5)
  229 FORMAT('*'/5(A,A) )
  228 FORMAT(5(A,1X))
      DO 16 JJJ=1,N
C Real part contains only N/2 data points.
      	IF (JJJ.LE. N/2) THEN
C Frequency of 0 is stored in data point 1. Hence subtract 1 from
C each point to get actual frequency.
      	WRITE (3,27) JJJ-1,TAB,BROAD(JJJ),TAB,FID(JJJ),TAB,XREAL(JJJ), 
     1	TAB, XIMAG(JJJ)
      	WRITE (*,27) JJJ-1,TAB,BROAD(JJJ),TAB,FID(JJJ),TAB,XREAL(JJJ), 
     1	TAB, XIMAG(JJJ)
      	ELSE
      	WRITE (3,28) JJJ, TAB, BROAD(JJJ), TAB, FID(JJJ)
      	ENDIF
 16   CONTINUE
 28   FORMAT(I4,A,F14.4,A,F14.4)
 27   FORMAT(I4,A,F14.4,A,F14.4,A,F14.4,A,F14.4)
 26   FORMAT(I14,1X,4(F14.4,1X))
      WRITE(*,*) 'File cricket.dat written. Transfer to Macintosh'
      GOTO 77
 299  STOP
      END
