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

using namespace std;
const double pi = acos(-1);

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

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))
		{
			//cerr << "x=" << x0 << "\tf=" << d(x0) <<"\td=" << f(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;
			
		}
		if (!maxTries) cerr << "maxTries reached in Newton Approximation" << endl;
		return xn;
	};
protected:
	Function& d;
};


//
//This implementation uses Miller's Algorithm to compute Bessel Function using it's
// the backward recursive algorithm Jn+1(x) + Jn-1(x) = (2n/x)Jn(x).  This recursion 
// will converge up to a constant coefficient given by the renormalization formula
// 1 = J0(x) + 2*J2(x) + 2*J4(x) ...
//
class BesselFunc : public Function
{
public:
	BesselFunc(int n, int maxN, double jInit=1) 
		: _n(n), _maxN(maxN), _jInit(jInit), j(0), jprime(0), d(0), sum(0)
	{
		//cout << "BesselFunc Constructor" << endl;
	};
	virtual ~BesselFunc()
	{
		if (j) 
		{
			free(j);
			free(jprime);
		}
		if (d)	delete d;
	};
	BesselFunc(const BesselFunc& bf)
	{
		_n = bf._n;
		_maxN = bf._maxN;
		_jInit = bf._jInit;
		sum = bf.sum;
		j = bf.j;
		jprime = bf.jprime;
		last_x = bf.last_x;
		d = bf.d;
	}
	virtual double operator()(double x)
	{
		if (j && (last_x==x)) return j[_n]/sum;

		last_x = x;
		sum = 0;
		if (!j)
		{
			j = (double*)malloc(sizeof(double)*(_maxN+2));
			jprime = (double*)malloc(sizeof(double)*(_maxN+2));
			memset(j, 0, sizeof(double)*(_maxN+2));
		};
		j[_maxN]=_jInit;
		for (int i=_maxN; i > 0;i--)
		{
			j[i-1] = (2*i/x)*j[i] -j[i+1];
			jprime[i] = (j[i-1]-j[i+1])/2;
			if (!(i%2)) sum+= j[i]; 
			if (sum > 1.0e+100) //for large _maxN parameters will run over maxDouble
			{
				sum /=1.0e+100;
				for (int k=_maxN;k>=i;k--)	
				{
					jprime[i] /=1.0e+100;
					j[k] /= 1.0e+100;
				}
				j[i-1] /=1.0e+100;
			};
			//cout << "j[" << i << "]=" << j[i] << "\tsum=" << sum << endl;
		};
		sum = 2*sum + j[0];	//normalize using equality 1 = j0(x) + 2*j2(x) + 2*j4(x) + 2*j6(x) ...
		jprime[0] = -j[1];
		//cout << "BesselFunc(" << _n << ", " << x << " ) returned : " << (j[_n]/sum) << endl;
		return j[_n]/sum;
	};

	friend class BesselDeriv;  //recursive derivative structure (should not be applied in this case)

	virtual Function& getDerivative();
	
protected:
	int _n;
	int _maxN;
	double _jInit;
	double sum;
	double *j;
	double *jprime;
	double last_x;
	Function * d;
};

//class used to demonstrate structure of recursive getDerivative function
//was not necessary for this example
class BesselDeriv  : public BesselFunc
	{
	public:
		BesselDeriv(BesselFunc* parent) :_parent(parent) ,BesselFunc(*parent) 
		{
			d = 0;
			if (parent->j)
			{
				_n=parent->_n;
				sum = parent->sum;
				j = parent->jprime;
				_maxN--;
				//calc_jprime();	//pro-active recursion
			}
		};
		virtual double operator()(double x)
		{
			if (_parent->j && (_parent->last_x==x)) 
			{
				j = _parent->jprime;
				sum = _parent->sum;
			}
			else
			{
				last_x=x;
				(*_parent)(x);
				_n = _parent->_n;
				sum = _parent->sum;
				j = _parent->jprime;
				//calc_jprime(); //recursion uneeded (only first derivative needed)
			}
			return j[_n]/sum;
		}
	protected: 
		void calc_jprime()
		{
			jprime[0] = -j[1];
			for (int i = 1;i<_maxN;i++)
				jprime[i] = (j[i+1]-j[i-1])/2;
		}
		BesselFunc* _parent;
	};

Function& BesselFunc::getDerivative()
{
		if (!d)
			d = new BesselDeriv(this);
		return *d;
};


int main(int argc, char* argv[])
{
	double maxN = 500;
	cout.precision(15);
	cout.width(15);

	cout << "# Welcome to the waveguide mode calculation using the MillerApproximation\n";
	cout << "# of the Bessel Recursion:" << endl << "#\n";
	cout << "# (2n/x)*J(n;x)-J(n+1;x) = J(n-1;x)\n" << endl;
	cout << "# The waveguide modes are defined by T(n;m;) where n is the order of the\n";
	cout << "# Bessel Function and the mth zero gives the solution to the longest\n";
	cout << "# wavelength (a cut-off wavelength lc) where J(2pi/lc) =0 where pi=" << pi << endl;
	cout << "# and the recursive relation was begun at " << maxN;
	cout << endl;

	BesselFunc millers0(0,maxN);
	BesselFunc millers1(1,maxN);
	BesselFunc millers2(2,maxN);
	NewtonApproximation m0(2,4, millers0, millers0.getDerivative()), 
		m1(3,4, millers1, millers1.getDerivative()),
		m2(4,6, millers2, millers2.getDerivative());

	cout << "# For 0th order mode, lc = " << (2*pi/(double)m0.calculate()) <<endl;
	cout << "# For 1st order mode, lc = " << (2*pi/(double)m1.calculate()) << endl;
	cout << "# For 2nd order mode, lc = " << (2*pi/(double)m2.calculate()) << endl;
	
	cout << "\n# The graphs of the Bessel functions and their derivatives are given by" << endl;
	double x;
	for (x=0.1;x<10;x+=0.1)
	{
		cout << x << "\t" << millers0(x) << "\t" << millers1(x) << "\t" << millers2(x) 
			<< "\t" << millers0.getDerivative()(x) << "\t" << millers1.getDerivative()(x) << "\t" << millers2.getDerivative()(x)<< endl;
	};
	
	for (int n=0;n<3;n++)
	{
		int m = 1;
		double b;
		BesselFunc millers(n, 500, 1);
		NewtonApproximation zeroOfMiller(1,2*sqrt(n)+6, millers, millers.getDerivative());
		b = zeroOfMiller.calculate();
		cout << "#Bessel = " << b << "\tT(" << n << ";" << m << ";) = " << 2*pi/b << endl;
	};
	return 0;
}