1 c
2 c http://www.netlib.org/fmm/spline.f
3 c
4 subroutine spline(n, x, y, b, c, d)
5 integer n
6 double precision x(n), y(n), b(n), c(n), d(n)
7 c
8 integer nmi, ib, i
9 double precision t
10 c
11 nmi = n-1
12 if ( n .lt. 2 ) return
13 if ( n .lt. 3 ) goto 50
14 c
15 d(1) = x(2) - x(1)
16 c(2) = ( y(2) - y(1) ) / d(1)
17 do 10 i = 2,nmi
18 d(i) = x(i+1) - x(i)
19 b(i) = 2.0*(d(i-1)+d(i))
20 c(i+1) = (y(i+1)-y(i))/d(i)
21 c(i) = c(i+1)-c(i)
22 10 continue
23 c
24 b(1) = -d(1)
25 b(n) = -d(n-1)
26 c(1) = 0.0
27 c(n) = 0.0
28 if ( n .eq. 3 ) goto 15
29 c(1) = c(3) / (x(4)-x(2)) - c(2)/(x(3)-x(1))
30 c(n) = c(n-1)/(x(n)-x(n-2)) - c(n-2)/(x(n-1)-x(n-3))
31 c(1) = c(1)*d(1)**2/(x(4)-x(1))
32 c(n) = -c(n)*d(n-1)**2/(x(n)-x(n-3))
33 c
34 15 do 20 i = 2, n
35 t = d(i-1)/b(i-1)
36 b(i) = b(i) - t*d(i-1)
37 c(i) = c(i) - t*c(i-1)
38 20 continue
39 c
40 c(n) = c(n)/b(n)
41 do 30 ib = 1, nmi
42 i = n - ib
43 c(i) = (c(i) - d(i)*c(i+1))/b(i)
44 30 continue
45 c
46 b(n) = (y(n) - y(nmi))/d(nmi) + d(nmi)*(c(nmi) + 2.0*c(n))
47 do 40 i = 1, nmi
48 b(i) = (y(i+1) - y(i))/d(i) - d(i)*(c(i+1) + 2.0*c(i))
49 d(i) = (c(i+1)-c(i))/d(i)
50 c(i) = 3.0*c(i)
51 40 continue
52 c(n) = 3.0*c(n)
53 d(n) = d(n-1)
54 return
55 c
56 50 b(i) = (y(2)-y(1))/(x(2)-x(1))
57 c(1) = 0.0
58 d(1) = 0.0
59 b(2) = b(1)
60 c(2) = 0.0
61 d(2) = 0.0
62 return
63 end
64 c
65 double precision function seval(n, u, x, y, b, c, d)
66 integer n
67 double precision u, x(n), y(n), b(n), c(n), d(n)
68 c
69 integer i, j, k
70 double precision dx
71 data i/1/
72 if ( i .ge. n ) i = 1
73 if ( u .lt. x(i) ) goto 10
74 if ( u .le. x(i+1) ) goto 30
75 10 i = 1
76 j = n + 1
77 20 k = (i+j)/2
78 if ( u .lt. x(k) ) j = k
79 if ( u .ge. x(k) ) i = k
80 if ( j .gt. i+1 ) goto 20
81 c
82 30 dx = u - x(i)
83 seval = y(i) + dx*(b(i) + dx*(c(i) + dx*d(i)))
84 return
85 end
86
87 c same but returns interpolation segment index
88 c side effect: f=f(u) and fp=f'(u)
89 integer function seval1(n, u, x, y, b, c, d, f, fp)
90 integer n
91 double precision u, x(n), y(n), b(n), c(n), d(n)
92 double precision f, fp
93 c
94 integer i, j, k
95 double precision dx
96 data i/1/
97 if ( i .ge. n ) i = 1
98 if ( u .lt. x(i) ) goto 10
99 if ( u .le. x(i+1) ) goto 30
100 10 i = 1
101 j = n + 1
102 20 k = (i+j)/2
103 if ( u .lt. x(k) ) j = k
104 if ( u .ge. x(k) ) i = k
105 if ( j .gt. i+1 ) goto 20
106 c
107 30 dx = u - x(i)
108 f = y(i) + dx*(b(i) + dx*(c(i) + dx*d(i)))
109 fp = b(i)+dx*(c(i)*2+dx*(d(i)*3))
110 seval1 = i
111 return
112 end