C    PROGRAM FFTBA2S BENSON ACKERMAN FINAL VERSION 03/22/74
C    SCALING ALGORITHM MODIFIED 12/31/81	
C      CONTAINS SUBROUTINES FFTBA2
C      UNSCRM
C      FOLDFT
C      COMPLX
      SUBROUTINE FFTBA2(XR,XI,N,L,JUNS)
      DOUBLE PRECISION TPI,DPCON,THETA,STHETA,CTHETA,SINTH,COSTH,
     1CHOLD,CSTAY
      DOUBLE PRECISION DATAN,DCOS,DSIN
      DIMENSION XR(N),XI(N)
      XN=N
      M=ALOG(XN)/ALOG(2.)+.5
      IF(L.EQ.0)GOTO21
      TPI=8.D+00*DATAN(1.D+00)
      IF(M.EQ.0)GOTO99
      THETA=TPI/N
      IF(L.LT.0) THETA=-THETA
      NSPAN=N
      NDIST=NSPAN/2
      CTHETA=DCOS(THETA)
      STHETA=DSIN(THETA)
      DO 2 I=1,M
      COSTH=1.D+00
      SINTH=0.
      DO 3 J=1,NDIST
      TR=COSTH
      TI=SINTH
      DO 5 K=J,N,NSPAN
      KANDS=K+NDIST
      HOLDR=XR(K)-XR(KANDS)
      HOLDI=XI(K)-XI(KANDS)
      XR(K)=XR(K)+XR(KANDS)
      XI(K)=XI(K)+XI(KANDS)
      XR(KANDS)=TR*HOLDR-TI*HOLDI
5     XI(KANDS)=TI*HOLDR+TR*HOLDI
      CHOLD=COSTH*CTHETA-SINTH*STHETA
      SINTH=SINTH*CTHETA+COSTH*STHETA
3     COSTH=CHOLD
      NSPAN=NDIST
      NDIST=NSPAN/2
      CSTAY=CTHETA*CTHETA-STHETA*STHETA
      STHETA=STHETA*CTHETA+CTHETA*STHETA
2     CTHETA=CSTAY
      IF(L.LT.0)GOTO4
      CON=1./N
      DO 1 J=1,N
      XR(J)=CON*XR(J)
1     XI(J)=CON*XI(J)
4     IF(M.EQ.1)GOTO99
21    IF(JUNS.EQ.0)GOTO99
      CALL UNSCRM(XR,XI,N,M)
99    RETURN
      END
      SUBROUTINE UNSCRM(XR,XI,N,M)
      DIMENSION XR(N),XI(N),IX(128)
      IF(M.EQ.1)GOTO99
      M2=M/2
      M3=M-2*M2
      K2M2=2**M2
      K2MX=K2M2
      IF(M3.EQ.1)K2MX=K2MX+K2MX
      IX(1)=0
      JLIM=1
      DO 3 I=1,M2
      DO 2 J=1,JLIM
      IX(J)=IX(J)+IX(J)
2     IX(JLIM+J)=IX(J)+K2MX
3     JLIM=JLIM+JLIM
      KUMX=K2M2-1
      IF(M3.EQ.1)GOTO10
      DO 11 KU=1,KUMX
      KVMN=KU+1
      DO 11 KV=KVMN,K2M2
      KP1=IX(KU)+KV
      KP2=IX(KV)+KU
      XSR=XR(KP1)
      XSI=XI(KP1)
      XR(KP1)=XR(KP2)
      XI(KP1)=XI(KP2)
      XR(KP2)=XSR
      XI(KP2)=XSI
11    CONTINUE
      GOTO99
