//This is a simple routinue to find the MINIMUM in a central force problem
//	using the Newton's Method
//	Since the minium is at a degenerate root, instead of finding the zero of f(x)
//		find the zero of f'(x)
//	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 MAXI 200
#define k 1.0
#define a 0.2
#define lcoef 0.01

//this problem is basically the same as central.c except the derivative of the function
//	is used instead
double findzero(double engy,double left,double right,int *counter);
double newtonstep(double eng,double rr);
double potrate(double EE,double rr);
double dpotrate(double rr);
double potential(double rr);

main()
	{
	double leftr,rightr,r;
	double energy;
	int count;

	FILE *outfile;

	printf("Finding turning points in central force\n");
	printf("Please input bracket\n");
	printf("\tLeft:\n");
	scanf("%lg",&leftr);
	printf("\tRight:\n");
	scanf("%lg",&rightr);

	r=findzero(0.0,leftr,rightr,&count);

	//output to file
	outfile=fopen("minimum.dat","w");
	fprintf(outfile,"Minimum for central force problem with U(r)=k(e^-r/a)/r\n");
	fprintf(outfile,"\tBracket for root: [%lg,%lg]\n",leftr,rightr);
	fprintf(outfile,"After %d Newton's Step, the estimated minimum is at %20.15lg\n",count,r);
	fprintf(outfile,"\tTotal Energy at minimum is: %20.15lg\n",potential(r));
	fprintf(outfile,"\n");
	fclose(outfile);

	return 0;
	}

double findzero(double engy,double left,double right,int *counter)
{
	double rr,tempr;
	double err;
	
	//initial counter
	*counter=0;

	rr=(right+left)/2.0;//midpoint of bracket as first guess

	do
		{

		tempr=newtonstep(engy,rr);
		
		//to check if newton's step is within bracket
		if(tempr>=left&&tempr<=right)
			{
			//error estimte
			err=tempr-rr;

			//if inside bracket, then take the Newton's step
			rr=tempr;

			}
		else
			{
			//error estimte
			err=right-left;

			//if outside bracket, take a bisection step
			rr=(right+left)/2.0;
			
			}
		//maintaining bracket
		if(potential(engy,left)*potential(engy,rr)<0)
			right=rr;
		else
			left=rr;
		
		*counter+=1;
		}
	while(fabs(err/rr)>EPS && *counter<MAXI);


	return rr;
}

double newtonstep(double eng,double rr)
{
	return rr-potrate(eng,rr)/dpotrate(rr);
}

double potrate(double EE,double rr)
//dU'/dr
{
	return k*exp(-rr/a)*(1.0/a+1.0/rr)/rr-2.0*lcoef/(rr*rr*rr);
}

double dpotrate(double rr)
//d^2U'/dr^2
{
	return -k*exp(-rr/a)*(1.0/(a*a)+2.0/(a*rr)+2.0/(rr*rr))/rr+6.0*lcoef/(rr*rr*rr*rr);
}

double potential(double rr)
//U'(r)
{
	return -k*exp(-rr/a)/rr+lcoef/(rr*rr);
}