C @(#)tsasin.for	17.1.1.1 (ES0-DMD) 01/25/02 17:20:37
C===========================================================================
C Copyright (C) 1995 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  @(#)tsasin.for	5.1 (ESO-IPG) 4/5/93 15:58:43
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT (c)  1992     European Southern Observatory and Warsaw Observatory
C.IDENT     tsasin.for
C.AUTHOR    Alex Schwarzenberg-Czerny, ESO and Warsaw Observatory
C.KEYWORD   MIDAS, time series, SINEFIT/TSA
C.LANGUAGE  FORTRAN 77
C.PURPOSE   Fit sine (Fourier) series
C.VERSION   0.0               June 1992
C.RETURNS   None
C.ENVIRON   TSA context
C-----------------------------------------------------------------------------
C
C
      INCLUDE 'MID_REL_INCL:TSA_DEF.INC'
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
C
      INTEGER      MXORDER
      PARAMETER   (MXORDER=20)          !  DEPENDS ON LENGTH OF OUTD
      REAL*8       FREQU                !  INITIAL PERIOD VALUE
      INTEGER      NORDER               !  NUMBER OF HARMONICS
      INTEGER      ITER                 !  NUMBER OF ITERATIONS
      CHARACTER*60 INAME                !  NAME OF INPUT TABLE
      CHARACTER*60 ONAME                !  NAME OF OUTPUT TABLE
      REAL*8       COEF(MXORDER+3)      !  NORDER COEFFICIENTS
C
      INTEGER      MODE,IACTS,KUN,KNUL
      INTEGER      TID,OTID,ITIME,IVAL
      INTEGER      NCOL,ICOL,NROW,IROW,ISOR
      INTEGER      LFIELD,TTYP,VTYP
      INTEGER      PTIME,PVAL,POTIME,POVAL
      CHARACTER*10 FORM
      CHARACTER*80 TEXT
      DOUBLE PRECISION E(MXORDER+3)           ! WORK TABLES
      DOUBLE PRECISION B(MXORDER+3,MXORDER+3) ! WORK TABLES
C
      INCLUDE 'MID_REL_INCL:TSA_DAT.INC'
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
      DATA MODE/1/
C
C   Get parameters
C
      CALL STSPRO ('tsasin')
      CALL STKRDC ('IN_A', 1,1,60,IACTS,INAME, KUN,KNUL,ISTAT)
      CALL STKRDC ('OUT_A',1,1,60,IACTS,ONAME, KUN,KNUL,ISTAT)
      CALL STKRDD ('STARTTSA', 1, 1,IACTS,FREQU, KUN,KNUL,ISTAT)
      CALL STKRDI ('ORDERTSA',  1, 1,IACTS,NORDER,KUN,KNUL,ISTAT)
      NORDER=2*NORDER+2
      IF (NORDER.GT.MXORDER+1) THEN
        NORDER=(MXORDER-1)/2*2+1
        WRITE(TEXT,'(A,I4)') 
     $   'Will use max. permitted ORDER of', 
     $    (MXORDER-1)/2
        CALL STTPUT(TEXT,ISTAT)
      ENDIF
      CALL STKRDI ('ITER',  1, 1, IACTS,ITER,KUN,KNUL,ISTAT)
      IF (ITER.LT.0) THEN
        CALL STKRDD ('OUTPUTD',1,NORDER-1,IACTS,COEF(1),
     $    KUN,KNUL,ISTAT)
      ENDIF
C
C   Map input/output data
C
      CALL TBTOPN (INAME,F_IO_MODE,TID,ISTAT)
      OTID=TID
      CALL TBIGET (TID,NCOL,NROW,ISOR,ICOL,IROW,ISTAT)
      CALL TBLSER (TID,'TIME',ITIME,ISTAT)
      IF (ITIME.LT.0) THEN
        CALL STETER(3,'Column :TIME not found')
      ENDIF
      CALL TBLSER (TID,'VALUE',IVAL,ISTAT)
      IF (IVAL.LT.0) THEN
        CALL STETER(4,'Column :VALUE not found')
      ENDIF
      CALL TBFGET (TID,ITIME,FORM,LFIELD,TTYP,ISTAT)
      CALL TBFGET (TID,IVAL, FORM,LFIELD,VTYP,ISTAT)
      CALL TBDGET (TID,ISTORE,ISTAT)
      IF (ISTORE.NE.F_TRANS) THEN
        TEXT='Input table '//INAME//' stored not transposed'
        CALL STETER(1,TEXT)
      ENDIF
      IF (TTYP.NE.D_R8_FORMAT.OR.VTYP.NE.D_R8_FORMAT) THEN
        CALL STETER(2,
     $    'Data column(s) must be of DOUBLE PRECISION type')
      ENDIF
      CALL TBCMAP (TID,ITIME,PTIME,ISTAT)
      POTIME=PTIME
      CALL TBCMAP (TID,IVAL, PVAL, ISTAT)
      POVAL=PVAL
C
C   Map output data
C
C M. Peron 22-02-93
C A. Schwarzenberg-Czerny 07-05-93
C
C   Perform single nonlinear fit iteration, store coeficients and residuals
C
      COEF(NORDER)=FREQU
      CALL TIMFIT(MADRID(PTIME),MADRID(PVAL),MADRID(POVAL),NROW,
     $   MODE,NORDER,ITER,MXORDER,COEF,E,B)
      FREQU=COEF(NORDER)
      CALL STKWRD ('STARTTSA',  FREQU,  1,       1,KUN,ISTAT)
      CALL STKWRD ('OUTPUTD',COEF(1),1,NORDER-1,KUN,ISTAT)
C
C   Wind-up
C
      CALL DSCUPT(OTID,OTID,' ',ISTAT)
      CALL STSEPI
C
      END
C
C
C
C
C
C
C                 Subroutine  FIT
C
C     Purpose: Given a trial period, fit a Fourier series by non-linear
C        least squares method.
C
C     Data: an ASCII file containing 3 columns (number of line(integer),
C        time, value (both real*4)) and less than maxm-2 rows.
C        A trial value of the period, an order of the series (<maxo/2-1)
C        and a name of the data file are prompted for.
C
C     Results: Correlation coefficients of the fitted parameters, their
C        values and errors are displayed on the screen after each
C        iteration. Finish iterations by typing 0 when changes in
C        the parameters and std.dev. are negligible. After the final
C        iteration values of the series are subtracted from the data and
C        output as fitsin.res file.
C
C     Method: The data are fitted by the following series:
C        f(1)+f(2)*sin(w*t)+f(3)*cos(w*t)+...+f(2*ord+1)*cos(2*ord*w*t)
C        by least squares. After the initial solution
C        the value of the frequency w is varied as well as the coefficients so
C        that due to nonlinearity of the problem several iterations are
C        required. To continue iterations type 1.
C
C     A. Schwarzenberg-Czerny                      June, 1987
C
      SUBROUTINE TIMFIT(T,X,OX,NX,MODE,NORD,NITER,MXORD,F,E,B)
C
      INCLUDE 'MID_REL_INCL:TSA_DEF.INC'
C

      INTEGER NX                 !  NUMBER OF DATA
      DOUBLE PRECISION  CRACOW             !  EXTERNAL FUNCTION
      DOUBLE PRECISION  T(NX)              !  TIME VECTOR
      DOUBLE PRECISION  X(NX)              !  INPUT DATA VECTOR
      DOUBLE PRECISION  OX(NX)             !  OUTPUT DATA VECTOR
      INTEGER MODE               !  MODE OF WORK
      INTEGER NORD               !  ORDER OF THE SERIES
      INTEGER MXORD              !  TABLE SIZE-3
      INTEGER NITER              !  NUMBER OF ITERATIONS
      DOUBLE PRECISION  F(MXORD+3)         !  FIT COEFFICIENTS (INCL. FREQ)
      DOUBLE PRECISION E(MXORD+3),B(MXORD+3,MXORD+3)   ! WORK TABLES
C
      DOUBLE PRECISION  FREQ
      LOGICAL START
      CHARACTER*80 MSG
      DOUBLE PRECISION TM,XM,TMX,TMI,DFREQ,PI2,CORRMX
      DOUBLE PRECISION PHMT,DPHMT,EPMAX, FREQSTEP
      INTEGER N,NP,I,N1,J,K,ITER
C
      INCLUDE 'MID_REL_INCL:TSA_CONST.INC'
C
C    get data
C
      PI2=TWO*TWO*TWO*ATAN(ONE)
      START=NITER.GT.0
      NITER=ABS(NITER)
      ITER=0
      TM=ZERO
      XM=ZERO
      TMX=T(1)
      TMI=TMX
      DO 1 I=1,NX
        TMX=MAX(TMX,T(I))
        TMI=MIN(TMI,T(I))
        TM=TM+T(I)
 1      XM=XM+X(I)
      TM=TM/NX
      XM=XM/NX
C
C    input parameters and output results of previous iteration
C
      IF (MODE.EQ.1) THEN
        IF (NORD.LT.1) CALL STETER(5,'Wrong parameter')
        F(NORD)=PI2*F(NORD)
        NP=2
        IF (START) THEN
           WRITE(MSG,'(a,1pg20.6,a,i4)')
     $         ' START: frequency=',F(NORD)/PI2,' order=',(NORD-2)/2
           CALL STTPUT(MSG,ISTAT)
           N=NORD-1
           DFREQ=ZERO
        ELSE
          N=NORD
        ENDIF
      ELSE
        CALL STETER(2,' Wrong mode')
      ENDIF
      IF (START) THEN
        DO 5 I=1,N
 5         F(I)=ZERO
      ENDIF
C
C    iterations for non-linear fit
C
 100  CONTINUE
C
C           freq - current frequency
C           dfreq - its error
C           tm - mean time of obserwation
C           phmt - phase of base sinusoid at mean time
C           dphmt - its error
C           epmax - epoch of maximum of the base sinusoid
C           The harmonics are ignored by phase calculation
C
C
C     set condition and normal equations
C
        N1=N+1
        ITER=ITER+1
        DO 60 I=1,N1
          DO 60 J=1,N1
 60         B(J,I)=ZERO
        DO 3 I=1,NX
          F(NORD+1)=T(I)-TM
          F(NORD+2)=X(I)-XM
          IF (MODE.EQ.1) THEN
            CALL FITSIN(E,N1,F,NP)
          ENDIF
C  *** prevent overwriting data before end of iterations, A.S-C, 07-05-93
          IF (ITER.GT.NITER) THEN
            OX(I)=E(N1)
          ENDIF
          DO 3 J=1,N1
            DO 3 K=1,N1
 3            B(K,J)=B(K,J)+E(J)*E(K)
C
C     solve normal equations and evaluate results
C
        IF (ITER.LE.NITER) THEN
          IF (CRACOW(B,N1,MXORD+3,NX).NE.ZERO) THEN
C
C   Update and output coefficients
C
            CORRMX=ZERO
            DO 7 J=1,N
              F(J)=F(J)+B(J,N1)
              DO 7 I=1,J-1
 7              CORRMX=MAX(CORRMX,ABS(B(I,J)/SQRT(B(J,J)*B(I,I))))
            IF (ITER.LE.1) THEN
              CALL STTPUT(
     $        '   iter.        freq.       std.dev.,     max. corr.',
     $          ISTAT)
            ELSE
              FREQSTEP=B(NORD,N1)/PI2
              IF (FREQSTEP*(TMX-TMI).GT.TWO) THEN
                CALL STTPUT(
     $            '*** Large frequency step, divergence?',ISTAT)
              ENDIF
            ENDIF
            WRITE(MSG,'(1x,i7,3(1pe15.7),'' '')')
     $         ITER,F(NORD)/PI2,SQRT(MAX(B(N1,N1),ZERO)),CORRMX
            CALL STTPUT(MSG,ISTAT)
            IF (ITER.EQ.NITER) THEN
              CALL STTPUT('   coefficients     errors',ISTAT)
              F(1)=F(1)+XM
              DO 6 J=1,N
                WRITE(MSG,'(1x,2(1pe15.7),'' '')')
     $                F(J),SQRT(MAX(B(J,J),ZERO))
 6            CALL STTPUT(MSG,ISTAT)
              F(1)=F(1)-XM
C
C   Output ephemeris for base sinusoid
C
              IF (NORD.GT.ZERO.AND.MODE.EQ.1) THEN
                FREQ=F(NORD)/PI2
                IF(.NOT.START) DFREQ=SQRT(MAX(B(N,N),ZERO))/PI2
                PHMT=ATAN2(F(3),F(2))/PI2
                DPHMT=SQRT(F(2)*F(2)*B(3,3)+F(3)*F(3)*B(2,2))/
     $             (F(2)*F(2)+F(3)*F(3))
                EPMAX=(HALF*HALF-PHMT)/FREQ+TM
                IF (PHMT.LT.ZERO) PHMT=PHMT+ONE
                MSG='   frequency        error         '//
     $            'phase at m.t.  error        '
                CALL STTPUT(MSG,ISTAT)
                WRITE(MSG,'(1x,5(1pe15.7))') FREQ,DFREQ,PHMT,DPHMT
                CALL STTPUT(MSG,ISTAT)
                CALL STTPUT(
     $          '   epoch of max.    mean time    ',ISTAT)
                WRITE(MSG,'(1x,5(1pe15.7))') EPMAX,TM
                CALL STTPUT(MSG,ISTAT)
              ENDIF
            ENDIF
C
C      next iterations with period correction
            N=NORD
            START=.FALSE.
            GOTO 100
          ELSE
            CALL STETER(4,
     $          ' Singular equations, or wrong dimensions ')
          ENDIF
        ENDIF
C        end loop 100
      IF (MODE.EQ.1) THEN
        F(NORD)=F(NORD)/PI2
      ENDIF
      END
C
C
C
      SUBROUTINE FITSIN(E,N1,P,NP)
C
C    condition equations for fit of Fourier series,
C    last term contains residuals
C
      INCLUDE 'MID_REL_INCL:TSA_DEF.INC'
C
      INTEGER N1,NP
      DOUBLE PRECISION E(N1),P(N1-1+NP)
C
      INTEGER J,NORD,NFR,L
      REAL SJ,CJ,T,X,FR,S,C,RHS,SAVE
C
      INCLUDE 'MID_REL_INCL:TSA_CONST.INC'
      NORD=(N1-2)/2
      NFR=2*(NORD+1)
      T=P(NFR+1)
      X=P(NFR+2)
      FR=P(NFR)
      S=SIN(T*FR)
      C=COS(T*FR)
      SJ=S
      CJ=C
      E(NFR)=ZERO
      E(1)=ONE
      RHS=X-P(1)
      DO 4 J=1,NORD
        L=2*J
        E(L)=SJ
        E(L+1)=CJ
        E(NFR)=E(NFR)+J*T*(P(L)*CJ-P(L+1)*SJ)
        RHS=RHS-P(L)*SJ-P(L+1)*CJ
        SAVE=CJ*C-SJ*S
        SJ=CJ*S+SJ*C
 4      CJ=SAVE
      E(N1)=RHS
      END
C
C
C
C
C
C
C
C
C