      SUBROUTINE FOURT2(DATA,N,NDIM,ISIGN,IFORM) 
C     COOLEY-TUKEY FAST FOURIER TRANSFORM IN ANS BASIC FORTRAN.
C     MULTI-DIMENSIONAL TRANSFORM, DIMENSIONS OF ARBITRARY SIZE, 
C     COMPLEX OR REAL DATA.  N POINTS CAN BE TRANSFORMED IN TIME 
C     PROPORTIONAL TO N*LOG(N), WHEREAS OTHER METHODS TAKE N**2 TIME.
C     ALSO, LESS ERROR IS BUILT UP. WRITTEN BY NORMAN BRENNER, MIT,
C     FEBRUARY 1969.  SEE-IEEE AUDIO AND ELECTROACOUSTICS TRANSACTIONS 
C     OF JUNE 1967 AMD JUNE 1969, SPECIAL FFT ISSUES.
C     DIMENSION DATA(N(1),N(2),...),TRANSFORM(N(1),N(2),...),N(NDIM) 
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), ETC., FOR ALL
C     NDIM SUBSCRIPTS.  NDIM MUST BE POSITIVE AND EACH N(IDIM) MAY BE
C     ANY POSITIVE INTEGER, NONE OF WHOSE PRIME FACTORS IS GREATER THAN
C     97.  DATA IS COMPLEX, REAL OR HALF-SIZE COMPLEX, AS IFORM = 1, 0 
C     OR -1.  IN THE FIRST TWO CASES, DATA IS DIMENSIONED N(1) BY N(2) 
C     BY ... BY N(NDIM) COMPLEX OR REAL, AS DESCRIBED IN THE N ARRAY.
C     IN THE LAST CASE, IT IS CONCEPTUALLY DIMENSIONED N(1) BY ... BY
C     N(NDIM) COMPLEX, BUT IN FACT IS ONLY N(1)/2+1 BY ... BY N(NDIM), 
C     SINCE THE MISSING VALUES ARE KNOWN TO BE COMPLEX CONJUGATES OF 
C     THOSE PRESENT.  SEE THE EXAMPLE.  ON RETURN, ARRAY DATA HAS BEEN 
C     REFILLED WITH THE TRANSFORM VALUES, WHICH ARE COMPLEX, HALF-SIZE 
C     COMPLEX OR REAL, RESPECTIVELY.  ISIGN IS +1 OR -1.  A +1 TRANS-
C     FORM FOLLOWED BY A -1 TRANSFORM (OR VICE VERSA) RETURNS NTOT 
C     TIMES THE ORIGINAL VALUES TO DATA, WHERE NTOT = N(1)* ...
C     *N(NDIM).  IF IFORM = 0 OR -1, N(1) MUST BE EVEN AND ENOUGH
C     STORAGE RESERVED FOR THE HALF-SIZE COMPLEX ARRAY.
C     EXAMPLE 1.  TRANSFORM OF A COMPLEX VECTOR 1960 LONG. 
C     DIMENSION DATA(1960) 
C     COMPLEX DATA 
C     CALL FOURT (DATA,1960,1,-1,+1) 
C
C     EXAMPLE 2.  TRANSFORM OF A REAL VECTOR 1960 LONG.
C     DIMENSION DATA(1960),TRAN(981) 
C     COMPLEX TRAN 
C     EQUIVALENCE (DATA,TRAN)
C     CALL FOURT (DATA,1960,1,-1,0)
C
C     DIRECTLY IMPLEMENTING THE SUMMATION TAKES A TIME PROPORTIONAL TO 
C     NTOT**2, WHILE FAST FOURIER TRANSFORM TAKES TIME PROPORTIONAL
C     TO NTOT*(SUM OF THE PRIME FACTORS OF NTOT).  THIS IS A TREMENDOUS
C     SPEEDUP FOR LARGE NTOT.  FACTORS OF TWO AND THREE ARE ESPECIALLY 
C     FAST.  IF THE NUMBER OF POINTS IS PRIME, SIMPLY ADD ZERO DATA TO 
C     PAD OUT TO A HIGHLY COMPOSITE LENGTH.  FINALLY, REAL OR HALF-SIZE
C     COMPLEX TRANSFORMS RUN TWICE AS FAST AS COMPLEX TRANSFORMS.
C
C     AN UPPER BOUND FOR THE RMS RELATIVE ERROR IS GIVEN BY GENTLEMAN
C     AND SANDE -- 3 * 2**(-B) * SUM(F**1.5), WHERE 2**(-B) IS THE 
C     SMALLEST BIT IN THE FLOATING POINT FRACTION AND THE SUM IS OVER
C     THE PRIME FACTORS OF NTOT. 
C
C     IF THE INPUT DATA ARE A TIME SERIES, WITH INDEX J REPRESENTING 
C     A TIME (J-1)*DELTAT, THEN THE CORRESPONDING INDEX K IN THE 
C     TRANSFORM REPRESENTS THE FREQUENCY (K-1)*2*PI/(N*DELTAT), WHICH
C     BY PERIODICITY, IS THE SAME AS FREQUENCY -(N-K+1)*2*PI/(N*DELTAT). 
C     THIS IS TRUE FOR N = EACH N(IDIM) INDEPENDENTLY. 
C
C     MACHINE/COMPILER DEPENDENCIES--
C     1.  REAL AND IMAGINARY PARTS MUST BE STORED ADJACENTLY (STANDARD). 
C     2.  ARRAYS ARE STORED LINEARLY WITH THE FIRST SUBSCRIPT INCREASING 
C     FASTEST (STANDARD).
C
      DIMENSION DATA(1),N(1),IFACT(32),IFSYM(32),IFCNT(10),WORK(2,97)
      NWORK=97 
