C dft - no number of points at beginning
c  Terence James Deeming  (1937-1992) Fourier Transform
C       
C ... AUTHOR: KEPLER
C
C ... prewheightening included on 05/17/86
c ... modified on 01.20.95 to normalize by average
c ...    includes dmod(twopi) factor

        real*8 t(130000)
        real*8 x(130000)
       character*20 infile, outfile
        real*8 f,fr,fi,a,c,s,ff,amplit,twopi,pi,tt,t1
        real*8 fbeg,fend,fstp,yaver,amp(20),fq(20),ph(20)  
        integer*4 i,num
C
      DATA N/0/,FRQ/0./,NPWT/0/,ITRIG/0/,i/0/ 
      pi=dacos(-1.0d0)
      TWOPI=2.0d0*PI 
      DELAY=0.0d0 
C
       ndata=130000
       write(*,600)
600       format(' Enter freq_beg, freq_end,freq_step: ',$)
      READ(*,*) fbeg,fend,fstp
       write(*,610)
610       format(' Enter name of data file (N<130000): ',$)
       read(*,620) infile
620       format(a20)
       open(1,file=infile,status='old')
       write(*,630)
630       format(' Enter name of amp spectra output file: ',$)
       read(*,620) outfile
       open(2,file=outfile,status='unknown')
	rewind(2)
C
C... asks for prewheitening
C
        write(*,308)
308     format(' how many sine curves do you want to remove from data?')
        read(*,*)itrig
        if(itrig.eq.0) goto 19
      DO 310 J=1,ITRIG
        write(*,309)j
309     format(1x,i3,' enter amplitude, period, tmax(sec)')
        READ (*,*) AMP(j),period,PH(j) 
310     fq(j)=1.0d0/period
C
C...  input file
19      READ(1,*)t1,X(1)
        t(1)=0.0d0
        do 26 i=2,ndata
        READ(1,*,end=28)tt,X(i)
        t(i)=tt-t1
26      continue
28        num=i-1
        CLOSE(1)
        write(*,31)NUM
31      format(' finished readind data, npoints= ',i6)

        numf=(fend-fbeg)/fstp
C     WRITE(2,*) numf, fstp
                             
C...here makes NORMALISES by average and make avergage = 0
      YAVER = 0.0d0
      DO 4 I=1,NUM
4        YAVER = YAVER + X(I)
      YAVER = YAVER / dble(NUM)
	write(*,41)yaver
41	format(' yaver=',e12.6)
      DO 5 I=1,NUM
5       X(I) = X(I)/YAVER - 1.0

C
C ... HERE IS WHERE IT PREWHITENS THE DATA.  
C
      IF(ITRIG.EQ.0)GOTO 46 
      DO 330 I=1,NUM
      DO 320 J=1,ITRIG
320   X(I)=X(I)-AMP(j)*dcos(TWOPI*FQ(j)*(t(i)-PH(j)))  
330     continue

46     CONTINUE 

36    F=FBEG
      N=0 
45    FR=0.0d0 
      FI=0.0d0 
        DO 50 I=1,NUM 
        A=TWOPI*F*T(I)
        a=dmod(a,twopi)
        C=DCOS(A)
        S=DSIN(A)
        FR=FR+X(I)*C
50      FI=FI+X(I)*S
        FR=FR/dble(NUM)
        FI=FI/dble(NUM)
        FF=FR*FR+FI*FI
c  amplitude spectra in file outfile
c  correction on amplitude of factor of 2
      AMPLIT = 2.0d0*DSQRT(FF)   
      WRITE(2,987)F,AMPLIT
987	format(f13.9,1x,1pe10.4)
C        write(*,51)n
51      format(' finished freq. step: ',i4)
52    N=N+1 
53    F=F+FSTP
C
C ... LOOP BACK UP AND DO IT AGAIN UNTIL WE'VE REACHED FEND=PMIN
C
      IF( (F - FEND) .LE. 0.0 ) GOTO 45 
        close(2)
      STOP
      END 
