C @(#)tsaort.for	17.1.1.1 (ESO-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 Massachusetts Ave, Cambridge, 
C MA 02139, USA.
C
C Correspondence 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  @(#)tsaort.for	5.1 (ESO-IPG) 4/5/93 15:58:45
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT (c)  1996     European Southern Observatory and Copernicus Center
C.IDENT     tsaort.for
C.AUTHOR    Alex Schwarzenberg-Czerny, ESO and Copernicus Center, Warsaw
C.KEYWORD   MIDAS, time series, ORT/TSA
C.LANGUAGE  FORTRAN 77
C.PURPOSE   Compute multiharmonic Fourier spectrum by orthogonal projections
C.VERSION   0.0               June 1996
C.RETURNS   None
C.ENVIRON   TIME/ANAL context
C-----------------------------------------------------------------------------
C
C
      INTEGER      NAXIS
      PARAMETER    (NAXIS=2)
      INTEGER      LENGTH(NAXIS)        !  NUMBER OF FREQUENCIES
      INTEGER      ORDER                !  NUMBER OF BINS
      INTEGER      MODE                
      INTEGER      IACTS,KUN,KNUL
      INTEGER      TID,ITIME,IVAL
      INTEGER      NCOL,ICOL,NROW,IROW,ISOR
      INTEGER      LFIELD,TTYP,VTYP
      INTEGER      PTIME,PVAL,PPER,IDPER
      INTEGER      ICZ,CZ,ICZN,CZN,ICCF,CCF,ICP,CP
      INTEGER      ICDZ,CDZ,ICDCF,CDCF,NOBS2,ASIZE
C 
      CHARACTER*10 FORM
      CHARACTER*80 TEXT
      CHARACTER*60 INAME                !  NAME OF INPUT TABLE
      CHARACTER*60 ONAME                !  NAME OF OUTPUT IMAGE
C
      DOUBLE PRECISION  START(NAXIS)         !  START FREQUENCY
      DOUBLE PRECISION  STEP(NAXIS)          !  STEP IN FREQUENCY
      DOUBLE PRECISION  OM
C
      INCLUDE 'MID_REL_INCL:TSA_DEF.INC'
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
C
      INCLUDE 'MID_REL_INCL:TSA_DAT.INC'
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
C   Get parameters
C
      CALL STSPRO ('tsaort')
      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,START(1) ,KUN,KNUL,ISTAT)
      CALL STKRDD ('STEPTSA',
     $     1, 1,IACTS,STEP(1)  ,KUN,KNUL,ISTAT)
      CALL STKRDI ('NSTEPS',  1, 1,IACTS,LENGTH(1),KUN,KNUL,ISTAT)
      CALL STKRDI ('ORDERTSA',
     $     1, 1,IACTS,ORDER, KUN,KNUL,ISTAT)
      MODE=1
      START(NAXIS)=  ZERO
      STEP(NAXIS)=   ONE
      LENGTH(NAXIS)= 2
C
C   Map input data
C
      CALL TBTOPN (INAME,F_I_MODE,TID,ISTAT)
      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(3,'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(2,TEXT)
      ENDIF
      IF (TTYP.NE.D_R8_FORMAT.OR.VTYP.NE.D_R8_FORMAT) THEN
        CALL STETER(1,'Column(s) must be of DOUBLE PRECISION type')
      ENDIF
      CALL TBCMAP (TID,ITIME,PTIME,ISTAT)
      CALL TBCMAP (TID,IVAL, PVAL, ISTAT)
C
C   Create and map work tables  Z,ZN,CF,P,DZ,DCF
C
      NOBS2=NROW+NROW
      CALL STFCRE('ZZMIDAWORK',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,
     $   NOBS2,ICZ,  ISTAT)
      CALL STFMAP(ICZ,  F_X_MODE,1,NOBS2,ASIZE,CZ,  ISTAT)
C
      CALL STFCRE('ZZMIDDWORK',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,
     $   NOBS2,ICZN, ISTAT)
      CALL STFMAP(ICZN, F_X_MODE,1,NOBS2,ASIZE,CZN, ISTAT)
C
      CALL STFCRE('ZZMIDEWORK',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,
     $   NOBS2,ICCF, ISTAT)
      CALL STFMAP(ICCF, F_X_MODE,1,NOBS2,ASIZE,CCF, ISTAT)
C
      CALL STFCRE('ZZMIDFWORK',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,
     $   NOBS2,ICP,  ISTAT)
      CALL STFMAP(ICP,  F_X_MODE,1,NOBS2,ASIZE,CP,  ISTAT)
C
      CALL STFCRE('ZZMIDGWORK',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,
     $   NOBS2,ICDZ, ISTAT)
      CALL STFMAP(ICDZ, F_X_MODE,1,NOBS2,ASIZE,CDZ, ISTAT)
C
      CALL STFCRE('ZZMIDHWORK',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,
     $   NOBS2,ICDCF,ISTAT)
      CALL STFMAP(ICDCF,F_X_MODE,1,NOBS2,ASIZE,CDCF,ISTAT)
C
C   Map output data
C
      CALL STIPUT (ONAME,D_R8_FORMAT,F_IO_MODE,
     $   F_IMA_TYPE,NAXIS,LENGTH,
     $   START,STEP,ONAME,'AXIS: 1/TIME DATA: UNITLESS',
     $   PPER,IDPER,ISTAT)
C
C   Compute multiharmonic Fourier periodogramme
C
      CALL TIMORT(NROW,MADRID(PTIME),MADRID(PVAL),LENGTH(1),START(1),
     $   STEP(1),ORDER,MODE,MADRID(PPER),OM,MADRID(CZ),MADRID(CZN),
     $   MADRID(CCF),MADRID(CP),MADRID(CDZ),MADRID(CDCF))
C
C   Wind-up
C
      CALL DSCUPT(IDPER,IDPER,' ',ISTAT)
      CALL STSEPI
C
      END
C
C
C
      subroutine TIMORT(kk,t,f,ll,fr0,frs,n2,mode,th,om,
     $    z,zn,cf,p,dz,dcf)
c
c   Fast Fit Periodogramme
c
c   Method: A. Schwarzenberg-Czerny, 
c       Ap.J. 460, L107, 1996 April 1
c
c   Input:
c     kk - number of observations
c     t(kk),f(kk) - times and values of observations
c     ll - number of frequencies
c     fr0 - first frequency
c     frs - frequency step
c     n2 - number of fitted harmonics
c     mode - mode of operation
c       1 - AOV periodogramme
c       2 - Chi2 periodogramme
c       3 - values in f(kk) replaced by residua, AOV periodogramme
c       4 - harmonics amplitudes and AOV periodogramme
c
c   Output
c     th(ll) - Periodogramme (type depends on 'mode')
c     om - frequency of maximum power
c
C
      INCLUDE  'MID_REL_INCL:TSA_DEF.INC'
C
      INTEGER  N2                       ! ORDER (NUMBER OF HARMONICS)
      INTEGER  KK                       ! NUMBER OF OBSERVATIONS
      INTEGER  LL                       ! NUMBER OF FREQUENCIES
      INTEGER  MODE                     ! MODE OF OPERATION

C is not really used...      INTEGER  NORDER

      integer k,n,nn,l,lmax,mtrig,ksign
C
      INTEGER  MN
      PARAMETER(MN=100)
C 
      DOUBLE PRECISION   T(KK),F(KK)    ! TIMES AND VALUES OF OBSERVATIONS
      DOUBLE PRECISION   TH(LL,2)       ! PERIODOGRAMME
      DOUBLE PRECISION   FR0,FRS,OM     ! FREQUENCY START AND INCREMENT
      DOUBLE PRECISION pi2,sav,df1,df2,df3,s2,thmax,tol,sn
C 
      COMPLEX*16 Z(KK),ZN(KK),CF(KK),P(KK),DZ(KK),DCF(KK) ! WORK TABLES
      complex*16 a,c,al,san,sad,sc,pc,amx,cmx
      dimension a(0:mn),c(0:mn),amx(0:mn),cmx(0:mn)
C
      CHARACTER*80 TEXT
C 
      logical test
c
      REAL   PROB
      REAL riembi
      EXTERNAL RIEMBI
C 
      parameter (mtrig=1,tol=1e-6)
C
      INCLUDE 'MID_REL_INCL:TSA_CONST.INC'
C
      CALL STTPUT('Mapped arrays',ISTAT)
C ...      IF (NORDER.GT.MN) THEN
C ...        WRITE(TEXT,'(A,I3)')
C ...     1     'WARNING: Used maximum allowed order of ',MN
C ...        CALL STTPUT(TEXT,ISTAT)
C ...        NORDER=MN
C ...      ENDIF

      nn=n2+n2
      if (nn.gt.mn) then
        CALL STTPUT(' TIMORT: Wrong value of ORDERTSA ',ISTAT)
        return
      endif
      pi2=two*atan2(zero,-one)
      df3=kk-1
      df1=nn+1
      df2=kk-df1
      thmax=zero
      s2=zero
      om=pi2*frs
      do 1 k=1,kk
        sav=om*t(k)
        dz(k)=cmplx(-two*sin(half*sav)**2,sin(sav))
        sav=n2*sav
        dcf(k)=cmplx(-two*sin(half*sav)**2,sin(sav))
        s2=s2+f(k)*f(k)
 1    continue

      do 10 l=1,ll
        th(l,1)=zero
        al=zero
        om=pi2*(fr0+(l-1)*frs)
        test=mod(l-1,mtrig).eq.0
        do 11 k=1,kk
          if(test) then
            sav=om*t(k)
            z(k)=cmplx(cos(sav),sin(sav))
            sav=n2*sav
            cf(k)=f(k)*cmplx(cos(sav),sin(sav))
          else
            z(k)=z(k)*dz(k)+z(k)
            cf(k)=cf(k)*dcf(k)+cf(k)
          endif
          zn(k)=conjg(z(k))
          p(k)=zn(k)
 11     continue

        do 12 n=0,nn
          san=zero
          sad=zero
          sn=zero
          sc=zero

          do 13 k=1,kk
            pc=p(k)
            pc=z(k)*pc-al*zn(k)*conjg(pc)
            p(k)=pc
            pc=conjg(pc)
            sav=pc*conjg(pc)
            sn=sn+sav
            sc=sc+cf(k)*pc
            zn(k)=zn(k)*z(k)
            sad=sad+zn(k)*pc*pc
            san=san+z(k)*sav
 13       continue
          if (min(abs(sad),abs(sn)).le.tol) then
            a(n)=zero
            c(n)=zero
          else
            a(n)=san/sad
            c(n)=sc/sn
          endif
          al=a(n)
          pc=c(n)*sqrt(sn)
          th(l,1)=th(l,1)+pc*conjg(pc)
 12     continue

        if (abs(a(nn)).eq.zero) then
          write(text,*) 'TIMORT: Singular at frequency',om/pi2,' '
          CALL STTPUT(TEXT,ISTAT)
        endif
        if (th(l,1).ge.thmax) then
          thmax=th(l,1)
          lmax=l
          do 14 n=0,nn
            amx(n)=a(n)
            cmx(n)=c(n)
 14       continue
        endif
        if (mode.eq.2) then
          th(l,1)=th(l,1)*df3/(df1*s2)
        else
          sav=(s2-th(l,1))/df2
          th(l,1)=th(l,1)/(df1*sav)
          th(l,2)=sav
        endif
 10   continue
      
      if (mode.eq.2) then
        prob=exp(-th(lmax,l))
      else 
        prob=riembi(REAL(half*df2),REAL(half*df1),
     $          REAL(df2/(df2+df1*th(lmax,1))))
      endif
      om=pi2*(fr0+(lmax-1)*frs)
 
      if (mode.eq.3) then
        ksign=0
        sn=zero
        do 21 k=1,kk
          sav=om*t(k)
          san=cmplx(cos(sav),sin(sav))
          sad=one
          pc=one
          sc=zero
          do 22 n=0,nn
            sc=cmx(n)*pc+sc
            pc=san*pc-amx(n)*sad*conjg(pc)
            sad=sad*san
 22       continue
          sav=-n2*sav
          f(k)=f(k)-sc*cmplx(cos(sav),sin(sav))
          sn=sn+f(k)*f(k)
 21     continue
        s2=abs(s2-sn-thmax)/s2
        if (s2.gt.tol) then
          write(text,*) 'TIMORT: Rounding error',s2,' '
          CALL STTPUT(TEXT,ISTAT)
        endif
      endif

      if (mode.eq.4) then
       
        CALL STTPUT('TIMORT: Power step expansion unimplemented yet '
     $     ,ISTAT)
      endif

      write(text,'(a,1pe15.6,2e12.3,0p)') 
     $   'TIMORT:max: f, th, prob=',om/pi2,th(lmax,1),prob
        CALL STTPUT(TEXT,ISTAT)

      end