#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 conduttori a potenziale +10 e -10  Volt
// (a1=0)<x<(a2=20)
// (b1=0)<y<(b2=10)
//  P1(5,5)=-10
//  P2(15,5)=10
//
void Esempio2(int iMethod=0) {

    gStyle->SetOptStat(0);
   // estremi della regione da esplorare
    double a1=0.,a2=20.,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]; TCanvas * c_array[11]; 
    for (int i=0; i<11; ++i)
      {
	std::string lab = "h_"+std::to_string(i);
	std::string tit = "Soluzione step = 100 * "+std::to_string(i+1);
	c_array[i] = new TCanvas(lab.c_str(),tit.c_str(),800,600);
	h2[i] = new TH2F(lab.c_str(),tit.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)/4,cx2=3*(nbinx+1)/4,cy1=(nbiny+1)/2,cy2=(nbiny+1)/2;
    cout<<" cx1, 2, cy "<<cx1<<" "<<cx2<<" "<<cy1<<" "<<cy2<<endl;
    h2tmp->SetBinContent(cx1,cy1,-5.);
    h2tmpold->SetBinContent(cx1,cy1,-5.);
    h2tmp->SetBinContent(cx2,cy2,5.);
    h2tmpold->SetBinContent(cx2,cy2,5.);

// Valore usato per la costante di rilassamento
    double w=0.5;
// procedura a numero di step fisso
    int NStep=1000;
    int ihisto= 0;
    //NStep = 1;
    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));
	}
      }
      var=0;
      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.;      
      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) continue;
             if (i==cx2 && j==cy2) continue;
             // algoritmo
             V_ij_old=h2tmpold->GetBinContent(i,j);
	     if (iMethod==0)
	       {
		 // Jakobi 
		 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);
	       }
	     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);
             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));
	  }
	}
      }
// salvo gli istogrammi ogni 100 step
      if ((k-1)%100==0)
	{
	  std::cout<<"Histogram ihisto ="<<ihisto<<" saved; for step k = "<<k<<std::endl;
	  c_array[ihisto]->cd();
	  h2[ihisto]->Draw("ZCOL");
	  //h2[ihisto]->Write();
	  ihisto++;
	}
    }

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

    f->Close();

    


}