C @(#)interpol.for	17.1.1.1 (ES0-DMD) 01/25/02 17:20:21
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  @(#)interpol.for	5.1 (ESO-IPG) 4/5/93 15:58:36
C
C
C
C   INTERPOLATE
C
C   Interpolate an unevensampled time series using its correlation matrix
C           Reference:  ApJ 385, 404
C
C   Input:
C   OBS1T(NOBS1)               - observation  matrix (time)
C   OBS1D(NOBS1)               - observation  matrix (value)
C   OBS1V(NOBS1)               - observation  matrix (variance)
C   OBS2T(NOBS2)               - interpolated matrix (time)     input
C   ACFUNC(TLAG,PARM)          - autocorrelation function
C   PARM(12)                   - its parameters
C   COBS(NOBS1+1,NOBS1+1)        - work file
C   NU(NOBS1)                 - work file
C   CT(NOBS1)                 - work file
C   Output:
C   OBS2D(NOBS2)               - interpolated matrix (value)    output
C   OBS2V(NOBS2)               - interpolated matrix (variance) output
C
      SUBROUTINE INTERPOL(OBS1T,OBS1D,OBS1V,OBS2T,OBS2D,OBS2V,
     $    ACFUNC,PARM,NOBS1,NOBS2,MODE,AVERT,VART,COBS,NU,CT)
      INCLUDE 'MID_REL_INCL:TSA_DEF.INC'
      INTEGER            NOBS1,NOBS2,MODE
      DOUBLE PRECISION   OBS1T(NOBS1),OBS1D(NOBS1),OBS1V(NOBS1)
      DOUBLE PRECISION   OBS2T(NOBS2),OBS2D(NOBS2),OBS2V(NOBS2)
      DOUBLE PRECISION   ACFUNC,CRACOW                    ! EXTERNAL FUNCTIONS
      DOUBLE PRECISION   PARM(12)
      DOUBLE PRECISION   AVERT,VART
      DOUBLE PRECISION   COBS(NOBS1+1,NOBS1+1)
      DOUBLE PRECISION   NU(NOBS1)
      DOUBLE PRECISION   CT(NOBS1)
      EXTERNAL ACFUNC
      EXTERNAL CRACOW
C
      INTEGER            IOBS,JOBS,KOBS
      DOUBLE PRECISION   AVER1,VAR1,ACFCOR
      DOUBLE PRECISION   SAV,WORK,VAR,COUNT,TINY
      INCLUDE 'MID_REL_INCL:TSA_CONST.INC'
      DATA TINY/1D-25/
C
C   Determine global characteristics of input data
C
      AVERT=ZERO
      VART= ZERO
      COUNT=ZERO
      DO 1 IOBS=1,NOBS1
        IF (OBS1T(IOBS).GE.OBS2T(1).AND.
     $      OBS1T(IOBS).LE.OBS2T(NOBS2)) THEN
          IF (MODE.EQ.1) THEN
            AVERT=AVERT+OBS1D(IOBS)
            VART=VART+OBS1D(IOBS)*OBS1D(IOBS)
            COUNT=COUNT+ONE
          ELSEIF (MODE.EQ.2) THEN
            AVERT=AVERT+OBS1D(IOBS)*OBS1V(IOBS)
            VART=VART+OBS1D(IOBS)*OBS1D(IOBS)*OBS1V(IOBS)
            COUNT=COUNT+OBS1V(IOBS)
          ENDIF
        ENDIF
 1    CONTINUE
      AVER1=AVERT/COUNT
      VAR1=VART/COUNT-AVER1*AVER1
      ACFCOR=VAR1-ACFUNC(TINY,PARM)
C
C   Build correlation matrix of observations C
C
      DO 9 IOBS=1,NOBS1
        DO 8 JOBS=IOBS,NOBS1
          SAV=ABS(OBS1T(IOBS)-OBS1T(JOBS))+TINY
C          COBS(JOBS,IOBS)=ACFUNC(SAV,PARM)+ACFCOR
          COBS(JOBS,IOBS)=ACFUNC(SAV,PARM)
          IF (IOBS.EQ.JOBS) THEN
            COBS(IOBS,IOBS)=COBS(IOBS,IOBS)+OBS1V(IOBS)
          ENDIF
 8      CONTINUE
 9    CONTINUE
C
C   Get its triangle inverse square root C**(-1/2)
C
      IF ( CRACOW(COBS,NOBS1+1,NOBS1+1,0).EQ.ZERO) THEN
        CALL STETER(110,
     $    'INTERPOL: Correlation matrix of observations singular')
      ENDIF
C
C   Compute the interpolation vector NU = C**(-1/2)*C**(-1/2)*OBS1D
C
      DO 11 IOBS=1,NOBS1
        SAV=ZERO
        DO 12 JOBS=IOBS,NOBS1
          SAV=SAV+COBS(IOBS,JOBS)*OBS1D(JOBS)
 12     CONTINUE
        NU(IOBS)=SAV
 11   CONTINUE
      DO 14 IOBS=1,NOBS1
        SAV=ZERO
        DO 15 JOBS=IOBS,NOBS1
          SAV=SAV+COBS(IOBS,JOBS)*NU(JOBS)
 15     CONTINUE
        NU(IOBS)=SAV
 14   CONTINUE
C
C   Do actual interpolation OBS2D = NU**T*CT  OBS2V=VAR2-(C**(-1/2)*CT)**2
C
      DO 21 IOBS=1,NOBS2
        SAV=ZERO
        WORK=OBS2T(IOBS)
        IF (WORK.GE.OBS1T(1).AND.WORK.LE.OBS1T(NOBS1)) THEN
          DO 22 JOBS=1,NOBS1
            CT(JOBS)=ACFUNC(OBS1T(JOBS)-WORK,PARM)
            SAV=SAV+NU(JOBS)*CT(JOBS)
 22       CONTINUE
C          OBS2D(IOBS)=SAV
          OBS2D(IOBS)=-SAV
C
C   Determine global characteristics of all data
C
          IF (MODE.EQ.1) THEN
            AVERT=AVERT+OBS2D(IOBS)
            VART= VART+ OBS2D(IOBS)*OBS2D(IOBS)
            COUNT=COUNT+ONE
          ELSEIF (MODE.EQ.2) THEN
            VAR=ZERO
            DO 31 KOBS=1,NOBS1
              SAV=ZERO
              DO 32 JOBS=1,KOBS
                SAV=SAV+COBS(JOBS,KOBS)*CT(JOBS)
 32           CONTINUE
              VAR=VAR+SAV*SAV
 31         CONTINUE
C            OBS2V(IOBS)=VAR1-VAR
            OBS2V(IOBS)=ACFUNC(TINY,PARM)-VAR
            AVERT=AVERT+OBS2D(IOBS)*OBS2V(IOBS)
            VART= VART+ OBS2D(IOBS)*OBS2D(IOBS)*OBS2V(IOBS)
            COUNT=COUNT+OBS2V(IOBS)
          ENDIF
        ENDIF
 21   CONTINUE
      AVERT=AVERT/COUNT
      VART=VART/COUNT-AVERT*AVERT
      END