home *** CD-ROM | disk | FTP | other *** search
/ minnie.tuhs.org / unixen.tar / unixen / PDP-11 / Trees / V7 / usr / src / cmd / spline.c < prev    next >
Encoding:
C/C++ Source or Header  |  1979-01-10  |  7.7 KB  |  334 lines
  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 hi=xi-xi-1
  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.     hiy"i-1+2(hi+hi+1)y"i+hi+1y"i+1
  24.     
  25.     = 6[(yi+1-yi)/hi+1-(yi-yi-1)/hi]   i=1,2,...,n
  26.  
  27. where y" = y"n+1 = 0
  28. This is a symmetric tridiagonal system of the form
  29.  
  30.     | a h               |  |y"|      |b|
  31.     | h a h            |  |y"|      |b|
  32.     |    h a h         |  |y"|  =   |b|
  33.     |         .           |  | .|      | .|
  34.     |            .        |  | .|      | .|
  35. It can be triangularized into
  36.     | d h               |  |y"|      |r|
  37.     |    d h            |  |y"|      |r|
  38.     |       d h         |  |y"|  =   |r|
  39.     |          .          |  | .|      | .|
  40.     |             .       |  | .|      | .|
  41. where
  42.     d = a
  43.  
  44.     r = 0
  45.  
  46.     di = ai - hi/di-1    1<i<_n
  47.  
  48.     ri = bi - hiri-1/di-1i    1<_i<_n
  49.  
  50. the back solution is
  51.     y"n = rn/dn
  52.  
  53.     y"i = (ri-hi+1y"i+1)/di    1<_i<n
  54.  
  55. superficially, di and ri don't have to be stored for they can be
  56. recalculated backward by the formulas
  57.  
  58.     di-1 = hi/(ai-di)    1<i<_n
  59.  
  60.     ri-1 = (bi-ri)di-1/hi    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.     y" = ky", yn+1"   = kyn"
  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 h to a and hn+1 to an.
  81.  
  82.  
  83. Periodic case_____________
  84.  
  85. To do this, add 1 more row and column thus
  86.  
  87.     | a h            h |  |y"|     |b|
  88.     | h a h            |  |y"|     |b|
  89.     |    h a h         |  |y"|     |b|
  90.     |                     |  | .|  =  | .|
  91.     |             .       |  | .|     | .|
  92.     | h            h a |  | .|     | .|
  93.  
  94. where h=_ hn+1
  95.  
  96. The same diagonalization procedure works, except for
  97. the effect of the 2 corner elements.  Let si be the part
  98. of the last element in the ith "diagonalized" row that
  99. arises from the extra top corner element.
  100.  
  101.         s = h
  102.  
  103.         si = -si-1hi/di-1    2<_i<_n+1
  104.  
  105. After "diagonalizing", the lower corner element remains.
  106. Call ti the bottom element that appears in the ith colomn
  107. as the bottom element to its left is eliminated
  108.  
  109.         t = h
  110.  
  111.         ti = -ti-1hi/di-1
  112.  
  113. Evidently ti = si.
  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.     u = v = 0
  120.  
  121.     ui = ui-1-si-1*ti-1/di-1
  122.  
  123.     vi = vi-1-ri-1*ti-1/di-1    2<_i<_n+1
  124.  
  125. The back solution is now obtained as follows
  126.  
  127.     y"n+1 = (rn+1+vn+1)/(dn+1+sn+1+tn+1+un+1)
  128.  
  129.     y"i = (ri-hi+1*yi+1-si*yn+1)/di    1<_i<_n
  130.  
  131. Interpolation in the interval xi<_x<_xi+1 is by the formula
  132.  
  133.     y = yix+ + yi+1x- -(hi+1/6)[y"i(x+-x+3)+y"i+1(x--x-)]
  134. where
  135.     x+ = xi+1-x
  136.  
  137.     x- = x-xi
  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); }
  333.  
  334.