#include <functional>
#include <iostream>
#include <math.h>

using namespace std;

class Function : public unary_function<double,double>
{
public:
	Function() {};
	virtual result_type operator()(argument_type x){return x;};
	virtual Function& getDerivative() { return *this;};
};


class ProtonNeutronFunc : public Function
{
public:
	ProtonNeutronFunc(double factor, double exponent_factor, 
		double angular, double energy) 
		:k(factor), a(exponent_factor), m(angular), E(energy), d(factor,exponent_factor,angular) 
	{};

	virtual result_type operator()(argument_type x)
	{
		return -1*k*exp(-1*x/a) + m/pow(x,2) - E;
	};

	virtual Function& getDerivative()
	{
		return d;
	};

protected:
	double k;
	double a;
	double m;
	double E;

	class ProtonNeutronDerivativeFunc : public Function
	{
	public:
		ProtonNeutronDerivativeFunc(double _k, double _a, double _m)
			:	k(_k),a(_a),m(_m), d(_k,_a,_m)
		{};
		result_type operator () (argument_type x) 
		{
			return (k/a)*exp(-1*x/a) - m*2/pow(x,3);
		};
		Function& getDerivative()
		{
			return d;
		};
	protected: 
		class PNddFunc: public Function
		{
		public:
			PNddFunc(double _k, double _a, double _m)
				:	k(_k),a(_a),m(_m)
			{};
			result_type operator () (argument_type x) 
			{
				return m*6/pow(x,4) - (k/pow(a,2))*exp(-1*x/a);
			};
		protected:
			double k;
			double a;
			double m;
		};
	protected:
		double k;
		double a;
		double m;
		PNddFunc d; 
	};
	
	ProtonNeutronDerivativeFunc d;
};

class Approximation
{
public:
	Approximation(double lowerBound, double upperBound, Function& function)
		:a(lowerBound), b(upperBound), f(function)
	{};

	virtual double calculate(double precision, int maxTries = 3000)
	{
		cout << "Default approximation not implemented" << endl;
		return -1;
	}

protected:
	double a;
	double b;
	Function& f;
};

class NewtonApproximation : public Approximation
{
public:
	NewtonApproximation(double lower, double upper, Function& f, Function& derivative)
		:Approximation(lower,upper,f), d(derivative)
	{};
	double calculate(double precision=0.0, int maxTries=3000)
	{
		double x0 = (a+b)/2;
		double xn;
		double p = b-a;
		while ((p > precision) && (--maxTries))
		{
			//cout << "x=" << x0 << "\tf=" << f(x0) <<"\td=" << d(x0) << endl;
			xn = x0 - f(x0)/d(x0);
			if (xn>b) xn = (x0+b)/2;
			if (xn<a) xn = (x0+a)/2;
			p = fabs(x0-xn);
			x0=xn;
			//cout << p << endl;
			
		}
		return xn;
	};
protected:
	Function& d;
};

int main(int argc, char* argv[])
{   
	double r;
	ProtonNeutronFunc func(70, 0.2, 0.01, 0);
	ProtonNeutronFunc func25(70, 0.2, 0.01, 25);
	ProtonNeutronFunc funcneg25(70, 0.2, 0.01, -25);

	cout << "Welcome to the Newton Approximation Program for the conservative" << endl;
	cout << "central force problem between a neutron and a proton where the " <<endl ;
	cout << "potential is given by -k*e^(-r/a)" << endl;
	cout << "For this example l^2/2m=0.01, k=70, a=0.2 " <<endl;
	cout << endl;


	NewtonApproximation newton25(0,0.04,func25, func25.getDerivative());
	cout << "For energy level 25 the zero found is at : " 
		<< newton25.calculate(0) <<endl;

	NewtonApproximation newtonneg25lower(0,0.04,funcneg25, func25.getDerivative());
	NewtonApproximation newtonneg25upper(0.05,1,funcneg25, func25.getDerivative());
	cout << "For energy level -25 the lower bound is found at : "
		<< newtonneg25lower.calculate(0) <<  endl
		<< " while the upper bound is found at :"
		<< newtonneg25upper.calculate(0) << endl;

	NewtonApproximation newtonmin(0,0.1, func.getDerivative(), func.getDerivative().getDerivative());
	cout << "The degenerate case in which the two bodies orbit around each other" << endl;
	cout << "is found when r = " << (r=newtonmin.calculate());
	cout << " and E = " << func(r) << endl;

	return 0;
}