//This is short program to find the perturbed energy levels by using the Romberg's Scheme
//	for a square well described in assignment #3, there are two even states and one odd state
//In this case, the first even and first odd eigenfunctions are hard wired into the program
//
//The Romberg's Scheme is to improve accuracy using extrapolation when h is halved
#include <stdio.h>
#include <math.h>

#define EPS 1e-10

#define ENGeven1 1.40721472477016
#define ENGeven2 9.99598073754667
#define ENGodd 5.37580591367022
#define Vp 0.01

double Tmat[2][10];//Tmat[][10] contains the nth row of the integrated values in Table 3-11
			  //Tmat[2][] alternate being the most recent and old values (nth and n-1th rows)

double alphaEven,alphaOdd;
double betaEven,betaOdd;
double normout;//normalization for part of wavefunction outside well

double rombint(double (*func)(double));
double romberg(double lastrow, double currentrow, int nstep);
double Trapz(double h,int Nbin,double (*fun)(double));
double powerOf2(int);
double eventop(double);
double evenbot(double);
double oddtop(double);
double oddbot(double);
double normouteval(double eng);

main()
{
	int choice;
	double perte;

	double (*topinteg)(double);
	double (*botinteg)(double);

	FILE *outfile;

	printf("Perturbation on which eigenstate:\n");
	printf("\t[0] 1st Odd\n\t[1] 1st Even\n\t[2] 2nd Even\n");
	scanf("%d",&choice);

	
	if(choice)
		{
		topinteg=&eventop;
		botinteg=&evenbot;
		if(choice==1)
			{
			alphaEven=sqrt(ENGeven1);
			betaEven=sqrt(10.0-ENGeven1);
			}
		if(choice==2)
			{
			alphaEven=sqrt(ENGeven2);
			betaEven=sqrt(10.0-ENGeven2);
			}
		normout=normouteval(betaEven);
		}
	else
		{
		topinteg=&oddtop;
		botinteg=&oddbot;
		alphaOdd=sqrt(ENGodd);
		betaOdd=sqrt(ENGodd);
		normout=normouteval(betaOdd);
		}
	

	perte=rombint(topinteg)/(rombint(botinteg)+normout);

	outfile=fopen("perteng.dat","a");
	fprintf(outfile,"Perturbation on which eigenstate:\n");
	fprintf(outfile,"\t[0] 1st Odd\n\t[1] 1st Even\n\t[2] 2nd Even\n");
	fprintf(outfile,"\tchoice: %d\n",choice);
	fprintf(outfile,"The perturbed energy is: %20.15lg\n",perte);

	return 0;
}

double normouteval(double beta)
{
	return exp(-2.0*beta)/beta;
}

double rombint(double (*func)(double))
{
	int i;
	double h;//stepsize
	int Numbin;
	int n;//iterative step
	int current,last;//keeping track of current row with respect to Tmat[2][]
	//Tmat[row][column]

	//the well is from -1 to 1 and total width is 2
	Numbin=4;//start with 2^2 bins
	h=2.0/Numbin;

	//Romberg iteration will start with h/2 row completed
	n=1;//matrix is 0 offset
	last=0;
	current=1;
	Tmat[last][n-1]=Trapz(h,Numbin,func);//Trapz integral with h
	Numbin*=2;//doubling bins
	h=2.0/Numbin;
	Tmat[current][n-1]=Trapz(h,Numbin,func);//Trapz integral with h/2
	Tmat[current][n]=romberg(Tmat[last][n-1],Tmat[current][n-1],n);
	
	while (fabs(Tmat[current][n]-Tmat[last][n-1])>EPS && n<8)
	{
		n++;
		Numbin*=2;
		h=2.0/Numbin;
		last=current;
		current=(current+1)%2;

		Tmat[current][0]=Trapz(h,Numbin,func);//first column in current row
		for(i=1;i<=n;i++)
			{
			Tmat[current][i]=romberg(Tmat[last][i-1],Tmat[current][i-1],i);
			}
	};

	return Tmat[current][n];
}

double powerOf2(int nn)
{
	double prod;
	int i;

	prod=4.0;
	for(i=1;i<nn;i++)
		prod*=4.0;

	return prod;
}

double romberg(double lastrow, double currentrow, int nstep)
{
	double coef;

	coef=powerOf2(nstep);
	return (coef*currentrow-lastrow)/(coef-1.0);
}

double Trapz(double stepp,int Nbin,double (*func)(double))
{
	double x;
	double intx;
	int i;

	x=-1.0;
	intx= func(x)/2.0;
	for(i=1;i<Nbin;i++)
	{
		x+=stepp;
		intx+= func(x);
	}
	x+=stepp;
	intx+= func(x)/2.0;

	return intx*stepp;
}

double eventop(double xx)
{
	double temp;

	temp=exp(-betaEven)*cos(alphaEven*xx)/cos(alphaEven);
	return temp*temp*Vp*exp(-xx*xx);
}

double evenbot(double xx)
{
	double temp;

	temp=exp(-betaEven)*cos(alphaEven*xx)/cos(alphaEven);
	return temp*temp;
}

double oddtop(double xx)
{
	double temp;

	temp=exp(-betaOdd)*sin(alphaOdd*xx)/sin(alphaOdd);
	return temp*temp*Vp*exp(-xx*xx);
}

double oddbot(double xx)
{
	double temp;

	temp=exp(-betaOdd)*sin(alphaOdd*xx)/sin(alphaOdd);
	return temp*temp;
}