c Replaces DATA by its discrete Fourier transform, if ISIGN is input as 1; 
c or replaces DATA by NN times its inverse discrete Fourier transform, 
c if ISIGN is input as -1.  DATA is a complex array of length NN or, 
c equivalently, a real array of length 2*NN 

      SUBROUTINE FOUR1(DATA,NN,ISIGN)
      
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER N,J,I,M
      REAL*8 WR,WI,WPR,WPI,WTEMP,THETA
      REAL*8 TEMPR
      DIMENSION DATA(2*NN)

      N = 2*NN
      J = 1
      DO I = 1, N, 2
         IF (J .GT. I) THEN
            TEMPR = DATA(J)
            TEMPI = DATA(J+1)
            DATA(J) = DATA(I)
            DATA(J+1) = DATA(I+1)
            DATA(I) = TEMPR
            DATA(I+1) = TEMPI
         ENDIF
         M = N/2
   1     IF ((M .GT. 2) .AND. (J .GT. M)) THEN
            J = J-M
            M = M/2
            GO TO 1
         ENDIF
         J = J+M
      ENDDO
      MMAX = 2
   2  IF (N .GT. MMAX) THEN  !Outer loop executed log2 NN times.
         ISTEP = 2*MMAX
         THETA = 6.28318530717959D0/(ISIGN*MMAX) !Initialize for the trigonometric recurrence.
         WPR = -2.D0*DSIN(0.5D0*THETA)**2
         WPI = DSIN(THETA)
         WR = 1.D0
         WI = 0.D0
         DO M=1, MMAX, 2
            DO I = M, N, ISTEP
               J = I+MMAX
               TEMPR = SNGL(WR)*DATA(J)-SNGL(WI)*DATA(J+1)
               TEMPI = SNGL(WR)*DATA(J+1)+SNGL(WI)*DATA(J)
               DATA(J) = DATA(I)-TEMPR
               DATA(J+1) = DATA(I+1)-TEMPI
               DATA(I) = DATA(I)+TEMPR
               DATA(I+1) = DATA(I+1)+TEMPI
            ENDDO
            WTEMP = WR
            WR = WR*WPR-WI*WPI+WR
            WI = WI*WPR+WTEMP*WPI+WI
         ENDDO
         MMAX = ISTEP
         GO TO 2         !Not yet done.
      ENDIF              !All done
      RETURN
      END
