///////////////////////////////////////////////////////////////////
//  What follows is that part of the GSL library which is used in the
//  core libraries in ODIN. It is used a replacement if libgsl is not installed.

#define  GSL_SUCCESS  0
#define  GSL_FAILURE  -1
#define  GSL_EINVAL   4

/* evaluation accelerator */
typedef struct {
  size_t  cache;        /* cache of index   */
  size_t  miss_count;   /* keep statistics  */
  size_t  hit_count;
}
gsl_interp_accel;

////////////////////////////////////////////////


/* interpolation object type */
typedef struct {
  const char * name;
  unsigned int min_size;
  void *  (*alloc) (size_t size);
  int     (*init)    (void *, const double xa[], const double ya[], size_t size);
  int     (*eval)    (const void *, const double xa[], const double ya[], size_t size, double x, gsl_interp_accel *, double * y);
  int     (*eval_deriv)  (const void *, const double xa[], const double ya[], size_t size, double x, gsl_interp_accel *, double * y_p);
  int     (*eval_deriv2) (const void *, const double xa[], const double ya[], size_t size, double x, gsl_interp_accel *, double * y_pp);
  int     (*eval_integ)  (const void *, const double xa[], const double ya[], size_t size, gsl_interp_accel *, double a, double b, double * result);
  void    (*free)         (void *);

} gsl_interp_type;

////////////////////////////////////////////////

/* general interpolation object */
typedef struct {
  const gsl_interp_type * type;
  double  xmin;
  double  xmax;
  size_t  size;
  void * state;
} gsl_interp;

////////////////////////////////////////////////

/* general interpolation object */
typedef struct {
  gsl_interp * interp;
  double  * x;
  double  * y;
  size_t  size;
} gsl_spline;

////////////////////////////////////////////////


gsl_interp_accel* gsl_interp_accel_alloc (void) {
  gsl_interp_accel *a = new gsl_interp_accel;
  a->cache = 0;
  a->hit_count = 0;
  a->miss_count = 0;
  return a;
}

gsl_interp* gsl_interp_alloc (const gsl_interp_type * T, size_t size) {
  if (size < T->min_size) {
   STD_cerr << "gsl_interp_alloc: insufficient number of points for interpolation type" << STD_endl;
    return 0;
  }
  gsl_interp * interp=new gsl_interp;
  interp->type = T;
  interp->size = size;

  if (interp->type->alloc == NULL) {
    interp->state = NULL;
    return interp;
  }

  interp->state = interp->type->alloc(size);
  if (interp->state == NULL) {
    delete interp;
    STD_cerr << "gsl_interp_alloc: failed to allocate space for interp state" << STD_endl;
    return 0;
  }

  return interp;
}

gsl_spline* gsl_spline_alloc (const gsl_interp_type * T, size_t size) {
  gsl_spline * spline = new gsl_spline;
  spline->interp = gsl_interp_alloc (T, size);
  if (spline->interp == NULL) {
    delete spline;
    STD_cerr << "gsl_interp_alloc: failed to allocate space for interp" << STD_endl;
    return 0;
  }
  spline->x = new double[size];
  spline->y = new double[size];
  spline->size = size;
  return spline;
}


int gsl_interp_init (gsl_interp * interp, const double x_array[], const double y_array[], size_t size) {
  if (size != interp->size) {
    STD_cerr << "gsl_interp_init: data must match size of interpolation object" << STD_endl;
  }

  interp->xmin = x_array[0];
  interp->xmax = x_array[size - 1];

  int status = interp->type->init(interp->state, x_array, y_array, size);
  return status;
}

int gsl_spline_init (gsl_spline * spline, const double x_array[], const double y_array[], size_t size) {
  if (size != spline->size) {
    STD_cerr << "gsl_spline_init: data must match size of spline object" << STD_endl;
  }

  memcpy (spline->x, x_array, size * sizeof(double));
  memcpy (spline->y, y_array, size * sizeof(double));

  int status = gsl_interp_init (spline->interp, x_array, y_array, size);
  return status;
}

void gsl_interp_accel_free(gsl_interp_accel * a) {delete a;}

void gsl_interp_free (gsl_interp * interp) {
  if (interp->type->free) interp->type->free (interp->state);
  delete interp;
}

void gsl_spline_free (gsl_spline * spline) {
  gsl_interp_free (spline->interp);
  delete[] spline->x;
  delete[] spline->y;
  delete spline;
}

///////////////////////////////////////////////////////////////////////

