# 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