1 /*-
2 * Copyright (c) 1991, 1993
3 * The Regents of the University of California. All rights reserved.
4 *
5 * %sccs.include.proprietary.c%
6 */
7
8 #ifndef lint
9 static char copyright[] =
10 "@(#) Copyright (c) 1991, 1993\n\
11 The Regents of the University of California. All rights reserved.\n";
12 #endif /* not lint */
13
14 #ifndef lint
15 static char sccsid[] = "@(#)spline.c 8.1 (Berkeley) 06/06/93";
16 #endif /* not lint */
17
18 #include <stdio.h>
19 #include <math.h>
20
21 #define NP 1000
22
23 struct proj { int lbf,ubf; float a,b,lb,ub,quant,mult,val[NP]; } x,y;
24 float *diag, *r;
25 float dx = 1.;
26 float ni = 100.;
27 int n;
28 int auta;
29 int periodic;
30 float konst = 0.0;
31 float zero = 0.;
32
33 /* Spline fit technique
34 let x,y be vectors of abscissas and ordinates
35 h be vector of differences hi=xi-xi-1
36 y" be vector of 2nd derivs of approx function
37 If the points are numbered 0,1,2,...,n+1 then y" satisfies
38 (R W Hamming, Numerical Methods for Engineers and Scientists,
39 2nd Ed, p349ff)
40 hiy"i-1+2(hi+hi+1)y"i+hi+1y"i+1
41
42 = 6[(yi+1-yi)/hi+1-(yi-yi-1)/hi] i=1,2,...,n
43
44 where y"0 = y"n+1 = 0
45 This is a symmetric tridiagonal system of the form
46
47 | a1 h2 | |y"1| |b1|
48 | h2 a2 h3 | |y"2| |b2|
49 | h3 a3 h4 | |y"3| = |b3|
50 | . | | .| | .|
51 | . | | .| | .|
52 It can be triangularized into
53 | d1 h2 | |y"1| |r1|
54 | d2 h3 | |y"2| |r2|
55 | d3 h4 | |y"3| = |r3|
56 | . | | .| | .|
57 | . | | .| | .|
58 where
59 d1 = a1
60
61 r0 = 0
62
63 di = ai - hi2/di-1 1<i<_n
64
65 ri = bi - hiri-1/di-1i 1<_i<_n
66
67 the back solution is
68 y"n = rn/dn
69
70 y"i = (ri-hi+1y"i+1)/di 1<_i<n
71
72 superficially, di and ri don't have to be stored for they can be
73 recalculated backward by the formulas
74
75 di-1 = hi2/(ai-di) 1<i<_n
76
77 ri-1 = (bi-ri)di-1/hi 1<i<_n
78
79 unhappily it turns out that the recursion forward for d
80 is quite strongly geometrically convergent--and is wildly
81 unstable going backward.
82 There's similar trouble with r, so the intermediate
83 results must be kept.
84
85 Note that n-1 in the program below plays the role of n+1 in the theory
86
87 Other boundary conditions_________________________
88
89 The boundary conditions are easily generalized to handle
90
91 y0" = ky1", yn+1" = kyn"
92
93 for some constant k. The above analysis was for k = 0;
94 k = 1 fits parabolas perfectly as well as stright lines;
95 k = 1/2 has been recommended as somehow pleasant.
96
97 All that is necessary is to add h1 to a1 and hn+1 to an.
98
99
100 Periodic case_____________
101
102 To do this, add 1 more row and column thus
103
104 | a1 h2 h1 | |y1"| |b1|
105 | h2 a2 h3 | |y2"| |b2|
106 | h3 a4 h4 | |y3"| |b3|
107 | | | .| = | .|
108 | . | | .| | .|
109 | h1 h0 a0 | | .| | .|
110
111 where h0=_ hn+1
112
113 The same diagonalization procedure works, except for
114 the effect of the 2 corner elements. Let si be the part
115 of the last element in the ith "diagonalized" row that
116 arises from the extra top corner element.
117
118 s1 = h1
119
120 si = -si-1hi/di-1 2<_i<_n+1
121
122 After "diagonalizing", the lower corner element remains.
123 Call ti the bottom element that appears in the ith colomn
124 as the bottom element to its left is eliminated
125
126 t1 = h1
127
128 ti = -ti-1hi/di-1
129
130 Evidently ti = si.
131 Elimination along the bottom row
132 introduces further corrections to the bottom right element
133 and to the last element of the right hand side.
134 Call these corrections u and v.
135
136 u1 = v1 = 0
137
138 ui = ui-1-si-1*ti-1/di-1
139
140 vi = vi-1-ri-1*ti-1/di-1 2<_i<_n+1
141
142 The back solution is now obtained as follows
143
144 y"n+1 = (rn+1+vn+1)/(dn+1+sn+1+tn+1+un+1)
145
146 y"i = (ri-hi+1*yi+1-si*yn+1)/di 1<_i<_n
147
148 Interpolation in the interval xi<_x<_xi+1 is by the formula
149
150 y = yix+ + yi+1x- -(h2i+1/6)[y"i(x+-x+3)+y"i+1(x--x-3)]
151 where
152 x+ = xi+1-x
153
154 x- = x-xi
155 */
156
157 float
rhs(i)158 rhs(i){
159 int i_;
160 double zz;
161 i_ = i==n-1?0:i;
162 zz = (y.val[i]-y.val[i-1])/(x.val[i]-x.val[i-1]);
163 return(6*((y.val[i_+1]-y.val[i_])/(x.val[i+1]-x.val[i]) - zz));
164 }
165
spline()166 spline(){
167 float d,s,u,v,hi,hi1;
168 float h;
169 float D2yi,D2yi1,D2yn1,x0,x1,yy,a;
170 int end;
171 float corr;
172 int i,j,m;
173 if(n<3) return(0);
174 if(periodic) konst = 0;
175 d = 1;
176 r[0] = 0;
177 s = periodic?-1:0;
178 for(i=0;++i<n-!periodic;){ /* triangularize */
179 hi = x.val[i]-x.val[i-1];
180 hi1 = i==n-1?x.val[1]-x.val[0]:
181 x.val[i+1]-x.val[i];
182 if(hi1*hi<=0) return(0);
183 u = i==1?zero:u-s*s/d;
184 v = i==1?zero:v-s*r[i-1]/d;
185 r[i] = rhs(i)-hi*r[i-1]/d;
186 s = -hi*s/d;
187 a = 2*(hi+hi1);
188 if(i==1) a += konst*hi;
189 if(i==n-2) a += konst*hi1;
190 diag[i] = d = i==1? a:
191 a - hi*hi/d;
192 }
193 D2yi = D2yn1 = 0;
194 for(i=n-!periodic;--i>=0;){ /* back substitute */
195 end = i==n-1;
196 hi1 = end?x.val[1]-x.val[0]:
197 x.val[i+1]-x.val[i];
198 D2yi1 = D2yi;
199 if(i>0){
200 hi = x.val[i]-x.val[i-1];
201 corr = end?2*s+u:zero;
202 D2yi = (end*v+r[i]-hi1*D2yi1-s*D2yn1)/
203 (diag[i]+corr);
204 if(end) D2yn1 = D2yi;
205 if(i>1){
206 a = 2*(hi+hi1);
207 if(i==1) a += konst*hi;
208 if(i==n-2) a += konst*hi1;
209 d = diag[i-1];
210 s = -s*d/hi;
211 }}
212 else D2yi = D2yn1;
213 if(!periodic) {
214 if(i==0) D2yi = konst*D2yi1;
215 if(i==n-2) D2yi1 = konst*D2yi;
216 }
217 if(end) continue;
218 m = hi1>0?ni:-ni;
219 m = 1.001*m*hi1/(x.ub-x.lb);
220 if(m<=0) m = 1;
221 h = hi1/m;
222 for(j=m;j>0||i==0&&j==0;j--){ /* interpolate */
223 x0 = (m-j)*h/hi1;
224 x1 = j*h/hi1;
225 yy = D2yi*(x0-x0*x0*x0)+D2yi1*(x1-x1*x1*x1);
226 yy = y.val[i]*x0+y.val[i+1]*x1 -hi1*hi1*yy/6;
227 printf("%f ",x.val[i]+j*h);
228 printf("%f\n",yy);
229 }
230 }
231 return(1);
232 }
readin()233 readin() {
234 for(n=0;n<NP;n++){
235 if(auta) x.val[n] = n*dx+x.lb;
236 else if(!getfloat(&x.val[n])) break;
237 if(!getfloat(&y.val[n])) break; } }
238
getfloat(p)239 getfloat(p)
240 float *p;{
241 char buf[30];
242 register c;
243 int i;
244 extern double atof();
245 for(;;){
246 c = getchar();
247 if (c==EOF) {
248 *buf = '\0';
249 return(0);
250 }
251 *buf = c;
252 switch(*buf){
253 case ' ':
254 case '\t':
255 case '\n':
256 continue;}
257 break;}
258 for(i=1;i<30;i++){
259 c = getchar();
260 if (c==EOF) {
261 buf[i] = '\0';
262 break;
263 }
264 buf[i] = c;
265 if('0'<=c && c<='9') continue;
266 switch(c) {
267 case '.':
268 case '+':
269 case '-':
270 case 'E':
271 case 'e':
272 continue;}
273 break; }
274 buf[i] = ' ';
275 *p = atof(buf);
276 return(1); }
277
278 getlim(p)
279 struct proj *p; {
280 register int i;
281 if (!p->lbf) {
282 p->lb = p->val[0];
283 for(i=1;i<n;i++) if (p->lb>p->val[i]) p->lb = p->val[i]; }
284 if (!p->ubf) {
285 p->ub = p->val[0];
286 for(i=1;i<n;i++) if (p->ub<p->val[i]) p->ub = p->val[i]; }
287 }
288
289
main(argc,argv)290 main(argc,argv)
291 char *argv[];{
292 extern char *malloc();
293 int i;
294 x.lbf = x.ubf = y.lbf = y.ubf = 0;
295 x.lb = x.ub = y.lb = y.ub = 0;
296 while(--argc > 0) {
297 argv++;
298 again: switch(argv[0][0]) {
299 case '-':
300 argv[0]++;
301 goto again;
302 case 'a':
303 auta = 1;
304 numb(&dx,&argc,&argv);
305 break;
306 case 'k':
307 numb(&konst,&argc,&argv);
308 break;
309 case 'n':
310 numb(&ni,&argc,&argv);
311 break;
312 case 'p':
313 periodic = 1;
314 break;
315 case 'x':
316 if(!numb(&x.lb,&argc,&argv)) break;
317 x.lbf = 1;
318 if(!numb(&x.ub,&argc,&argv)) break;
319 x.ubf = 1;
320 break;
321 default:
322 (void)fprintf(stderr, "spline: illegal option -- %c\n",
323 argv[0][0]);
324 (void)fprintf(stderr, "usage: spline [-aknpx]\n");
325 exit(1);
326 }
327 }
328 if(auta&&!x.lbf) x.lb = 0;
329 readin();
330 getlim(&x);
331 getlim(&y);
332 i = (n+1)*sizeof(dx);
333 diag = (float *)malloc((unsigned)i);
334 r = (float *)malloc((unsigned)i);
335 if(r==NULL||!spline()) for(i=0;i<n;i++){
336 printf("%f ",x.val[i]);
337 printf("%f\n",y.val[i]); }
338 exit(0);
339 }
numb(np,argcp,argvp)340 numb(np,argcp,argvp)
341 int *argcp;
342 float *np;
343 char ***argvp;{
344 double atof();
345 char c;
346 if(*argcp<=1) return(0);
347 c = (*argvp)[1][0];
348 if(!('0'<=c&&c<='9' || c=='-' || c== '.' )) return(0);
349 *np = atof((*argvp)[1]);
350 (*argcp)--;
351 (*argvp)++;
352 return(1); }
353