//This program will use the Neville's Algorithm to interpolate
//	an resonant energy spectrum
//There will be three outputs for this program
//	the output file name will be resonantE.dat
//	-interpolated result with all 9 pts
//	-interpolated result with every three points
//	-data with theoretical curve
#include <stdio.h>
#include <math.h>
#include "numcip.h"

double lagran(double e,double **dlist,double *cc,double *dd, int order);
double crosssection(double e);
void putindata(double **);

main()
{
	double **datapt;
	double *c,*d;
	double eng,sig9,sig3,theo;
	int n;

	FILE *outfile;

	//initialize program
		//vectors and matrix allocations
	datapt=dmatrix(1,9,1,2);//this will contain the data points
		//c[] and d[] will be updated at each 
	c=dvector(1,9);//downward difference
	d=dvector(1,9);//upward difference

		//put data in datapt
	putindata(datapt);
		//set up output file
	outfile=fopen("resonantE.dat","w");
	fprintf(outfile,"eng\tint9\tint3\ttheo\n");
	


	//main loop

	eng=0.0;
	for(n=1;n<40;n++)//this loop over the whole energy range
		{
		sig9=lagran(eng,datapt,c,d,9);
		sig3=lagran(eng,datapt,c,d,3);
		theo=crosssection(eng);

		fprintf(outfile,"%lg\t%lg\t%lg\t%lg\n",
			eng,sig9,sig3,theo);
		
		eng+=5.0;//increment eng by 5Mev
		}
	//output for final point at e=200
	fprintf(outfile,"%lg\t%lg\t%lg\t%lg\n",
		200.0,datapt[9][2],datapt[9][2],crosssection(200.0));

	return 0;
}

double lagran(double e,double **dlist,double *cc,double *dd, int order)
{
	double **currentpt;//will point to first data point in chunk
	int left;//index for first point in chunk
	int i,k;
	double y;
	double nec,ned,den,cminusd;

	//put e within data range
	if(order==3)
		{
		left=(int)(e/50);//this give the ith data chunck(starting with 0th) 
		left=2*left;//index needed to shifted to get to current data chunck
		}
	//if order = 9, it is the full range
	if(order==9)
		left=0;

	currentpt=dlist+left;//currentpt point to first in chunck
	//initialize c's and d's
	for(i=1;i<=order;i++)
		cc[i]=dd[i]=currentpt[i][2];

	//going down the tree column by column
	//	we start at the top left and move y down the top edge
	y=currentpt[1][2];
	for(k=1;k<order;k++)//k is the column number
		{
		//update c's and d's
		for(i=1;i<=order-k;i++)
			{
			ned=currentpt[i+k][1]-e;//top x factor for d
			nec=currentpt[i][1]-e;//top x factor for c
			den=currentpt[i][1]-currentpt[i+k][1];//bottom x factor for both
			cminusd=cc[i+1]-dd[i];//top factor with c-d
			
			//new c and d
			dd[i]=ned*cminusd/den;
			cc[i]=nec*cminusd/den;
			}
		//move y one column forward
		y+=cc[1];
		}
	return y;
}

#define sig0 6.389e4
#define engr 78.0
#define halfgamma 27.5

double crosssection(double e)
{
	return sig0/((e-engr)*(e-engr)+halfgamma*halfgamma);
}

void putindata(double **datalist)
//[][1]-energy; [][2]-cross section
{
	datalist[1][1]=0.0;
	datalist[1][2]=10.6;
	datalist[2][1]=25.0;
	datalist[2][2]=16.0;
	datalist[3][1]=50.0;
	datalist[3][2]=45.0;
	datalist[4][1]=75.0;
	datalist[4][2]=83.5;
	datalist[5][1]=100.0;
	datalist[5][2]=52.8;
	datalist[6][1]=125.0;
	datalist[6][2]=19.9;
	datalist[7][1]=150.0;
	datalist[7][2]=10.8;
	datalist[8][1]=175.0;
	datalist[8][2]=8.25;
	datalist[9][1]=200.0;
	datalist[9][2]=4.7;
}