SUBROUTINE four1(data,nn,isign)
c     Fast Fourier Tranform (FFT) from Numerical recepies, with some
c     changes for real*8 precision
      implicit none
c     nn .. number of complex numbers (2*nn real ones)
c     data(1) Real f0
c     data(2) Im f0
c     data(3) Real f1
c     data(4) Im f1   etc
c      return the same ordering. f0 zero-frequency 
c     Replaces data(1:2*nn) by its discrete Fourier transform, if isign is 
c     input as 1; or replaces data(1:2*nn) by nn times its inverse discrete
c     Fourier transform, if isign is input as -1.
c     data is a complex array of length nn or, equivalently, a real array 
c     of length 2*nn. nn MUST be an integer power of 2 
c     (this is not checked for!).
      INTEGER*4 isign,nn
      REAL*8 data(2*nn)
      INTEGER*4 i,istep,j,m,mmax,n
      REAL*8 tempi,tempr,theta,wi,wpi,wpr,wr,wtemp 
c
      n=2*nn
      j=1
c
c     bit reversal section:
      do 11 i=1,n,2 
      if(j.gt.i)then
c      exchange two complex numbers:
       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=nn
1     if ((m.ge.2).and.(j.gt.m)) then
       j=j-m
       m=m/2
       goto 1
      endif
11    j=j+m
c
c     Danielson-Lanczos section:
      mmax=2 
c     Outer loop executed log2 nn times.
2     if (n.gt.mmax) then
       istep=2*mmax
c      Initialize for the trigonometric recurrence:
       theta=6.28318530717959d0/(isign*mmax) 
       wpr=-2.d0*sin(0.5d0*theta)**2 
       wpi=sin(theta)
       wr=1.d0
       wi=0.d0
c      Two nested inner loops:
       do 13 m=1,mmax,2
       do 12 i=m,n,istep
       j=i+mmax
c      Danielson-Lanczos formula:
       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
12     data(i+1)=data(i+1)+tempi
c      Trigonometric recurrence:
       wtemp=wr
       wr=wr*wpr-wi*wpi+wr
13     wi=wi*wpr+wtemp*wpi+wi
       mmax=istep
c      Not yet done:
       goto 2 
      endif 
c     All done:
      return
      END