//risolvere problema elettrostatica condensatore piano facce paralle distanti 0.1 e di lunghezza 1, una faccia carica 1Kv e l'altra -1Kv

#define Nx 201
#define Ny 111

void Condensatore()
{
	double Soluzione[Nx][Ny]={0};
	double h=0.01;
	double xGriglia[Nx]={0}, yGriglia[Ny]={0}; //li riempiamo con i valori per capire come graficare il condensatore nella griglia
	double SolGS[Nx][Ny]={0}; //Gauss-S
	double SolRil[Nx][Ny]={0}; //Rilassamento

	for(int i=0;i<Nx;i++)
	{
		xGriglia[i]=h*i;
	}

	for(int j=0;j<Ny;j++)
	{
		yGriglia[j]=h*j;
	}

	printf ("xMax %f  yMax %f",xGriglia[Nx-1],yGriglia[Ny-1]);

//condizioni iniziali
	for (int i=50;i<=150;i++)
	{
		Soluzione[i][50]=-1000;
		Soluzione[i][60]=1000;

	}

	for (int j=0;j<Nx;j++)
	{
		Soluzione[j][0]=0;
		Soluzione[j][Ny-1]=0;
	}

	for (int k=0;k<Ny;k++)
	{
		Soluzione[0][k]=0;
		Soluzione[Nx-1][k]=0;
	}

	//controllo condizione iniziale attraverso un istogramma con lo stesso numero di bin della matrice in modo tale che il cetro di ogni bin coincida con il valore 

	TH2F* hCondIn= new TH2F("hCondIn","hCondIn",Nx,-h/2,xGriglia[Nx-1]+h/2,Ny,-h/2,yGriglia[Ny-1]+h/2); //-h/2 per avere valore 0 centrato

	for(int i=0;i<Nx;i++)
	{
		for(int j=0;j<Ny;j++)
		{
			hCondIn->SetBinContent(i+1,j+1,Soluzione[i][j]); //aggiungiamo +1 perchè root usa il bin 0 per gli underflow
		}
	}

	//per jacobi ci serve una matrice temporanea che conservi la soluzione al passo precedente
	double MatrTemp[Nx][Ny];
	double TempGS[Nx][Ny]; //matrice tempr per G-S
	double TempRil[Nx][Ny];
	

	//inizializziamo la matrice con la condiz iniz
	for(int j=0;j<Ny;j++)
	{
		for(int i=0;i<Nx;i++)
		{ //matrici uguali con stesse cond iniz
			MatrTemp[i][j]=Soluzione[i][j];
			SolGS[i][j]=Soluzione[i][j];
			TempGS[i][j]=Soluzione[i][j];
			SolRil[i][j]=Soluzione[i][j];
			TempRil[i][j]=Soluzione[i][j];
		}
	}

	double err=0, errmax=1, errloc=0; //Jacobi
	double errGS=0, errlocGS=0; //G-S 
	double errRil=0, errlocRil=0; //Rilassamento
	int k=0, kGS=0, kRil=0;
	double w=0.5; //omega per rilass
	double besterrGS=10, besterrJ=10, besterrRil=10; //devono esser> errmax

	do
	//usiamo ciclo do while perche il metodo di Jacobi termina una volta superata una certa soglia, usiamo il do while così possiamo fare una prima iterazione, dato che abbiamo err=0
	{
		err=0;
		errGS=0;
		errRil=0;
		//aggiorniamo i valori della matrice
		printf("besterrJ: %f   besterrGS:%f  erroreJacobi:%f  erroreGS:%f  k:%d  kGS:%d \n ",besterrJ, besterrGS, err,errGS,k,kGS);
		for (int i = 1; i < Nx-1; ++i) //partiamo da 1 fino a Nx-1 per non modificare la cond iniz
		{
			for (int j = 1; j< Ny-1; ++j)
			{
				if(j==60||j==50) //per non mod la cond iniz
				{
					if(i>=50&&i<=150)
						continue; //continue fa passare all'iterazione successiva
				}

				
				//JACOBI (serve matrice temporanea)
				if(besterrJ>errmax)
				{
					Soluzione[i][j]=(MatrTemp[i][j-1]+MatrTemp[i][j+1]+MatrTemp[i-1][j]+MatrTemp[i+1][j])/4;
					errloc=fabs(Soluzione[i][j]-MatrTemp[i][j]);//errore locale sul singolo punto. ci serve l'errore globale, cioe l'errore massimo su tutta la matrice
					if(errloc>err)
						err=errloc;
				}
				

				//GS
				if(besterrGS>errmax)
				{
					
					SolGS[i][j]=(SolGS[i][j-1]+TempGS[i][j+1]+SolGS[i-1][j]+TempGS[i+1][j])/4;
					errlocGS=fabs(SolGS[i][j]-TempGS[i][j]);//errore locale sul singolo punto. ci serve l'errore globale, cioe l'errore massimo su tutta la matrice
					if(errlocGS>errGS)
						errGS=errlocGS;
				}
				
				//Ril
				if(besterrRil>errmax)
				{
					
					SolRil[i][j]=(1-w)*TempRil[i][j]+w*(SolRil[i][j-1]+TempRil[i][j+1]+SolRil[i-1][j]+TempRil[i+1][j])/4;
					errlocRil=fabs(SolRil[i][j]-TempRil[i][j]);//errore locale sul singolo punto. ci serve l'errore globale, cioe l'errore massimo su tutta la matrice
					if(errlocRil>errRil)
						errRil=errlocRil;
				}
			}
		}
		for(int j=0;j<Ny;j++)
		{
			for(int i=0;i<Nx;i++)
			{
				MatrTemp[i][j]=Soluzione[i][j];
				TempGS[i][j]=SolGS[i][j];
				TempRil[i][j]=SolRil[i][j];
			}
		}
		if(besterrJ>errmax)
			k++;
		if(besterrGS>errmax)
			kGS++;
		if(besterrRil>errmax)
			kRil++;
		if(err>0) //aggiorniamo i besterr solo se err>0
			besterrJ=err;
		if(errGS>0)
			besterrGS=errGS; 
		if(errRil>0)
			besterrRil=errRil;

		if(k==kGS)
			printf("\nIterazione num %d    Errore massimo Jacobi %f   Errore massimo GS %f\n",k,err,errGS);
		if(kGS<k) 
			printf ("\nIterazione num %d    Errore massimo Jacobi %f\n",k,err);
		if(kGS>k)
			printf ("\nIterazione num %d    Errore massimo GS %f\n",kGS,errGS);
	}
	while (err>errmax||errGS>errmax||errRil>errmax);//continua fino a quando non sono soddisfatte entrabe le condiz

	printf("k: %d  kGS:%d  kRil:%d besterrJ:%f  besterrGS:%f  besterrRil:%f",k,kGS,kRil,besterrJ,besterrGS,besterrRil);
		
 //disegnare la soluzione con istogramma, sia per ultima matrice che per una intermedia
		//impleentare Gauss sidel e rilassamento
		TH2F* hFinale= new TH2F("hFinale","hFinale",Nx,-h/2,xGriglia[Nx-1]+h/2,Ny,-h/2,yGriglia[Ny-1]+h/2);
		TH2F* hGS=new TH2F("hGS","hGS",Nx,-h/2,xGriglia[Nx-1]+h/2,Ny,-h/2,yGriglia[Ny-1]+h/2);
		TH2F* hRil=new TH2F("hRil","hRil",Nx,-h/2,xGriglia[Nx-1]+h/2,Ny,-h/2,yGriglia[Ny-1]+h/2);

	for(int i=0;i<Nx;i++)
	{
		for(int j=0;j<Ny;j++)
		{
			hFinale->SetBinContent(i+1,j+1,Soluzione[i][j]); //aggiungiamo +1 perchè root usa il bin 0 per gli underflow
			hGS->SetBinContent(i+1,j+1,SolGS[i][j]);
			hRil->SetBinContent(i+1,j+1,SolRil[i][j]);
		}
	}


}	

//confrontare i metodi
//Jacobi potrebbe essere il più veloce (distanza tra due iterazioni minore), ma non è detto che sia il più vicino alla soluzione reale.
//bisogna confrontare i metodi con errmax=0.1(più vicina alla sol reale) e con errmax=1, conservare le soluzion in un istogramma e fare la differenza