double gsl_interp_eval (const gsl_interp * interp, const double xa[], const double ya[], double x, gsl_interp_accel * a) {
  double y;
  int status = interp->type->eval (interp->state, xa, ya, interp->size, x, a, &y);
//  DISCARD_STATUS(status);
  return y;
}

double gsl_spline_eval (const gsl_spline * spline, double x, gsl_interp_accel * a) {
  return gsl_interp_eval (spline->interp,
                          spline->x, spline->y,
                          x, a);
}

///////////////////////////////////////////////////////////////////////

size_t gsl_interp_bsearch (
  const double x_array[], double x,
  size_t index_lo,
  size_t index_hi
  )
{
  size_t ilo = index_lo;
  size_t ihi = index_hi;
  while (ihi > ilo + 1)
    {
      size_t i = (ihi + ilo) / 2;
      if (x_array[i] > x)
        ihi = i;
      else
        ilo = i;
    }

  return ilo;
}

///////////////////////////////////////////////////////////////////////

size_t
gsl_interp_accel_find (gsl_interp_accel * a, const double xa[], size_t len, double x)
{
  size_t x_index = a->cache;

  if (x < xa[x_index])
    {
      a->miss_count++;
      a->cache = gsl_interp_bsearch (xa, x, 0, x_index);
    }
  else if (x > xa[x_index + 1])
    {
      a->miss_count++;
      a->cache = gsl_interp_bsearch (xa, x, x_index, len - 1);
    }
  else
    {
      a->hit_count++;
    }

  return a->cache;
}

///////////////////////////////////////////////////////////////////////

static inline double
integ_eval (double ai, double bi, double ci, double di, double xi, double a,
	    double b)
{
  const double t0 = b + a;
  const double t1 = a * a + a * b + b * b;
  const double t2 = a * a * a + a * a * b + b * b * a + b * b * b;
  const double bterm = 0.5 * bi * (t0 - 2.0 * xi);
  const double cterm = ci / 3.0 * (t1 - 3.0 * xi * (t0 - xi));
  const double dterm =
    di / 4.0 * (t2 - 2.0 * xi * (2.0 * t1 - xi * (3.0 * t0 - 2.0 * xi)));
  return (b - a) * (ai + bterm + cterm + dterm);
}

///////////////////////////////////////////////////////////////////////

typedef struct
{
  double * b;
  double * c;
  double * d;
  double * _m;
} akima_state_t;


/* common creation */
static void *
akima_alloc (size_t size)
{
  akima_state_t *state = new akima_state_t;
  state->b = new double[size];
  state->c = new double[size];
  state->d = new double[size];
  state->_m = new double[size + 4];
  return state;
}


/* common calculation */
static void
akima_calc (const double x_array[], double b[],  double c[],  double d[], size_t size, double m[])
{
  size_t i;

  for (i = 0; i < (size - 1); i++)
    {
      const double NE = fabs (m[i + 1] - m[i]) + fabs (m[i - 1] - m[i - 2]);
      if (NE == 0.0)
	{
	  b[i] = m[i];
	  c[i] = 0.0;
	  d[i] = 0.0;
	}
      else
	{
	  const double h_i = x_array[i + 1] - x_array[i];
	  const double NE_next = fabs (m[i + 2] - m[i + 1]) + fabs (m[i] - m[i - 1]);
	  const double alpha_i = fabs (m[i - 1] - m[i - 2]) / NE;
	  double alpha_ip1;
	  double tL_ip1;
	  if (NE_next == 0.0)
	    {
	      tL_ip1 = m[i];
	    }
	  else
	    {
	      alpha_ip1 = fabs (m[i] - m[i - 1]) / NE_next;
	      tL_ip1 = (1.0 - alpha_ip1) * m[i] + alpha_ip1 * m[i + 1];
	    }
	  b[i] = (1.0 - alpha_i) * m[i - 1] + alpha_i * m[i];
	  c[i] = (3.0 * m[i] - 2.0 * b[i] - tL_ip1) / h_i;
	  d[i] = (b[i] + tL_ip1 - 2.0 * m[i]) / (h_i * h_i);
	}
    }
}


static int
akima_init (void * vstate, const double x_array[], const double y_array[],
            size_t size)
{
  akima_state_t *state = (akima_state_t *) vstate;

  double * m = state->_m + 2; /* offset so we can address the -1,-2
                                 components */

  size_t i;
  for (i = 0; i <= size - 2; i++)
    {
      m[i] = (y_array[i + 1] - y_array[i]) / (x_array[i + 1] - x_array[i]);
    }

  /* non-periodic boundary conditions */
  m[-2] = 3.0 * m[0] - 2.0 * m[1];
  m[-1] = 2.0 * m[0] - m[1];
  m[size - 1] = 2.0 * m[size - 2] - m[size - 3];
  m[size] = 3.0 * m[size - 2] - 2.0 * m[size - 3];

  akima_calc (x_array, state->b, state->c, state->d, size, m);

  return GSL_SUCCESS;
}


