/***************************************************************************
    gmf_bessel.cpp - wrapper for gsl bessel 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 */
    
    
#include "gmf_bessel.h" 


double BesselK(const double& Nu, const double& x)
{
	gsl_sf_result gsl_result;
	int gsl_status = 0;
	
	// avoid gsl errors causing an abort
	gsl_set_error_handler_off ();
	
	// The meat!
	gsl_status = gsl_sf_bessel_Knu_e (Nu, x, &gsl_result);
	if (gsl_status == GSL_SUCCESS)
	{
		return gsl_result.val;
	}
	else if (gsl_status == GSL_EUNDRFLW)  // check for obvious errors
	{
		//dout(WARN) << "Underflow detected in gsl_sf_bessel_Knu_e()" << endl;
		//dout(WARN) << "Maybe we should be avoiding this!"
      //           << endl;
		return 0.0;
	}
	else 
	{
		dout(ERR) << "Unexpected return from gsl_sf_bessel_Knu_e() : " 
		          << gsl_status << endl;
	}
	// FIXME - should really be throwing some exception here I guess
	dout(ERR) << "Reached unplanned code block in BesselK()" << endl;
	return 0.0;
}

double BesselK_gsl_form (double x, void * params) 
{
	double Nu = *(double *) params;
	return BesselK(Nu, x);
}