-typedef struct {
- double a, b, c, d;
-} spline_coeff_t;
-
-void compute_spline(int n, double *x, double *f, spline_coeff_t *p)
-{
- double *h, *alpha, *B, *m;
- int i;
- int orig_n = n;
- double new_x[3], new_f[3];
-
- if (n==0) return;
- if (n==1) {
- new_x[0] = x[0];
- new_f[0] = f[0];
- new_x[1] = x[0]+0.00001;
- new_f[1] = f[0];
- x = new_x;
- f = new_f;
- n = 3;
- }
- if (n==2) {
- new_x[0] = x[0];
- new_f[0] = f[0];
- new_x[1] = x[1];
- new_f[1] = f[1];
- new_x[2] = x[1] + x[1]-x[0];
- new_f[2] = f[1] + f[1]-f[0];
- x = new_x;
- f = new_f;
- n = 3;
- }
-
- /* we're solving a linear system of equations of the form Ax = B.
- * The matrix a is tridiagonal and consists of coefficients in
- * the h[] and alpha[] arrays.
- */
-
- h = (double *)malloc(sizeof(double) * (n-1));
- for (i=0; i<n-1; i++) {
- h[i] = x[i+1]-x[i];
- }
-
- alpha = (double *)malloc(sizeof(double) * (n-2));
- for (i=0; i<n-2; i++) {
- alpha[i] = 2 * (h[i] + h[i+1]);
- }
-
- /* B[] is the vector on the right hand side of the equation */
- B = (double *)malloc(sizeof(double) * (n-2));
- for (i=0; i<n-2; i++) {
- B[i] = 6 * ((f[i+2] - f[i+1])/h[i+1] - (f[i+1] - f[i])/h[i]);
- }
-
- /* Now solve the n-2 by n-2 system */
- m = (double *)malloc(sizeof(double) * (n-2));
- for (i=1; i<=n-3; i++) {
- /*
- d0 = alpha 0
- a0 = h1
- c0 = h1
-
- di = di - (ai-1 / di-1) * ci-1
- bi = bi - (ai-1 / di-1) * bi-1
- ;
- */
- alpha[i] = alpha[i] - (h[i]/alpha[i-1]) * h[i];
- B[i] = B[i] - (h[i]/alpha[i-1]) * B[i-1];
- }
- /* xn-3 = bn-3 / dn-3; */
- m[n-3] = B[n-3]/alpha[n-3];
- for (i=n-4; i>=0; i--) {
- m[i] = (B[i]-h[i+1]*m[i+1])/alpha[i];
- }
-
- for (i=0; i<orig_n-1; i++) {
- double mi, mi1;
- mi = (i==(n-2)) ? 0 : m[i];
- mi1 = (i==0) ? 0 : m[i-1];
-
- p[i].a = f[i+1];
- p[i].b = (f[i+1] - f[i]) / h[i] + h[i] * (2*mi + mi1) / 6;
- p[i].c = mi/2;
- p[i].d = (mi-mi1)/(6*h[i]);
- }
-
- free(alpha);
- free(B);
- free(h);
- free(m);
-}