static void
akima_free (void * vstate)
{
  akima_state_t *state = (akima_state_t *) vstate;

  delete[] (state->b);
  delete[] (state->c);
  delete[] (state->d);
  delete[] (state->_m);
  delete (state);
}


static
int
akima_eval (const void * vstate,
            const double x_array[], const double y_array[], size_t size,
            double x,
            gsl_interp_accel * a,
            double *y)
{
  const akima_state_t *state = (const akima_state_t *) vstate;

  size_t index;

  if (a != 0)
    {
      index = gsl_interp_accel_find (a, x_array, size, x);
    }
  else
    {
      index = gsl_interp_bsearch (x_array, x, 0, size - 1);
    }

  /* evaluate */
  {
    const double x_lo = x_array[index];
    const double delx = x - x_lo;
    const double b = state->b[index];
    const double c = state->c[index];
    const double d = state->d[index];
    *y = y_array[index] + delx * (b + delx * (c + d * delx));
    return GSL_SUCCESS;
  }
}


static int
akima_eval_deriv (const void * vstate,
                  const double x_array[], const double y_array[], size_t size,
                  double x,
                  gsl_interp_accel * a,
                  double *dydx)
{
  const akima_state_t *state = (const akima_state_t *) vstate;

  size_t index;

//  DISCARD_POINTER(y_array); /* prevent warning about unused parameter */

  if (a != 0)
    {
      index = gsl_interp_accel_find (a, x_array, size, x);
    }
  else
    {
      index = gsl_interp_bsearch (x_array, x, 0, size - 1);
    }

  /* evaluate */
  {
    double x_lo = x_array[index];
    double delx = x - x_lo;
    double b = state->b[index];
    double c = state->c[index];
    double d = state->d[index];
    *dydx = b + delx * (2.0 * c + 3.0 * d * delx);
    return GSL_SUCCESS;
  }
}


static
int
akima_eval_deriv2 (const void * vstate,
                   const double x_array[], const double y_array[], size_t size,
                   double x,
                   gsl_interp_accel * a,
                   double *y_pp)
{
  const akima_state_t *state = (const akima_state_t *) vstate;

  size_t index;

//  DISCARD_POINTER(y_array); /* prevent warning about unused parameter */

  if (a != 0)
    {
      index = gsl_interp_accel_find (a, x_array, size, x);
    }
  else
    {
      index = gsl_interp_bsearch (x_array, x, 0, size - 1);
    }

  /* evaluate */
  {
    const double x_lo = x_array[index];
    const double delx = x - x_lo;
    const double c = state->c[index];
    const double d = state->d[index];
    *y_pp = 2.0 * c + 6.0 * d * delx;
    return GSL_SUCCESS;
  }
}


static
int
akima_eval_integ (const void * vstate,
                  const double x_array[], const double y_array[], size_t size,
                  gsl_interp_accel * acc,
                  double a, double b,
                  double * result)
{
  const akima_state_t *state = (const akima_state_t *) vstate;

  size_t i, index_a, index_b;

  if (acc != 0)
    {
      index_a = gsl_interp_accel_find (acc, x_array, size, a);
      index_b = gsl_interp_accel_find (acc, x_array, size, b);
    }
  else
    {
      index_a = gsl_interp_bsearch (x_array, a, 0, size - 1);
      index_b = gsl_interp_bsearch (x_array, b, 0, size - 1);
    }

  *result = 0.0;

  /* interior intervals */

  for(i=index_a; i<=index_b; i++) {
    const double x_hi = x_array[i + 1];
    const double x_lo = x_array[i];
    const double y_lo = y_array[i];
    const double dx = x_hi - x_lo;
    if(dx != 0.0) {

      if (i == index_a || i == index_b)
        {
          double x1 = (i == index_a) ? a : x_lo;
          double x2 = (i == index_b) ? b : x_hi;
          *result += integ_eval (y_lo, state->b[i], state->c[i], state->d[i],
                                 x_lo, x1, x2);
        }
      else
        {
          *result += dx * (y_lo
                           + dx*(0.5*state->b[i]
                                 + dx*(state->c[i]/3.0
                                       + 0.25*state->d[i]*dx)));
        }
    }
    else {
      *result = 0.0;
      return GSL_FAILURE;
    }
  }

  return GSL_SUCCESS;
}


