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 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 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 } 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 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 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 } 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