~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina4/src/math/splines/ics.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /**@(#)Interpolation Cubic Spline. Functions converted from Numeric
  2  * @(#)Recipies to perform 1d cubic spline interpolation
  3  */
  4 
  5 #include <tina/sys.h>
  6 #include <tina/sysfuncs.h>
  7 
  8 void    spl_ci_gen(float *x, float *y, int n, double *ypf, double *ypl, float *y2)
  9 {
 10     int     i, k;
 11     float   p, qn, sig, un, *u;
 12 
 13     u = fvector_alloc(0, n - 1);/* intermediate values */
 14 
 15     if (ypf == NULL)            /* natural splines: 0 second derivative
 16                                  * at ends */
 17         y2[0] = u[0] = 0.0;
 18     else
 19     {
 20         y2[0] = -0.5;
 21         u[0] = (3.0 / (x[1] - x[0])) * ((y[1] - y[0]) / (x[1] - x[0]) - *ypf);
 22     }
 23 
 24     for (i = 1; i <= n - 2; i++)
 25     {
 26         sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]);
 27         p = sig * y2[i - 1] + 2.0;
 28         y2[i] = (sig - 1.0) / p;
 29         u[i] = (y[i + 1] - y[i]) / (x[i + 1] - x[i]) - (y[i] - y[i - 1]) / (x[i] - x[i - 1]);
 30         u[i] = (6.0 * u[i] / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p;
 31     }
 32 
 33     if (ypl == NULL)
 34         qn = un = 0.0;
 35     else
 36     {
 37         qn = 0.5;
 38         un = (3.0 / (x[n - 1] - x[n - 2])) * (*ypl - (y[n - 1] - y[n - 2]) / (x[n - 1] - x[n - 2]));
 39     }
 40     y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0);
 41 
 42     for (k = n - 2; k >= 0; k--)
 43         y2[k] = y2[k] * y2[k + 1] + u[k];
 44 
 45     fvector_free((void *) u, 0);
 46 }
 47 
 48 /* ARGSUSED Quieten Lint */
 49 void    spl_ci_sbp(float *x, float *y, int n, double *ypf, double *ypl, float *y2)
 50 {
 51     int     i;
 52     float  *u, *g;
 53 
 54     u = fvector_alloc(0, n);    /* intermediate values */
 55     g = fvector_alloc(0, n);    /* intermediate values */
 56 
 57     u[0] = g[0] = 0.0;          /* start condition */
 58 
 59     for (i = 1; i <= n - 2; i++)
 60     {
 61         double  a, b, c, d, temp;
 62 
 63         a = x[i] - x[i - 1];
 64         b = (x[i + 1] - x[i - 1]) / 3;
 65         c = x[i + 1] - x[i];
 66         d = (y[i + 1] - y[i]) / c - (y[i] - y[i - 1]) / a;
 67         a /= 6;
 68         c /= 6;
 69 
 70         temp = b - a * g[i - 1];
 71         u[i] = (d - a * u[i - 1]) / temp;
 72         g[i] = c / temp;
 73     }
 74 
 75     u[n - 1] = g[n - 1] = 0.0;  /* end condition */
 76 
 77     y2[n - 1] = u[n - 1];
 78     for (i = n - 2; i >= 0; i--)
 79         y2[i] = u[i] + y2[i + 1] * g[i];
 80 
 81     fvector_free((void *) u, 0);
 82     fvector_free((void *) g, 0);
 83 }
 84 
 85 double  spl_ci_val(float *xa, float *ya, float *y2a, int n, double x)
 86 {
 87     int     klo, khi, k;
 88     float   h, b, a, y;
 89 
 90     klo = 0;
 91     khi = n - 1;
 92     while (khi - klo > 1)
 93     {                           /* binary search */
 94         k = (khi + klo) >> 1;   /* mid point */
 95         if (xa[k] > x)
 96             khi = k;
 97         else
 98             klo = k;
 99     }
100     h = xa[khi] - xa[klo];      /* should be strictly increasing */
101     if (h == 0.0)
102     {
103         error("spl_ci_val: non monotonic absisa", warning);
104         return (0.0);
105     }
106     a = (xa[khi] - x) / h;
107     b = (x - xa[klo]) / h;
108     y = a * ya[klo] + b * ya[khi] + ((a * a * a - a) * y2a[klo] + (b * b * b - b) * y2a[khi]) * (h * h) / 6.0;
109     return (y);
110 }
111 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.