//****************************************************************************

//! \file Besneum.hpp
//! \date 07-11-2003

//! \class Bessel
//! \brief The class representing the Bessel functions

//! \class Neumann
//! \brief The class representing the Neumann functions

//*****************************************************************************

//     This file is part of the streuung.
//
//     streuung 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.
//
//     streuung 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 streuung; if not, write to the Free Software
//     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

//*****************************************************************************

#ifndef BESNEUM_HPP
#define BESNEUM_HPP

#include <cmath>
#include <iostream>
#include <gsl/gsl_sf_bessel.h>

//*****************************************************************************
class Bessel
//-----------------------------------------------------------------------------
//! This is a class designed just for one purpose: calculating the Bessel functions
//! /attention this works just for l<10 for l>=10 it will return -1
//*****************************************************************************
{

protected:

	double si; //!< this will contain the current s value (l) for the recursion
	double si1; //!< this will contain the s value for l-1 for the recursion
	double buf; //!< buffer for the s values during calculation

	virtual inline const double s0( const double x) const;//!< calculates the Bessel func for l=0
	virtual inline const double s1( const double x) const;//!< calculates the Bessel func for l=1

public:

	//! \brief constructor doing nothing for now
	Bessel() {};

	//! \brief destructor doing nothing for now
	~Bessel() {};

	//! \brief the main function doing the actual work
	//! \param l the l Value for the Calculation of the Bessel function
	//! \param x the x Value for the Calculation of the Bessel function
	//! \return the value of the Bessel function for this x and l Value
	inline double operator()( const int l, const double x);
};

//*****************************************************************************
const double Bessel::s0( const double x) const
//*****************************************************************************
{
	return sin( x)/x;
}

//*****************************************************************************
const double Bessel::s1( const double x) const
//*****************************************************************************
{
	return ( sin( x)/( x*x) - cos( x)/x);
}

//*****************************************************************************
double Bessel::operator()( const int l, const double x)
//*****************************************************************************
{
	if( l==0) return s0( x);
	if( l==1) return s1( x);
	if( l>=10) return -1.;

	si = s1( x);
	si1 = s0( x);
	for( int i = 1; i<l; ++i)
	{
		buf = si;
		si = si*( 2.*i /*+ 1.*/)/x - si1;
		si1 = buf;
	}

	return si;
}

//*****************************************************************************
class Neumann : public Bessel
//-----------------------------------------------------------------------------
//! This is a class designed just for one purpose: calculating the Neumann functions
//! derived from Bessel for my convenience!
//! /attention this works just for l<10 for l>=10 it will return -1
//*****************************************************************************
{

protected:

	virtual inline const double s0( const double x) const;//!< calculates the Neumann func for l=0
	virtual inline const double s1( const double x) const;//!< calculates the Neumann func for l=1

public:

	//! \brief constructor doing nothing for now
	Neumann() : Bessel() {};

	//! \brief destructor doing nothing for now
	~Neumann() {};

};

//*****************************************************************************
const double Neumann::s0( const double x) const
//*****************************************************************************
{
	return (-1.)*cos( x)/x;
}

//*****************************************************************************
const double Neumann::s1( const double x) const
//*****************************************************************************
{
	return ( ( -1.)*cos( x)/( x*x) - sin( x)/x);
}


#endif