//This prgram will evaluate the electric potential of a conducting sphere with V(theta) given at r=1
//	-the solution is expanded into a series of Legendre polynomials
//	-the coefficients to the Legendre series is by integrating P(theta) with V(theta)
//	-the number of terms is truncated so that the boundary condition is satified up to single precision
//NOTE: With the change of variable x=cos(theta), the integral can be simplified
//		I[theta]->I(x)=I[e^(-x^2)P_n(x)dx, {x,-1,1}]
//		Since e^(-x^2) is symmetric, I[x]=0 for n=odd
//			this gives ai[odd]=0 
#include <stdio.h>
#include <math.h>

//accuracy for Simpson's Rule
#define EPS 1.0e-7 //relative accuracy
#define MaxM 100000 //maximum number of intervals

double ai[100];//coefficients to the Legendre series

double legendre(int,double);
double integrateSR(int pindex,double tol);
double Ekernal(int pindex,double xx);
int getcoefficient(void);


main()
{
	FILE *outfile;
	//plotting variables
	int Nr,Nx;
	double rstep,xstep;
	double r,x;
	int i,j;
	double powr;
	double potential;

	//number of terms in Legendre series needed
	int maxN;
	int n;

	maxN=getcoefficient();

	//
	//plotting section
	//	inside of the sphere: [r:0,1]x[t:0,pi]->[r:0,1]x[x:-1,1]
	//
	Nr=200;
	Nx=200;
	rstep=1.0/Nr;
	xstep=2.0/Nx;//[0,pi]->[1,-1]
	
	outfile=fopen("EMsphere.txt","w");
	fprintf(outfile,"rval\txval\tpo\n");

	r=0.0;
	for(i=0;i<=Nr;i++)
	{
		x=-1.0;
		for(j=0;j<=Nx;j++)
		{
			potential=ai[0];// +ai[1]*r*x;//ai[1]=0 by symmetry
			powr=r;//keeping track of the power of r
			for(n=2;n<=maxN;n+=2)//ai[odd]=0 by symmetry
				{
				powr*=r;
				potential+=(ai[n]*powr*legendre(n,x));
				}

			fprintf(outfile,"%lg\t%lg\t%lg\n",r,x,potential);
			x+=xstep;
		}
		r+=rstep;
	}
	fclose(outfile);

	return 0;
}


int getcoefficient(void)
{ 
	int n;
	double x;
	double rhs,lhs;

	//pick a random point on the sphere r=1 for checking b.c.
	x=0.5467812;
	lhs=exp(-x*x);

	//calculate ai[0] and ai[1] first
	ai[0]=integrateSR(0,EPS)/2.0;
	// ai[1]=3*integrateSR(1,EPS)/2.0; - ai[1]=0 by symmetry
	rhs=ai[0];// ai[1]=0 by symmetry

	n=0;
	do
	{
		n+=2;
		//generating coefficients
		ai[n]=(2.0*n+1.0)*integrateSR(n,EPS)/2.0;

		//check sum with boundary condition
		rhs+=(ai[n]*legendre(n,x));

	}
	while(fabs((rhs-lhs)/rhs)>EPS && n<99);

	if(n>=100)
		printf("Too many coefficients but with not enough tolerance\n");

	return n;
}


// Simpson's Rule 

double integrateSR(int nord,double tol)
{
	int i;
	int Ny;
	double hy,xmin;
	double Sye,Syo,Syd;
	double Iy,Iyold;

	Ny=6;
	//the integration range is remapped from [0,pi] to [-1,1]
	xmin=-1.0;
	hy=2.0/Ny;

	Syd=Ekernal(nord,-1.0)+Ekernal(nord,1.0);
	Syo=Sye=0.0;
	for(i=1;i<Ny;i+=2)
		Syo+=Ekernal(nord,i*hy+xmin);
	for(i=2;i<Ny;i+=2)
		Sye+=Ekernal(nord,i*hy+xmin);
	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(nord,i*hy+xmin);
			Iy=hy/3.0*(Syd+2.0*Sye+4.0*Syo);

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

//Evaluation function of the integrant
//a change of variable [x=cos(theta)] is used
//with x, the integrant is -e^(-x^2)*P(x)dx
//	it will call legendre()
double Ekernal(int porder, double xx)
{
	if(porder==0)
		return exp(-xx*xx);
	if(porder==1)
		return exp(-xx*xx)*xx;
	else
		return exp(-xx*xx)*legendre(porder,xx);
}

//	Legendre Polynomials

double legendre(int n,double xx)
//this evaluate the nth order Legendre polynomials at x
{
	int i;
	double pn0,pn1,pn;

	pn0=1;
	pn1=xx;

	for(i=2;i<=n;i++)
	{
		pn=2.0*xx*pn1-pn0-(xx*pn1-pn0)/i;

		pn0=pn1;
		pn1=pn;
	}

	return pn;
}