10    DO 12 KU=1,KUMX
      KVMN=KU+1
      DO 12 KV=KVMN,K2M2
      KP1=IX(KU)+KV
      KP2=IX(KV)+KU
      XSR=XR(KP1)
      XSI=XI(KP1)
      XR(KP1)=XR(KP2)
      XI(KP1)=XI(KP2)
      XR(KP2)=XSR
      XI(KP2)=XSI
      KP1=KP1+K2M2
      KP2=KP2+K2M2
      XSR=XR(KP1)
      XSI=XI(KP1)
      XR(KP1)=XR(KP2)
      XI(KP1)=XI(KP2)
      XR(KP2)=XSR
      XI(KP2)=XSI
12    CONTINUE
99    RETURN
      END
      SUBROUTINE FOLDFT(XR,XI,N,ITYPE)
      DIMENSION XR(N),XI(N)
C      IF ITYPE < 0:
C      VECTORS XR AND XI ARE FOLDED ABOUT THE N/2+1 TERM. REAL FOLDED
C      VALUES ARE ADDED, IMAGINARY VALUES ARE SUBTRACTED. THE
C      RESULTING VALUES ARE STORED IN THE FIRST N/2 LOCATIONS.
C
C      IF ITYPE > 0: VECTORS XR AND XI ARE UNFOLDED ABOUT THE N/2+1 TER
C     M.
C      ONLY THE FIRST N/2 TERMS ARE USED. ALL VALUES ARE DIVIDED BY 2 A
C     ND
C      IN ADDITION THE IMAGINARY TERMS ARE NEGATED. THIS PROCEDURE IS T
C     HE
C      INVERSE OF THAT WHEN ITYPE < 0.ANY ORIGINAL N/2+2 TO N VALUES AR
C     E
C      DESTROYED.
C
C      IT ITYPE = 0: THE FIRST TO N-1 ALTERNATE TERMS OF XR ARE STORED
C      IN THE FIRST N/2 TERMS OF XR. THE SECOND TO N ALTERNATE TERMS OF
C      XR ARE STORED IN THE FIRST N/2 TERMS OF XI.ANY ORIGIONAL VALUED
C      STORES IN THE FIRST N/2 TERMS OF XI ARE DESTROYED. THE N/2+1 TO
C     N
C      TERMS OF BOTH XR AND XI ARE NOT USED.
      NHALF=N/2
      NP2=N+2
      IF(ITYPE)21,51,41
21    XR(1)=2.*XR(1)
      DO 1 I=2,NHALF
      IDEX=NP2-I
      XR(I)=XR(I)+XR(IDEX)
      XI(I)=XI(I)-XI(IDEX)
      XR(IDEX)=0.
1     XI(IDEX)=0.
      GOTO99
41    XR(1)=.5*XR(1)
      XI(1)=0.
      DO 2 I=2,NHALF
      IDEX=NP2-I
      XR(IDEX)=.5*XR(I)
      XR(I)=XR(IDEX)
      XI(IDEX)=-.5*XI(I)
2     XI(I)=-XI(IDEX)
      GOTO99
51    DO 3 I=1,NHALF
      I2=I+I
      XR(I)=XR(I2-1)
3     XI(I)=XR(I2)
99    RETURN
      END
      SUBROUTINE COMPLX(XR,XI,XRN,XIN,NHALF)
      DIMENSION XR(NHALF),XI(NHALF),XRN(NHALF),XIN(NHALF)
      FACT=4.*ATAN(1.)/NHALF
      XRN(1)=(XR(1)+XI(1))*.5
      XIN(1)=(XI(1)-XR(1))*.5
      IF(NHALF.LT.2)GOTO99
      ARG=FACT
      DO 1 I=2,NHALF
      IDEX=NHALF+2-I
      XRD=XR(I)-XR(IDEX)
      XIS=XI(I)+XI(IDEX)
      CFACT=COS(ARG)
      SFACT=SIN(ARG)
      ARG=FACT*I
      XRN(I)=(XR(I)+XR(IDEX)+CFACT*XIS-SFACT*XRD)*.5
1     XIN(I)=(XI(I)-XI(IDEX)-SFACT*XIS-CFACT*XRD)*.5
99    RETURN
      END