      SUBROUTINE COOL (DATA,NPREV,N,NREM,ISIGN,IFACT,WORK)
C     FOURIER TRANSFORM OF LENGTH N.  IN PLACE COOLEY-TUKEY METHOD, 
C     DIGIT-REVERSED TO NORMAL ORDER, SANDE-TUKEY FACTORING (2).
C     DIMENSION DATA(NPREV,N,NREM)
C     COMPLEX DATA
C     DATA(I1,J2,I3) = SUM(DATA(I1,I2,I3)*EXP(ISIGN*2*PI*I*((I2-1)* 
C     (J2-1)/N))), SUMMED OVER I2 = 1 TO N FOR ALL I1 FROM 1 TO NPREV,
C     J2 FROM 1 TO N AND I3 FROM 1 TO NREM.  THE FACTORS OF N ARE GIVEN 
C     IN ANY ORDER IN ARRAY IFACT.  FACTORS OF TWO ARE DONE IN PAIRS
C     AS MUCH AS POSSIBLE (FOURIER TRANSFORM OF LENGTH FOUR), FACTORS OF
C     THREE ARE DONE SEPARATELY, AND ALL FACTORS FIVE OR HIGHER 
C     ARE DONE BY GOERTZEL'S ALGORITHM (4). 
      IMPLICIT REAL (A-H,O-Z),
     +INTEGER (I-N)
      DIMENSION DATA(1), WORK(1), IFACT(1)
      TWOPI=6.283185307E0*FLOAT(ISIGN)
      IP0=2 
      IP1=IP0*NPREV 
      IP4=IP1*N 
      IP5=IP4*NREM
      IF=0
      IP2=IP1 
 10   IF (IP2-IP4) 20,240,240 
 20   IF=IF+1 
      IFCUR=IFACT(IF) 
      IF (IFCUR-2) 60,30,60 
 30   IF (4*IP2-IP4) 40,40,60 
 40   IF (IFACT(IF+1)-2) 60,50,60 
 50   IF=IF+1 
      IFCUR=4 
 60   IP3=IP2*IFCUR 
      THETA=TWOPI/FLOAT(IFCUR)
      SINTH=SIN(THETA/2.E0) 
      ROOTR=-2.E0*SINTH*SINTH 
C     COS(THETA)-1, FOR ACCURACY. 
      ROOTI=SIN(THETA)
      THETA=TWOPI/FLOAT(IP3/IP1)
      SINTH=SIN(THETA/2.E0) 
      WSTPR=-2.E0*SINTH*SINTH 
      WSTPI=SIN(THETA)
      WR=1.E0 
      WI=0.E0 
      DO 230 I2=1,IP2,IP1 
      IF (IFCUR-4) 70,70,210
 70   IF ((I2-1)*(IFCUR-2)) 240,90,80 
 80   W2R=WR*WR-WI*WI 
      W2I=2.E0*WR*WI
      W3R=W2R*WR-W2I*WI 
      W3I=W2R*WI+W2I*WR 
 90   I1MAX=I2+IP1-IP0
      DO 200 I1=I2,I1MAX,IP0
      DO 200 I5=I1,IP5,IP3
      J0=I5 
      J1=J0+IP2 
      J2=J1+IP2 
      J3=J2+IP2 
      IF (I2-1) 140,140,100 
 100  IF (IFCUR-3) 130,120,110
C     APPLY THE PHASE SHIFT FACTORS 
 110  TEMPR=DATA(J3)
      DATA(J3)=W3R*TEMPR-W3I*DATA(J3+1) 
      DATA(J3+1)=W3R*DATA(J3+1)+W3I*TEMPR 
      TEMPR=DATA(J2)
      DATA(J2)=WR*TEMPR-WI*DATA(J2+1) 
      DATA(J2+1)=WR*DATA(J2+1)+WI*TEMPR 
      TEMPR=DATA(J1)
      DATA(J1)=W2R*TEMPR-W2I*DATA(J1+1) 
      DATA(J1+1)=W2R*DATA(J1+1)+W2I*TEMPR 
      GO TO 140 
 120  TEMPR=DATA(J2)
      DATA(J2)=W2R*TEMPR-W2I*DATA(J2+1) 
      DATA(J2+1)=W2R*DATA(J2+1)+W2I*TEMPR 
 130  TEMPR=DATA(J1)
      DATA(J1)=WR*TEMPR-WI*DATA(J1+1) 
      DATA(J1+1)=WR*DATA(J1+1)+WI*TEMPR 
 140  IF (IFCUR-3) 150,160,170
C     DO A FOURIER TRANSFORM OF LENGTH TWO
 150  TEMPR=DATA(J1)
      TEMPI=DATA(J1+1)
      DATA(J1)=DATA(J0)-TEMPR 
      DATA(J1+1)=DATA(J0+1)-TEMPI 
      DATA(J0)=DATA(J0)+TEMPR 
      DATA(J0+1)=DATA(J0+1)+TEMPI 
      GO TO 200 
C     DO A FOURIER TRANSFORM OF LENGTH THREE
 160  SUMR=DATA(J1)+DATA(J2)
      SUMI=DATA(J1+1)+DATA(J2+1)
      TEMPR=DATA(J0)-0.5E0*SUMR
      TEMPI=DATA(J0+1)-0.5E0*SUMI
      DATA(J0)=DATA(J0)+SUMR
      DATA(J0+1)=DATA(J0+1)+SUMI
      DIFR=ROOTI*(DATA(J2+1)-DATA(J1+1))
      DIFI=ROOTI*(DATA(J1)-DATA(J2))
      DATA(J1)=TEMPR+DIFR 
      DATA(J1+1)=TEMPI+DIFI 
      DATA(J2)=TEMPR-DIFR 
      DATA(J2+1)=TEMPI-DIFI 
      GO TO 200 
C     DO A FOURIER TRANSFORM OF LENGTH FOUR (FROM BIT REVERSED ORDER) 
 170  T0R=DATA(J0)+DATA(J1) 
      T0I=DATA(J0+1)+DATA(J1+1) 
      T1R=DATA(J0)-DATA(J1) 
      T1I=DATA(J0+1)-DATA(J1+1) 
      T2R=DATA(J2)+DATA(J3) 
      T2I=DATA(J2+1)+DATA(J3+1) 
      T3R=DATA(J2)-DATA(J3) 
      T3I=DATA(J2+1)-DATA(J3+1) 
      DATA(J0)=T0R+T2R
      DATA(J0+1)=T0I+T2I
      DATA(J2)=T0R-T2R
      DATA(J2+1)=T0I-T2I
      IF (ISIGN) 180,180,190
 180  T3R=-T3R
      T3I=-T3I
 190  DATA(J1)=T1R-T3I
      DATA(J1+1)=T1I+T3R
      DATA(J3)=T1R+T3I
      DATA(J3+1)=T1I-T3R
 200  CONTINUE
      GO TO 220 
C     DO A FOURIER TRANSFORM OF LENGTH FIVE OR MORE 
 210  CALL GOERT (DATA(I2),NPREV,IP2/IP1,IFCUR,IP5/IP3,WORK,WR,WI,ROOTR,
     XROOTI)
 220  TEMPR=WR
      WR=WSTPR*TEMPR-WSTPI*WI+TEMPR 
 230  WI=WSTPR*WI+WSTPI*TEMPR+WI
      IP2=IP3 
      GO TO 10
 240  RETURN
      END 
