include "register.h"

#       Given arrays XA and YA, each of length N, and given a value X, this
#       routine returns a value Y, and an estimated error DY.  If P(x) is
#       the polynomial of degree N-1 such that P(XAi) = YAi, i = 1,...,N
#       then the returned value Y = P(X).  From Numerical Recipes, p82,
#       Press et al.  This version will interpolate for polynomials up to
#       order 20.

procedure polint (xa, ya, n, x, y, dy, status)

int n
real xa[n], ya[n]
real x, y, dy
real c[21], d[21]
int i, ns, m
real dif, dift, ho, hp, w, den
int status
char context[SZ_LINE]

string zderr "Adjacent Xs are equal: zero divide error"

begin 

        status = OK

        ns = 1
        dif = abs(x - xa[1])

        do i = 1, n {

#          Here we find the index of the closest table entry

           dift = abs(x - xa[i])
           if (dift <  dif) {
              ns = i
              dif = dift
           }

#          and initialize the tableau of C's and D's

           c[i] = ya[i]
           d[i] = ya[i]

        }

#       This is the initial approximation to Y

        y = ya[ns]
        ns = ns - 1

        do m = 1, n-1 {
           do i = 1, n-m {

#             For each column of the tableau, we loop over
#             the current C's and D's and update them

              ho = xa[i] - x
              hp = xa[i+m] - x
              w = c[i+1] -d[i]
              den = ho - hp
              if (den == 0.0) {
              call strcpy(zderr, context, SZ_LINE)
	      call eprintf(context)
              status = FAIL
              }
              den = w/den

#             Update the C's and D's

              d[i] = hp*den
              c[i] = ho*den

           }

#          After each column in the tableau is completed, we decide
#          which correction, C or D we want to add to our accumulating
#          value of Y, i.e. which path to take through the tableau--
#          forking up or down.  We do this in such a way as to take the
#          most "straight line" route through the tableau to its apex,
#          updating ns accordingly to keep track of where we are.
#          This route keeps the partial approximations centered
#          (insofar as possible) on the target X.  The last DY added
#          is thus the error indication.

           if (2*ns < n-m) {
              dy = c[ns+1]
           } else {
              dy = d[ns]
              ns = ns - 1
           }
           y = y + dy

        }


end