C     MAXIMUM OF 32 PRIME FACTORS PER N(I), MAXIMUM OF 10 UNPAIRED 
C     FACTORS, MAXIMUM FACTOR OF NWORK (USUALLY 97)
      IF (IFORM) 10,10,40
 10   IF (N(1)-2*(N(1)/2)) 20,40,20
 20   WRITE (6,30) IFORM,N(1)
 30   FORMAT (26H0ERROR IN FOURT.  IFORM = ,I2,41H (REAL OR HALF-SIZE CO 
     *MPLEX), BUT N(1) = ,I5,37H IS NOT EVEN.  NO TRANSFORM WAS DONE.) 
      RETURN 
 40   NTOT=1 
      DO 50 IDIM=1,NDIM
 50   NTOT=NTOT*N(IDIM)
      NREM=NTOT
      IF (IFORM) 60,70,70
 60   NREM=1 
      NTOT=(NTOT/N(1))*(N(1)/2+1)
C     LOOP OVER ALL DIMENSIONS.
 70   DO 230 JDIM=1,NDIM 
      IF (IFORM) 80,90,90
 80   IDIM=NDIM+1-JDIM 
      GO TO 100
 90   IDIM=JDIM
      NREM=NREM/N(IDIM)
 100  NCURR=N(IDIM)
      IF (IDIM-1) 110,110,140
 110  IF (IFORM) 120,130,140 
C     PRE-ADJUST A HALF-SIZE COMPLEX ARRAY 
 120  CALL FIXRL (DATA,N(1),NREM,ISIGN,IFORM)
      NTOT=(NTOT/(N(1)/2+1))*N(1)
 130  NCURR=NCURR/2
 140  IF (NCURR-1) 190,190,150 
C     FACTOR N(IDIM), AS 1960 = 2*2*2*3*7*7. 
 150  CALL FACTR (NCURR,IFACT,NFACT) 
      IF (IFACT(NFACT)-NWORK) 180,180,160
 160  WRITE (6,170) NWORK,IFACT(NFACT),IDIM,N(IDIM)
 170  FORMAT (44H0ERROR IN FOURT.  THE WORK ARRAY, OF LENGTH ,I3,
     *  19H, IS SMALLER THAN ,I5,32H, THE LARGEST PRIME FACTOR OF N(,I1, 
     *  4H) = ,I5,1H./23H NO TRANSFORM WAS DONE.)
      RETURN 
