//this program integrates the square charge distribution using two different methods:
//	Monte Carlo and 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>
#include "numcip.h"

//constants for ran1 and 2(int *)
#define M1 259200
#define IA1 7141
#define IC1 54773
#define RM1 (1.0/M1)
#define M2 124456
#define IA2 8121
#define IC2 28411
#define RM2 (1.0/M2)
#define M3 243000
#define IA3 4561
#define IC3 51349

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

double pi;


//dig variables
FILE *outfile;
double *hist;

double integrateMC(double xxp, double yyp, double tol, long *count);
double Ekernal(double xx,double yy,double xxp,double yyp);
float rann1(int *idum);
float rann2(int *idum);


double integrateSR(double xxp,double yyp,double tol,long *count);
double innerSR(double xx,double xxp,double yyp,double tol,long *count);

main()
	{
	double xp,yp;
	double phi,phi2;
	long counter,counter2;
	int choice;

	pi=4.0*atan(1.0);
	//userio
	printf("Picking a random observation point?\n");
	printf("\t[0]outside of charge distribution\n");
	printf("\t[1]inside of charge distribution\n");
	scanf("%d",&choice);
	if(choice)
	{
		//inside charge distribution
		xp=pi/5.3;
		yp=pi/4.2;
	}
	else
	{
		//inside charge distribution
		xp=pi/2.0;
		yp=pi/3.0;
	}
	printf("\tx_p:%lg\n",xp);
	printf("\ty_p:%lg\n",yp);


	//integration
	phi=integrateMC(xp,yp,EPS,&counter);
	phi2=integrateSR(xp,yp,EPS,&counter2);


	printf("\nFor the accuracy of %g, the electric potential is:\n",EPS);
	printf("\tMonte Carlo method:%lg with %ld evaluation pts\n",phi,counter);
	printf("\tSimpson's Rule:%lg with %ld evaluation pts\n",phi2,counter2);

	outfile=fopen("efficiency_test.dat","w");
	fprintf(outfile,"\nFor the accuracy of %g, the electric potential is:\n",EPS);
	fprintf(outfile,"Observation point: [%lg,%lg]\n\n",xp,yp);
	fprintf(outfile,"\tMonte Carlo method:%lg with %ld evaluation pts\n",phi,counter);
	fprintf(outfile,"\tSimpson's Rule:%lg with %ld evaluation pts\n",phi2,counter2);

	return 0;
	}

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


// Monte Carlo Integrator
double integrateMC(double xxp, double yyp, double tol, long *count)
//integral = sum( f(x) del_vol) and x is sample at random location within [-1,1]X[-1,1]
// del_vol - differential vol = total vol/number of evaluation points.
//adapted from Box 2-5
	{
	int seed1,seed2;//seeds for random generator with uniform distribution
	long Num,i;
	double xx,yy;//integration variables
	double totalvol;
	double val1,val2;
	int inflag;

	FILE *digfile;//outfput file for MC convergence rate plot
	
	seed1=-1;//initialize random generator1
	seed2=-3;//initialize random generator2

	Num=10;
	totalvol=4.0;//total volume of charge distribution
	val1=val2=0.0;

	digfile=fopen("mc_rate.dat","w");
	//calculate val1 with Num
	for(i=1;i<=Num;i++)
		{
		//find a random points with [-1,1]X[-1,1]
		xx=(rann1(&seed1)-0.5)*2.0;
		yy=(rann2(&seed2)-0.5)*2.0;

		//summing integrant at xx,yy
		val1+=Ekernal(xx,yy,xxp,yyp);
		}
	val1*=(totalvol/(double)Num);

	fprintf(digfile,"%lg\t%lg\n",log10((double)Num),val1);
	
	//calculate val2 with a different set of Num
	inflag=1;
	do
		{
		for(i=1;i<=Num;i++)
			{
			//find a random points with [-1,1]X[-1,1]
			xx=(rann1(&seed1)-0.5)*2.0;
			yy=(rann2(&seed2)-0.5)*2.0;

			//summing integrant at xx,yy
			val2+=Ekernal(xx,yy,xxp,yyp);
			}
		val2*=(totalvol/(double)Num);

		//check for convergence
		if(fabs((val1-val2)/val1)<tol)
			{
			inflag=0;
			*count=2*Num;
			}
		else
			{
			val1=(val1+val2)/2.0;
			
			val2=0.0;
			Num+=Num;

			fprintf(digfile,"%lg\t%lg\n",log10((double)Num),val1);
			}
		}
	while (inflag && Num<MaxN);

	
	if(Num<MaxN)
	{
		val1=(val1+val2)/2.0;
		fprintf(digfile,"%lg\t%lg\n",log10((double)*count),val1);

		return val1;
	}
	else
		{
		*count=Num;
		fprintf(digfile,"%lg\t%lg\n",log10((double)*count),val1);

		return val1;
		}


	}

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

	Nx=6;
	hx=2.0/Nx;
	Sxd=innerSR(-1.0,xxp,yyp,tol,&incount)+innerSR(1.0,xxp,yyp,tol,&incount2);
	*count=incount+incount2;
	Sxo=Sxe=0.0;
	for(i=1;i<Nx;i+=2)
		{
		Sxo+=innerSR(-1.0+i*hx,xxp,yyp,tol,&incount);
		*count+=incount;
		}
	for(i=2;i<Nx;i+=2)
		{
		Sxe+=innerSR(-1.0+i*hx,xxp,yyp,tol,&incount);
		*count+=incount;
		}
	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,&incount);
				*count+=incount;
				}
			Ix=hx/3.0*(Sxd+2.0*Sxe+4.0*Sxo);

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

	return Ix;
	}

