#include <cstdlib>
#include <iostream> //std::cout, std::fixed
#include <fstream>

#include <TApplication.h>
#include <TCanvas.h>
#include <TH1F.h>

#include "Garfield/ComponentComsol.hh"
#include "Garfield/ComponentGrid.hh"
#include "Garfield/ComponentAnalyticField.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/ViewFEMesh.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/DriftLineRKF.hh"
#include "Garfield/AvalancheMC.hh"
#include "Garfield/Random.hh"

#include <cmath>
#include <numeric>
#include <string>
#include "Garfield/FundamentalConstants.hh"
#include "Garfield/GarfieldConstants.hh"
#include "Garfield/Numerics.hh"

#include <vector>
#include <iomanip> //std::setprecision

using namespace Garfield;

int main(int argc, char * argv[]) {

  TApplication app("app", &argc, argv);

  // Load the field map.
  ComponentComsol fm;
  //fm.Initialise("mesh.mphtxt", "dielectrics.dat", "field.txt", "mm");
  fm.Initialise("a_mesh.mphtxt", "a_dielectrics.txt", "a_potentialfield.txt", "mm");
  fm.PrintRange();
  
  ViewField fieldView;
  constexpr bool plotField = true;
  if (plotField) {
    fieldView.SetComponent(&fm);
    // Set the normal vector of the viewing plane (xz plane).
    fieldView.SetPlane(0, -1, 0, 0, 0, 0);

    // Set the plot limits in the current viewing plane.
    fieldView.SetArea(-67.5, -67.5, 67.5, 67.5);
    fieldView.SetVoltageRange(0., 2030.);
    TCanvas* cf = new TCanvas("cf", "", 600, 600);
    cf->SetLeftMargin(0.16);
    fieldView.SetCanvas(cf);
    fieldView.PlotContour();
  }
  

  // Setup the gas.
  MediumMagboltz gas;
  gas.SetComposition("ch4", 100.);  
  gas.SetTemperature(293.); // Kelvin
  gas.SetPressure(101.258);
  
  gas.Initialise(true);  
  gas.LoadGasFile("ch4_100.gas");
  
  // Associate the gas with the corresponding field map material. 
  const unsigned int nMaterials = fm.GetNumberOfMaterials();
  for (unsigned int i = 0; i < nMaterials; ++i) {
    const double eps = fm.GetPermittivity(i);
    if (eps == 1.) fm.SetMedium(i, &gas);
  }
  fm.PrintMaterials();

  // Create the sensor.
  Sensor sensor;
  sensor.AddComponent(&fm);
  
  AvalancheMicroscopic aval;
  aval.SetSensor(&sensor);

  double x0 = 0., y0 = 0., z0 = 5, t0 = 0.;   // position: [cm] 
  double e0 = 0.; //Set the initial energy [eV]
  double dx0 = 0., dy0 = 0., dz0 = -1.;
  //Calculate an electron avalanche.
  aval.AvalancheElectron(x0, y0, z0, t0, e0, dx0, dy0, dz0);

  int ne,ni;
  //Get the number of electrons and ions in the avalanche.
  aval.GetAvalancheSize(ne, ni);
  //Get the number of electron drift lines.
  int np = aval.GetNumberOfElectronEndpoints();
  //Initial postion and time
  //double x1, y1, z1, t1;
  //Final position and time
  //double x2, y2, z2, t2;
  //Initial and final energy
  //double e1, e2;
  //Flag indicating why the tracking of an electron was stopped.
  //int status;
  
  std::cout<<"ne, ni, np: "<<ne<<" "<<ni<<" "<<np<<std::endl;

  ViewDrift driftView;
  constexpr bool plotDrift = true;
  if (plotDrift) {
    aval.EnablePlotting(&driftView);
    //drift.EnablePlotting(&driftView);
  }
  if (plotDrift) {
    TCanvas* cd = new TCanvas();
    constexpr bool plotMesh = true;
    if (plotMesh) {
      ViewFEMesh* meshView = new ViewFEMesh();
      //meshView->SetArea(-67.5,-67.5,-67.5,
      //		         67.5, 67.5, 67.5);
      meshView->SetArea(-0.05,-0.02, 0,
                         0.05, 0.02, 1.6);
      //meshView->SetArea(-0.05, 0,
      //		         0.05, 1.6);
      meshView->SetCanvas(cd);
      meshView->SetComponent(&fm);
      // x-z projection.
      meshView->SetPlane(0, -1, 0, 0, 0, 0);
      meshView->SetFillMesh(true);
      meshView->SetColor(0, kYellow + 3);
      meshView->SetColor(1, kGray);
      //meshView->SetColor(2, kGreen);
      meshView->SetColor(3, kRed);
      meshView->EnableAxes();
      meshView->SetViewDrift(&driftView);
      meshView->Plot();
    } else {
      driftView.SetPlane(0, -1, 0, 0, 0, 0);
      //driftView.SetArea(-6.75, -6.75, 6.75, 6.75);
      //driftView.SetArea(-67.5, -67.5, 67.5, 67.5);
      driftView.SetArea(-0.05, 0,
                         0.05, 1.6);
      driftView.SetCanvas(cd);
      constexpr bool twod = true;
      driftView.Plot2d(twod);
    }
  }
  
  

  app.Run(kTRUE);
}