Return to spline.c CVS log

Up to [CSRG BSD Unix] / 42BSD / usr.bin

Annotation of 42BSD/usr.bin/spline.c, revision 1.1

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

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.