+ } else
+ *up = *down = VIK_DEFAULT_ALTITUDE;
+}
+
+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];