double spline_eval(spline *s, double t)
{
  int    i, lo, mid, hi;
  double h, value;

  if ( s == NULL )
    spline_error("spline_eval: null pointer passed");

  /* if length == 1 we return ylist[0] */
  /* if length == 0 we raise an error */
  if ( s->length <= 0 )
    spline_error("spline_eval: length <= 0");
  if ( s->length == 1 )
    return s->ylist[0];

  /* we compute i so that xlist[i] <= t <= xlist[i+1], 
     or set i=0 if t < xlist[0], 
     or set i=length-2 if t > xlist[length-1] */
  lo = 0;
  hi = s->length - 1;

  if ( t <= s->xlist[lo] )
    i = 0;
  else if ( t >= s->xlist[hi] )
    i = s->length - 2;
  else
  {
    /* use binary search */
    while ( hi > lo + 1 )
    { /* loop invariant: 
         s->xlist[lo] <= t <= s->xlist[hi] */
      mid = (lo+hi)/2;
      if ( t <= s->xlist[mid] )
        hi = mid;
      else /* t > x->xlist[mid] */
        lo = mid;
    }
    i = lo;
  }

  /* now use the standard formula for 
     interval [xlist[i], xlist[i+1]] */
  h = s->xlist[i+1] - s->xlist[i];
  diff1 = t - s->xlist[i];
  diff2 = s->xlist[i+1] - t;
  value = diff1*(s->ylist[i+1]/h - s->Mlist[i+1]*h/6.0
            + diff1*diff1*s->Mlist[i+1]/(6.0*h))
        + diff2*(s->ylist[i  ]/h - s->Mlist[i  ]*h/6.0
            + diff2*diff2*s->Mlist[i  ]/(6.0*h));

  return value;
}