C @(#)tsadel.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  @(#)tsadel.for	5.1 (ESO-IPG) 4/5/93 15:58:47
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT (c)  1992     European Southern Observatory and Warsaw Observatory
C.IDENT     tsadel.for
C.AUTHOR    Alex Schwarzenberg-Czerny, ESO and Warsaw Observatory
C.KEYWORD   MIDAS, time analysis, DELAY/TIME

C.LANGUAGE  FORTRAN 77
C.PURPOSE   Compute chi2-time lag function
C           Reference:  ApJ 385, 404
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 CHI2 TABLE
      REAL*8              TLAG1          !  START TIME LAG
      REAL*8              DTLAG          !  TIME LAG INCREMENT
      INTEGER             NLAG           !  NUMBER OF TIME LAGS
      CHARACTER*3         CFUNC          !  TYPE OF ACF
      CHARACTER*3         CMODE          !  MODE
      DOUBLE PRECISION    PARM(12)       !  PARAMETERS OF ACF
      INTEGER             MODE
      DOUBLE PRECISION ACFLIN,ACFPOL,ACFPOW,ACFEXP,
     $                 ACFLOG,ACFIPO
      EXTERNAL         ACFLIN,ACFPOL,ACFPOW,ACFEXP,
     $                 ACFLOG,ACFIPO
      DOUBLE PRECISION TSADELUR0,TSADELUR1,TSADELUR2,TSADELUR3,
     $                 TSADELUR4,TSADELUR5,TSADELUR6,TSADELUR7,
     $                 TSADELUR8,TSADELUR9
      EXTERNAL         TSADELUR0,TSADELUR1,TSADELUR2,TSADELUR3,
     $                 TSADELUR4,TSADELUR5,TSADELUR6,TSADELUR7,
     $                 TSADELUR8,TSADELUR9
C
      INTEGER      NOBS1, NOBS2, NOBS,  NOBSP1
      INTEGER      NPAR,  ISIZE, ASIZE
      INTEGER      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      IOBS
      INTEGER      POBS
      INTEGER      ICOBS, ICOND, IPAR
      INTEGER      PCOBS, PCOND, PPAR
      INTEGER      ICHI2L,ICHI2D
      INTEGER      PCHI2L,PCHI2D
      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 ('tsadel')
      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 ('CFUNC',  1,1, 3,IACTS,CFUNC, KUN,KNUL,ISTAT)
      CALL STKRDC ('CMODE',  1,1, 3,IACTS,CMODE, KUN,KNUL,ISTAT)
      CALL STKRDD ('INPUTD',   1,12,IACTS,PARM,  KUN,KNUL,ISTAT)
      CALL STKRDI ('ORDERTSA',
     $    1, 1,IACTS,MODE,  KUN,KNUL,ISTAT)
      IF     (LLE(CMODE,'NOR').AND.LGE(CMODE,'NOR')) THEN
        MODE=1
        NPAR=2
      ELSEIF (LLE(CMODE,'ORI').AND.LGE(CMODE,'ORI')) THEN
        MODE=2
        NPAR=1
      ELSE
        CALL STETER(10,'Wrong MODE')
      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(31,'Column :TIME not found in 1st table')
      ENDIF
      CALL TBLSER (TID1,'VALUE',IDAT ,ISTAT)
      IF (IDAT.LT.0) THEN
        CALL STETER(32,'Column :VALUE not found in 1st table')
      ENDIF
      CALL TBLSER (TID1,'VAR'  ,IVAR, ISTAT)
      IF (IVAR.LT.0) THEN
        CALL STETER(33,'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(2,
     $     '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
      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(34,'Column :TIME not found in 2nd table')
      ENDIF
      CALL TBLSER (TID2,'VALUE',IDAT ,ISTAT)
      IF (IDAT.LT.0) THEN
        CALL STETER(35,'Column :VALUE not found in 2nd table')
      ENDIF
      CALL TBLSER (TID2,'VAR'  ,IVAR, ISTAT)
      IF (IVAR.LT.0) THEN
        CALL STETER(36,'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)
C
C   Create and map work tables  OBS,OBS1,OBS2,COBS,COND,CPAR
C
      NOBS=NOBS1+NOBS2
      NOBSP1=NOBS+1
      ISIZE=NOBS*3
      CALL STFCRE('ZZMIDAWORK',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,
     $   ISIZE,IOBS, ISTAT)
      CALL STFMAP(IOBS, F_X_MODE,1,ISIZE,ASIZE,POBS, ISTAT)
C
      ISIZE=NOBSP1*NOBSP1
      CALL STFCRE('ZZMIDDWORK',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,
     $   ISIZE,ICOBS,ISTAT)
      CALL STFMAP(ICOBS,F_X_MODE,1,ISIZE,ASIZE,PCOBS,ISTAT)
C
      ISIZE=(NPAR+1)*NOBS
      CALL STFCRE('ZZMIDEWORK',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,
     $   ISIZE,ICOND,ISTAT)
      CALL STFMAP(ICOND,F_X_MODE,1,ISIZE,ASIZE,PCOND,ISTAT)
C
      ISIZE=(NPAR+1)*(NPAR+1)
      CALL STFCRE('ZZMIDFWORK',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,
     $   ISIZE,IPAR, ISTAT)
      CALL STFMAP(IPAR, F_X_MODE,1,ISIZE,ASIZE,PPAR, ISTAT)
C
C   Create and map output table
C
      NCOL=2
      CALL TBTINI (ONAME,F_TRANS,F_IO_MODE,NCOL,NLAG,OTID,ISTAT)
      CALL TBCINI (OTID,D_R8_FORMAT,1,'E15.6',' ','LAG', ICHI2L,ISTAT)
      CALL TBCINI (OTID,D_R8_FORMAT,1,'E15.6',' ','CHI2',ICHI2D,ISTAT)
      CALL TBCMAP (OTID,ICHI2L,PCHI2L,ISTAT)
      CALL TBCMAP (OTID,ICHI2D,PCHI2D,ISTAT)
C
C   Compute covariance function
C
C
      IF (LLE(CFUNC,'LIN').AND.LGE(CFUNC,'LIN')) THEN
        CALL DELAY(
     $    MADRID(POBS1T), MADRID(POBS1D), MADRID(POBS1V),
     $    MADRID(POBS2T), MADRID(POBS2D), MADRID(POBS2V),
     $    MADRID(POBS),
     $    ACFLIN,PARM,    MADRID(PCOBS),  MADRID(PCOND),
     $    MADRID(PPAR),   MADRID(PCHI2L), MADRID(PCHI2D),
     $    NOBS1,NOBS2,NPAR+1,MODE,NLAG,TLAG1,DTLAG)
      ELSEIF (LLE(CFUNC,'POL').AND.LGE(CFUNC,'POL')) THEN
        CALL DELAY(
     $    MADRID(POBS1T), MADRID(POBS1D), MADRID(POBS1V),
     $    MADRID(POBS2T), MADRID(POBS2D), MADRID(POBS2V),
     $    MADRID(POBS),
     $    ACFPOL,PARM,    MADRID(PCOBS),  MADRID(PCOND),
     $    MADRID(PPAR),   MADRID(PCHI2L), MADRID(PCHI2D),
     $    NOBS1,NOBS2,NPAR+1,MODE,NLAG,TLAG1,DTLAG)
      ELSEIF (LLE(CFUNC,'POW').AND.LGE(CFUNC,'POW')) THEN
        CALL DELAY(
     $    MADRID(POBS1T), MADRID(POBS1D), MADRID(POBS1V),
     $    MADRID(POBS2T), MADRID(POBS2D), MADRID(POBS2V),
     $    MADRID(POBS),
     $    ACFPOW,PARM,    MADRID(PCOBS),  MADRID(PCOND),
     $    MADRID(PPAR),   MADRID(PCHI2L), MADRID(PCHI2D),
     $    NOBS1,NOBS2,NPAR+1,MODE,NLAG,TLAG1,DTLAG)
      ELSEIF (LLE(CFUNC,'EXP').AND.LGE(CFUNC,'EXP')) THEN
        CALL DELAY(
     $    MADRID(POBS1T), MADRID(POBS1D), MADRID(POBS1V),
     $    MADRID(POBS2T), MADRID(POBS2D), MADRID(POBS2V),
     $    MADRID(POBS),
     $    ACFEXP,PARM,    MADRID(PCOBS),  MADRID(PCOND),
     $    MADRID(PPAR),   MADRID(PCHI2L), MADRID(PCHI2D),
     $    NOBS1,NOBS2,NPAR+1,MODE,NLAG,TLAG1,DTLAG)
      ELSEIF (LLE(CFUNC,'LOG').AND.LGE(CFUNC,'LOG')) THEN
        CALL DELAY(
     $    MADRID(POBS1T), MADRID(POBS1D), MADRID(POBS1V),
     $    MADRID(POBS2T), MADRID(POBS2D), MADRID(POBS2V),
     $    MADRID(POBS),
     $    ACFLOG,PARM,    MADRID(PCOBS),  MADRID(PCOND),
     $    MADRID(PPAR),   MADRID(PCHI2L), MADRID(PCHI2D),
     $    NOBS1,NOBS2,NPAR+1,MODE,NLAG,TLAG1,DTLAG)
      ELSEIF (LLE(CFUNC,'IPO').AND.LGE(CFUNC,'IPO')) THEN
        CALL DELAY(
     $    MADRID(POBS1T), MADRID(POBS1D), MADRID(POBS1V),
     $    MADRID(POBS2T), MADRID(POBS2D), MADRID(POBS2V),
     $    MADRID(POBS),
     $    ACFIPO,PARM,    MADRID(PCOBS),  MADRID(PCOND),
     $    MADRID(PPAR),   MADRID(PCHI2L), MADRID(PCHI2D),
     $    NOBS1,NOBS2,NPAR+1,MODE,NLAG,TLAG1,DTLAG)
      ELSEIF (LLE(CFUNC,'UR0').AND.LGE(CFUNC,'UR0')) THEN
        CALL DELAY(
     $    MADRID(POBS1T), MADRID(POBS1D), MADRID(POBS1V),
     $    MADRID(POBS2T), MADRID(POBS2D), MADRID(POBS2V),
     $    MADRID(POBS),
     $    TSADELUR0,PARM, MADRID(PCOBS),  MADRID(PCOND),
     $    MADRID(PPAR),   MADRID(PCHI2L), MADRID(PCHI2D),
     $    NOBS1,NOBS2,NPAR+1,MODE,NLAG,TLAG1,DTLAG)
      ELSEIF (LLE(CFUNC,'UR1').AND.LGE(CFUNC,'UR1')) THEN
        CALL DELAY(
     $    MADRID(POBS1T), MADRID(POBS1D), MADRID(POBS1V),
     $    MADRID(POBS2T), MADRID(POBS2D), MADRID(POBS2V),
     $    MADRID(POBS),
     $    TSADELUR1,PARM, MADRID(PCOBS),  MADRID(PCOND),
     $    MADRID(PPAR),   MADRID(PCHI2L), MADRID(PCHI2D),
     $    NOBS1,NOBS2,NPAR+1,MODE,NLAG,TLAG1,DTLAG)
      ELSEIF (LLE(CFUNC,'UR2').AND.LGE(CFUNC,'UR2')) THEN
        CALL DELAY(
     $    MADRID(POBS1T), MADRID(POBS1D), MADRID(POBS1V),
     $    MADRID(POBS2T), MADRID(POBS2D), MADRID(POBS2V),
     $    MADRID(POBS),
     $    TSADELUR2,PARM, MADRID(PCOBS),  MADRID(PCOND),
     $    MADRID(PPAR),   MADRID(PCHI2L), MADRID(PCHI2D),
     $    NOBS1,NOBS2,NPAR+1,MODE,NLAG,TLAG1,DTLAG)
      ELSEIF (LLE(CFUNC,'UR3').AND.LGE(CFUNC,'UR3')) THEN
        CALL DELAY(
     $    MADRID(POBS1T), MADRID(POBS1D), MADRID(POBS1V),
     $    MADRID(POBS2T), MADRID(POBS2D), MADRID(POBS2V),
     $    MADRID(POBS),
     $    TSADELUR3,PARM, MADRID(PCOBS),  MADRID(PCOND),
     $    MADRID(PPAR),   MADRID(PCHI2L), MADRID(PCHI2D),
     $    NOBS1,NOBS2,NPAR+1,MODE,NLAG,TLAG1,DTLAG)
      ELSEIF (LLE(CFUNC,'UR4').AND.LGE(CFUNC,'UR4')) THEN
        CALL DELAY(
     $    MADRID(POBS1T), MADRID(POBS1D), MADRID(POBS1V),
     $    MADRID(POBS2T), MADRID(POBS2D), MADRID(POBS2V),
     $    MADRID(POBS),
     $    TSADELUR4,PARM, MADRID(PCOBS),  MADRID(PCOND),
     $    MADRID(PPAR),   MADRID(PCHI2L), MADRID(PCHI2D),
     $    NOBS1,NOBS2,NPAR+1,MODE,NLAG,TLAG1,DTLAG)
      ELSEIF (LLE(CFUNC,'UR5').AND.LGE(CFUNC,'UR5')) THEN
        CALL DELAY(
     $    MADRID(POBS1T), MADRID(POBS1D), MADRID(POBS1V),
     $    MADRID(POBS2T), MADRID(POBS2D), MADRID(POBS2V),
     $    MADRID(POBS),
     $    TSADELUR5,PARM, MADRID(PCOBS),  MADRID(PCOND),
     $    MADRID(PPAR),   MADRID(PCHI2L), MADRID(PCHI2D),
     $    NOBS1,NOBS2,NPAR+1,MODE,NLAG,TLAG1,DTLAG)
      ELSEIF (LLE(CFUNC,'UR6').AND.LGE(CFUNC,'UR6')) THEN
        CALL DELAY(
     $    MADRID(POBS1T), MADRID(POBS1D), MADRID(POBS1V),
     $    MADRID(POBS2T), MADRID(POBS2D), MADRID(POBS2V),
     $    MADRID(POBS),
     $    TSADELUR6,PARM, MADRID(PCOBS),  MADRID(PCOND),
     $    MADRID(PPAR),   MADRID(PCHI2L), MADRID(PCHI2D),
     $    NOBS1,NOBS2,NPAR+1,MODE,NLAG,TLAG1,DTLAG)
      ELSEIF (LLE(CFUNC,'UR7').AND.LGE(CFUNC,'UR7')) THEN
        CALL DELAY(
     $    MADRID(POBS1T), MADRID(POBS1D), MADRID(POBS1V),
     $    MADRID(POBS2T), MADRID(POBS2D), MADRID(POBS2V),
     $    MADRID(POBS),
     $    TSADELUR7,PARM, MADRID(PCOBS),  MADRID(PCOND),
     $    MADRID(PPAR),   MADRID(PCHI2L), MADRID(PCHI2D),
     $    NOBS1,NOBS2,NPAR+1,MODE,NLAG,TLAG1,DTLAG)
      ELSEIF (LLE(CFUNC,'UR8').AND.LGE(CFUNC,'UR8')) THEN
        CALL DELAY(
     $    MADRID(POBS1T), MADRID(POBS1D), MADRID(POBS1V),
     $    MADRID(POBS2T), MADRID(POBS2D), MADRID(POBS2V),
     $    MADRID(POBS),
     $    TSADELUR8,PARM, MADRID(PCOBS),  MADRID(PCOND),
     $    MADRID(PPAR),   MADRID(PCHI2L), MADRID(PCHI2D),
     $    NOBS1,NOBS2,NPAR+1,MODE,NLAG,TLAG1,DTLAG)
      ELSEIF (LLE(CFUNC,'UR9').AND.LGE(CFUNC,'UR9')) THEN
        CALL DELAY(
     $    MADRID(POBS1T), MADRID(POBS1D), MADRID(POBS1V),
     $    MADRID(POBS2T), MADRID(POBS2D), MADRID(POBS2V),
     $    MADRID(POBS),
     $    TSADELUR9,PARM, MADRID(PCOBS),  MADRID(PCOND),
     $    MADRID(PPAR),   MADRID(PCHI2L), MADRID(PCHI2D),
     $    NOBS1,NOBS2,NPAR+1,MODE,NLAG,TLAG1,DTLAG)
      ELSE
        CALL STETER(5,'Wrong name for ACF type')
      ENDIF
      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
C   DELAY
C
C   Determine joint chi2 of two observation samples as a function of time
C   and signal shift.
C
C   Input:
C   NOBS                      - number of observations in 1+2 sample
C   NOBSP1                    - number of observations in 1+2 sample +1
C   NOBS1                     - number of observations in 1 sample
C   NOBS2                     - number of observations in 2 sample
C   OBS1T,OBS1D,OBS1V(NOBS1)  - 1st observation set (time, value, variance)
C   OBS2T,OBS2D,OBS2V(NOBS2)  - 2nd observation set (time, value, variance)
C   ACFUNC(TLAG,PARM)         - autocorrelation function (ACF)
C   PARM(12)                  - parameters of ACF
C   NLAG                      - number of time lags to be considered
C   TLAG1                     - value of the 1st lag
C   DLAG                      - lag increment
C
C   Work tables:
C   OBS(NOBS,3)
C   COBS(NOBSP1,NOBSP1)
C   COND(NPAR1,NOBS)
C   CPAR(NPAR1,NPAR1)
C
C   Output:
C   CHI2L(NLAG)                - lags
C   CHI2D(NLAG)                - chi2 for given lags
C
      SUBROUTINE DELAY(
     $   OBS1T,OBS1D,OBS1V,OBS2T,OBS2D,OBS2V,OBS,
     $   ACFUNC,PARM,COBS,COND,CPAR,CHI2L,CHI2D,
     $   NOBS1,NOBS2,NPAR1,MODE,NLAG,TLAG1,DLAG)
C
      INCLUDE 'MID_REL_INCL:TSA_DEF.INC'
      EXTERNAL ACFUNC
      INTEGER          NOBS1, NOBS2
      INTEGER          NPAR1, MODE,  NLAG
      DOUBLE PRECISION OBS1T(NOBS1), OBS1D(NOBS1),  OBS1V(NOBS1)
      DOUBLE PRECISION OBS2T(NOBS2), OBS2D(NOBS2),  OBS2V(NOBS2)
      DOUBLE PRECISION OBS(NOBS1+NOBS2,3)
      DOUBLE PRECISION ACFUNC                    ! EXTERNAL FUNCTION
      DOUBLE PRECISION PARM(12)
      DOUBLE PRECISION COBS(NOBS1+NOBS2+1,NOBS1+NOBS2+1)
      DOUBLE PRECISION COND(NPAR1,NOBS1+NOBS2)
      DOUBLE PRECISION CPAR(NPAR1,NPAR1)
      DOUBLE PRECISION CHI2L(NLAG),  CHI2D(NLAG)
      REAL*8           TLAG1, DLAG
C
      INTEGER          NOBS,  NOBSP1
      INTEGER          NFIRST,IOBS,  JOBS,  ILAG,  I
      REAL*8           TLAG,  SAV
      DOUBLE PRECISION WGHT1, WGHT2, AVRG1, AVRG2, VAR1, VAR2
      CHARACTER*80     TEXT
      INCLUDE 'MID_REL_INCL:TSA_CONST.INC'
C
      NOBS=NOBS1+NOBS2
      NOBSP1=NOBS+1
      NFIRST=NOBS1
      DO 1 ILAG=1,NLAG
        TLAG=TLAG1+DLAG*(ILAG-1)
C
C   Copy time and apply time lag
C
        CALL INPCPY(OBS1T, OBS(1,1),       NOBS1)
        CALL INPCPY(OBS2T, OBS(NOBS1+1,1), NOBS2)
        DO 2 IOBS=NOBS1+1,NOBS
           OBS(IOBS,1)=OBS(IOBS,1)-TLAG
 2      CONTINUE
C
C   Get parameters for condition equations
C
        CALL INPCPY(OBS1D, OBS(1,2),       NOBS1)
        CALL INPCPY(OBS1V, OBS(1,3),       NOBS1)
        CALL INPCPY(OBS2D, OBS(NOBS1+1,2), NOBS2)
        CALL INPCPY(OBS2V, OBS(NOBS1+1,3), NOBS2)
        IF (MODE.EQ.1) THEN
          WGHT1=-HALF
          WGHT2=+HALF
        ELSEIF (MODE.EQ.2) THEN
          CALL INTERPOL(
     $       OBS(1,1),       OBS(1,2),       OBS(1,3),
     $       OBS(NOBS1+1,1), OBS(NOBS1+1,2), OBS(NOBS1+1,3),
     $       ACFUNC,PARM,NOBS1,NOBS2,2,AVRG1,VAR1,
     $       COBS,COND(1,1),COND(2,1))
          CALL INPCPY(OBS2D, OBS(NOBS1+1,2), NOBS2)
          CALL INPCPY(OBS2V, OBS(NOBS1+1,3), NOBS2)
          CALL INTERPOL(
     $       OBS(NOBS1+1,1), OBS(NOBS1+1,2), OBS(NOBS1+1,3),
     $       OBS(1,1),       OBS(1,2),       OBS(1,3),
     $       ACFUNC,PARM,NOBS2,NOBS1,2,AVRG2,VAR2,
     $       COBS,COND(1,1),COND(2,1))
          CALL INPCPY(OBS1D, OBS(1,2),       NOBS1)
          CALL INPCPY(OBS1V, OBS(1,3),       NOBS1)
          DO 3 IOBS=1,NOBS1
            OBS(IOBS,2)=OBS(IOBS,2)-AVRG1
 3        CONTINUE
          DO 4 IOBS=NOBS1+1,NOBS
            OBS(IOBS,2)=OBS(IOBS,2)-AVRG2
 4        CONTINUE
        ELSEIF (MODE.EQ.3) THEN
          AVRG1=ZERO
          AVRG2=ZERO
          VAR1=1E-32
          VAR2=1E-32
        ENDIF
C
C   Set condition equations
C
        DO 10 IOBS=1,NOBS
          IF (MODE.EQ.1) THEN
            IF (IOBS.GT.NFIRST) THEN
              COND(2,IOBS)=WGHT2
            ELSE
              COND(2,IOBS)=WGHT1
            ENDIF
          ELSEIF (MODE.EQ.2.OR.MODE.EQ.3) THEN
          ENDIF
          COND(1,IOBS)=ONE
          COND(NPAR1,IOBS)=OBS(IOBS,2)
 10     CONTINUE
C
C   Set covariance matrix
C
        DO 20 IOBS=1,NOBS
          DO 21 JOBS=IOBS,NOBS
            SAV=OBS(IOBS,1)-OBS(JOBS,1)
            COBS(JOBS,IOBS)=ACFUNC(SAV,PARM)
            IF (IOBS.EQ.JOBS) THEN
              COBS(IOBS,IOBS)=COBS(IOBS,IOBS)+OBS(IOBS,3)
            ENDIF
 21       CONTINUE
 20     CONTINUE
        CALL FITCOR(COND,COBS,CPAR,NOBS,NPAR1,CHI2D(ILAG))
        CHI2D(ILAG)=CHI2D(ILAG)*(NOBS-2)
        CHI2L(ILAG)=TLAG
C
C      Print progress report
C
        IF (MOD(ILAG,NLAG/10+1).EQ.0) THEN
          IF (ILAG.EQ.NLAG/10+1) THEN
            TEXT='    Lag     Std.Dev.'//
     $        '    Mean      Error   Difference  Error'
            CALL STTPUT(TEXT,ISTAT)
          ENDIF
          SAV=SQRT(CPAR(NPAR1,NPAR1))
          IF (MODE.EQ.1.OR.MODE.EQ.3) THEN
            WRITE(TEXT,'(1PE10.3,1PE10.3,3(1PE12.4,1PE8.1))')
     $        TLAG,SAV,(CPAR(I,NPAR1),SQRT(CPAR(I,I)),I=1,NPAR1-1)
          ELSEIF (MODE.EQ.2) THEN
            WRITE(TEXT,'(1PE10.3,1PE10.3,3(1PE12.4,1PE8.1))')
     $        TLAG,SAV,(CPAR(I,NPAR1)+HALF*(AVRG1+AVRG2),
     $        SQRT(CPAR(I,I)),I=1,NPAR1-1),
     $        AVRG2-AVRG1,SQRT((VAR1/(NOBS1-1)+VAR2/(NOBS2-1)))
          ENDIF
          CALL STTPUT(TEXT,ISTAT)
        ENDIF
 1    CONTINUE
      END
C
C
C
C
C    FITCOR
C
C   General purpose least squares routine for correlated observations.
C   The correlation function must be given analyticaly
C   and observation equations in the matrix form.
C
C   A. Schwarzenberg-Czerny, ESO & Warsaw Observatory, June 1992
C
C   Input:
C   COBS(NOBS,NOBS)     - observation correlation matrix (to be destroyed)
C   COND(NPAR1,NPAR1)   - normal equations
C   NOBS                - number of observations
C   NPAR1               - number of parameters+1 (1 for chi2 only)
C   MODE                - returned RHS are:
C                       - 0 - none
C                       - 1 - restored by least squares
C                       - 2 - detrended
C                       - 3 - restored and detrended
C
C   Output:
C   CHI2DF              - chi2 per degree of freedom
C   CPAR(NPAR1,NPAR1)   - fitted parameters and their covariance matrix
C                         (See CRACOW for details)
C
C   Method:
C   The AC matrix is inverted first. This requires its storage.
C   The inversion can be done with a call of CRACOW. Both mean that
C   both observation equations and the AC matrix have to be stored.
C   Algorithm:
C   Find triangle root of CO, invert it and multiply CND by the triangle
C   inverse. Square the result to get the normal equations. Solve them and
C   return values of the parameters and their correlation matrix.
C
      SUBROUTINE FITCOR(COND,COBS,CPAR,NOBS,NPAR1,CHI2DF)
C
      INCLUDE 'MID_REL_INCL:TSA_DEF.INC'
      DOUBLE PRECISION CRACOW
      INTEGER NOBS,NPAR1
      DOUBLE PRECISION COND(NPAR1,NOBS)
      DOUBLE PRECISION COBS(NOBS+1,NOBS+1)
      DOUBLE PRECISION CPAR(NPAR1,NPAR1)
      DOUBLE PRECISION CHI2DF
C
      INTEGER I,J,K,NOBSP1
      DOUBLE PRECISION SUM
      INCLUDE 'MID_REL_INCL:TSA_CONST.INC'
C
      NOBSP1=NOBS+1
      CHI2DF=ZERO
C   Find inverse triangle root and get normal equations
      IF (CRACOW(COBS,NOBSP1,NOBSP1,0).EQ.ZERO) THEN
        CALL STETER(20,
     $   ' Covariance matrix of observations is singular (empty?)')
        RETURN
      ENDIF
      DO 1 I=NOBS,1,-1
        DO 1 J=1,NPAR1
          SUM=ZERO
          DO 2 K=1,I
            SUM=SUM+COBS(K,I)*COND(J,K)
 2        CONTINUE
          COND(J,I)=SUM
 1    CONTINUE
      DO 10 I=1,NPAR1
        DO 10 J=I,NPAR1
          SUM=ZERO
          DO 11 K=1,NOBS
            SUM=SUM+COND(I,K)*COND(J,K)
 11       CONTINUE
          CPAR(J,I)=SUM
 10   CONTINUE
C
C   Fit parameters by least squares and return chi2
C
      IF (CRACOW(CPAR,NPAR1,NPAR1,NOBS).EQ.ZERO) THEN
        CALL STETER(21,
     $   ' Normal equations are singular, modify trend functions')
        RETURN
      ENDIF
      CHI2DF=CPAR(NPAR1,NPAR1)
      END
C
C
C