SUBROUTINE SINFT(Y,N)
C *** Calculates the sine transform of a set of N real-valued points
C *** stored in Y, N must be a power of 2.  Also calculates N/2 times
C *** the inverse transform.  F(f)=Sum (j=1 to N-1) f(j)sin(pi*j*f/N)         The essential definition
      REAL*8 WR,WI,WPR,WPI,WTEMP,THETA
      DIMENSION Y(N+2)
      THETA=3.14159265358979D0/N
      WR=1
      WI=0
      WPR=-2*DSIN(.5D0*THETA)**2
      WPI=DSIN(THETA)
      Y(1)=0.0
      M=N/2
      DO J=1,M
        WTEMP=WR
        WR=WR*WPR-WI*WPI+WR
        WI=WI*WPR+WTEMP*WPI+WI
        Y1=WI*(Y(J+1)+Y(N-J+1))
        Y2=0.5D0*(Y(J+1)-Y(N-J+1))
        Y(J+1)=Y1+Y2
        Y(N-J+1)=Y1-Y2
      ENDDO
      CALL REALFT(Y,M,+1)
      SUM=0
      Y(1)=0.5D0*Y(1)
      Y(2)=0
      DO J=1,N-1,2
        SUM=SUM+Y(J)
        Y(J)=Y(J+1)
        Y(J+1)=SUM
      ENDDO
      RETURN
      END

      SUBROUTINE REALFT(DATA,N,ISIGN)
C *** calculates Fourier Transform of 2N real valued data points
C *** N must be a power of 2 (with ISIGN=-1  gives N times inverse
C *** transform)
      REAL*8 WR,WI,WPR,WPI,WTEMP,THETA
      DIMENSION DATA(2*N+2)
      THETA=6.28318530717959D0/(2*N)
      WR=1
      WI=0
      C1=0.5D0
      IF(ISIGN.EQ.1)THEN
        C2=-C1
        CALL FFT(DATA,N,+1)
        DATA(2*N+1)=DATA(1)
        DATA(2*N+2)=DATA(2)
      ELSE
        C2=C1
        THETA=-THETA
        DATA(2*N+1)=DATA(2)
        DATA(2*N+2)=0
        DATA(2)=0
      ENDIF
      WPR=-2*DSIN(0.5D0*THETA)**2
      WPI=DSIN(THETA)
      N2P3=N*2+3
      DO 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
      ENDDO
      IF(ISIGN.EQ.1)THEN
        DATA(2)=DATA(2*N+1)
      ELSE
        CALL FFT(DATA,N,-1)
      ENDIF
      RETURN
      END