4.3BSD - /usr/src/usr.bin/spline.c

   1: #ifndef lint
   2: static char *sccsid = "@(#)spline.c	4.3 (Berkeley) 9/21/85";
   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: }
 325: numb(np,argcp,argvp)
 326:     int *argcp;
 327:     float *np;
 328:     char ***argvp;{
 329:     double atof();
 330:     char c;
 331:     if(*argcp<=1) return(0);
 332:     c = (*argvp)[1][0];
 333:     if(!('0'<=c&&c<='9' || c=='-' || c== '.' )) return(0);
 334:     *np = atof((*argvp)[1]);
 335:     (*argcp)--;
 336:     (*argvp)++;
 337:     return(1); }

Defined functions

Defined variables

Defined struct's

Defined macros