# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.

.help
		Procedure iipol_terp

A polynomial interpolator with x and y arrays given.
This algorithm is based on the Newton form as described in
de Boor's book, A Practical Guide to Splines, 1978.
There is no error checking - this is meant to be used only by calls
from more complete routines that take care of such things.

Maximum number of terms is 6.
.endhelp

real procedure iipol_terp(x,y,n,x0)

real x[ARB],y[ARB]	# x and y array
real x0			# desired x
int n			# number of points in x and y = number of
			# terms in polynomial = order + 1

int k,i
real d[6]

begin

	# Fill in entries for divided difference table.
	do i = 1,n
	    d[i] = y[i]

	# Generate divided differences
	do k = 1,n-1
	    do i = 1,n-k
		d[i] = (d[i+1] - d[i])/(x[i+k] - x[i])

	# Shift divided difference table to center on x0
	do i = 2,n
	    d[i] = d[i] + d[i-1] * (x0 - x[i])

	return(d[n])
end