SUBROUTINE DREALFT(DATA,NIN,ISIGN)
*     ================================
* Double precision version
* Data in FOUR2 format i.e. complex array runs from 0->N/2, thus 2 extra
* storage elements required.
* 2-May-1991 DFB Change normalisation so it's the same as for FOUR2 i.e.
*                a forward transform followed by a backward transform
*                gives original data multiplied by NIN, where NIN is the
*                length of the real dataset or equivalently the length
*                of the complex dataset when expanded to include the
*                symmetric half.
      IMPLICIT REAL*8 (A-H, O-Z)
      REAL*8 WR,WI,WPR,WPI,WTEMP,THETA
      DIMENSION DATA(*)

      N=NIN/2
      THETA=6.28318530717959D0/2.0D0/DBLE(N)
      WR=1.0D0
      WI=0.0D0
      IF (ISIGN.EQ.1) THEN
        C1=0.5D0
        C2=-C1
        CALL DFOUR1(DATA,N,+1)
        DATA(2*N+1)=DATA(1)
        DATA(2*N+2)=DATA(2)
      ELSE
        C1=1D0
        C2=C1
        THETA=-THETA
C - Change to FOUR2 format
C        DATA(2*N+1)=DATA(2)
C        DATA(2*N+2)=0.0
        DATA(2)=0.0
      ENDIF
      WPR=-2.0D0*DSIN(0.5D0*THETA)**2
      WPI=DSIN(THETA)
      N2P3=2*N+3
      DO 11 I=1,N/2+1
        I1=2*I-1
        I2=I1+1
        I3=N2P3-I2
        I4=I3+1
        H1R=C1*(DATA(I1)+DATA(I3))
        H1I=C1*(DATA(I2)-DATA(I4))
        H2R=-C2*(DATA(I2)+DATA(I4))
        H2I=C2*(DATA(I1)-DATA(I3))
        DATA(I1)=H1R+WR*H2R-WI*H2I
        DATA(I2)=H1I+WR*H2I+WI*H2R
        DATA(I3)=H1R-WR*H2R+WI*H2I
        DATA(I4)=-H1I+WR*H2I+WI*H2R
        WTEMP=WR
        WR=WR*WPR-WI*WPI+WR
        WI=WI*WPR+WTEMP*WPI+WI
11    CONTINUE
      IF (ISIGN.EQ.1) THEN
C - Change to FOUR2 format
C        DATA(2)=DATA(2*N+1)
      ELSE
        CALL DFOUR1(DATA,N,-1)
      ENDIF
      RETURN
      END


      SUBROUTINE DFOUR1(DATA,NN,ISIGN)
* Double precision version of Num. Recipes routine.
* Recursion for COS and SIN removed.
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 WR,WI,WPR,WPI,WTEMP,THETA
      DIMENSION DATA(*)
      N=2*NN
      J=1
      DO 11 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.GE.2).AND.(J.GT.M)) THEN
          J=J-M
          M=M/2
        GO TO 1
        ENDIF
        J=J+M
11    CONTINUE
      MMAX=2
2     IF (N.GT.MMAX) THEN
        ISTEP=2*MMAX
        THETA=6.28318530717958647692D0/(ISIGN*MMAX)
c        WPR=-2.D0*DSIN(0.5D0*THETA)**2
c        WPI=DSIN(THETA)
        WR=1.D0
        WI=0.D0
        icntr=1
        DO 13 M=1,MMAX,2
          DO 12 I=M,N,ISTEP
            J=I+MMAX
            TEMPR=WR*DATA(J)-WI*DATA(J+1)
            TEMPI=WR*DATA(J+1)+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
12        CONTINUE
c          WTEMP=WR
c          WR=WR*WPR-WI*WPI+WR
c          WI=WI*WPR+WTEMP*WPI+WI
          wr=dcos(icntr*theta)
          wi=dsin(icntr*theta)
          icntr=icntr+1
13      CONTINUE
        MMAX=ISTEP
      GO TO 2
      ENDIF
      RETURN
      END