#include <TFile.h>
// determiniamo la soluzione dell'equazione di Laplace
// con il metodo del rilassamento (o di Gauss-Seidel)
// all'interno di un sistema di conduttori quadrato a potenziali 0., 10.  Volt
// (a1=0)<x<(a2=10)
// (b1=0)<y<(b2=10)
//
#include <string>

void Esempio1(int iMethod=0) {

  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
  double varrel=0;  // variazione relativa massima rispetto allo step precedente
  double localvar=0;  // variazione locale rispetto allo step precedente
  double localvarrel=0;  // variazione relativa locale rispetto allo step precedente
  
// costruiamo un istogramma che istrogramma conterrą ad ogni passo k l'evoluzione del potenziale
// i sui estremi e il numero di bin dipendono dalla regione da analizzare
// questo istogramma mi serve per aggiornare la mappa del potenziale bin per bin ad ogni iterazione (Gauss Seidel)
    int nbinx=1+(a2-a1)/h,nbiny=1+(b2-b1)/h;
    std::cout<<"estremi della regione in x "<<a1<<" "<<a2<<" --- numero di bin in x ="<<nbinx<<std::endl; 
    std::cout<<"estremi della regione in y "<<b1<<" "<<b2<<" --- numero di bin in y ="<<nbiny<<std::endl; 
    TH2D * h2tmp    = new TH2D("h2","Soluzione temporanea",nbinx,a1,a2,nbiny,b1,b2);
// questo istogramma mi serve per avere a disposizione ad ogni iterazione la soluzione completa
// all'iterazione precedente (Jacobi)
    TH2D * h2tmpold = new TH2D("h2_tmpold","Soluzione old",nbinx,a1,a2,nbiny,b1,b2);

// conserviamo in questi 11 istogrammi le soluzioni ottenute ogni 100 iterazioni    
    TH2F * h2[11]; TCanvas * c_array[11]; 
    for (int i=0; i<11; ++i)
      {
	std::string lab = "h_"+std::to_string(i+1);
	std::string tit = "Soluzione step = 1 + 100 * "+std::to_string(i);

	//std::string lab = label[i];
	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;
      }

//
// Metodo del rilassamento a due variabili
//
// inserisco le condizioni al contorno alla Dirichlet:
//
// pareti verticali a V=0    
    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.);
    }
// pareti orizzontali a V=0    
    for (int i=1; i<=nbinx;i++) {
      h2tmp->SetBinContent(i,1,0.);
      h2tmpold->SetBinContent(i,1,0.);
      h2tmp->SetBinContent(i,nbiny,0.);
      h2tmpold->SetBinContent(i,nbiny,0.);
    }
// filo centrale a V=10 V     
    int cx=(nbinx+1)/2,cy=(nbiny+1)/2;
    cout<<"Bin del filo "<<cx<<" "<<cy<<endl;
    h2tmp->SetBinContent(cx,cy,10.);
    h2tmpold->SetBinContent(cx,cy,10.);
    
// Valore usato per la costante di rilassamento
    //double w=0.5;
    double w=1.;
// procedura a numero di step fisso
    int NStep=5000;
    int ihisto= 0;

    /// iterazioni ...     
    for (k=1; k<NStep+2; k++) {

      /// inizializziamo l'istogramma che descrive la siluzione completa all'iterazione precedente (k-1)
      // allo stato dell'istogramma temporaneo 
      for (int j=1;j<nbiny; j++) {
	for (int i=1;i<nbinx;i++) {
	  h2tmp->SetBinContent(i,j,h2tmpold->GetBinContent(i,j));
	}
      }
      var=0;  /// variazione assoluta massima sulla griglia di discretiozzazioma tra iterazione k e k-1
      localvar=0.;
      varrel=0;  /// variazione relativa massima sulla griglia di discretiozzazioma tra iterazione k e k-1
      localvarrel=0.;
      for (int j=2;j<nbiny; j++) {
          for (int i=2;i<nbinx;i++) {
            // per non alterare il potenziale nel punto centrale.
             if (i==cx && j==cy) continue;
             // algoritmo  
             double V_ij_old=h2tmpold->GetBinContent(i,j);
             double V_imj   =h2tmp->GetBinContent(i-1,j);
             double V_ipj   =h2tmp->GetBinContent(i+1,j);
             double V_ijm   =h2tmp->GetBinContent(i,j-1);
             double V_ijp   =h2tmp->GetBinContent(i,j+1);
             double V_ij    =(V_imj+V_ipj+V_ijm+V_ijp)/4.;
	     //Rilassamento ????
             V_ij = V_ij*w + (1.-w)*V_ij_old;
	     h2tmp->SetBinContent(i,j,V_ij);
	     localvar=fabs(V_ij-V_ij_old);
	     if (fabs(V_ij)>1.E-6) localvarrel=fabs((V_ij-V_ij_old)/V_ij);
             if (localvar   >var   )  var   =localvar   ;
             if (localvarrel>varrel)  varrel=localvarrel;
          }
      }
      
      if ((k-1)%100==0) std::cout<<"iterazione k = "<<k<<" massima variazione (relativa) rispetto a iter k-1 = "
	       <<var<<" ("<<varrel<<")"
	       <<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 && ihisto<11) {
	    h2[ihisto]->SetBinContent(i,j,h2tmp->GetBinContent(i,j));
	  }
	}
      }
// salvo gli istogrammi ogni 100 step
      if ((k-1)%100==0 && ihisto<11) {
	std::cout<<"Histogram ihisto ="<<ihisto<<" saved; for step k = "<<k<<std::endl;
	c_array[ihisto]->cd();
	h2[ihisto]->Draw("ZCOL");
        //h2[ihisto]->Write();
	ihisto++;
      }
    }

    // soluzione finale
    TCanvas* Csoluzione = new TCanvas("Csol","soluzione",800,600);
    h2tmp->Draw("ZCOL");
    
    
//
    // Apro file
    TFile * f=new TFile("Esempio1.root","RECREATE");
    for (int i=0; i<11; ++i)
      {
	h2[i]->Write();
	std::cout<<"Histogram named "<<h2[i]->GetName()<<" saved on file"<<std::endl;
      }
    h2tmp->Write();
    std::cout<<"Histogram named "<<h2tmp->GetName()<<" saved on file"<<std::endl;

    f->Close();





}