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

double nevilles(double x, double *data, double * f, int lower, int upper) 
{
	int a=lower-1;
	int b=upper-1;
	//cout << a << "," << b << endl;
	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);
}
int main (int argc,char* argv[])
{	
	const int num_of_values = 9;
	double data[num_of_values] = {0,25,50,75, 100,125,150,175,200};
	double  values[num_of_values] = {10.6, 16.0, 45.0, 83.5, 52.8, 19.9, 10.8, 8.25, 4.7};
	int inc =3;
	int i, j, k,m;
	double x, y, average_x=0, average_y=0;
	//double x ;
    
	//while (1)
	//{	
		cout << "# Welcome to the Nevilles Interpolation Alogrithm program" << endl;
		//cout << "# These are the data points" << endl;
		//for (int i=0;i<num_of_values;i++) 
		//	cout << "# " << data[i] << "\t" << values[i] << endl;
		//cout << "# Please enter an x value to interpolate : ";
		//cin >> x;
		//cout << "# The nevilles algorithm interpolated value is " << nevilles(x, data, values, 1, num_of_values) << endl;
	//}

	for (i = 0; i <= 200; i+=5)
		cout << i << "\t" << nevilles(i, data, values, 1, num_of_values) << endl;
	
	cout << endl;

	for (i=0, j = inc, k=0,m=0 ;i <= 200; i+=5, k++,m++)
	{
		cout << i << "\t" << nevilles(i, data, values, j -inc+1, j) << endl;
		if (i >= data[j-1] )  j+=inc-1;
	}

	for (i=0, j = inc, k=0;i <= 200; i+=5)
	{
		double value = nevilles(i, data, values, j -inc+1, j);
		for (int i2=i+5, j2 = j; ((i2 <= 200) && (j2<=num_of_values)); i2+=5)
		{
			double  inner_value = nevilles(i2, data, values, j2 -inc+1, j2);
		
			x = (63890*(1/value - 1/inner_value) - (i)*(i) +i2*i2)/(2*(i2-i));
			if ((63890/value - (x-i)*(x-i)) > 0) 
			{	
				y = 2*sqrt(63890/value - (x-i)*(x-i));
				average_x +=  x;
				average_y += y;
				k++;
				//cout << x  << "\t" << y << endl;
			}
			if (i2 >=data[j2-1]) j2+=inc-1;
		}
		if (i >= data[j-1] )  j+=inc-1;
	}

	cout << "# Average Energy : " << average_x/(k) << endl;
	cout << "# Average cross section: " << average_y/(k) << endl;

	return 0;
}