C     ARRANGE THE FACTORS SYMMETRICALLY FOR SIMPLER DIGIT REVERSAL.
 180  CALL SMFAC (IFACT,NFACT,ISYM,IFSYM,NFSYM,ICENT,IFCNT,NFCNT)
      NPREV=NTOT/(N(IDIM)*NREM)
C     DIGIT REVERSE ON SYMMETRIC FACTORS, FOR EXAMPLE 2*7*6*7*2. 
      CALL SYMRV (DATA,NPREV,NCURR,NREM,IFSYM,NFSYM) 
C     DIGIT REVERSE THE ASYMMETRIC CENTER, FOR EXAMPLE, ON 6 =2*3. 
      CALL ASMRV(DATA,NPREV*ISYM,ICENT,ISYM*NREM,IFCNT,NFCNT,WORK,NWORK) 
C     FOURIER TRANSFORM ON EACH FACTOR, FOR EXAMPLE, ON 2,7,2,3,7 AND 2. 
      CALL COOL(DATA,NPREV,NCURR,NREM,ISIGN,IFACT,WORK)
 190  IF (IFORM) 200,210,230 
 200  NREM=NREM*N(IDIM)
      GO TO 230
 210  IF (IDIM-1) 220,220,230
C     POST-ADJUST THE TRANSFORM OF A REAL ARRAY, DISGUISED AS COMPLEX
 220  CALL FIXRL (DATA,N(1),NREM,ISIGN,IFORM)
      NTOT=NTOT/N(1)*(N(1)/2+1)
 230  CONTINUE 
      RETURN 
      END
      SUBROUTINE ASMRV (DATA,NPREV,N,NREM,IFACT,NFACT,IWORK,NWORK) 
C     PERMUTE THE DATA ARRAY BY REVERSING THE DIGITS OF ONE INDEX. 
C     THE OPERATION IS THE SAME AS IN SYMRV, EXCEPT THAT THE FACTORS 
C     NEED NOT BE SYMMETRICALLY ARRANGED, I.E., GENERALLY IFACT (IF) 
C     NOT = IFACT(NFACT+1-IF).  SANDE"S METHOD OF CYCLE TRACING IS USED, 
C     WHICH TRANSFORMS THE PERMUTATION INTO A SET OF PAIR INTERCHANGES.
C     A WORK ARRAY IS USED TO KEEP TRACK OF EXCHANGES.  THE LARGER 
C     THE WORK ARRAY, THE FASTER THIS SUBROUTINE RUNS, THO IT IS 
C     NEVER AS FAST AS SYMRV.  NFACT MUST BE LESS THAN 11.  JULY 1970. 
      DIMENSION DATA(1), IFACT(1), IWORK(1), MOD(10) 
      IF (NFACT-1) 190,190,10
 10   IP0=2
      IP1=IP0*NPREV
      IP2=IP1*N
      IP3=IP2*NREM 
      MULT=N/IFACT(1)
C     COMPUTE THE MODULOES.
      IPROD=MULT*IP1 
      INVPR=IPROD
      DO 20 IF=2,NFACT 
      IPROD=IPROD*IFACT(IF-1)
      INVPR=INVPR/IFACT(IF)
 20   MOD(IF)=IPROD-INVPR
      IWMAX=2*NWORK
C     IWORK IS REALLY A COMPLEX ARRAY DISGUISED AS INTEGER.
      I2NXT=1
      DO 30 IW=1,IWMAX 
      IWORK(IW)=I2NXT
 30   I2NXT=I2NXT+IP1
      I2WMX=1+IP1*(IWMAX-1)
      DO 180 I2=1,IP2,IP1
      I2NXT=I2 
