blob: c04a94deae6292215bfc387e5fc66b7459e4a080 (
plain) (
blame)
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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
# II_SPLINE -- This procedure fits uniformly spaced data with a cubic
# spline. The spline is given as basis-spline coefficients that replace
# the data values.
#
# Storage at call time:
#
# bcoeff[1] = second derivative at x = 1
# bcoeff[2] = first data point y[1]
# bcoeff[3] = y[2]
#
# bcoeff[n+1] = y[n]
# bcoeff[n+2] = second derivative at x = n
#
# Storage after call:
#
# bcoeff[1] ... bcoeff[n+2] = the n + 2 basis-spline coefficients for the
# basis splines as defined in P.M. Prenter's book "Splines and Variational
# Methods", Wiley, 1975.
procedure ii_spline (bcoeff, diag, npts)
real bcoeff[ARB] # data in and also bspline coefficients out
real diag[ARB] # needed for offdiagnol matrix elements
int npts # number of data points
int i
begin
diag[1] = -2.
bcoeff[1] = bcoeff[1] / 6.
diag[2] = 0.
bcoeff[2] = (bcoeff[2] - bcoeff[1]) / 6.
# Gaussian elimination - diagnol below main is made zero
# and main diagnol is made all 1's
do i = 3, npts + 1 {
diag[i] = 1. / (4. - diag[i-1])
bcoeff[i] = diag[i] * (bcoeff[i] - bcoeff[i-1])
}
# Find last b spline coefficient first - overlay r.h.s.'s
bcoeff[npts+2] = ((diag[npts] + 2.) * bcoeff[npts+1] - bcoeff[npts] +
bcoeff[npts+2] / 6.) / (1. + diag[npts+1] * (diag[npts] + 2.))
# back substitute filling in coefficients for b splines
# note bcoeff[npts+1] is evaluated correctly as can be checked
# bcoeff[2] is already set since offdiagnol is 0.
do i = npts + 1, 3, -1
bcoeff[i] = bcoeff[i] - diag[i] * bcoeff[i+1]
# evaluate bcoeff[1]
bcoeff[1] = bcoeff[1] + 2. * bcoeff[2] - bcoeff[3]
end
|
