C
C $Id: fourier.F,v 1.4 1998/07/16 16:39:47 jjv5 Exp arjan $
C
C------------------------------------------------------------------------
C     Fourier transformations
      SUBROUTINE FOURN(DATA,NN,NDIM,ISIGN)
      INTEGER ISIGN,NDIM,NN(NDIM)
      DOUBLE PRECISION DATA(*)
      INTEGER I1,I2,I2REV,I3,I3REV,IBIT,IDIM,IFP1,IFP2,IIP1,IIP2,
     *      IIP3,K1,K2,N,NPREV,NREM,NTOT
      DOUBLE PRECISION TEMPI,TEMPR
      DOUBLE PRECISION THETA,WI,WPI,WPR,WR,WTEMP
      NTOT=1
      DO 11 IDIM=1,NDIM
         NTOT=NTOT*NN(IDIM)
11    CONTINUE
      NPREV=1
      DO 18 IDIM=1,NDIM
         N=NN(IDIM)
         NREM=NTOT/(N*NPREV)
         IIP1=2*NPREV
         IIP2=IIP1*N
         IIP3=IIP2*NREM
         I2REV=1
         DO 14 I2=1,IIP2,IIP1
            IF(I2.LT.I2REV)THEN
               DO 13 I1=I2,I2+IIP1-2,2
                  DO 12 I3=I1,IIP3,IIP2
                     I3REV=I2REV+I3-I2
                     TEMPR=DATA(I3)
                     TEMPI=DATA(I3+1)
                     DATA(I3)=DATA(I3REV)
                     DATA(I3+1)=DATA(I3REV+1)
                     DATA(I3REV)=TEMPR
                     DATA(I3REV+1)=TEMPI
12                CONTINUE
13             CONTINUE
            ENDIF
            IBIT=IIP2/2
1           IF ((IBIT.GE.IIP1).AND.(I2REV.GT.IBIT)) THEN
               I2REV=I2REV-IBIT
               IBIT=IBIT/2
            GO TO 1
            ENDIF
            I2REV=I2REV+IBIT
14       CONTINUE
         IFP1=IIP1
2        IF(IFP1.LT.IIP2)THEN
            IFP2=2*IFP1
            THETA=ISIGN*6.28318530717959D0/(IFP2/IIP1)
            WPR=-2.D0*DSIN(0.5D0*THETA)**2
            WPI=DSIN(THETA)
            WR=1.D0
            WI=0.D0
            DO 17 I3=1,IFP1,IIP1
               DO 16 I1=I3,I3+IIP1-2,2
                  DO 15 I2=I1,IIP3,IFP2
                     K1=I2
                     K2=K1+IFP1
                     TEMPR=SNGL(WR)*DATA(K2)-SNGL(WI)*DATA(K2+1)
                     TEMPI=SNGL(WR)*DATA(K2+1)+SNGL(WI)*DATA(K2)
                     DATA(K2)=DATA(K1)-TEMPR
                     DATA(K2+1)=DATA(K1+1)-TEMPI
                     DATA(K1)=DATA(K1)+TEMPR
                     DATA(K1+1)=DATA(K1+1)+TEMPI
15                CONTINUE
16             CONTINUE
               WTEMP=WR
               WR=WR*WPR-WI*WPI+WR
               WI=WI*WPR+WTEMP*WPI+WI
17          CONTINUE
            IFP1=IFP2
         GO TO 2
         ENDIF
         NPREV=N*NPREV
18    CONTINUE
      RETURN
      END