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