1 files changed, 137 insertions, 0 deletions
diff --git a/math/interp/iieval.x b/math/interp/iieval.x
new file mode 100644
index 00000000..87538525
--- /dev/null
+++ b/math/interp/iieval.x
@@ -0,0 +1,137 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+.help
+		procedures to evaluate interpolants
+
+     These procedures are seperated from the other programs in order
+to make it easier to optimise these in case it becomes necessary to
+obtain the fastest possible interpolants.  The interpolation package
+spends most of its time in these routines.
+.endhelp
+
+procedure iievne(x,y,n,a)   # nearest neighbor
+
+real	x[ARB]	# x values, must be within [1,n]
+real	y[ARB]	# interpolated values returned to user
+int	n	# number of x values
+real	a[ARB]	# data to be interpolated
+
+int	i
+
+begin
+	do i = 1, n
+	    y[i] = a[int(x[i] + 0.5)]
+	return
+end
+
+procedure iievli(x,y,n,a)   # linear
+
+real	x[ARB]	# x values, must be within [1,n]
+real	y[ARB]	# interpolated values returned to user
+int	n	# number of x values
+real	a[ARB]	# data to be interpolated
+
+int	i,nx
+
+begin
+	
+	do i = 1,n {
+	    nx = x[i]
+	    y[i] =  (x[i] - nx) * a[nx + 1] + (nx + 1 - x[i]) * a[nx]
+	}
+	return
+end
+
+procedure iievp3(x,y,n,a)   # interior third order polynomial
+
+real	x[ARB]	# x values, must be within [1,n]
+real	y[ARB]	# interpolated values returned to user
+int	n	# number of x values
+real	a[ARB]	# data to be interpolated from a[0] to a[n+2]
+
+int	i,nx,nxold
+real	s,t,cd20,cd21
+
+begin
+	nxold = -1
+	do i = 1,n {
+	    nx = x[i]
+	    s = x[i] - nx
+	    t = 1. - s
+	    if (nx != nxold) {
+
+		# second central differences:
+		cd20 = 1./6. * (a[nx+1] - 2. * a[nx] + a[nx-1])
+		cd21 = 1./6. * (a[nx+2] - 2. * a[nx+1] + a[nx])
+		nxold = nx
+	    }
+	    y[i] = s * (a[nx+1] + (s * s - 1.) * cd21) +
+		    t * (a[nx] + (t * t - 1.) * cd20)
+
+	}
+	return
+end
+
+procedure iievp5(x,y,n,a)   # interior fifth order polynomial
+
+real	x[ARB]	# x values, must be within [1,n]
+real	y[ARB]	# interpolated values returned to user
+int	n	# number of x values
+real	a[ARB]	# data to be interpolated - from a[-1] to a[n+3]
+
+int	i,nx,nxold
+real	s,t,cd20,cd21,cd40,cd41
+
+begin
+	nxold = -1
+	do i = 1,n {
+	    nx = x[i]
+	    s = x[i] - nx
+	    t = 1. - s
+	    if (nx != nxold) {
+		cd20 = 1./6. * (a[nx+1] - 2. * a[nx] + a[nx-1])
+		cd21 = 1./6. * (a[nx+2] - 2. * a[nx+1] + a[nx])
+
+		# fourth central differences
+		cd40 = 1./120. * (a[nx-2] - 4. * a[nx-1] +
+			6. * a[nx] - 4. * a[nx+1] + a[nx+2])
+		cd41 = 1./120. * (a[nx-1] - 4. * a[nx] +
+			6. * a[nx+1] - 4. * a[nx+2] + a[nx+3])
+		nxold = nx
+	    }
+	    y[i] = s * (a[nx+1] + (s * s - 1.) *
+			 (cd21 + (s * s - 4.) * cd41)) +
+		   t * (a[nx] + (t * t - 1.) *
+			 (cd20 + (t * t - 4.) * cd40)) 
+	}
+	return
+end
+
+procedure iievs3(x,y,n,a)   # cubic spline evaluator
+
+real	x[ARB]	# x values, must be within [1,n]
+real	y[ARB]	# interpolated values returned to user
+int	n	# number of x values
+real	a[ARB]	# basis spline coefficients - from a[0] to a[n+1]
+
+int	i,nx,nxold
+real	s,c0,c1,c2,c3
+
+begin
+	nxold = -1
+	do i = 1,n {
+	    nx = x[i]
+	    s = x[i] - nx
+	    if (nx != nxold) {
+
+		# convert b-spline coeff's to poly. coeff's
+		c0 = a[nx-1] + 4. * a[nx] + a[nx+1]
+		c1 = 3. * (a[nx+1] - a[nx-1])
+		c2 = 3. * (a[nx-1] - 2. * a[nx] + a[nx+1])
+		c3 = -a[nx-1] + 3. * a[nx] - 3. * a[nx+1] + a[nx+2]
+		nxold = nx
+	    }
+	    y[i] = c0 + s * (c1 + s * (c2 + s * c3) )
+	}
+	return
+end