UNIX 7th Ed. - /usr/src/cmd/spline.c

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

Defined functions

Defined variables

Defined struct's

Defined macros