C @(#)tsacov.for	19.1 (ESO-DMD) 02/25/03 13:33:25
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 Correspence 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+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
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.RETURNS   None
C.ENVIRON   TSA context
C.VERSION   0.0               June 1992
C 
C 021031	last modif
C 
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      ICOVN, ICOVC, ICOVV
C 
      INTEGER*8      POBS1T,POBS1D,POBS1V
      INTEGER*8      POBS2T,POBS2D,POBS2V
      INTEGER*8      PCOVN, PCOVC, PCOVV
C 
      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

