C @(#)timdft.for	17.1.1.1 (ES0-DMD) 01/25/02 17:20:21
C===========================================================================
C Copyright (C) 1996 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
C  @(#)timdft.for	5.1 (ESO-IPG) 28/07/96 15:58:35
C
C   Fourier and normalized Scargle power spectrum
C   (see ApJ 338,277 & ApJ 263,835)
C                               alex schwarzenberg-czerny, JULY 1996
c Tested OK 28.07.1996
C
      SUBROUTINE TIMDFT(T,X,N,P,NFREQ,CMS,CPS,CD,SD,
     $   FREQIN,FSTEP,SCARGL)
C
      INCLUDE 'MID_REL_INCL:TSA_DEF.INC'
      LOGICAL SCARGL
      INTEGER N,NFREQ
      DOUBLE PRECISION  T(N),X(N),CMS(N),CPS(N),CD(N),SD(N),
     $   P(NFREQ,2),FREQIN,FSTEP
C
      LOGICAL NEXT
      INTEGER I,IF
      DOUBLE PRECISION PI2,START,STEP,EN,ANORM,
     $   C,S,XC,XS,A,B,E,F,TERM,EPS
      PARAMETER (EPS=1D-8)
      INCLUDE 'MID_REL_INCL:TSA_CONST.INC'
C  Set up trig functions recurrence
      PI2=   ABS(ATAN2(ZERO,-ONE))*TWO
      START= FREQIN*PI2
      STEP=  FSTEP*PI2
      EN=    N
      XC=    ZERO
      XS=    ZERO
C
      DO 1 I=1,N
        TERM=  T(I)*START
        C=     COS(TERM)
        S=     SIN(TERM)
        CMS(I)=C-S
        CPS(I)=C+S
        TERM=  T(I)*STEP
        S=     SIN(TERM*HALF)
        CD(I)=-TWO*S*S
        SD(I)= SIN(TERM)
        TERM=X(I)
        XC=    TERM+XC
        XS=    TERM*TERM+XS
 1    CONTINUE
C  Normalization
      XC= XC/EN
      XS= (XS-EN*XC*XC)/(EN-ONE)
      IF (SCARGL) THEN
        ANORM= HALF*HALF/XS/EN
      ELSE
        ANORM= HALF/EN/EN
        TERM=  ONE
      ENDIF
C  Compute periodogramme
      DO 2 IF=1,NFREQ
        NEXT=.TRUE.
 4      XC=ZERO
        XS=ZERO
        C= ZERO
        S= ZERO
        DO 3 I=1,N
          A=CMS(I)
          B=CPS(I)
          IF (NEXT) THEN
C Next frequency or ...
            CMS(I)=A*CD(I)-B*SD(I)+A
            CPS(I)=B*CD(I)+A*SD(I)+B
          ELSE
C ... recompute the spectrum at the frequency near the singularity
            E=0.9999*STEP*T(I)
            F=SIN(E)
            E=SIN(E*HALF)
            E=-TWO*E*E
            TERM=A*E+B*F+A
            B=B*E-A*F+B
            A=TERM
          ENDIF
          XC=X(I)*A+XC
          XS=X(I)*B+XS
          C=A*B+C
          S=B*B-ONE+S
 3      CONTINUE          
C  Lomb-Scargle statistics
        IF (SCARGL) THEN
          TERM=SQRT(C*C+S*S)
          C=C/TERM
          S=S/TERM
          E=SQRT((ABS(C)+ONE)*HALF)
          IF (C.GE.ZERO) THEN
            S=S/E*HALF
            C=E
          ELSE
            C=S/E*HALF
            S=E
          ENDIF
          TERM=ONE-TERM/EN
C
          IF (NEXT.AND.TERM.LT.EPS) THEN
            NEXT=.FALSE.
            GOTO 4
          ENDIF
          A=XS+XC
          B=XS-XC
          XC=A*C+B*S
          XS=B*C-A*S
        ENDIF
        P(IF,1)=(XS*XS/TERM+XC*XC/(TWO-TERM))*ANORM
 2    CONTINUE
C  Spectral window
      IF (.NOT.SCARGL) THEN
        DO 11 I=1,N
          CMS(I)=ONE
          CPS(I)=ZERO
 11     CONTINUE
        DO 12 IF=1,NFREQ
          C=ZERO
          S=ZERO
          DO 13 I=1,N
            A=      CMS(I)
            B=      CPS(I)
            CMS(I)= A*CD(I)-B*SD(I)+A
            CPS(I)= B*CD(I)+A*SD(I)+B
            C=      A+C
            S=      B+S
 13       CONTINUE          
          P(IF,2)=  (C*C+S*S)/EN
 12     CONTINUE
      ENDIF
      END
c