1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
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
|
