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
|
subroutine splev(t,n,c,k,x,y,m,e,ier)
c subroutine splev evaluates in a number of points x(i),i=1,2,...,m
c a spline s(x) of degree k, given in its b-spline representation.
c
c calling sequence:
c call splev(t,n,c,k,x,y,m,e,ier)
c
c input parameters:
c t : array,length n, which contains the position of the knots.
c n : integer, giving the total number of knots of s(x).
c c : array,length n, which contains the b-spline coefficients.
c k : integer, giving the degree of s(x).
c x : array,length m, which contains the points where s(x) must
c be evaluated.
c m : integer, giving the number of points where s(x) must be
c evaluated.
c e : integer, if 0 the spline is extrapolated from the end
c spans for points not in the support, if 1 the spline
c evaluates to zero for those points, if 2 ier is set to
c 1 and the subroutine returns, and if 3 the spline evaluates
c to the value of the nearest boundary point.
c
c output parameter:
c y : array,length m, giving the value of s(x) at the different
c points.
c ier : error flag
c ier = 0 : normal return
c ier = 1 : argument out of bounds and e == 2
c ier =10 : invalid input data (see restrictions)
c
c restrictions:
c m >= 1
c-- t(k+1) <= x(i) <= x(i+1) <= t(n-k) , i=1,2,...,m-1.
c
c other subroutines required: fpbspl.
c
c references :
c de boor c : on calculating with b-splines, j. approximation theory
c 6 (1972) 50-62.
c cox m.g. : the numerical evaluation of b-splines, j. inst. maths
c applics 10 (1972) 134-149.
c dierckx p. : curve and surface fitting with splines, monographs on
c numerical analysis, oxford university press, 1993.
c
c author :
c p.dierckx
c dept. computer science, k.u.leuven
c celestijnenlaan 200a, b-3001 heverlee, belgium.
c e-mail : Paul.Dierckx@cs.kuleuven.ac.be
c
c latest update : march 1987
c
c++ pearu: 11 aug 2003
c++ - disabled cliping x values to interval [min(t),max(t)]
c++ - removed the restriction of the orderness of x values
c++ - fixed initialization of sp to double precision value
c
c ..scalar arguments..
integer n, k, m, e, ier
c ..array arguments..
real*8 t(n), c(n), x(m), y(m)
c ..local scalars..
integer i, j, k1, l, ll, l1, nk1
c++..
integer k2
c..++
real*8 arg, sp, tb, te
c ..local array..
real*8 h(20)
c ..
c before starting computations a data check is made. if the input data
c are invalid control is immediately repassed to the calling program.
ier = 10
c-- if(m-1) 100,30,10
c++..
if (m .lt. 1) go to 100
c..++
c-- 10 do 20 i=2,m
c-- if(x(i).lt.x(i-1)) go to 100
c-- 20 continue
ier = 0
c fetch tb and te, the boundaries of the approximation interval.
k1 = k + 1
c++..
k2 = k1 + 1
c..++
nk1 = n - k1
tb = t(k1)
te = t(nk1 + 1)
l = k1
l1 = l + 1
c main loop for the different points.
do 80 i = 1, m
c fetch a new x-value arg.
arg = x(i)
c check if arg is in the support
if (arg .lt. tb .or. arg .gt. te) then
if (e .eq. 0) then
goto 35
else if (e .eq. 1) then
y(i) = 0
goto 80
else if (e .eq. 2) then
ier = 1
goto 100
else if (e .eq. 3) then
if (arg .lt. tb) then
arg = tb
else
arg = te
endif
endif
endif
c search for knot interval t(l) <= arg < t(l+1)
c++..
35 if (arg .ge. t(l) .or. l1 .eq. k2) go to 40
l1 = l
l = l - 1
go to 35
c..++
40 if(arg .lt. t(l1) .or. l .eq. nk1) go to 50
l = l1
l1 = l + 1
go to 40
c evaluate the non-zero b-splines at arg.
50 call fpbspl(t, n, k, arg, l, h)
c find the value of s(x) at x=arg.
sp = 0.0d0
ll = l - k1
do 60 j = 1, k1
ll = ll + 1
sp = sp + c(ll)*h(j)
60 continue
y(i) = sp
80 continue
100 return
end
|
