C  
C-----------------------------------------------------------------------
C SUBROUTINE: FOUREA
C PERFORMS COOLEY-TUKEY FAST FOURIER TRANSFORM
C-----------------------------------------------------------------------
C  
      SUBROUTINE FOUREA(ID1,ID2,IC,IR,N,ISI,IRSLT)
C  
C THE COOLEY-TUKEY FAST FOURIER TRANSFORM IN ANSI FORTRAN
C  
C DATA IS A ONE-DIMENSIONAL COMPLEX ARRAY WHOSE LENGTH, N, IS A
C POWER OF TWO.  ISI IS +1 FOR AN INVERSE TRANSFORM AND -1 FOR A
C FORWARD TRANSFORM.  TRANSFORM VALUES ARE RETURNED IN THE INPUT
C ARRAY, REPLACING THE INPUT.
C TRANSFORM(J)=SUM(DATA(I)*W**((I-1)*(J-1))), WHERE I AND J RUN
C FROM 1 TO N AND W = EXP (ISI*2*PI*SQRT(-1)/N).  PROGRAM ALSO
C COMPUTES INVERSE TRANSFORM, FOR WHICH THE DEFINING EXPRESSION
C IS INVTR(J)=(1/N)*SUM(DATA(I)*W**((I-1)*(J-1))).
C RUNNING TIME IS PROPORTIONAL TO N*LOG2(N), RATHER THAN TO THE
C CLASSICAL N**2.
C AFTER PROGRAM BY BRENNER, JUNE 1967. THIS IS A VERY SHORT VERSION
C OF THE FFT AND IS INTENDED MAINLY FOR DEMONSTRATION. PROGRAMS
C ARE AVAILABLE IN THIS COLLECTION WHICH RUN FASTER AND ARE NOT
C RESTRICTED TO POWERS OF 2 OR TO ONE-DIMENSIONAL ARRAYS.
C SEE -- IEEE TRANS AUDIO (JUNE 1967), SPECIAL ISSUE ON FFT.
C  
C      COMPLEX DATA(1)
C      COMPLEX TEMP, W
C MAKE THIS A REAL FFT, NOT COMPLEX...
	REAL*8 DATA(1),TEMP,W
	INTEGER*2 ID1,ID2,IC,IR
C  
C CHECK FOR POWER OF TWO UP TO 15
C  
	IRSLT=1
C INITIALLY SAY ALL OK
      NN = 1
      DO 10 I=1,15
        M = I
        NN = NN*2
        IF (NN.EQ.N) GO TO 20
  10  CONTINUE
	IRSLT=3
C 3 = BAD ARGS
	RETURN
  20  CONTINUE
C  
      PI = 4.*ATAN(1.)
      FN = N
C  
C THIS SECTION PUTS DATA IN BIT-REVERSED ORDER
C  
      J = 1
      DO 80 I=1,N
C  
C AT THIS POINT, I AND J ARE A BIT REVERSED PAIR (EXCEPT FOR THE
C DISPLACEMENT OF +1)
C  
	IF(I.GE.J)GOTO 40
C  
C EXCHANGE DATA(I) WITH DATA(J) IF I.LT.J.
C  
  30    TEMP = DATA(J)
        DATA(J) = DATA(I)
        DATA(I) = TEMP
C  
C IMPLEMENT J=J+1, BIT-REVERSED COUNTER
C  
  40    M = N/2
  50    IF (J.LE.M) GOTO 70
  60    J = J - M
        M = (M+1)/2
        GO TO 50
  70    J = J + M
  80  CONTINUE
C  
C NOW COMPUTE THE BUTTERFLIES
C  
      MMAX = 1
  90  IF (MMAX.GE.N)GOTO 130
 100  ISTEP = 2*MMAX
      DO 120 M=1,MMAX
        THETA = PI*FLOAT(ISI*(M-1))/FLOAT(MMAX)
	 W = COS(THETA)
C        W = CMPLX(COS(THETA),SIN(THETA))
        DO 110 I=M,N,ISTEP
          J = I + MMAX
          TEMP = W*DATA(J)
          DATA(J) = DATA(I) - TEMP
          DATA(I) = DATA(I) + TEMP
 110    CONTINUE
 120  CONTINUE
      MMAX = ISTEP
      GO TO 90
  130  IF (ISI.LT.0) GOTO 160
C  
C FOR INV TRANS -- ISI=1 -- MULTIPLY OUTPUT BY 1/N
C  
 140  DO 150 I=1,N
        DATA(I) = DATA(I)/FN
 150  CONTINUE
 160  RETURN
      END