#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>

// soluzione analitica
double EAnalitico(double x, double y, double dphi, double r_c, double r_w, double xc, double yc)
{
  return  (dphi / log(r_c/r_w))/sqrt( (x-xc)*(x-xc)+(y-yc)*(y-yc) );
}
double phiAnalitico(double x, double y, double dphi, double r_c, double r_w, double xc, double yc)
{
  return  (500. - (dphi / log(r_c/r_w))*log( sqrt( (x-xc)*(x-xc)+(y-yc)*(y-yc) )/r_w ));
}
  
void CondCilindrico(int iMethod=-1, double omega=1.) {

  if (iMethod==-1)
    {
      std::cout<<"Argomento di ingresso non specificato. "<<std::endl;
      std::cout<<"Esegui la macro con: "<<std::endl;
      std::cout<<".x CondCilindrico.C([metodo], [omega])          [metodo=0 per Jacobi, 1 per Gauss-Seidel] "<<std::endl;
      std::cout<<"                                                [omega (tra 0 e 1[default])] parametro di rilassamento; corregge Gauss Seidel se < 1]"<<std::endl;
      return;
    }
  // Jacobi       0.05 e 50000 ==>>> errore max 185 
  // Gauss Seidel 0.05 e 50000 ==>>> errore max 175
  // Gauss Seidel 0.05 e 50000, omega = 0.9 ==>>> errore max 176.5
  
  
  gStyle->SetOptStat(0);
  // estremi della regione da esplorare
  double a1=0.,a2=22.,b1=0.,b2=22.;
  // Passo per la discretizzazione
  //double h = 0.025;
  double h = 0.05;
  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;
  int nbiny=1+(b2-b1)/h;
  bool enableMatrix[1000][1000] = {false};
  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,800);
	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:
//
    double RCilindro = 10.;
    double RCilindro2 = RCilindro*RCilindro;
    double wparete2 = h*h*2;
    double Vparete = -1000.; 
    double Vfilo =     500.;
    //double Vparete =     0.; 
    //double Vfilo =    1500.;
    
// 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,Vfilo);
    h2tmpold->SetBinContent(cx,cy,Vfilo);

    TH2D* hsa = (TH2D*)h2tmp->Clone();
    hsa->SetTitle("Soluzione Analitica");
    double sa=0;

   double xcenter = h2tmp->GetXaxis()->GetBinCenter(cx);
   double ycenter = h2tmp->GetYaxis()->GetBinCenter(cy);
   std::cout<<" bin del filo centrale, cx, cy = "<<cx<<" "<<cy<<" x,y center = "<<xcenter<<" "<<ycenter<<std::endl;
   std::cout<<" wparete, al quadrato2 = "<<sqrt(wparete2)<<" "<<wparete2<<std::endl;
   double xbin=0.;
   double ybin=0.;
   double dist2 = 0;
   bool firstInRow=true;
// pareti verticali a V=0    
    for (int j=1; j<=nbiny;j++) {
      firstInRow=true;
      for (int k=1; k<=nbinx;k++) {
 	xbin = h2tmp->GetXaxis()->GetBinCenter(k);
 	ybin = h2tmp->GetYaxis()->GetBinCenter(j);

	dist2 = (xbin-xcenter)*(xbin-xcenter)+(ybin-ycenter)*(ybin-ycenter);
	//if (j==300) std::cout<<"k,j "<<k<<" "<<j<<" xbin, ybin, dist2, dist  "<<xbin<<" "<<ybin<<" "<<dist2<<" "<<sqrt(dist2)<<std::endl; 
	if (fabs(sqrt(dist2)-RCilindro) < h/2. ) {
	  // parete catodica
	  if (firstInRow) 
	    std::cout<<" al catodo - riga j = "<<j<<" colonne "<<k;
	  else std::cout<<" "<<k;
	  h2tmp->SetBinContent(j,k,Vparete);
	  h2tmpold->SetBinContent(j,k,Vparete);
	  firstInRow=false;
	}
	else if (sqrt(dist2)<RCilindro) 
	  {
	    if (k==cx && j==cy) continue;
	    enableMatrix[k][j]=true;
	    sa = phiAnalitico(xbin, ybin, Vfilo-Vparete, RCilindro, h/2., xcenter, ycenter);
	    hsa->SetBinContent(k,j,sa);
	    h2tmp->SetBinContent(j,k,Vparete);
	    h2tmpold->SetBinContent(j,k,Vparete);
	  }
      }
      if (!firstInRow)std::cout<<std::endl;
    }

    // disegno le condizioni al contorno
    TCanvas* CCC = new TCanvas("Ccc","Condizioni Al contorno",800,800);
    TH2D* hcc = (TH2D*)h2tmp->Clone();
    hcc->Draw("zcol");
    // disegno la soluzione analitica
    TCanvas* Csa = new TCanvas("Csa","Soluzione analitica",800,800);
    hsa->Draw("zcol");
    //return;
    
// Valore usato per la costante di rilassamento
    double w=1.;
    if (iMethod==1) w=omega; // rilassamento applicato solo in caso di metodo Gauss-Seidel
// procedura a numero di step fisso
    int NStep=50000;
    int ihisto= 0;

    double V_imj   =0.;
    double V_ipj   =0.;
    double V_ijm   =0.;
    double V_ijp   =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.;
      double errore_max = 0.;
      double errore_loc = 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;
	     if (!enableMatrix[i][j]) continue;
             // algoritmo  
             double 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);
	     }
	     else{
	       // 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);
	     }
             double V_ij    =(V_imj+V_ipj+V_ijm+V_ijp)/4.;
	     //Gauss-Seidel + Rilassamento 
             if (iMethod==1) V_ij = w*V_ij + (1.-w)*V_ij_old;
	     h2tmp->SetBinContent(i,j,V_ij);
	     localvar=fabs(V_ij-V_ij_old);
	     errore_loc = fabs(V_ij - hsa->GetBinContent(i,j));
	     if (errore_loc>errore_max) errore_max = errore_loc;
	     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<<" max Delta (relativa) rispetto a iter k-1 = "
	       <<var<<" ("<<varrel<<"); errore globale = "<< errore_max 
	       <<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,800);
    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();





}