C @(#)tsacov.for	17.1.1.1 (ES0-DMD) 01/25/02 17:20:36
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  @(#)tsacov.for	5.1 (ESO-IPG) 4/5/93 15:58:50
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT (c)  1992     European Southern Observatory and Warsaw Observatory
C.IDENT     tsacov.for
C.AUTHOR    Alex Schwarzenberg-Czerny, ESO and Warsaw Observatory
C.KEYWORD   MIDAS, time series, COVAR/TSA
C.LANGUAGE  FORTRAN 77
C.PURPOSE   Compute discrete covariance function
C           Reference: Edelson R.A. and Krolik, J.H., 1988, ApJ 333, 646
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
      CHARACTER*60 INAME1               !  NAME OF 1ST OBSERVATION TABLE
      CHARACTER*60 INAME2               !  NAME OF 2ND OBSERVATION TABLE
      CHARACTER*60 ONAME                !  NAME OF OUTPUT COVARIANCE TABLE
      REAL*8       TLAG1                !  START TIME LAG
      REAL*8       DTLAG                !  TIME LAG INCREMENT
      INTEGER      NLAG                 !  NUMBER OF TIME LAGS
C
      INTEGER      NOBS1, NOBS2
      INTEGER      ISCALE,IFUNCT,IACTS, KUN,   KNUL
      INTEGER      TID1,  TID2,  OTID,  ITIME, IDAT, IVAR
      INTEGER      NCOL,  ICOL,  NROW,  ISOR
      INTEGER      LFIELD,TTYP,  DTYP,  VTYP
      INTEGER      POBS1T,POBS1D,POBS1V
      INTEGER      POBS2T,POBS2D,POBS2V
      INTEGER      ICOVN, ICOVC, ICOVV
      INTEGER      PCOVN, PCOVC, PCOVV
      DOUBLE PRECISION       AV1,VAR1,AV2,VAR2
      CHARACTER*3  CSCALE
      CHARACTER*1  CFUNCT
      CHARACTER*10 FORM
      CHARACTER*80 TEXT
C
      INCLUDE 'MID_REL_INCL:TSA_DAT.INC'
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
C
C   Get parameters
C
      CALL STSPRO ('tsacov')
      CALL STKRDC ('IN_A',   1,1,60,IACTS,INAME1,KUN,KNUL,ISTAT)
      CALL STKRDC ('IN_B',   1,1,60,IACTS,INAME2,KUN,KNUL,ISTAT)
      CALL STKRDC ('OUT_A',  1,1,60,IACTS,ONAME, KUN,KNUL,ISTAT)
      CALL STKRDD ('STARTTSA',    1, 1,IACTS,TLAG1, KUN,KNUL,ISTAT)
      CALL STKRDD ('STEPTSA',     1, 1,IACTS,DTLAG, KUN,KNUL,ISTAT)
      CALL STKRDI ('NSTEPS',   1, 1,IACTS,NLAG,  KUN,KNUL,ISTAT)
      CALL STKRDC ('CSCALE', 1,1, 3,IACTS,CSCALE,KUN,KNUL,ISTAT)
      IF     (LLE(CSCALE,'LIN').AND.LGE(CSCALE,'LIN')) THEN
        ISCALE=1
      ELSEIF (LLE(CSCALE,'LOG').AND.LGE(CSCALE,'LOG')) THEN
        ISCALE=2
      ELSE
        CALL STETER(6,'SCALE must be LIN or LOG')
      ENDIF
      CALL STKRDC ('CFUNCT', 1,1, 1,IACTS,CFUNCT,KUN,KNUL,ISTAT)
      IF     (LLE(CFUNCT,'C').AND.LGE(CFUNCT,'C')) THEN
        IFUNCT=1
      ELSEIF (LLE(CFUNCT,'S').AND.LGE(CFUNCT,'S')) THEN
        IFUNCT=2
      ELSE
        CALL STETER(6,'FUNCT must be C or S')
      ENDIF

C
C   Map input data
C
      CALL TBTOPN (INAME1,F_I_MODE,TID1,ISTAT)
      CALL TBIGET (TID1,NCOL,NOBS1,ISOR,ICOL,NROW,ISTAT)
      CALL TBLSER (TID1,'TIME' ,ITIME,ISTAT)
      IF (ITIME.LT.0) THEN
        CALL STETER(5,'Column :TIME not found in 1st table')
      ENDIF
      CALL TBLSER (TID1,'VALUE',IDAT ,ISTAT)
      IF (IDAT.LT.0) THEN
        CALL STETER(6,'Column :VALUE not found in 1st table')
      ENDIF
      CALL TBLSER (TID1,'VAR'  ,IVAR, ISTAT)
      IF (IVAR.LT.0) THEN
        CALL STETER(7,'Column :VAR not found in 1st table')
      ENDIF
      CALL TBFGET (TID1,ITIME,FORM,LFIELD,TTYP,ISTAT)
      CALL TBFGET (TID1,IDAT, FORM,LFIELD,DTYP,ISTAT)
      CALL TBFGET (TID1,IVAR, FORM,LFIELD,VTYP,ISTAT)
      CALL TBDGET (TID1,ISTORE,ISTAT)
      IF (ISTORE.NE.F_TRANS) THEN
        TEXT='Input table '//INAME1//' stored not transposed'
        CALL STETER(1,TEXT)
      ENDIF
      IF (TTYP.NE.D_R8_FORMAT.OR.DTYP.NE.D_R8_FORMAT.OR.
     $     VTYP.NE.D_R8_FORMAT) THEN
        CALL STETER(1,
     $    'Column(s) in 1st table must be in DOUBLE PRECISION')
      ENDIF
      CALL TBCMAP (TID1,ITIME,POBS1T,ISTAT)
      CALL TBCMAP (TID1,IDAT, POBS1D,ISTAT)
      CALL TBCMAP (TID1,IVAR, POBS1V,ISTAT)
C
      IF (LLE(INAME2,INAME1).AND.LGE(INAME2,INAME1)) THEN
        POBS2T=POBS1T
        POBS2D=POBS1D
        POBS2V=POBS1V
      ELSE
        CALL TBTOPN (INAME2,F_I_MODE,TID2,ISTAT)
        CALL TBIGET (TID2,NCOL,NOBS2,ISOR,ICOL,NROW,ISTAT)
        CALL TBLSER (TID2,'TIME' ,ITIME,ISTAT)
        IF (ITIME.LT.0) THEN
          CALL STETER(8,'Column :TIME not found in 2nd table')
        ENDIF
        CALL TBLSER (TID2,'VALUE',IDAT ,ISTAT)
        IF (IDAT.LT.0) THEN
          CALL STETER(9,'Column :VALUE not found in 2nd table')
        ENDIF
        CALL TBLSER (TID2,'VAR'  ,IVAR, ISTAT)
        IF (IVAR.LT.0) THEN
          CALL STETER(10,'Column :VAR not found in 2nd table')
        ENDIF
        CALL TBFGET (TID2,ITIME,FORM,LFIELD,TTYP,ISTAT)
        CALL TBFGET (TID2,IDAT, FORM,LFIELD,DTYP,ISTAT)
        CALL TBFGET (TID2,IVAR, FORM,LFIELD,VTYP,ISTAT)
        CALL TBDGET (TID2,ISTORE,ISTAT)
        IF (ISTORE.NE.F_TRANS) THEN
          TEXT='Input table '//INAME2//' stored not transposed'
          CALL STETER(3,TEXT)
        ENDIF
        IF (TTYP.NE.D_R8_FORMAT.OR.DTYP.NE.D_R8_FORMAT.OR.
     $       VTYP.NE.D_R8_FORMAT) THEN
          CALL STETER(4,
     $      'Column(s) in 2nd table must be in DOUBLE PRECISION')
        ENDIF
        CALL TBCMAP (TID2,ITIME,POBS2T,ISTAT)
        CALL TBCMAP (TID2,IDAT, POBS2D,ISTAT)
        CALL TBCMAP (TID2,IVAR, POBS2V,ISTAT)
      ENDIF
C
C   Create and map output table
C
      NCOL=3
      CALL TBTINI (ONAME,F_TRANS,F_IO_MODE,NCOL,NLAG,OTID,ISTAT)
      CALL TBCINI (OTID,D_R8_FORMAT,1,'E15.6',' ','LAG',  ICOVN,ISTAT)
      CALL TBCINI (OTID,D_R8_FORMAT,1,'E15.6',' ','COVAR',ICOVC,ISTAT)
      CALL TBCINI (OTID,D_R8_FORMAT,1,'E15.6',' ','ERROR',ICOVV,ISTAT)
      CALL TBCMAP (OTID,ICOVN,PCOVN,ISTAT)
      CALL TBCMAP (OTID,ICOVC,PCOVC,ISTAT)
      CALL TBCMAP (OTID,ICOVV,PCOVV,ISTAT)
C
C   Compute covariance function
C
      CALL COVAR(
     $  MADRID(POBS1T),MADRID(POBS1D),MADRID(POBS1V),
     $  MADRID(POBS2T),MADRID(POBS2D),MADRID(POBS2V),
     $  MADRID(PCOVC), MADRID(PCOVV), MADRID(PCOVN),
     $  NOBS1,NOBS2,NLAG,ISCALE,IFUNCT,TLAG1,DTLAG,AV1,VAR1,AV2,VAR2)
      CALL TBIPUT (OTID,NCOL,NLAG,ISTAT)
C
C   Wind-up
C
      CALL DSCUPT(OTID,OTID,' ',ISTAT)
      CALL STSEPI
C
      END
C
C
C
C
C   COVAR
C
C   Compute discrete covariance function for unequaly sampled signal
C   Reference: Edelson R.A. and Krolik, J.H., 1988, ApJ 333, 646
C
C   Input:
C   OBS1T(NOBS1) - 1st set of observations: time
C   OBS1D(NOBS1) - 1st set of observations:     value
C   OBS1V(NOBS1) - 1st set of observations:          variance
C   OBS2T(NOBS2) - 2nd set of observations: time
C   OBS2D(NOBS2) - 2nd set of observations:     value
C   OBS2V(NOBS2) - 2nd set of observations:          variance
C   COVN (NLAG)  - time lag,...
C   COVC (NLAG)  -       ... covariance, ...
C   COVV (NLAG)  -                   ... and its error
C   ISCALE          - mode of output scale (linear of logarythmic)
C   IFUNCT          - type of output function (correlation or structure)
C   NOBS1         - number of observations in 1st set
C   NOBS2         - number of observations in 2nd set
C   TLAG1         - first time lag
C   DTLAG         - lag increment
C   NLAG          - number of lags
C
C   Output:
C   COV(NLAG)     - structure function
C
      SUBROUTINE COVAR(
     $  OBS1T,OBS1D,OBS1V,OBS2T,OBS2D,OBS2V,
     $  COVC,COVV,COVN,NOBS1,NOBS2,NLAG,ISCALE,IFUNCT,TLAG1,DTLAG,
     $  AV1,VAR1,AV2,VAR2)
C
      INCLUDE 'MID_REL_INCL:TSA_DEF.INC'
C
      INTEGER NOBS1,NOBS2,NLAG,ISCALE,IFUNCT
      DOUBLE PRECISION  OBS1T(NOBS1),OBS1D(NOBS1),OBS1V(NOBS1)
      DOUBLE PRECISION  OBS2T(NOBS2),OBS2D(NOBS2),OBS2V(NOBS2)
      DOUBLE PRECISION  COVC(NLAG),  COVV(NLAG),  COVN(NLAG)
      DOUBLE PRECISION  TLAG1,DTLAG
      DOUBLE PRECISION  AV1,VAR1,AV2,VAR2,VAR12
C
      INTEGER IOBS1,IOBS2,ILAG
      DOUBLE PRECISION  CT1,CT2,COUNT,OLDAVER,AVER
      DOUBLE PRECISION  V,VI,EI
      DOUBLE PRECISION  T,TLAG0,TI
C
      INCLUDE 'MID_REL_INCL:TSA_CONST.INC'
C
C   Compute weighted averages and variances
C
      AV1=ZERO
      AV2=ZERO
      VAR1=ZERO
      VAR2=ZERO
      CT1=ZERO
      CT2=ZERO
      DO 3 IOBS1=1,NOBS1
        AV1=AV1+OBS1D(IOBS1)/OBS1V(IOBS1)
 3      CT1=CT1+ONE/OBS1V(IOBS1)
      AV1=AV1/CT1
      DO 31 IOBS1=1,NOBS1
        V=OBS1D(IOBS1)-AV1
 31     VAR1=VAR1+V*V/OBS1V(IOBS1)
      VAR1=VAR1/(NOBS1-ONE)*NOBS1/CT1
      DO 4 IOBS2=1,NOBS2
        AV2=AV2+OBS2D(IOBS2)/OBS2V(IOBS2)
 4      CT2=CT2+ONE/OBS2V(IOBS2)
      AV2=AV2/CT2
      DO 41 IOBS2=1,NOBS2
        V=OBS2D(IOBS2)-AV2
 41     VAR2=VAR2+V*V/OBS2V(IOBS2)
      VAR2=VAR2/(NOBS2-ONE)*NOBS2/CT2
      VAR12=SQRT(VAR1*VAR2)
C
C   Compute binned discrete covariance function
C
      DO 5 ILAG=1,NLAG
        COVC(ILAG)=ZERO
        COVV(ILAG)=ZERO
 5      COVN(ILAG)=ZERO
      TLAG0=TLAG1
      DO 1 IOBS1=1,NOBS1
        TI=OBS1T(IOBS1)
        VI=OBS1D(IOBS1)-AV1
        EI=OBS1V(IOBS1)
        DO 1 IOBS2=1,NOBS2
          T=TI-OBS2T(IOBS2)
          IF (T.GT.ZERO) THEN
            IF (ISCALE.EQ.2) THEN
              T=LOG10(T)
            ENDIF
            ILAG=(T-TLAG0)/DTLAG
            IF (ILAG.GT.0.AND.ILAG.LE.NLAG) THEN
              V=AV2-OBS2D(IOBS2)+VI
              V=V*V-EI-OBS2V(IOBS2)
C
C   Use rounding proof algorithm for accumulation of means and variances
C
              COUNT=COVN(ILAG)+ONE
              OLDAVER=COVC(ILAG)
              AVER=OLDAVER+(V-OLDAVER)/COUNT
              COVC(ILAG)=AVER
              COVV(ILAG)=COVV(ILAG)+(V-OLDAVER)*(V-AVER)
              COVN(ILAG)=COUNT
            ENDIF
          ENDIF
 1    CONTINUE
C
C   Get errors and time lag grid
C
      DO 6 ILAG=1,NLAG
        IF (COVN(ILAG).GT.1.5) THEN
          COVV(ILAG)=SQRT(COVV(ILAG)/(COVN(ILAG)-ONE))
          IF (IFUNCT.EQ.1) THEN
            COVV(ILAG)=VAR12-COVV(ILAG)
          ENDIF
        ENDIF
        IF (IFUNCT.EQ.1) THEN
          COVC(ILAG)=VAR12-COVC(ILAG)
        ENDIF
        COVN(ILAG)=TLAG1+DTLAG*(ILAG-1)
 6    CONTINUE
      END