#include <functional>  //for degenerate bisection problem
#include <iostream>
#include <math.h>

using namespace std;

double rvalues[11] = {0,1,2,3,4,5,6,7,8,9,10};

double bessel0data[11] = {1,0.7651976866,0.2238907791,-0.2600519549, -0.3971498099, -0.1775967713,0.1506452573,0.3000792705,0.1716508071,-0.0903336112,-0.2459357645};
double bessel1data[11] = {0,0.4400505857,0.5767248078,0.3390589585,-0.0660433280,-0.3275791376,-0.2766838581,-0.0046828235,0.2346363469,0.2453117866,0.0434727462};
double bessel2data[11] = {0,0.1149034849,0.3528340286,0.4860912606,0.3641281459,0.0465651163,-0.2428732100,-0.3014172201,-0.1129917204,0.1448473415,0.2546303137};


//recursive function calculatiing nevilles approximation
//
//	input: point to approximate, x values, y values, and the bounds
//  output: approximation for the value of 
//                 f (x) 
//
//  calulation is based on traveling from the value of the of the root
//  to the nodes and then performing the equation up the tree,  
// recusively representing the (upper-lower)th degree polynomial
//
double nevilles(double x, double *data, double * f, int lower, int upper) 
{
	int a=lower-1;
	int b=upper-1;
	if (a==b) return f[a];
	return (x-data[b])/(data[a]-data[b]) *nevilles(x, data, f, lower, upper-1) +
		(x-data[a])/(data[b]-data[a]) *nevilles(x, data, f, lower+1, upper);
}

// Implementation of a generic soultion to choosing bounds
// has not been verified
// 
// input: point given to interpolate,  increment of the interpolation (2 = quadratic, 3=cubic)
//
// output: sets the value of the lower and upper bound within the array 
//
// calculation: chooses lowest bracket for all numbers less than half the increment
//                        chooses highest bracket for all numbers in the upper half increment
//                        uses the bracket closes to the center of the graph for other points
//
void setBounds(double r, int inc, int& lower, int& upper)
{
	const int max = 11;
    int i = inc;
	
	if (r < rvalues[inc/2+1])
	{
		lower = 1;
		upper = 1+inc;
	}
	else
	{
		if (r > rvalues[max-inc/2-1])
		{
			lower = max-inc;
			upper = max;
		}
		else
		{
			while (r > rvalues[i]) i++;
			if (i>max/2) 
			{
				lower = i-inc+inc/2+1;
				upper = i+inc-inc/2;
			}
			else
			{
				lower = i-inc+inc/2+1;
				upper = i+inc-inc/2;
			}
		}
	}
}

//function that calculates the airy equation using nevilles approximation of the bessel function
//
double airy(double r)
{
	int lower, upper;
	
	setBounds(r,3,lower,upper);   //uses cubic approximation
	
	//cout << "# lower = " << lower << "\tupper = " << upper << endl;
    //cout << nevilles(r, rvalues, bessel1data, lower, upper) << endl;;
	return pow((2*nevilles(r, rvalues, bessel1data, lower, upper)/r), 2); 
}

//function to calculate the bessel using a simplified cubic Hermite interpolation
//
double bessel(double r)
{
	double val=0;
	int lower, upper;

	/*if (r<=1) 
	{
		lower = 0;
		upper = 1;
	}
	else
	{	
		lower = r-0.0000000001;
		upper = r-0.0000000001+1;
	}*/
	setBounds(r, 1, lower,upper);
	lower-=1;
	upper-=1;
	return  ((1+2*(r-rvalues[lower]))*pow(r-rvalues[upper],2) * bessel1data[lower] +
		(1-2*(r-rvalues[upper]))*pow(r-rvalues[lower],2) * bessel1data[upper] +
		(r-rvalues[lower])*pow(r-rvalues[upper],2) * (bessel0data[lower] - bessel2data[lower])/2  +
		(r-rvalues[upper])*pow(r-rvalues[lower],2) * (bessel0data[upper] - bessel2data[upper])/2);		
}

