$STORAGE: 2
$NOFLOATCALLS
C  
C-----------------------------------------------------------------------
C SUBROUTINE: FOUREA
C PERFORMS COOLEY-TUKEY FAST FOURIER TRANSFORM
C-----------------------------------------------------------------------
C  
      SUBROUTINE FOUREA(ID1,ID2,IC,IR,IVN,ISI)
C ID1,ID2 = COORDS OF FIRST CELL. IC AND IR ARE 0, OR 1
C ONLY ONE OF IC, IR MAY BE NONZERO. (FLAGS HORIZ/VERTICAL
C DATA AREA)
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 ASSUMES THAT FIRST N/2 ELEMENTS ARE REAL, SECOND COMPLEX...
C STORES DATA THAT WAY ALSO...
C
C      COMPLEX DATA(1)
C      COMPLEX TEMP, W
C MAKE THIS A REAL FFT, NOT COMPLEX...
	REAL*8 DATA(1),TEMP,W,TEMP2,TEMPI,WI
	INTEGER*2 ID1,ID2,IC,IR,IRX,IRXX,IVN,N
C SET UP STMT FUNCTIONS...
	ID1F(K)=ID1+IC*(K-1)
	ID2F(K)=ID2+IR*(K-1)
	N=IVN
C  
C CHECK FOR POWER OF TWO UP TO 14
C  
C INITIALLY SAY ALL OK
      NN = 1
      DO 10 I=1,14
        M = I
        NN = NN*2
        IF (NN.EQ.N) GO TO 20
	IF(NN.GT.N)GOTO 11
  10  CONTINUE
11	CONTINUE
	N=NN/2
C USE NEXT SMALLER POWER OF 2 ARRAY...
C	RETURN
C HERE BEGINNETH ACTUAL WORK.
C SET UP DATA COORDS ON THE FLY. NORMALLY I,J RUN IN RANGE 1 TO N
C SO WHERE K=(I OR J) (I.E., ONE OF THE TWO) WE USE A RELATION
C ID1V=ID1+IC*(K-1) AND ID2V=ID2+IR*(K-1). WE USE STMT FUNCTIONS
C ID1F AND ID2F FOR THIS.
  20  CONTINUE
	NOV2=N/2
C  
C      PI = 4.*ATAN(1.)
	PI=3.14159265358979323846264
      FN = NOV2
C  
C THIS SECTION PUTS DATA IN BIT-REVERSED ORDER
C  
      J = 1
      DO 80 I=1,NOV2
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	CONTINUE
C EXCHANGE DATA(J), DATA(I)
	CALL XVBLGT(ID1F(J),ID2F(J),TEMP)
	CALL XVBLGT(ID1F(I),ID2F(I),DATA(1))
	CALL XVBLST(ID1F(J),ID2F(J),DATA(1))
	CALL XVBLST(ID1F(I),ID2F(I),TEMP)
C FLIP BOTH REAL AND COMPLEX PARTS OF DATA
	CALL XVBLGT(ID1F(J+NOV2),ID2F(J+NOV2),TEMP)
	CALL XVBLGT(ID1F(I+NOV2),ID2F(I+NOV2),DATA(1))
	CALL XVBLST(ID1F(J+NOV2),ID2F(J+NOV2),DATA(1))
	CALL XVBLST(ID1F(I+NOV2),ID2F(I+NOV2),TEMP)
C  30    TEMP = DATA(J)
C        DATA(J) = DATA(I)
C        DATA(I) = TEMP
C  
C IMPLEMENT J=J+1, BIT-REVERSED COUNTER
C  
  40    M = NOV2/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.NOV2)GOTO 130
 100  ISTEP = 2*MMAX
      DO 120 M=1,MMAX
        THETA = PI*FLOAT(ISI*(M-1))/FLOAT(MMAX)
	 W = COS(THETA)
        WI = SIN(THETA)
C        W = CMPLX(COS(THETA),SIN(THETA))
        DO 110 I=M,NOV2,ISTEP
          J = I + MMAX
C GET REAL AND IMAG HALVES OF NUMBER...
	  CALL XVBLGT(ID1F(J),ID2F(J),TEMP)
	  CALL XVBLGT(ID1F(J+NOV2),ID2F(J+NOV2),TEMPI)
C DO COMPLEX MULTIPLICATION BY HAND TO AVOID LARGE RUNTIME SYSTEM
C ROUTINE INCLUSION.
	  TEMP2=W*TEMP-WI*TEMPI
	  TEMPI=WI*TEMP+W*TEMPI
	TEMP=TEMP2
C          TEMP = W*DATA(J)
C          DATA(J) = DATA(I) - TEMP
C          DATA(I) = DATA(I) + TEMP
	   CALL XVBLGT(ID1F(I),ID2F(I),DATA(1))
	   TEMP2=DATA(1)+TEMP
	   DATA(1)=DATA(1) - TEMP
	   CALL XVBLST(ID1F(J),ID2F(J),DATA(1))
	   CALL XVBLST(ID1F(I),ID2F(I),TEMP2)
C COMPLEX PART
	   CALL XVBLGT(ID1F(I+NOV2),ID2F(I+NOV2),DATA(1))
	   TEMP2=DATA(1)+TEMPI
	   DATA(1)=DATA(1) - TEMPI
	   CALL XVBLST(ID1F(J+NOV2),ID2F(J+NOV2),DATA(1))
	   CALL XVBLST(ID1F(I+NOV2),ID2F(I+NOV2),TEMP2)
 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
C        DATA(I) = DATA(I)/FN
	CALL XVBLGT(ID1F(I),ID2F(I),TEMP)
	TEMP=TEMP/FN
	CALL XVBLST(ID1F(I),ID2F(I),TEMP)
 150  CONTINUE
 160  RETURN
      END


