Initial commit

1 files changed, 39 insertions, 0 deletions

diff --git a/math/iminterp/ii_polterp.x b/math/iminterp/ii_polterp.x
new file mode 100644
index 00000000..24216751
--- /dev/null
+++ b/math/iminterp/ii_polterp.x
@@ -0,0 +1,39 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include "im1interpdef.h"
+
+# II_POLTERP -- polynomial interpolator with x and y arrays given.
+# This algorithm is based on the Newton form as described in
+# C. 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 MAX_NDERIVS.
+
+real procedure ii_polterp (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[MAX_NDERIVS]
+
+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