# include "data.h"

/* construct the natural spline for x, y in y2 */
void spline (float *x, float *y, int N, float *y2) {

  int i;
  float dy, dx, *tmp;
  
  ALLOCATE (tmp, float, N);

  y2[0] = tmp[0] = 0.0;
  
  for (i = 1; i < N-1; i++) {
    dx = (x[i+0] - x[i-1]) / (x[i+1] - x[i-1]);
    dy = dx * y2[i-1] + 2.0;
    y2[i] = (dx - 1.0) / dy;
    tmp[i] = (y[i+1] - y[i+0]) / (x[i+1] - x[i+0]) - (y[i+0] - y[i-1]) / (x[i+0] - x[i-1]);
    tmp[i] = (6.0 * tmp[i] / (x[i+1] - x[i-1]) - dx*tmp[i-1]) / dy;
  }
  
  y2[N-1] = 0;
  for (i = N-2; i >= 1; i--)
    y2[i] = y2[i]*y2[i+1] + tmp[i];

  free (tmp);
}

/* evaluate spline for x, y, y2 at X */
float splint (float *x, float *y, float *y2, int N, float X) {

  int i, lo, hi;
  float dx, a, b, value;
  
  /* find correct element in array (x must be sorted) */
  lo = 0;
  hi = N-1;
  while (hi - lo > 1) {
    i = (hi+lo) >> 1;
    if (x[i] > X) 
      hi = i;
    else 
      lo = i;
  }

  /* error condition: duplicate abssisca */
  dx = x[hi] - x[lo];
  if (dx == 0.0) {
    return (HUGE_VAL);
  }

  /* evaluate spline */
  a = (x[hi] - X) / dx;
  b = (X - x[lo]) / dx;

  value = a*y[lo] + b*y[hi] + ((a*a*a - a)*y2[lo] + (b*b*b - b)*y2[hi])*(dx*dx) / 6.0;
  return (value);

}