C     COMPUTE THE NEXT ELEMENT IN THE PERMUTATION SUBCYCLE.
 40   I2NXT=(I2NXT-1)*MULT+1 
      IF=NFACT 
      GO TO 60 
 50   I2NXT=I2NXT-MOD(IF)*((I2NXT-1)/MOD(IF))
      IF=IF-1
 60   IF (IF-2) 70,70,50 
 70   IQUOT=(I2NXT-1)/MOD(2) 
      IF (IQUOT-IFACT(2)) 90,80,80 
 80   IQUOT=IQUOT-1
 90   I2NXT=I2NXT-MOD(2)*IQUOT 
C     EXCHANGE ELEMENTS (IF I2 IS LESS THAN I2NXT), OR 
C     EXIT IF THE SUBCYCLE HAS CLOSED (IF I2=I2NXT), OR
C     GET THE NEXT SUBCYCLE ELEMENT (IF I2 IS GREATER THAN I2NXT). 
 100  IF (I2-I2NXT) 130,180,110
C     SAVE SOME TIME IF THE NEXT ELEMENT NUMBER IS IN THE IWORK TABLE. 
 110  IF (I2NXT-I2WMX) 120,120,40
 120  IWN=1+(I2NXT-1)/IP1
      I2NXT=IWORK(IWN) 
      GO TO 100
 130  IF (I2-I2WMX) 140,140,160
C     UPDATE THE IWORK TABLE.
 140  IW2=1+(I2-1)/IP1 
      ITEMP=IWORK(IW2) 
      IWT=1+(ITEMP-1)/IP1
      IWORK(IWT)=I2NXT 
      IF (I2NXT-I2WMX) 150,150,160 
 150  IWN=1+(I2NXT-1)/IP1
      IWORK(IWN)=ITEMP 
 160  I1MAX=I2+IP1-IP0 
      DO 170 I1=I2,I1MAX,IP0 
      DO 170 I3=I1,IP3,IP2 
      I3NXT=I3+I2NXT-I2
      TEMPR=DATA(I3) 
      TEMPI=DATA(I3+1) 
      DATA(I3)=DATA(I3NXT) 
      DATA(I3+1)=DATA(I3NXT+1) 
      DATA(I3NXT)=TEMPR
 170  DATA(I3NXT+1)=TEMPI
 180  CONTINUE 
 190  RETURN 
      END
      SUBROUTINE SYMRV (DATA,NPREV,N,NREM,IFACT,NFACT) 
C     PERMUTE THE DATA ARRAY BY REVERSING THE DIGITS OF ONE INDEX. 
C     DIMENSION DATA(NPREV,N,NREM) 
C     REPLACE DATA(I1,I2,I3) BY DATA(I1,I2REV,I3) FOR ALL I1 FROM 1 TO 
C     NPREV, I2 FROM 1 TO N AND I3 FROM 1 TO NREM.  I2REV-1 IS THE 
C     INTEGER WHOSE DIGIT REPRESENTATION IN THE MULTI-RADIX NOTATION 
C     OF FACTORS IFACT(IF) IS THE REVERSE OF THE REPRESENTATION OF I2-1. 
C     FOR EXAMPLE, IF ALL IFACT(IF) = 2, I2-1 = 11001, I2REV-1 = 10011.
C     THE FACTORS MUST BE SYMMETRICALLY ARRANGED, I.E., IFACT(IF) =
C     IFACT(NFACT+1-IF).  THEN, PAIRWISE EXCHANGE CAN BE DONE. 
      DIMENSION DATA(1), IFACT(1)
      IF (NFACT-1) 80,80,10
 10   IP0=2
      IP1=IP0*NPREV
      IP4=IP1*N
      IP5=IP4*NREM 
      I4REV=1
      DO 70 I4=1,IP4,IP1 
      IF (I4-I4REV) 20,40,40 
 20   I1MAX=I4+IP1-IP0 
      DO 30 I1=I4,I1MAX,IP0
      DO 30 I5=I1,IP5,IP4
      I5REV=I4REV+I5-I4
      TEMPR=DATA(I5) 
      TEMPI=DATA(I5+1) 
      DATA(I5)=DATA(I5REV) 
      DATA(I5+1)=DATA(I5REV+1) 
      DATA(I5REV)=TEMPR
 30   DATA(I5REV+1)=TEMPI
 40   IP3=IP4
      DO 60 IF=1,NFACT 
      IP2=IP3/IFACT(IF)
      I4REV=I4REV+IP2
      IF (I4REV-IP3) 70,70,50
 50   I4REV=I4REV-IP3
 60   IP3=IP2
 70   CONTINUE 
 80   RETURN 
      END
      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. 
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).
      DIMENSION DATA(1), WORK(1), IFACT(1) 
      DOUBLE PRECISION TWOPI,THETA,SIN,SINTH,ROOTR,ROOTI,WSTPR,WSTPI,
     *  WR,WI,TEMPR
      SIN(THETA)=DSIN(THETA) 
      TWOPI=6.2831853071795865*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.)
      ROOTR=-2.*SINTH*SINTH
