//This is a simple routinue to find the eigenenergies of a square well by
//	Newton's Method
//	NOTE:  The initial guess is important for Newton Method
//		The initial guess should be fairly close to the root and
//		Newton's Method will be very efficient in converging to it
//		The initial guess can be effectively done by graphical method
#include <stdio.h>
#include <math.h>

#define EPS 1e-12
#define MaxVal 1e7
#define MAXI 200
#define V0 10.0

double findeven(double engy,int *counter);
double findodd(double engy,int *counter);
double evenf(double x,double y);
double oddf(double x,double y);
double devenf(double x,double y);
double doddf(double x,double y);

main()
	{
	double e;
	int choice;
	int count;


	printf("Choose one: [0] Even state [1] Odd State\n");
	scanf("%d",&choice);
	printf("Please input initial guess for eigen-energy\n");
	scanf("%lg", &e);

	if(choice)
		e=findodd(e,&count);
	else
		e=findeven(e,&count);

	printf("After %d Newton's Step, the estimated eigenenergy is %20.15lg\n",count,e);

	return 0;
	}

double findeven(double newengy,int *counter)
{
	double engy;
	double alpha, beta;

	*counter=0;
	do
		{
		engy=newengy;
		alpha=sqrt(engy);
		beta=sqrt(V0-engy);

		newengy=engy-evenf(alpha,beta)/devenf(alpha,beta);
		if(fabs(newengy)>MaxVal)
			{
			printf("Newton Blowup... Try New Guess\n");
			*counter=MAXI;
			}
		*counter+=1;
		}
	while(fabs(newengy-engy)>EPS && *counter<MAXI);

	return newengy;

}

double findodd(double newengy,int *counter)
{
	double engy;
	double alpha, beta;

	*counter=0;
	do
		{
		engy=newengy;
		alpha=sqrt(engy);
		beta=sqrt(V0-engy);

		newengy=engy-oddf(alpha,beta)/doddf(alpha,beta);
		if(fabs(newengy)>MaxVal)
			{
			printf("Newton Blowup... Try New Guess\n");
			*counter=MAXI;
			}
		*counter+=1;
		}
	while(fabs(newengy-engy)>EPS&& *counter<MAXI);

	return newengy;

}

double evenf(double aa,double bb)
{
	return aa*tan(aa)-bb;
}

double oddf(double aa,double bb)
{
	return -aa/tan(aa)-bb;
}

double devenf(double aa,double bb)
{
	double coseng;

	coseng=cos(aa);
	return ((1.0/coseng+sin(aa)/aa)/coseng+1.0/bb)/2.0;
}

double doddf(double aa,double bb)
{
	double sineng;

	sineng=sin(aa);
	return ((1.0/sineng-cos(aa)/aa)/sineng+1.0/bb)/2.0;
}