static const gsl_interp_type akima_type =
{
  "akima",
  5,
  &akima_alloc,
  &akima_init,
  &akima_eval,
  &akima_eval_deriv,
  &akima_eval_deriv2,
  &akima_eval_integ,
  &akima_free
};

const gsl_interp_type * gsl_interp_akima = &akima_type;


///////////////////////////////////////////////////////////////////

static int
linear_init (void * vstate,
             const double x_array[],
             const double y_array[],
             size_t size)
{
  return GSL_SUCCESS;
}

static
int
linear_eval (const void * vstate,
             const double x_array[], const double y_array[], size_t size,
             double x,
             gsl_interp_accel * a,
             double *y)
{
  double x_lo, x_hi;
  double y_lo, y_hi;
  double dx;
  size_t index;
  
  if (a != 0)
    {
      index = gsl_interp_accel_find (a, x_array, size, x);
    }
  else
    {
      index = gsl_interp_bsearch (x_array, x, 0, size - 1);
    }
  
  /* evaluate */
  x_lo = x_array[index];
  x_hi = x_array[index + 1];
  y_lo = y_array[index];
  y_hi = y_array[index + 1];
  dx = x_hi - x_lo;
  if (dx > 0.0)
    {
      *y = y_lo + (x - x_lo) / dx * (y_hi - y_lo);
      return GSL_SUCCESS;
    }
  else
    {
      *y = 0.0;
      return GSL_EINVAL;
    }
}


static
int
linear_eval_deriv (const void * vstate,
                   const double x_array[], const double y_array[], size_t size,
                   double x,
                   gsl_interp_accel * a,
                   double *dydx)
{
  double x_lo, x_hi;
  double y_lo, y_hi;
  double dx;
  double dy;
  size_t index;
  
  if (a != 0)
    {
      index = gsl_interp_accel_find (a, x_array, size, x);
    }
  else
    {
      index = gsl_interp_bsearch (x_array, x, 0, size - 1);
    }
  
  /* evaluate */
  x_lo = x_array[index];
  x_hi = x_array[index + 1];
  y_lo = y_array[index];
  y_hi = y_array[index + 1];
  dx = x_hi - x_lo;
  dy = y_hi - y_lo;
  if (dx > 0.0)
    {
      *dydx = dy / dx;;
      return GSL_SUCCESS;
    }
  else
    {
      *dydx = 0.0;
      return GSL_EINVAL;
    }
}


static
int
linear_eval_deriv2 (const void * vstate,
                    const double x_array[], const double y_array[], size_t size,
                    double x,
                    gsl_interp_accel * a,
                    double *y_pp)
{
  *y_pp = 0.0;

  return GSL_SUCCESS;
}


static
int
linear_eval_integ (const void * vstate,
                   const double x_array[], const double y_array[], size_t size,
                   gsl_interp_accel * acc,
                   double a, double b,
                   double * result)
{
  size_t i, index_a, index_b;
  
  if (acc != 0)
    {
      index_a = gsl_interp_accel_find (acc, x_array, size, a);
      index_b = gsl_interp_accel_find (acc, x_array, size, b);
    }
  else
    {
      index_a = gsl_interp_bsearch (x_array, a, 0, size - 1);
      index_b = gsl_interp_bsearch (x_array, b, 0, size - 1);
    }
  
    /* endpoints span more than one interval */

  *result = 0.0;
  
  /* interior intervals */
  for(i=index_a; i<=index_b; i++) {
    const double x_hi = x_array[i + 1];
    const double x_lo = x_array[i];
    const double y_lo = y_array[i];
    const double y_hi = y_array[i + 1];
    const double dx = x_hi - x_lo;

    if(dx != 0.0) {
      if (i == index_a || i == index_b)
        {
          double x1 = (i == index_a) ? a : x_lo;
          double x2 = (i == index_b) ? b : x_hi;
          const double D = (y_hi-y_lo)/dx;
          *result += (x2-x1) * (y_lo + 0.5*D*((x2-x_lo)+(x1-x_lo)));
        }
      else
        {
          *result += 0.5 * dx * (y_lo + y_hi);
        }
    }
  }
    
  return GSL_SUCCESS;
}

static const gsl_interp_type linear_type = 
{
  "linear", 
  2,
  NULL, /* alloc, not applicable */
  &linear_init,
  &linear_eval,
  &linear_eval_deriv,
  &linear_eval_deriv2,
  &linear_eval_integ,
  NULL, /* free, not applicable */
};

const gsl_interp_type * gsl_interp_linear = &linear_type;



///////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////