C     COS(THETA)-1, FOR ACCURACY.  USE DOUBLE PRECISION, IF AVAILABLE. 
      ROOTI=SIN(THETA) 
      THETA=TWOPI/FLOAT(IP3/IP1) 
      SINTH=SIN(THETA/2.)
      WSTPR=-2.*SINTH*SINTH
      WSTPI=SIN(THETA) 
      WR=1.
      WI=0.
      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.*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)-.5*SUMR 
      TEMPI=DATA(J0+1)-.5*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, 
     $ROOTI) 
 220  TEMPR=WR 
      WR=WSTPR*TEMPR-WSTPI*WI+TEMPR
  230 WI=WSTPR*WI+WSTPI*TEMPR+WI 
      IP2=IP3
      GO TO 10 
 240  RETURN 
      END
      SUBROUTINE SMFAC (IFACT,NFACT,ISYM,IFSYM,NFSYM,ICENT,IFCNT,NFCNT)
C     REARRANGE THE PRIME FACTORS OF N INTO A SQUARE AND A NON-
C     SQUARE.  N = ISYM*ICENT*ISYM, WHERE ICENT IS SQUARE-FREE.
C     ISYM = IFSYM(1)*...*IFSYM(NFSYM), EACH A PRIME FACTOR. 
C     ICENT = IFCNT(1)*...*IFCNT(NFCNT), EACH A PRIME FACTOR.
C     FOR EXAMPLE, N = 1960 = 14*10*14.  THEN ISYM = 14, ICENT = 10, 
C     NFSYM = 2, NFCNT = 2, NFACT = 6, IFSYM(IFS) = 2, 7, IFCNT(IFC) = 
C     2, 5 AND IFACT(IF) = 2, 7, 2, 5, 7, 2. 
      DIMENSION IFSYM(1), IFCNT(1), IFACT(1) 
      ISYM=1 
      ICENT=1
      IFS=0
      IFC=0
      IF=1 
 10   IF (IF-NFACT) 20,40,50 
 20   IF (IFACT(IF)-IFACT(IF+1)) 40,30,40
 30   IFS=IFS+1
      IFSYM(IFS)=IFACT(IF) 
      ISYM=IFACT(IF)*ISYM
      IF=IF+2
      GO TO 10 
 40   IFC=IFC+1
      IFCNT(IFC)=IFACT(IF) 
      ICENT=IFACT(IF)*ICENT
      IF=IF+1
      GO TO 10 
 50   NFSYM=IFS
      NFCNT=IFC
      NFSM2=2*NFSYM
      NFACT=2*NFSYM+NFCNT
      IF (NFCNT) 80,80,60
 60   NFSM2=NFSM2+1
      IFSYM(NFSYM+1)=ICENT 
      DO 70 IFC=1,NFCNT
      IF=NFSYM+IFC 
 70   IFACT(IF)=IFCNT(IFC) 
 80   IF (NFSYM) 110,110,90
 90   DO 100 IFS=1,NFSYM 
      IFSCJ=NFSM2+1-IFS
      IFSYM(IFSCJ)=IFSYM(IFS)
      IFACT(IFS)=IFSYM(IFS)
      IFCNJ=NFACT+1-IFS
 100  IFACT(IFCNJ)=IFSYM(IFS)
 110  NFSYM=NFSM2
      RETURN 
      END
      SUBROUTINE FIXRL (DATA,N,NREM,ISIGN,IFORM) 
