      SUBROUTINE FOUR2 (DATA,N,NDIM,ISIGN,IFORM)
C     COOLEY-TUKEY FAST FOURIER TRANSFORM IN USASI BASIC FORTRAN. 
C     MULTI-DIMENSIONAL TRANSFORM, EACH DIMENSION A POWER OF TWO, 
C     COMPLEX OR REAL DATA. 
C     TRANSFORM(K1,K2,...) = SUM(DATA(J1,J2,...)*EXP(ISIGN*2*PI*SQRT(-1)
C     *((J1-1)*(K1-1)/N(1)+(J2-1)*(K2-1)/N(2)+...))), SUMMED FOR ALL
C     J1 AND K1 FROM 1 TO N(1), J2 AND K2 FROM 1 TO N(2), 
C     ETC. FOR ALL NDIM SUBSCRIPTS.  NDIM MUST BE POSITIVE AND
C     EACH N(IDIM) MUST BE A POWER OF TWO.  ISIGN IS +1 OR -1.
C     LET NTOT = N(1)*N(2)*...*N(NDIM).  THEN A -1 TRANSFORM
C     FOLLOWED BY A +1 ONE (OR VICE VERSA) RETURNS NTOT 
C     TIMES THE ORIGINAL DATA.  IFORM = 1, 0 OR -1, AS DATA IS
C     COMPLEX, REAL OR THE FIRST HALF OF A COMPLEX ARRAY.  TRANSFORM
C     VALUES ARE RETURNED TO ARRAY DATA.  THEY ARE COMPLEX, REAL OR 
C     THE FIRST HALF OF A COMPLEX ARRAY, AS IFORM = 1, -1 OR 0. 
C     THE TRANSFORM OF A REAL ARRAY (IFORM = 0) DIMENSIONED N(1) BY N(2)
C     BY ... WILL BE RETURNED IN THE SAME ARRAY, NOW CONSIDERED TO
C     BE COMPLEX OF DIMENSIONS N(1)/2+1 BY N(2) BY ....  NOTE THAT IF 
C     IFORM = 0 OR -1, N(1) MUST BE EVEN, AND ENOUGH ROOM MUST BE 
C     RESERVED.  THE MISSING VALUES MAY BE OBTAINED BY COMPLEX CONJUGA- 
C     TION.  THE REVERSE TRANSFORMATION, OF A HALF COMPLEX ARRAY DIMEN- 
C     SIONED N(1)/2+1 BY N(2) BY ..., IS ACCOMPLISHED BY SETTING IFORM
C     TO -1.  IN THE N ARRAY, N(1) MUST BE THE TRUE N(1), NOT N(1)/2+1. 
C     THE TRANSFORM WILL BE REAL AND RETURNED TO THE INPUT ARRAY. 
C     RUNNING TIME IS PROPORTIONAL TO NTOT*LOG2(NTOT), RATHER THAN
C     THE NAIVE NTOT**2.  FURTHERMORE, LESS ERROR IS BUILT UP.
C     WRITTEN BY NORMAN BRENNER OF MIT LINCOLN LABORATORY, JANUARY 1969.
C     SEE-- IEEE AUDIO TRANSACTIONS (JUNE 1967), SPECIAL ISSUE ON FFT.
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
      DIMENSION DATA(1), N(1) 
      NTOT=1
      DO 10 IDIM=1,NDIM 
 10   NTOT=NTOT*N(IDIM) 
      IF (IFORM) 70,20,20 
 20   NREM=NTOT 
      DO 60 IDIM=1,NDIM 
      NREM=NREM/N(IDIM) 
      NPREV=NTOT/(N(IDIM)*NREM) 
      NCURR=N(IDIM) 
      IF (IDIM-1+IFORM) 30,30,40
 30   NCURR=NCURR/2 
 40   CALL BITRV (DATA,NPREV,NCURR,NREM)
      CALL COOL2 (DATA,NPREV,NCURR,NREM,ISIGN)
      IF (IDIM-1+IFORM) 50,50,60
 50   CALL FIXRL (DATA,N(1),NREM,ISIGN,IFORM) 
      NTOT=(NTOT/N(1))*(N(1)/2+1) 
 60   CONTINUE
      RETURN
 70   NTOT=(NTOT/N(1))*(N(1)/2+1) 
      NREM=1
      DO 100 JDIM=1,NDIM
      IDIM=NDIM+1-JDIM
      NCURR=N(IDIM) 
      IF (IDIM-1) 80,80,90
 80   NCURR=NCURR/2 
      CALL FIXRL (DATA,N(1),NREM,ISIGN,IFORM) 
      NTOT=NTOT/(N(1)/2+1)*N(1) 
 90   NPREV=NTOT/(N(IDIM)*NREM) 
      CALL BITRV (DATA,NPREV,NCURR,NREM)
      CALL COOL2 (DATA,NPREV,NCURR,NREM,ISIGN)
 100  NREM=NREM*N(IDIM) 
      RETURN
      END 
