SUBROUTINE VXYZIN (
*
*  inputs
*
     :                  MJD, XYZ, NPTS, NLAG, EPOCH,
*
*  outputs
*
     :                  PXYZ0, STATUS)
*
*  Module number:
*
*  Module name:
*
*  Keyphrase:
*  ----------
*  state vector interpolation for a specified epoch
*
*  Description:
*  ------------
*  Using Lagrange's formula of polynomial interpolation, obtain the state
*  vector of any epoch from the input state vector table.
*  The input state vectors must be tabulated in the order of increasing epoch.
*
*  FORTRAN name: VXYZIN.FOR
*
*  Keywords of accessed files and tables:
*  --------------------------------------
*  Name                         I/O     Description / Comments
*
*  Subroutines Called:
*  -------------------
*  CDBS:
*       VPOLIN
*  SDAS:
*       UMSPUT
*  Others:
*       None
*
*  History:
*  --------
*  Version      Date        Author          Description
*     1       02-28-88     J.-C. HSU        coding
*
*-------------------------------------------------------------------------------
*
*== input:
*                                       --number of entries of the input state
*                                       --vector table
        INTEGER                 NPTS,
*                                       --number of points used in Lagrange's
*                                       --polynomial interpolation fornula
     :                          NLAG
*                                       --epoch array from the state vector
*                                       --table (in MJD)
        DOUBLE PRECISION        MJD(1),
*                                       --input state vectors
     :                          XYZ(3, 1),
*                                       --desired epoch
     :                          EPOCH
*
*== output:
*                                       --state vector of the desired epoch
        DOUBLE PRECISION        PXYZ0(3)
*                                       --error status
        INTEGER                 STATUS
*
*== local:
*
        INTEGER                 I, J, ISTART, ISTOP, DIM, STATOK
        PARAMETER               (DIM = 20)
        DOUBLE PRECISION        TA(DIM), XA(DIM), YA(DIM), ZA(DIM),
*                                       --working array
     :                          C(DIM), D(DIM),
     :                          DX, DY, DZ
        CHARACTER*130           CONTXT, MESS
*
*------------------------------------------------------------------------------
*=========================begin hsp.inc=========================================
*                                       --status return code
        INTEGER         OK, ERRNUM(20)
        INTEGER         DEST, PRIO
        DATA OK /0/
        DATA ERRNUM /701, 702, 703, 704, 705, 706, 707, 708, 709, 710,
     :               711, 712, 713, 714, 715, 716, 717, 718, 719, 720/
*                                       --message destination and priority
        DATA DEST, PRIO /1, 0/
*=========================end hsp.inc===========================================
*
*  check the number points is no larger than the size of working areas
*
        IF (NLAG .GT. DIM) THEN
            STATUS = ERRNUM(1)
            CONTXT = 'number of interpolation points is too large'
            GO TO 999
        END IF
*
*  find the index of the desired epoch
*
        DO 10 I = 1, NPTS-1
            IF (EPOCH .LT. MJD(I+1) .AND. EPOCH .GE. MJD(I)) GO TO 20
   10   CONTINUE
*
*  if no index is found, something is wrong, exit
*
        STATUS = ERRNUM(1)
        CONTXT = 'cannot find proper epoch range in the input state '
     :                  // 'vector table'
        GO TO 999
*
*  try to center the desired epoch, but if it is at either ends of the table
*  have to shift the usable array
*
   20   ISTART = I - NLAG / 2
        IF (ISTART .LT. 1) ISTART = 1
*
        ISTOP = ISTART + NLAG - 1
        IF (ISTOP .GT. NPTS) ISTART = NPTS - NLAG + 1
*
        DO 30 J = 1, NLAG
            TA(J) = MJD(ISTART+J-1)
            XA(J) = XYZ(1, ISTART+J-1)
            YA(J) = XYZ(2, ISTART+J-1)
            ZA(J) = XYZ(3, ISTART+J-1)
   30   CONTINUE
*
*  perform the interpolation
*
        CALL VPOLIN (TA, XA, NLAG, C, D, EPOCH,    PXYZ0(1), DX)
        CALL VPOLIN (TA, YA, NLAG, C, D, EPOCH,    PXYZ0(2), DY)
        CALL VPOLIN (TA, ZA, NLAG, C, D, EPOCH,    PXYZ0(3), DZ)
*
        GO TO 1000
*
  999   MESS = 'VXYZIN: ' // CONTXT
        CALL UMSPUT (MESS, DEST, PRIO,    STATOK)
*
 1000   RETURN
        END