C     FOR IFORM = 0, CONVERT THE TRANSFORM OF A DOUBLED-UP REAL ARRAY, 
C     CONSIDERED COMPLEX, INTO ITS TRUE TRANSFORM.  SUPPLY ONLY THE
C     FIRST HALF OF THE COMPLEX TRANSFORM, AS THE SECOND HALF HAS
C     CONJUGATE SYMMETRY.  FOR IFORM = -1, CONVERT THE FIRST HALF
C     OF THE TRUE TRANSFORM INTO THE TRANSFORM OF A DOUBLED-UP REAL
C     ARRAY.  N MUST BE EVEN.
C     USING COMPLEX NOTATION AND SUBSCRIPTS STARTING AT ZERO, THE
C     TRANSFORMATION IS--
C     DIMENSION DATA(N/2,NREM) 
C     ZSTP = EXP(ISIGN*2*PI*I/N) 
C     DO 10 I2=0,NREM-1
C     DATA(0,I2) = CONJ(DATA(0,I2))*(1+I)
C     DO 10 I1=1,N/4 
C     Z = (1+(2*IFORM+1)*I*ZSTP**I1)/2 
C     I1CNJ = N/2-I1 
C     DIF = DATA(I1,I2)-CONJ(DATA(I1CNJ,I2)) 
C     TEMP = Z*DIF 
C     DATA(I1,I2) = (DATA(I1,I2)-TEMP)*(1-IFORM) 
C 10  DATA(I1CNJ,I2) = (DATA(I1CNJ,I2)+CONJ(TEMP))*(1-IFORM) 
C     IF I1=I1CNJ, THE CALCULATION FOR THAT VALUE COLLAPSES INTO 
C     A SIMPLE CONJUGATION OF DATA(I1,I2). 
      DIMENSION DATA(1)
      DOUBLE PRECISION TWOPI,THETA,SIN,SINTH,ZSTPR,ZSTPI,ZR,ZI,TEMPR 
      SIN(THETA)=DSIN(THETA) 
      TWOPI=6.2831853071795865*FLOAT(ISIGN)
      IP0=2
      IP1=IP0*(N/2)
      IP2=IP1*NREM 
      IF (IFORM) 10,70,70
