diff options
Diffstat (limited to 'math/deboor/seval.x')
| -rw-r--r-- | math/deboor/seval.x | 159 |
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 |