//value of the airy function using bessel hermite interpolation
double value(double i)
{
	if (i)
		return pow(2*bessel(i)/i,2);
	else
		return 1;
}

//derivative of the airy function using bessel data tables
double deriv(int i)
{
	if (i)
		return 8*bessel1data[i]*(i*(bessel0data[i]-bessel2data[i])/2 - bessel1data[i])/(i*i*i);
	else
		return 0;
}


//calculates value of airy function at a given point using the herimte interpolation
double hermite(double r)
{
/*	double val=0;
	double airydata[11];
	double j1;
	int i;

    /*
	//attempt to use nevilles approximation 
	//
	for (i = 1; i < 11 ; i++)
		airydata[i] = pow(2*bessel1data[i]/i,2);
	j1 = nevilles(r,rvalues, airydata, 3, 11);
	for ( i=2;i<11;i++)
	{
		val += (1-2*(r-rvalues[i]))*pow(j1,2)* value(i) + (r-rvalues[i])*pow(j1,2)*deriv(i);
	}
	return val;
	
	//lower-=1;
	//upper-=1;	

	if (r<=2) 
	{
		lower = 1;
		upper = 2;
	}
	else
	{	
		lower = r-0.0000000001;
		upper = r-0.0000000001+1;
	}
	*/

	int lower, upper;
	setBounds(r, 1, lower,upper);
	lower-=1;
	upper-=1;
	//cout << "# lower = " << lower << "\tupper = " << upper << endl;
	
	return  ((1+2*(r-rvalues[lower]))*pow(r-rvalues[upper],2) * value(lower) +
		(1-2*(r-rvalues[upper]))*pow(r-rvalues[lower],2) * value(upper)  +
		(r-rvalues[lower])*pow(r-rvalues[upper],2) * deriv(lower)  +
		(r-rvalues[upper])*pow(r-rvalues[lower],2) * deriv(upper) );
		
}

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

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;
};

//Bisection Approximation for the degenerate case
//
// bisects twice and chooses 3/4 of the interval around the least value
//
// has not been proven
//
class BisectionApproximation : public Approximation
{	
public:
	BisectionApproximation(double lower, double upper, Function& f)
		: Approximation(lower,upper,f)
	{};
	double calculate(double precision=0.0, int maxTries=2000)
	{
		double p=precision+1;
		while ((p > precision) && (--maxTries))
		{
			if ((f((a+(b-a)/4)))< (f(a+3*(b-a)/4)))
				b = a+3*(b-a)/4;
			else 
				a = a+(b-a)/4;
			p = fabs(a-b);
		}
		return (a+b)/2;
	};
};

int main (int argc, char* argv[])
{
	double r;  

	cout << "# Welcome to the Sir George Airy diffraction Pattern Bessel function interpolation problem" << endl;
	cout << "#  Data from the bessel function table given below will be used to interpolate the Airy function, I = I0 [2*J1(r)/r]^2" << endl;
	cout << "#          r\t        J0(r)\t         J1(r)\t         J2(r)" << endl;
	for (int i = 0; i < 11; i++)
		cout << "# " << i << "\t" << bessel0data[i] << "\t" << bessel1data[i] << "\t" << bessel2data[i] << endl;
	cout  << endl;

	cout << "# The relative intensity (I/I0) as a function of r across the airy diffraction pattern is given by the following data points" << endl;

	for (r = 0;r<=10; r+=0.1)
		cout << r << "\t" << airy(r) << endl;

	cout << endl <<  "# The relative intensity calculated using a cubic Hermite interpolation algorithm" << endl;
	for (r = 0;r<=10;r+=0.1)
		cout << r << "\t" << hermite(r) << endl;

	Hermite h;
	BisectionApproximation b(3,5,h);
	cout << endl << "# The cross section of the scattering has a zero at : " << b.calculate() <<endl;

	return 0;
};