C     PACK THE REAL INPUT VALUES (TWO PER COLUMN)
 10   J1=IP1+1 
      DATA(2)=DATA(J1) 
      IF (NREM-1) 70,70,20 
 20   J1=J1+IP0
      I2MIN=IP1+1
      DO 60 I2=I2MIN,IP2,IP1 
      DATA(I2)=DATA(J1)
      J1=J1+IP0
      IF (N-2) 50,50,30
 30   I1MIN=I2+IP0 
      I1MAX=I2+IP1-IP0 
      DO 40 I1=I1MIN,I1MAX,IP0 
      DATA(I1)=DATA(J1)
      DATA(I1+1)=DATA(J1+1)
 40   J1=J1+IP0
 50   DATA(I2+1)=DATA(J1)
 60   J1=J1+IP0
 70   DO 80 I2=1,IP2,IP1 
      TEMPR=DATA(I2) 
      DATA(I2)=DATA(I2)+DATA(I2+1) 
 80   DATA(I2+1)=TEMPR-DATA(I2+1)
      IF (N-2) 200,200,90
 90   THETA=TWOPI/FLOAT(N) 
      SINTH=SIN(THETA/2.)
      ZSTPR=-2.*SINTH*SINTH
      ZSTPI=SIN(THETA) 
      ZR=(1.-ZSTPI)/2. 
      ZI=(1.+ZSTPR)/2. 
      IF (IFORM) 100,110,110 
 100  ZR=1.-ZR 
      ZI=-ZI 
 110  I1MIN=IP0+1
      I1MAX=IP0*(N/4)+1
      DO 190 I1=I1MIN,I1MAX,IP0
      DO 180 I2=I1,IP2,IP1 
      I2CNJ=I2+IP0*(N/2)-2*(I1-1)
      IF (I2-I2CNJ) 150,120,120
 120  IF (ISIGN*(2*IFORM+1)) 130,140,140 
 130  DATA(I2+1)=-DATA(I2+1) 
 140  IF (IFORM) 170,180,180 
 150  DIFR=DATA(I2)-DATA(I2CNJ)
      DIFI=DATA(I2+1)+DATA(I2CNJ+1)
      TEMPR=DIFR*ZR-DIFI*ZI
      TEMPI=DIFR*ZI+DIFI*ZR
      DATA(I2)=DATA(I2)-TEMPR
      DATA(I2+1)=DATA(I2+1)-TEMPI
      DATA(I2CNJ)=DATA(I2CNJ)+TEMPR
      DATA(I2CNJ+1)=DATA(I2CNJ+1)-TEMPI
      IF (IFORM) 160,180,180 
 160  DATA(I2CNJ)=DATA(I2CNJ)+DATA(I2CNJ)
      DATA(I2CNJ+1)=DATA(I2CNJ+1)+DATA(I2CNJ+1)
 170  DATA(I2)=DATA(I2)+DATA(I2) 
      DATA(I2+1)=DATA(I2+1)+DATA(I2+1) 
 180  CONTINUE 
      TEMPR=ZR-.5
      ZR=ZSTPR*TEMPR-ZSTPI*ZI+ZR 
 190  ZI=ZSTPR*ZI+ZSTPI*TEMPR+ZI 
C     RECURSION SAVES TIME, AT A SLIGHT LOSS IN ACCURACY.  IF AVAILABLE, 
C     USE DOUBLE PRECISION TO COMPUTE ZR AND ZI. 
 200  IF (IFORM) 270,210,210 
C     UNPACK THE REAL TRANSFORM VALUES (TWO PER COLUMN)
 210  I2=IP2+1 
      I1=I2
      J1=IP0*(N/2+1)*NREM+1
      GO TO 250
 220  DATA(J1)=DATA(I1)
      DATA(J1+1)=DATA(I1+1)
      I1=I1-IP0
      J1=J1-IP0
 230  IF (I2-I1) 220,240,240 
 240  DATA(J1)=DATA(I1)
      DATA(J1+1)=0.
 250  I2=I2-IP1
      J1=J1-IP0
      DATA(J1)=DATA(I2+1)
      DATA(J1+1)=0.
      I1=I1-IP0
      J1=J1-IP0
      IF (I2-1) 260,260,230
 260  DATA(2)=0. 
 270  RETURN 
      END
      SUBROUTINE GOERT(DATA,NPREV,IPROD,IFACT,IREM,WORK,WMINR,WMINI, 
     $ ROOTR,ROOTI)
      DOUBLE PRECISION WMINR,WMINI,ROOTR,ROOTI 
C     PHASE-SHIFTED FOURIER TRANSFORM OF LENGTH IFACT BY THE GOERTZEL
C     ALGORITHM.  IFACT MUST BE ODD AND AT LEAST FIVE.  FURTHER SPEED
C     IS GAINED BY COMPUTING TWO TRANSFORM VALUES AT THE SAME TIME.
C     DIMENSION DATA(NPREV,IPROD,IFACT,IREM) 
C     DATA(I1,1,J3,I5) = SUM(DATA(I1,1,I3,I5) * W**(I3-1)), SUMMED 
C     OVER I3 = 1 TO IFACT FOR ALL I1 FROM 1 TO NPREV, J3 FROM 1 TO
C     IFACT AND I5 FROM 1 TO IREM. 
C     W = WMIN * EXP(ISIGN*2*PI*I*(J3-1)/IFACT). 
      DIMENSION DATA(1), WORK(1) 
      IP0=2
      IP1=IP0*NPREV
      IP2=IP1*IPROD
      IP3=IP2*IFACT
      IP5=IP3*IREM 
      IF (WMINI) 10,40,10
