#include <TFile.h>
// determiniamo la soluzione dell'equazione di Laplace
// con il metodo del rilassamento (o di Gauss-Seidel)
// all'interno di una gabbia a potenziale 0 con all'interno due piani a potenziale +10 e -10  Volt
// (a1=0)<x<(a2=10)
// (b1=0)<y<(b2=10)
//  P1(4.5,2.5-7.5)=-10
//  P2(5.5,2.5-7.5)=10
//
void Esempio4(int iMethod=0, float w=1.) {

  gStyle->SetOptStat(0);
   // estremi della regione da esplorare
    double a1=0.,a2=10.,b1=0.,b2=10.;
   // Passo per la discretizzazione
   double h=0.1;

   int k=0;  // step
    double var=0;  // variazione massima rispetto allo step precedente
// questo istrogramma sarą salvato nel TTree e conterrą ad ogni passo
// l'evoluzione del campo potenziale
// i sui estremi e il numero di bin dipendono dalla regione da analizzare
    int nbinx=(a2-a1)/h,nbiny=(b2-b1)/h;

    TH2D * h2tmp    = new TH2D("h2_old","Soluzione old",nbinx,a1,a2,nbiny,b1,b2);
    TH2F * h2[11]; TH2F * hE[11]; TH2F * hEx[11]; TH2F * hEy[11]; TCanvas * c_array[11]; 
    for (int i=0; i<11; ++i)
      {
	std::string labV  = "hV_"+std::to_string(i);
	std::string labE  = "hE_"+std::to_string(i);
	std::string labEx = "hEx_"+std::to_string(i);
	std::string labEy = "hEy_"+std::to_string(i);
	std::string titV  = "V Soluzione step = 100 * "+std::to_string(i+1);
	std::string titE  = "E Soluzione step = 100 * "+std::to_string(i+1);
	std::string titEx = "Ex Soluzione step = 100 * "+std::to_string(i+1);
	std::string titEy = "Ey Soluzione step = 100 * "+std::to_string(i+1);
	c_array[i] = new TCanvas(labV.c_str(),titV.c_str(),800,600);
	h2[i] = new TH2F(labV.c_str(),titV.c_str(),nbinx,a1,a2,nbiny,b1,b2);
	hE[i] = new TH2F(labE.c_str(),titE.c_str(),nbinx,a1,a2,nbiny,b1,b2);
	hEx[i] = new TH2F(labEx.c_str(),titEx.c_str(),nbinx,a1,a2,nbiny,b1,b2);
	hEy[i] = new TH2F(labEy.c_str(),titEy.c_str(),nbinx,a1,a2,nbiny,b1,b2);
	std::cout<<"Histogram named "<<h2[i]->GetName()<<" created"<<std::endl;
      }
    
    TH2D * h2tmpold = new TH2D("h2_old","Soluzione old",nbinx,a1,a2,nbiny,b1,b2);

// Metodo del rilassamento a due variabili
//
// inserisco le condizioni al contorno alla Dirichlet:
//
//

   for (int j=1; j<=nbiny;j++) {
     h2tmp->SetBinContent(1,j,0.);
     h2tmpold->SetBinContent(1,j,0.);
     h2tmp->SetBinContent(nbinx,j,0.);
     h2tmpold->SetBinContent(nbinx,j,0.);
   }
   for (int i=2; i<nbinx;i++) {
     h2tmp->SetBinContent(i,1,0.);
     h2tmpold->SetBinContent(i,1,0.);
     h2tmp->SetBinContent(i,nbiny,0.);
     h2tmpold->SetBinContent(i,nbiny,0.);
   }

   int cx1=(nbinx+1)/2-10,cx2=(nbinx+1)/2+10;
   int cy1=(nbiny+1)/4,cy2=(nbiny+1)*3/4;
   cout<<cx1<<" "<<cx2<<endl;
   cout<<cy1<<" "<<cy2<<endl;
   for (int j=cy1;j<=cy2;j++) {
     h2tmp->SetBinContent(cx1,j,-10.);
     h2tmpold->SetBinContent(cx1,j,-10.);
     h2tmp->SetBinContent(cx2,j,10.);
     h2tmpold->SetBinContent(cx2,j,10.);
   }
//
// Valore usato per la costante di rilassamento
   //double w=0.5;
// procedura a numero di step fisso
    int ihisto= 0;
    int NStep=1000;
    for (k=1; k<NStep+2; k++) {

      for (int j=1;j<nbiny; j++) {
	for (int i=1;i<nbinx;i++) {
	  h2tmp->SetBinContent(i,j,h2tmpold->GetBinContent(i,j));
	}
      }
      double V_ij_old=0.;
      double V_imj   =0.;
      double V_ipj   =0.;
      double V_ijm   =0.;
      double V_ijp   =0.;
      double V_ij    =0.;
      var=0;
      for (int j=2;j<nbiny; j++) {
          for (int i=2;i<nbinx;i++) {
            // per non alterare il potenziale nei punti a potenziale fisso.
             if (i==cx1 && j>=cy1&& j<=cy2) continue;
             if (i==cx2 && j>=cy1&& j<=cy2) continue;
             // algoritmo
	     V_ij_old=h2tmpold->GetBinContent(i,j);
	     if (iMethod==0)
	       {
		 // Jacobi
		 V_imj   =h2tmpold->GetBinContent(i-1,j);
		 V_ipj   =h2tmpold->GetBinContent(i+1,j);
		 V_ijm   =h2tmpold->GetBinContent(i,j-1);
		 V_ijp   =h2tmpold->GetBinContent(i,j+1);
		 V_ij    =1/4.*(V_imj+V_ipj+V_ijm+V_ijp);
	       }
	     else if (iMethod==1)
	       {
		 // Gauss Seidel
		 V_imj   =h2tmp->GetBinContent(i-1,j);
		 V_ipj   =h2tmp->GetBinContent(i+1,j);
		 V_ijm   =h2tmp->GetBinContent(i,j-1);
		 V_ijp   =h2tmp->GetBinContent(i,j+1);
		 V_ij    =1/4.*(V_imj+V_ipj+V_ijm+V_ijp);
		 // rilassamento
		 V_ij=V_ij*w+(1.-w)*V_ij_old;
	       }
             h2tmp->SetBinContent(i,j,V_ij);
             double localvar=fabs(V_ij-V_ij_old);
             if (localvar>var)var=localvar;
          }
      }
      std::cout<<"iterazione k = "<<k<<" massima variazione rispetto a iter k-1 = "
	       <<var
	       <<std::endl;
            // salvo la matrice (soluzione numerica sulla griglia discretizzata) all'iterazione corrente per il calcolo al prossimo passo 
      for (int j=1;j<nbiny; j++) {
	for (int i=1;i<nbinx;i++) {
	  h2tmpold->SetBinContent(i,j,h2tmp->GetBinContent(i,j));
	  if ((k-1)%100==0) {
	    h2[ihisto]->SetBinContent(i,j,h2tmp->GetBinContent(i,j));

             double Ex_ij,Ey_ij;
	     double V_ij=h2tmp->GetBinContent(i,j);
             if (i==1) {
                double V_ipj=h2tmp->GetBinContent(i+1,j);
                Ex_ij=(V_ij-V_ipj)/h;
             } else if (i==nbinx) {
                double V_imj=h2tmp->GetBinContent(i-1,j);
                Ex_ij=(V_ij-V_imj)/h;
             } else {
               double V_imj=h2tmp->GetBinContent(i-1,j);
               double V_ipj=h2tmp->GetBinContent(i+1,j);
               Ex_ij=(V_ipj-V_imj)/2./h;

             }
             if (j==1) {
                double V_ijp=h2tmp->GetBinContent(i,j+1);
                Ey_ij=(V_ij-V_ijp)/h;
             } else if (j==nbiny) {
                double V_ijm=h2tmp->GetBinContent(i,j-1);
                Ey_ij=(V_ij-V_ijm)/h;
             } else {
               double V_ijm=h2tmp->GetBinContent(i,j-1);
               double V_ijp=h2tmp->GetBinContent(i,j+1);
               Ey_ij=(V_ijp-V_ijm)/2./h;

             }
             double E=sqrt(Ex_ij*Ex_ij +Ey_ij*Ey_ij );
             hE[ihisto]->SetBinContent(i,j,E);
             hEx[ihisto]->SetBinContent(i,j,Ex_ij);
             hEy[ihisto]->SetBinContent(i,j,Ey_ij);


	    
	  }
	}
      }

// riempio TTree ogni 100 step
      if ((k-1)%100==0)
	{
	  std::cout<<"Histogram ihisto ="<<ihisto<<" saved; for step k = "<<k<<std::endl;
	  c_array[ihisto]->Divide(2,2);
	  c_array[ihisto]->cd(1);
	  h2[ihisto]->Draw("ZCOL");
	  c_array[ihisto]->cd(2);
	  hE[ihisto]->Draw("ZCOL");
	  c_array[ihisto]->cd(3);
	  hEx[ihisto]->Draw("ZCOL");
	  c_array[ihisto]->cd(4);
	  hEy[ihisto]->Draw("ZCOL");
	  //h2[ihisto]->Write();
	  ihisto++;
	}
    }

//

    // Apro file
    TFile * f=new TFile("Esempio4.root","RECREATE");
    for (int i=0; i<11; ++i)
      {
	h2[i]->Write();
	hE[i]->Write();
	hEx[i]->Write();
	hEy[i]->Write();
	std::cout<<"Histogram named "<<h2[i]->GetName()<<" saved on file"<<std::endl;
      }

    f->Close();



}