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