C     APPLY THE PHASE SHIFT FACTORS
 10   WR=WMINR 
      WI=WMINI 
      I3MIN=1+IP2
      DO 30 I3=I3MIN,IP3,IP2 
      I1MAX=I3+IP1-IP0 
      DO 20 I1=I3,I1MAX,IP0
      DO 20 I5=I1,IP5,IP3
      TEMPR=DATA(I5) 
      DATA(I5)=WR*TEMPR-WI*DATA(I5+1)
 20   DATA(I5+1)=WR*DATA(I5+1)+WI*TEMPR
      TEMPR=WR 
      WR=WMINR*TEMPR-WMINI*WI
 30   WI=WMINR*WI+WMINI*TEMPR
 40   DO 90 I1=1,IP1,IP0 
      DO 90 I5=I1,IP5,IP3
C     STRAIGHT SUMMATION FOR THE FIRST TERM
      SUMR=0.
      SUMI=0.
      I3MAX=I5+IP3-IP2 
      DO 50 I3=I5,I3MAX,IP2
      SUMR=SUMR+DATA(I3) 
 50   SUMI=SUMI+DATA(I3+1) 
      WORK(1)=SUMR 
      WORK(2)=SUMI 
      WR=ROOTR+1.
      WI=ROOTI 
      IWMIN=1+IP0
      IWMAX=IP0*((IFACT+1)/2)-1
      DO 80 IWORK=IWMIN,IWMAX,IP0
      TWOWR=WR+WR
      I3=I3MAX 
      OLDSR=0. 
      OLDSI=0. 
      SUMR=DATA(I3)
      SUMI=DATA(I3+1)
      I3=I3-IP2
 60   TEMPR=SUMR 
      TEMPI=SUMI 
      SUMR=TWOWR*SUMR-OLDSR+DATA(I3) 
      SUMI=TWOWR*SUMI-OLDSI+DATA(I3+1) 
      OLDSR=TEMPR
      OLDSI=TEMPI
      I3=I3-IP2
      IF (I3-I5) 70,70,60
C     IN A FOURIER TRANSFORM THE W CORRESPONDING TO THE POINT AT K 
C     IS THE CONJUGATE OF THAT AT IFACT-K (THAT IS, EXP(TWOPI*I* 
C     K/IFACT) = CONJ(EXP(TWOPI*I*(IFACT-K)/IFACT))).  SINCE THE 
C     MAIN LOOP OF GOERTZELS ALGORITHM IS INDIFFERENT TO THE IMAGINARY 
C     PART OF W, IT NEED BE SUPPLIED ONLY AT THE END.
 70   TEMPR=-WI*SUMI 
      TEMPI=WI*SUMR
      SUMR=WR*SUMR-OLDSR+DATA(I3)
      SUMI=WR*SUMI-OLDSI+DATA(I3+1)
      WORK(IWORK)=SUMR+TEMPR 
      WORK(IWORK+1)=SUMI+TEMPI 
      IWCNJ=IP0*(IFACT+1)-IWORK
      WORK(IWCNJ)=SUMR-TEMPR 
      WORK(IWCNJ+1)=SUMI-TEMPI 
C     SINGLETON"S RECURSION, FOR ACCURACY AND SPEED (5). 
      TEMPR=WR 
      WR=WR*ROOTR-WI*ROOTI+WR
 80   WI=TEMPR*ROOTI+WI*ROOTR+WI 
      IWORK=1
      DO 90 I3=I5,I3MAX,IP2
      DATA(I3)=WORK(IWORK) 
      DATA(I3+1)=WORK(IWORK+1) 
 90   IWORK=IWORK+IP0
      RETURN 
      END
