C
C********************************************************************
      SUBROUTINE FFT1(A,M,N)
C***********************************************************************
C
C *** DECIMATION IN TIME FFT (COOLEY-TUKEY ALGORITHM) FOR NATURALLY
C *** ORDERED INPUTS
C
C *** FFT1 COMPUTES FORWARD DFT.  INVERSE DFT CAN BE COMPUTED WITH 
C *** FFT1 BY FOLLOWING SEQUENCE:
C ***     * COMPLEX CONJUGATE INPUT A
C ***     * CALL FFT1
C ***     * COMPLEX CONJUGATE RESULT A
C ***     * NORMALIZE BY 1/N IF DESIRED
C
C ***   INPUTS:
C ***     A          COMPLEX VECTOR OF LENGTH N CONTAINING NATURALLY
C ***                ORDERED VALUES
C ***     N          FFT SIZE, MUST BE POWER OF 2, 2**M
C ***     M          LOG2(N)
C
C ***   OUTPUTS:
C ***     A          FORWARD N-POINT DISCRETE TRANSFORM
C
      COMPLEX A,U,W,T
C
      DIMENSION A(N)
C
      DATA PI/3.14159265/
C
      N=2**M
C
C *** SCRAMBLE INPUT ARRAY (BIT REVERSE) SO THAT OUTPUTS WILL BE
C *** NATURALLY ORDERED
      NV2=N/2
      NM1=N-1
      J=1
      DO 40 I=1,NM1
        IF(I.LT.J) THEN
          T=A(J)
          A(J)=A(I)
          A(I)=T
        ENDIF
        K=NV2
   20   CONTINUE
        IF(K.LT.J) THEN
          J=J-K
          K=K/2
          GO TO 20
        ENDIF
        J=J+K
   40 CONTINUE
C
C *** PERFORM DECIMATION IN TIME FFT
      DO 60 L=1,M
        LE=2**L
        LE1=LE/2
        U=CMPLX(1.,0.)
        W=CMPLX(COS(PI/LE1),-SIN(PI/LE1))
        DO 60 J=1,LE1
          DO 50 I=J,N,LE
            IP=I+LE1
            T=A(IP)*U
            A(IP)=A(I)-T
            A(I)=A(I)+T
   50     CONTINUE
          
      END
