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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include "bspln.h"
.help spline 2 "math library"
.ih
NAME
spline -- generalized spline interpolation with b-splines.
.ih
USAGE
spline (x, y, n, bspln, q, k, ier)
.ih
PARAMETERS
.ls x,y
(real[n]). Abcissae and ordinates of the N data points.
.le
.ls n
The number of input data points, and the number of spline coefficients
to be generated.
.le
.ls bspln
(real[2*n+30]). Output array of N b-spline coefficients, N+K knots,
and the spline descriptor array.
.le
.ls q
(real[(2*K-1)*N]). On output, will contain the
triangular factorization of the coefficient matrix of the linear
system for the b-coefficients of the spline interpolant. If a new
data set differs only in the Y values, Q may be input to FSPLIN to
efficiently solve for the b-spline coefficients of the new data set.
.le
.ls k
The order of the spline (2-20). Must be EVEN at present
(k=4 is the cubic spline).
.le
.ls ier
Zero if ok, nonzero if error (invalid input parameters).
.le
.ih
DESCRIPTION
General spline interpolation. SPLINE is a thinly disguised version of SPLINT
(DeBoor, pg. 204). SPLINE differs from SPLINT in that the knot spacing T is
calculated automatically, using the not-a-knot boundary conditions. Only even
K are permitted at present (K = 2, 4, 6,...). The cubic spline interpolant
corresponds to K = 4. SPLINE saves a complete description of the spline in the
BSPLN array, for use later by SEVAL (used to evaluate the fitted spline) and
FSPLIN (used after a call to SPLINE to more efficiently fit subsequent data
sets).
.ih
SOURCE
See the listing of SPLINT, or Chapter XIII of DeBoors book.
.ih
SEE ALSO
seval(2), fsplin(2).
.endhelp ______________________________________________________________
procedure spline (x, y, n, bspln, q, k, ier)
real x[n], y[n]
int n, k, ier
real q[ARB], bspln[ARB] #q[(2*k-1)*n], bspln[2*n+30]
int m, i, knot
begin
if (k < 2 || mod (k, 2) != 0) {
ier = 1
return
}
NCOEF = n #set up spline descriptor
ORDER = k
XMIN = x[1]
XMAX = x[n]
knot = int (KNOT1) - 1 #offset to knot array in bspln
KINDEX = knot + k #initial posn for SEVAL
m = knot + n
do i = 1, k { #not-a-knot knot boundary conditions
bspln[knot+i] = XMIN
bspln[m+i] = XMAX
}
m = k / 2 #use data x-values inside
do i = k+1, n
bspln[knot+i] = x[i+m-k]
call splint (x, y, bspln[knot+1], n, k, q, bspln[COEF1], ier)
if (ier == 1) #1=success, 2==failure
ier = 0
end
|
