SUBROUTINE VPOLIN (
*
*  inputs
*
     :                  XA, YA, N, C, D, X,
*
*  outputs
*
     :                  Y, DY)
*
*  Module number:
*
*  Module name:
*
*  Keyphrase:
*  ----------
*  Lagrange's formula of polynomial interpolation
*
*  Description:
*  ------------
*  Neville's algorithm of Lagrange's formula of polynomial interpolation.
*  Adapted from the subroutine POLINT in p. 82 of Numerical Recipes by
*  Press et. al. (1986),
*
*  Given arrays XA and YA, each of length N, and given a value X, this routine
*  returns a value Y, and an error estimate DY.
*
*  FORTRAN name: VPOLIN.FOR
*
*  Keywords of accessed files and tables:
*  --------------------------------------
*  Name                         I/O     Description / Comments
*
*  Subroutines Called:
*  -------------------
*  CDBS:
*       None
*  SDAS:
*       None
*  Others:
*       None
*
*  History:
*  --------
*  Version      Date        Author          Description
*     1       02-28-88     J.-C. HSU        coding
*
*-------------------------------------------------------------------------------
*
*== input:
*                                       --size of arrays XA and YA
        INTEGER                 N
*                                       --input arrays
        DOUBLE PRECISION        XA(1), YA(1),
*                                       --working array
     :                          C(1), D(1),
     :                          X
*
*== output:
*                                       --interpolated value and its error
        DOUBLE PRECISION        Y, DY
*
*== local:
*
        INTEGER                 M, NS, I
        DOUBLE PRECISION        HO, HP, W, DEN
*
*------------------------------------------------------------------------------
*
*  initialize the tableau C and D
*
        DO 10 I = 1, N
            C(I) = YA(I)
            D(I) = YA(I)
   10   CONTINUE
*
        NS = N / 2 + 1
*
*  initial approximate value of Y
*
        Y = YA(NS)
*
        NS = NS-1
*
*  for each column of the tableau, loop over the current C's and D's and update
*  them
*
        DO 30 M = 1, N-1
            DO 20 I = 1, N-M
                HO = XA(I) - X
                HP = XA(I+M) - X
                W = C(I+1) - D(I)
                DEN = HO - HP
*
                DEN = W / DEN
                D(I) = HP * DEN
                C(I) = HO * DEN
   20       CONTINUE
*
            IF (2*NS .LT. N-M) THEN
                DY = C(NS+1)
            ELSE
                DY = D(NS)
                NS = NS - 1
            ENDIF
*
            Y = Y + DY
   30   CONTINUE
*
        RETURN
        END