/***************************************************************************
    aon_interpolation.cpp - wrappers for GNU Scientific Library 
                            interpolation functions 

                             -------------------
    begin				: 2004-02-17
    copyright			: (C) 2004 by Andrew O'Neill, HiRes Grad Student
    email				: oneill@phys.columbia.edu
 ***************************************************************************/
 /* 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 2 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */
    
using namespace std;
   
#include "aon_interpolation.h"

// Constructor
// Set up the size of the interpolation and create a spline object
Interpolation::Interpolation(const int NumPoints, 
                             const gsl_interp_type *InterpType)
: itsNumPoints (NumPoints)
{
	itsSpline = gsl_spline_alloc (InterpType, NumPoints);
	itsAccel = gsl_interp_accel_alloc ();
	itsInitializedFlag = 0;
	dout(MEM) << "Constructor Interpolation::Interpolation()" << endl;
}

// also, use cspline interpolation by default
Interpolation::Interpolation(const int NumPoints)
: itsNumPoints (NumPoints)
{
	itsSpline = gsl_spline_alloc (gsl_interp_cspline, NumPoints);
	itsAccel = gsl_interp_accel_alloc ();
	itsInitializedFlag = 0;
}	
	
// Destructor
// Free the spline object
Interpolation::~Interpolation()
{
	dout(MEM) << "Destructor Interpolation::~Interpolation()" << endl;
	if (itsSpline != NULL)
	{
		gsl_spline_free (itsSpline);
		itsSpline = NULL;
	}
	else 
	{
		dout(WARN) << "In Interpolation destructor I found a NULL pointer to "
		           << "spline." << endl;
	}
	if (itsAccel != NULL)
	{
		gsl_interp_accel_free (itsAccel);
		itsAccel = NULL;
	}
	else 
	{
		dout(WARN) << "In Interpolation destructor I found a NULL pointer to "
		           << "accel." << endl;
	}
}


/// Evaluate(x)
/// uses the interpolated values we calculated to obtain the value of the 
/// function at position x
double Interpolation::Evaluate (const double& Arg)
{	
	if (itsInitializedFlag ==0)
	{
		dout(ERR) << "Trying to evaluate interpolation "
					 <<"when it hasn't been initialised" << endl;
	}
	double answer = 0.0;
	int status = gsl_spline_eval_e(itsSpline, Arg, itsAccel, &answer);
	switch (status)
	{
		case GSL_SUCCESS:
			break;
		case GSL_EDOM:
			dout(ERR) << "Domain error tring to evaluate an interpolation at"
			          << Arg << endl;
			break;
		default:
			dout(ERR) << "Unexpected error in gsl_spine_eval_e() called from " 
			          << "Evaluate()" << endl;
			dout(ERR) << "Trying to evaluate argument " << Arg << endl;
			dout(ERR) << "Error encountered: " << status << endl;	
			break;
	}
	return answer;
}

/// EvaluateIntegral(a, b)
/// evaluates our interpolated function's numerical integral between a & b
double Interpolation::EvaluateIntegral (const double& a, const double& b)
{
	if (itsInitializedFlag ==0)
	{
		dout(ERR) << "Trying to integrate interp when it hasn't been initialised"
		          << endl;
	}
	double answer = 0.0;
	int status = gsl_spline_eval_integ_e(itsSpline, a, b, itsAccel, &answer);
	switch (status)
	{
		case GSL_SUCCESS:
			break;
		case GSL_EDOM:
			dout(ERR) << "Domain error tring to evaluate the integral between "
			          << a << " and " << b << endl;
			break;
		default:
			dout(ERR) << "Unexpected error in gsl_spine_eval_integ_e() called from " 
			          << "EvaluateIntegral()" << endl;
      	dout(ERR) << "Trying to evaluate between" << a << " and " << b << endl;
      	dout(ERR) << "Error encountered: " << status << endl;	
      	break;
	}
	return answer;
}