# 2.11BSD#431 - /usr/src/usr.bin/spline.c

##### Functions

    getfloat          X
    getlim            X
    main              X
    numb              X
    readin            X
    rhs               X
    spline            X

##### Variables

    auta              X
    diag              X
    dx                X
    konst             X
    n                 X
    ni                X
    periodic          X
    r                 X
    x                 X
    y                 X
    zero              X

##### Macros

    INF               X
    NP                X

##### struct's

    proj              X

       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

getfloat defined in line 227; used 2 times

- in line 224-225(2)

getlim defined in line 266; used 2 times

- in line 316-317(2)

main defined in line 275; never used
numb defined in line 325; used 5 times

- in line 292-298(3), 304-306(2)

readin defined in line 221; used 1 times

- in line 315

rhs defined in line 145; used 1 times

- in line 173

spline defined in line 154; used 1 times

- in line 321

#### Defined variables

auta defined in line 16; used 3 times

- in line 223, 291, 314

diag defined in line 12; used 4 times

- in line 178, 191, 197, 319

dx defined in line 13; used 3 times

- in line 223, 292, 318

konst defined in line 18; used 8 times

- in line 162, 176-177(2), 195-196(2), 202-203(2), 295

n defined in line 15; used 19 times

- in line 149, 161-168(3), 177-183(3), 196, 203, 222-225(7), 269, 318-321(2)

ni defined in line 14; used 3 times

- in line 206(2), 298

periodic defined in line 17; used 6 times

- in line 162-166(3), 182, 201, 301

r defined in line 12; used 7 times

- in line 164, 172-173(3), 190, 320-321(2)

sccsid defined in line 2; never used
x defined in line 11; used 34 times

- in line 150-151(4), 167-169(6), 184-188(6), 207(2), 215, 223-224(3), 279-281(4), 304-307(4), 314-316(3), 322

y defined in line 11; used 13 times

- in line 150-151(4), 214(2), 225, 279-283(4), 317, 323

zero defined in line 19; used 3 times

- in line 171-172(2), 189

#### Defined struct's

proj defined in line 11; used 2 times

- in line 267(2)

#### Defined macros

INF defined in line 9; used 4 times

- in line 280-283(4)

NP defined in line 8; used 2 times

- in line 11, 222

Last modified: 1987-02-17 Generated: 2016-12-26 | Generated by src2html V0.67 page hit count: 1357 |