double innerSR(double xx,double xxp,double yyp,double tol,long *count)
	{
	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;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);
	
	*count=Ny+1;
	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));
	}


//
//
//Random Number generators
//
//
//to generate uniform distribution in 2D, we used two identical
//	INDEPENDENT generators ran1 and ran2
float rann1(int *idum)
{
	static long ix1,ix2,ix3;
	static float r[98];
	float temp;
	static int iff=0;
	int j;
	
	extern void nrerror(char []);

	if (*idum < 0 || iff == 0)
		{
		iff=1;
		ix1=(IC1-(*idum)) % M1;
		ix1=(IA1*ix1+IC1) % M1;
		ix2=ix1 % M2;
		ix1=(IA1*ix1+IC1) % M1;
		ix3=ix1 % M3;
		for (j=1;j<=97;j++)
			{
			ix1=(IA1*ix1+IC1) % M1;
			ix2=(IA2*ix2+IC2) % M2;
			r[j]=(ix1+ix2*RM2)*RM1;
			}
		*idum=1;
		}
	ix1=(IA1*ix1+IC1) % M1;
	ix2=(IA2*ix2+IC2) % M2;
	ix3=(IA3*ix3+IC3) % M3;
	j=1 + ((97*ix3)/M3);
	if (j > 97 || j < 1)
		nrerror("ran1: This cannot happen.");
	temp=r[j];
	r[j]=(ix1+ix2*RM2)*RM1;
	return temp;
}


float rann2(int *idum)
{
	static long ix1,ix2,ix3;
	static float r[98];
	float temp;
	static int iff=0;
	int j;
	
	extern void nrerror(char []);

	if (*idum < 0 || iff == 0)
		{
		iff=1;
		ix1=(IC1-(*idum)) % M1;
		ix1=(IA1*ix1+IC1) % M1;
		ix2=ix1 % M2;
		ix1=(IA1*ix1+IC1) % M1;
		ix3=ix1 % M3;
		for (j=1;j<=97;j++)
			{
			ix1=(IA1*ix1+IC1) % M1;
			ix2=(IA2*ix2+IC2) % M2;
			r[j]=(ix1+ix2*RM2)*RM1;
			}
		*idum=1;
		}
	ix1=(IA1*ix1+IC1) % M1;
	ix2=(IA2*ix2+IC2) % M2;
	ix3=(IA3*ix3+IC3) % M3;
	j=1 + ((97*ix3)/M3);
	if (j > 97 || j < 1)
		nrerror("ran1: This cannot happen.");
	temp=r[j];
	r[j]=(ix1+ix2*RM2)*RM1;
	return temp;
}