/* interpolation/spline.c
 * 
 * Copyright (C) 2001, 2007 Brian Gough
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <string.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_spline.h>

gsl_spline *
gsl_spline_alloc (const gsl_interp_type * T, size_t size)
{
  gsl_spline * spline = (gsl_spline *) malloc (sizeof(gsl_spline));
  
  if (spline == NULL)
    {
      GSL_ERROR_NULL ("failed to allocate space for spline struct", 
                      GSL_ENOMEM);
    }
  
  spline->interp = gsl_interp_alloc (T, size);
  
  if (spline->interp == NULL)
    {
      free (spline);          
      GSL_ERROR_NULL ("failed to allocate space for interp", GSL_ENOMEM);
    };
    
  spline->x = (double *) malloc (size * sizeof(double));

  if (spline->x == NULL)
    {
      gsl_interp_free(spline->interp);
      free(spline);
      GSL_ERROR_NULL ("failed to allocate space for x", GSL_ENOMEM);
    }

  spline->y = (double *) malloc (size * sizeof(double));

  if (spline->y == NULL)
    {
      free(spline->x);
      gsl_interp_free(spline->interp);
      free(spline);
      GSL_ERROR_NULL ("failed to allocate space for y", GSL_ENOMEM);
    }
  
  spline->size = size;

  return spline;
}

int
gsl_spline_init (gsl_spline * spline, const double x_array[], const double y_array[], size_t size)
{
  if (size != spline->size)
    {
      GSL_ERROR ("data must match size of spline object", GSL_EINVAL);
    }
  
  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;
  }
}

const char *
gsl_spline_name(const gsl_spline * spline)
{
  return gsl_interp_name(spline->interp);
}

unsigned int
gsl_spline_min_size(const gsl_spline * spline)
{
  return gsl_interp_min_size(spline->interp);
}

void
gsl_spline_free (gsl_spline * spline)
{
  RETURN_IF_NULL (spline);
  gsl_interp_free (spline->interp);
  free (spline->x);
  free (spline->y);
  free (spline);
}

int
gsl_spline_eval_e (const gsl_spline * spline, 
                   double x,
                   gsl_interp_accel * a, double *y)
{
  return gsl_interp_eval_e (spline->interp, 
                            spline->x, spline->y,
                            x, a, 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);
}


int
gsl_spline_eval_deriv_e (const gsl_spline * spline,
                         double x,
                         gsl_interp_accel * a,
                         double *dydx)
{
  return gsl_interp_eval_deriv_e (spline->interp, 
                                  spline->x, spline->y,
                                  x, a, dydx);
}

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


int
gsl_spline_eval_deriv2_e (const gsl_spline * spline,
                          double x,
                          gsl_interp_accel * a,
                          double * d2)
{
  return gsl_interp_eval_deriv2_e (spline->interp, 
                                   spline->x, spline->y,
                                   x, a, d2);
}

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


int
gsl_spline_eval_integ_e (const gsl_spline * spline,
                         double a, double b,
                         gsl_interp_accel * acc,
                         double * result)
{
  return gsl_interp_eval_integ_e (spline->interp, 
                                  spline->x, spline->y,
                                  a, b, acc, result);
}


double
gsl_spline_eval_integ (const gsl_spline * spline,
                       double a, double b,
                       gsl_interp_accel * acc)
{
  return gsl_interp_eval_integ (spline->interp, 
                                spline->x, spline->y,
                                a, b, acc);
}