#include <stdio.h>
#include <stdlib.h>
#include <math.h>

//accuracy for Simpson's Rule
#define EPS 1.0e-7
#define MaxN 1000000

double pi;
double eng;
double rstart,rstop;
double **orbit;

double integrateSR(double rstar,double rend,double tol,int *count);
double Ekernal(double xx,double angle);

double *dvector(int nl,int nh);
double **dmatrix(int nrl,int nrh,int ncl,int nch);
void nrerror(char error_text[]);


main()
	{
	int i;
	int n,N;
	int choice;
	int Nloop;
	int count;
	double midpt;
	double r;
	double kstep;
	double **flip;

	FILE *outfile;

	pi=4.0*atan(1.0);

	//user io
	printf("Please choose a type of orbit:\n");
	printf("\t[0]Unbounded orbits\n");
	printf("\t[1]Bounded orbits\n");
	scanf("%d",&choice);
	if(!choice)
		{
		printf("Starting r [1.0]\n");
		rstart=1.0;
		printf("Stopping r [0.00953478]\n");
		rstop=0.00953478;
		eng=10.0;
		Nloop=1;//number of extensions to be made
		outfile=fopen("orb_ub.dat","w");
		}
	else
		{
		printf("Starting r [0.01221621]\n");
		rstart=0.01221621;
		printf("Stopping r [0.0576140]\n");
		rstop=0.0576140;
		eng=-10.0;
		Nloop=20;//number of extensions to be made
		outfile=fopen("orb_bd.dat","w");
		}

	printf("The number of subdivision for plotting? [200]\n");
	scanf("%d",&N);
	kstep=(rstop-rstart)/N;

	//get storage for basic orbit
	orbit=dmatrix(0,N,1,2);
	flip=dmatrix(0,N,1,2);

	fprintf(outfile,"rr\ttheta\txx\tyy\n");
 
	//generating the basic orbit
	if(!choice)
		{
		r=rstart;
		for(i=0;i<=N;i++)
			{
			orbit[i][1]=r;//storing radius
			orbit[i][2]=integrateSR(rstart,r,EPS,&count);//storing angle

			printf("Counter: %d\n",count);//check for number of steps in SR

			r+=kstep;
			}
		}
	else
		{
		//the routine is more efficient if we start at a midpoint
		//	and integrate out to the turning pts separately
		midpt=rstart+N/2*kstep;
		r=midpt;
		for(i=0;i<=N/2;i++)//toward rstart
			{
			orbit[N/2-i][1]=r;//storing radius
			orbit[N/2-i][2]=integrateSR(midpt,r,EPS,&count);//storing angle

			printf("Counter: %d\n",count);

			r-=kstep;
			}
		r=midpt+kstep;
		for(i=N/2+1;i<=N;i++)//toward rstop
			{
			orbit[i][1]=r;//storing radius
			orbit[i][2]=integrateSR(midpt,r,EPS,&count);//storing angle

			printf("Counter: %d\n",count);

			r+=kstep;
			}
		}

	//print the basic orbit to data file
	for(i=0;i<=N;i++)
		fprintf(outfile,"%lg\t%lg\t%lg\t%lg\n",orbit[i][1],orbit[i][2],orbit[i][1]*cos(orbit[i][2]),
			orbit[i][1]*sin(orbit[i][2]));


	//getting the other extensions and print out data to file
	for(n=1;n<=Nloop;n++)
		{
		printf("Extension %d\n",n);
		
		for(i=0;i<=N;i++)
		{
			flip[i][1]=orbit[N-i][1];
			flip[i][2]=2.0*orbit[N][2]-orbit[N-i][2];
		}

		for(i=0;i<=N;i++)
			{
			orbit[i][1]=flip[i][1];
			orbit[i][2]=flip[i][2];
			fprintf(outfile,"%lg\t%lg\t%lg\t%lg\n",orbit[i][1],orbit[i][2],orbit[i][1]*cos(orbit[i][2]),
				orbit[i][1]*sin(orbit[i][2]));
			}
		}

	return 0;
	}
		
//
//
//Integration Routinues
//
//

// Simpson's Rule

double integrateSR(double rstar,double rend,double tol,int *counter)
	{
	int i;
	int Ny;
	double hy;
	double Sye,Syo,Syd;
	double Iy,Iyold;

	*counter=0;
	Ny=6;
	hy=(rend-rstar)/Ny;
	Syd=Ekernal(rstar,eng)+Ekernal(rend,eng);
	Syo=Sye=0.0;
	for(i=1;i<Ny;i+=2)
		Syo+=Ekernal(rstar+i*hy,eng);
	for(i=2;i<Ny;i+=2)
		Sye+=Ekernal(rstar+i*hy,eng);
	Iy=hy/3.0*(Syd+2.0*Sye+4.0*Syo);

	do
		{
			Iyold=Iy;
			Ny*=2;
			hy/=2.0;
			Sye+=Syo;
			Syo=0.0;
			for(i=1;i<Ny;i+=2)
				Syo+=Ekernal(rstar+i*hy,eng);
			Iy=hy/3.0*(Syd+2.0*Sye+4.0*Syo);
			*counter+=1;

		}while(fabs((Iy-Iyold)/Iyold)>tol && Ny<MaxN);
	
	return Iy;
	}


//Evaluation function of the integrant
double Ekernal(double rr,double eeng)
	{
	double temp;

	temp=0.01/(eeng+1.0*exp(-rr/0.2)/rr-0.01/(rr*rr));
	if(temp>0)
		return sqrt(temp)/(rr*rr);
	else
		{
		printf("r value in forbidden region ... stop ...\n");
		exit(1);
		}
	}

//
//
// vector and matrix allocations
//
//

double *dvector(int nl,int nh)
{
	double *v;
	
	v=(double *) calloc((unsigned) (nh-nl+1), (unsigned) sizeof(double));
	if (!v) nrerror("allocation failure in dvector()");
	return v-nl;
}

double **dmatrix(int nrl,int nrh,int ncl,int nch)
{
	int i;
	double **m;
	
	m=(double **) calloc((unsigned) (nrh-nrl+1), (unsigned) sizeof(double*));
	if (!m) 
		nrerror("allocation failure 1 in dmatrix()");
	m -= nrl;
	for (i=nrl; i<=nrh;i++)
		{
		m[i]=(double *) calloc((unsigned) (nch-ncl+1), (unsigned) sizeof(double));
		if (!m[i])
			nrerror ("allocation failure 2 in dmatrix()");
		m[i] -= ncl;
		}
	return m;
}

void nrerror(char error_text[])
{	
	fprintf(stderr,"Numerical Recipes run-time error...\n");
	fprintf(stderr,"%s\n",error_text);
	fprintf(stderr,"...now exiting to system...\n");
	exit(1);
}