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

//! \file Numerov.hpp
//! \date 04-11-2003

//! \class Numerov
//! \brief The class representing Numerov's Method of solving Differential Equations

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

//     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 NUMEROV_HPP
#define NUMEROV_HPP

#include <cmath>
#include <iostream>

const double eps = 5.9;//!< the epsilon from the L-J Potential
const double sig = 3.57; //!< the sigma from the L-J Potential
const double lam = 2.099; //!< the lambda for the calculation of f

//*****************************************************************************
class Numerov
//-----------------------------------------------------------------------------
//! This is the class encapsulating all Numerov related stuff
//*****************************************************************************
{
	double x0; //!< the start value for x
	double w0; //!< the numerov start value
	double wb0; //!< the numerov value before the start value (w before 0)
	double h; //!< the stepsize
	double r; //!< the function parameter
	double h212; //!< h*h/12
	double h2; //!< h*h
	int l; //!< the l for the calculation of f
	double E; //!< the E for the calculation of f
	double f; //!< the "Energy" value
	double buf; //!< buffer for the old w value

	//! \brief calculate the "Energy" for the Numerov formula
	//! \param t the parameter for the Energy
	//! \return the f value
	inline double ff( const double t) const;

	//! \brief calculates w from x
	//! \param x the x value to be transformed to w
	//! \return the w value
	inline const double w( const double x) const;

	//! \brief calculates x from w
	//! \param w the w value to be transformed to x
	//! \return the x value
	inline const double x( const double w) const;

	//! \brief calculates the L-J Pot
	//! \param r the distance for the L-J calculation
	inline const double V( const double r) const;

public:
	
	//! \brief the constructor
	//! \param inh the stepsize
	//! \param inl the l Value for f
	//! \param ine the given Energy
	//! \param startr the r start value
	//! \param x0 the start value
	//! \param xb0 the value one step before the  start value
	inline Numerov( const double inh, const int inl = 0, const double ine = 0, const double startr = 0, 
		const double x0 = 0, const double xb0 = 0);

	//! \brief the destructor
	~Numerov(){};

	//! \brief resets all calculation parameters
	//! \param inl the l Value for f
	//! \param ine the given Energy
	//! \param startr the r start value
	//! \param inx0 the start value
	//! \param inxb0 the value one step before the  start value
	inline void Init( const int inl = 0, const double ine = 0, const double startr = 0, 
		const double inx0 = 0, const double inxb0 = 0);

	//! \brief the main function doing the actual work
	//! \return the new x value
	inline double step();

	//! \brief doing N steps
	//! \param N the number of steps
	//! \return the x value of the last step
	inline double loopstep( const int N);

	//! \brief returns the current r value
	//! \return the r value
	inline const double getR() const { return r;}
};

//*****************************************************************************
Numerov::Numerov( const double inh, const int inl, const double inE, const double startr,
	const double inx0, const double inxb0)
: h( inh), h2( inh*inh), h212( inh*inh/12.), l( inl), E( inE), r( startr), x0( inx0)
//*****************************************************************************
{
	if( r > 0.)
	{
		f = ff( r - h);
		wb0 = w( inxb0);
		f = ff( r);
		w0 = w( x0);
	}
}

//*****************************************************************************
void Numerov::Init( const int inl, const double inE, const double startr, const double inx0, const double inxb0)
//*****************************************************************************
{
	E = inE;
	l = inl;
	r = startr;
	x0 = inx0;
	f = ff( r - h);
	wb0 = w( inxb0);
	f = ff( r);
	w0 = w( x0);
}

//*****************************************************************************
double Numerov::step()
//*****************************************************************************
{
	buf = w0;
	w0 = h2*x0*f+2*w0-wb0;
	wb0 = buf;
	r += h;
	f = ff(r);
	x0 = x(w0);
	return x0;
}

//*****************************************************************************
double Numerov::loopstep( const int N)
//*****************************************************************************
{
	for( int i = 0; i< N-1; ++i) step();

	return step();
}

//*****************************************************************************
double Numerov::ff( const double t) const
//*****************************************************************************
{
// 	if( ( l==4)&&( E==0.01))std::cout<<t<<" "<<l*( l + 1.)/( t*t)+V(t)/lam<<std::endl;
	return ( l*( l + 1.)/( t*t) + ( V( t) - E)/lam);
}

//*****************************************************************************
const double Numerov::w( const double x) const
//*****************************************************************************
{
	return x*( 1. - h212*f);
}

//*****************************************************************************
const double Numerov::x( const double w) const
//*****************************************************************************
{
	return w/( 1. - h212*f);
}

//*****************************************************************************
const double Numerov::V( const double r) const
//*****************************************************************************
{

	return eps*( pow( sig/r, 12.) - 2.*pow( sig/r, 6.));
}

#endif