/* interpolation.cpp */

#include "interpolation.h"


/* default constructor */
interpolation::interpolation():
_x(), _y(),
_x_min(0.0), _x_max(0.0),
_acc(nullptr,gsl_interp_accel_free),
_spline(nullptr,gsl_spline_free)
{}

/* constructor */
interpolation::interpolation(std::vector<double> x, std::vector<double> y):
_x( x ), _y( y ),
_x_min( _x.front() ), _x_max( _x.back() ),
_acc(gsl_interp_accel_alloc(),gsl_interp_accel_free),
_spline(gsl_spline_alloc( gsl_interp_cspline, _x.size()),gsl_spline_free)
{
	gsl_spline_init( _spline.get(), _x.data(), _y.data(), _x.size() );
}

/* set */
auto interpolation::set(std::vector<double> x, std::vector<double> y) -> void
{
	_x = x;
	_y = y;
	_x_min = _x.front();
	_x_max = _x.back();

	_acc = std::unique_ptr<gsl_interp_accel,void(*)(gsl_interp_accel*)>
		(gsl_interp_accel_alloc(),gsl_interp_accel_free);

	_spline = std::unique_ptr<gsl_spline,void(*)(gsl_spline*)>
		(gsl_spline_alloc( gsl_interp_cspline, _x.size()),gsl_spline_free);

	gsl_spline_init( _spline.get(), _x.data(), _y.data(), _x.size() );
}

/* operator */
auto interpolation::operator()(double x) -> double
{
	if( x < _x.front() )
		return _y.front();

	if( x > _x.back() )
		return _y.back();

	return gsl_spline_eval( _spline.get(), x , _acc.get() );
}