Initial commit

1 files changed, 41 insertions, 0 deletions

diff --git a/math/interp/iipol_terp.x b/math/interp/iipol_terp.x
new file mode 100644
index 00000000..030d1964
--- /dev/null
+++ b/math/interp/iipol_terp.x
@@ -0,0 +1,41 @@
+# 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