PROGRAM FFT
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 PI
      PARAMETER (N0=4)
      PARAMETER (N1=2**N0,PI=3.1415926535898D0)
C The N1 by 2 array P contains the complex input data
C P(K,1) is the real part of data point K
C P(K,2) is the imaginary part of data point K
      DIMENSION P(N1,2)
C Use ISW to select fourier transformation or back transformation
      ISW=1
C ISW=1  FFT, ISW=-1 BACK TRANSFORMATION
      R2=PI/N1
      IF(ISW.GT.0)GOTO 2120
      R2=-R2
 2120 N6=N1+1
      N7=2*N1
      DO 2350 N2=1,N0
      N7=N7/2
      N8=N7/2
      R2=R2*2
      NN=N6-N7
      N3=1
 2181 R3=-R2
      N5=N3-1
      DO 2330 N4=1,N8
      R3=R3+R2
      N9=N5+N4
      N10=N9+N8
C In a real application, these sin and cos values come from a table.
      R8=DCOS(R3)
      R6=DSIN(R3)
      R7=R8*(P(N9,1)-P(N10,1))-R6*(P(N9,2)-P(N10,2))
      R6=R6*(P(N9,1)-P(N10,1))+R8*(P(N9,2)-P(N10,2))
      P(N9,1)=P(N9,1)+P(N10,1)
      P(N9,2)=P(N9,2)+P(N10,2)
      P(N10,1)=R7
      P(N10,2)=R6
 2330 CONTINUE
      N3=N3+N7
      IF(N3.LT.NN+N7)GOTO 2181
 2340 CONTINUE
 2350 CONTINUE
C The actual FFT is already computed, but the output is in a
C scrambled order. Unscramble:
      DO 2540 N2=1,N1
      N9=N2-1
      N10=0
      DO 2450 N3=1,N0
      N10=2*N10
      R1=N9
      R1=R1/2
      N5=R1
      R1=R1-N5
      IF(ABS(R1).LT..1/N1)GOTO 2240
      N10=N10+1
 2240 N9=N5
 2450 CONTINUE
      N10=N10+1
      IF(N10.LE.N2)GOTO 2540
      R7=P(N2,1)
      R6=P(N2,2)
      P(N2,1)=P(N10,1)
      P(N2,2)=P(N10,2)
      P(N10,1)=R7
      P(N10,2)=R6
 2540 CONTINUE
C The conventional FFT stops here, but if we want frequencies to
C start with the most negative in the first point, then going downwards
C through zero and stop at the highest positive frequency we have
C to rearrange a bit more.
      JJ=N1/2
      N3=JJ
      DO 2700 N2=1,JJ
      N3=N3+1
      R4=P(N2,1)
      R5=P(N2,2)
      P(N2,1)=P(N3,1)
      P(N2,2)=P(N3,2)
      P(N3,1)=R4
      P(N3,2)=R5
 2700 CONTINUE
      R4=P(1,1)
      R5=P(1,2)
      DO 2800 N2=2,N1
      N3=N2-1
      P(N3,1)=P(N2,1)
 2800 P(N3,2)=P(N2,2)
      P(N1,1)=R4
      P(N1,2)=R5
C Now the FFT is complete and in correct order.
C The power spectrum is an array of length N1.
C POWER(K)=P(K,1)**2+P(K,2)**2
C This is the quantity to average over time when looking for weak signals.
      STOP
      END