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