Initial commit

1 files changed, 159 insertions, 0 deletions

diff --git a/math/deboor/seval.x b/math/deboor/seval.x
new file mode 100644
index 00000000..ee9cd1e9
--- /dev/null
+++ b/math/deboor/seval.x
@@ -0,0 +1,159 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+
+include <iraf.h>
+include	"bspln.h"
+
+.help seval 2 "math library"
+.ih
+NAME
+seval -- evaluate the Dth derivative of a Kth order b-spline at X.
+.ih
+USAGE
+y = seval (x, deriv, bspln)
+.ih
+PARAMETERS
+The single precision real value returned by SEVAL is the value of the D-th
+derivative of the b-spline at X.  INDEF is returned if X lies outside the
+domain of the spline.
+.ls x
+(real).  Logical x-value at which the spline is to be evaluated.
+.le
+.ls deriv
+The derivative to be evaluated.  The zeroth derivative is the function value.
+.le
+.ls bspln
+(real[2*n+30]).  B-spline descriptor structure, generated by a previous
+call to SPLINE, SPLLSQ, etc.
+.le
+.ih
+DESCRIPTION
+Calculate the Dth derivative of a Kth order B-spline.  Based on BVALUE, from
+"A Practical Guide to Splines", by C. DeBoor.  The main changes incorporated
+in SEVAL involve the use of the BSPLN array to convey all of the information
+needed to describe the spline.  This simplifies the calling sequences, and
+improves the communication between routines.  SEVAL is designed to be easy to
+use, and reasonably efficient in the case where a N point spline is to be
+evaluated N times or less.  If a spline is to be fitted and then evaluated
+many times, conversion to PP representation may be warranted to gain maximum
+efficiency.
+.ih
+METHOD
+See DeBoor, pg. 144, or the listing of BVALUE, for a detailed explanation
+of the mathematics involved.  In short, we find the knot interval containing
+X, choosing the next interval to the left if X happens to fall square on a
+knot.  The distances of X from the K-1 knots to the left and right are then
+calculated, and the K b-spline coefficients contributing to the interval are
+extracted.  Calculation of divided differences follows, if a derivative is
+to be calculated.  Finally the value of the b-spline at X is evaluated and
+returned.
+.ih
+BUGS
+The special case ND == 0 and ORDER == 4 (cubic spline) should be recognized,
+and specially optimized code used to evaluate the spline in that case.
+.ih
+SEE ALSO
+spline(2), fsplin(2), spllsq(2).
+.endhelp _____________________________________________________________
+
+
+real procedure seval (x, deriv, bspln)
+
+real	x			# logical x value
+real	bspln[ARB]		# ARB = [2 * NCOEF + 30]
+int	deriv			# derivative = 0,1,2,...,k-1
+
+int	i, j, k, n, jj, ilo, kmj, km1, offset, knots, coefs
+real	aj[KMAX], dl[KMAX], dr[KMAX], fkmj
+
+begin
+	if (x < XMIN || x > XMAX)	# x out of range?
+	    return (INDEF)
+
+	k = ORDER			# fetch k from bspln header
+	if (deriv >= k)			# Kth deriv K degree polyn is zero
+	    return (0.)
+ 
+# Find i such that X is in the interval T[i] <= X < T[i+1].  First we fetch
+# KINDEX from the bspln header.  This is the index into the knot array from
+# the last call (many, if not most, calls to SEVAL are sequential in nature).
+# Usually, KINDEX will already be the right interval, or will be one off.  If
+# more than one off, DeBoors INTERV is called to find the correct index via a
+# binary search.  INTERV handles the special cases (x == XMAX). 
+
+	i = KINDEX			# previous index for this spline
+	n = NCOEF
+	knots = KNOT1
+
+	if (x >= bspln[i+1]) {			# go right
+	    i = i + 1
+	    if (x >= bspln[i+1]) {
+		call interv (bspln[knots], n, x, i, j)
+		i = i + knots - 1
+	    }
+
+	} else if (x < bspln[i]) {		# go left
+	    i = i - 1
+	    if (x < bspln[i]) {
+		call interv (bspln[knots], n, x, i, j)
+		i = i + knots - 1
+	    }
+	}
+	KINDEX = i
+
+
+	km1 = k - 1
+	coefs = i - knots - k + COEF1
+
+	if (km1 <= 0)		# if order=1 (& deriv=0), bvalue = bspln[i].
+	    return (bspln[coefs+1])
+ 
+
+# Store the ORDER b-spline coefficients relevant for the knot interval
+# (T[i],T[i+1]) in aj[1],...,aj[k] and compute dl[j] = x - t[i+k-j],
+# dr[j] = t[i+k-1+j] - x, j=1,...,k-1.  Note that the K knots at each endpoint
+# all have the same T-value, hence we use the lookup table T, rather than
+# calculate the DL and DR from the knot spacing.
+
+	offset = i + 1
+	do j = 1, km1
+	    dl[j] = x - bspln[offset-j]
+
+	offset = offset - 1
+	do j = 1, km1
+	    dr[j] = bspln[i+j] - x
+
+	do j = 1, k		  		# fetch K coefficients
+	    aj[j] = bspln[coefs+j]
+
+
+# Difference the coefficients deriv times.
+
+	if (deriv != 0) {			# only if taking a derivative
+	    do j = 1, deriv {
+		kmj = k - j
+		fkmj = real (kmj)
+		ilo = kmj
+		do jj = 1, kmj {
+		    aj[jj] = ((aj[jj+1] - aj[jj]) / (dl[ilo] + dr[jj])) * fkmj
+		    ilo = ilo - 1
+		}
+	    }
+	}
+
+
+# Compute value at X in (t[i],t[i+])) of deriv-th derivative, given its
+# relevant b-spline coeffs in aj[1],...,aj[k-deriv].
+
+	if (deriv != km1) {
+	    do j = deriv + 1, km1 {
+		kmj = k - j
+		ilo = kmj
+		do jj = 1, kmj {
+		    aj[jj] = (aj[jj+1] * dl[ilo] + aj[jj] * dr[jj]) / (dl[ilo] + dr[jj])
+		    ilo = ilo - 1
+		}
+	    }
+	}
+	return (aj[1])
+end