//this program calculate the electric potential for a uniform square at [-1,1]X[-1,1]
//	this is based on assg2team1.c - integration routinue used is Simpson's rule
//NOTE:  if the observation point (x_p,y_p) is within the charge distribution, the
//	integrant is singular at that point but the integral remains integrable.

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

//accuracy for integrator
#define EPS 1.0e-4
#define MaxN 100000000

double pi;

double Ekernal(double xx,double yy,double xxp,double yyp);
double integrateSR(double xxp,double yyp,double tol);
double innerSR(double xx,double xxp,double yyp,double tol);

main()
	{
	int i,j;
	double xp,yp;
	int NNx,NNy;
	double stepx,stepy;
	double offset;
	double phi;

	FILE *outfile;

	pi=4.0*atan(1.0);
	offset=1e-3*pi;

	outfile=fopen("charge.dat","w");
	
	printf("Number of pixel within [-4,4]X[-4,4]\n");
	printf("\tNx:\n");
	scanf("%d",&NNx);
	printf("\tNy:\n");
	scanf("%d",&NNy);
	stepx=8.0/NNx;
	stepy=8.0/NNy;

	xp=-4.0+offset;
	yp=-4.0-offset;
	for(i=0;i<NNx;i++)
		{
		for(j=0;j<NNy;j++)
			{
			phi=integrateSR(xp,yp,EPS);
			fprintf(outfile,"%lg\t",phi);
			xp+=stepx;
			}
		fprintf(outfile,"\n");
		fflush(outfile);
		xp=-4.0+offset;		
		yp+=stepy;
		}
	
	return 0;
	}

//
//
//Integrating Routinues
//
//

// Simpson's Rule
double integrateSR(double xxp,double yyp,double tol)
	{
	int i;
	int Nx;
	double hx;
	double Sxe,Sxo,Sxd;
	double Ix,Ixold;

	Nx=6;
	hx=2.0/Nx;
	Sxd=innerSR(-1.0,xxp,yyp,tol)+innerSR(1.0,xxp,yyp,tol);
	Sxo=Sxe=0.0;
	for(i=1;i<Nx;i+=2)
		{
		Sxo+=innerSR(-1.0+i*hx,xxp,yyp,tol);
		}
	for(i=2;i<Nx-1;i+=2)
		{
		Sxe+=innerSR(-1.0+i*hx,xxp,yyp,tol);
		}
	Ix=hx/3.0*(Sxd+2.0*Sxe+4.0*Sxo);
	
	//
	do	
		{
			Ixold=Ix;
			Nx*=2;
			hx/=2.0;
			Sxe+=Sxo;
			Sxo=0.0;
			for(i=1;i<=Nx;i+=2)
				{
				Sxo+=innerSR(-1.0+i*hx,xxp,yyp,tol);
				}
			Ix=hx/3.0*(Sxd+2.0*Sxe+4.0*Sxo);

		}while(fabs(Ix-Ixold)>tol && Nx<MaxN);

	return Ix;
	}

double innerSR(double xx,double xxp,double yyp,double tol)
	{
	int i;
	int Ny;
	double hy;
	double Sye,Syo,Syd;
	double Iy,Iyold;

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

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

//Evaluation function of the integrant
double Ekernal(double xx,double yy,double xxp,double yyp)
	{
	return 1.0/sqrt((xx-xxp)*(xx-xxp)+(yy-yyp)*(yy-yyp));
	}