#include <iostream>
#include <fstream>
#include <string>
#include <ctime>
#include <vector>
#include <dirent.h>
#include <math.h>

#include <stdio.h>
#include <unistd.h>

#include "Garfield/ComponentAnsys123.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/ViewMedium.hh"
#include "Garfield/ViewSignal.hh"
#include "Garfield/ViewFEMesh.hh"
#include "Garfield/Sensor.hh"
#include "Garfield/AvalancheMicroscopic.hh"
#include "Garfield/AvalancheMC.hh"
#include "Garfield/Random.hh"
#include "Garfield/TrackHeed.hh"
#include "Garfield/Plotting.hh"


using namespace Garfield;

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

double percent = 70.;
  double ppm = 0.;
  double penning = 0.57;
  const char *noble = "ar";
  double temp = 293.15;
  double press = 760;

  ComponentAnsys123 *fm = new ComponentAnsys123();
  fm->Initialise("ELIST.lis",
                 "NLIST.lis",
                 "MPLIST.lis",
                 "PRNSOL.lis", "mm");
  fm->SetWeightingField("PRNSOL.lis", "WT");
  fm->EnableMirrorPeriodicityX();
  fm->EnableMirrorPeriodicityY();
  fm->PrintRange();

  double co2 = 100. - percent;
  double ar = percent;

MediumMagboltz *gas = new MediumMagboltz();
  gas->LoadGasFile("ar_70_co2_30.gas");

  gas->EnableDrift();
  gas->SetMaxElectronEnergy(18000);
  gas->Initialise();

  if (penning != 0.)
  {
    gas->EnablePenningTransfer(penning, 0., noble);
  }
  else
  {
    gas->DisablePenningTransfer();
  }

  gas->LoadIonMobility("IonMobility_Ar+_Ar.txt");
  int nMaterials = fm->GetNumberOfMaterials();
  for (int i = 0; i < nMaterials; ++i)
  {
    const double eps = fm->GetPermittivity(i);
    if (fabs(eps - 1.) < 1.e-3)
      fm->SetMedium(i, gas);
  }

  fm->PrintMaterials();

  double zmax = 0.3;
  double zmin = -0.1;
  double driftzone = 0.0001;
  int n_events = 1;
  const double t0 = 0., tEnd = 1.E3;
  const int nTimeBins = 100;

  Sensor *sensor = new Sensor();
  sensor->AddComponent(fm);
  sensor->SetArea(-10, -10, zmin, 10, 10, zmax);
  sensor->AddElectrode(fm, "WT");
  sensor->SetTimeWindow(t0, (tEnd - t0) / nTimeBins, nTimeBins);

ViewDrift *vDrift = new ViewDrift();
  vDrift->SetArea(-0.5, -0.5, zmin, 0.5, 0.5, zmax);
  TrackHeed *track = new TrackHeed();
  track->SetSensor(sensor);
  track->SetParticle("mu-"); 
  track->SetEnergy(4E9);     
  track->EnableElectricField();
  track->EnablePlotting(vDrift);

   AvalancheMicroscopic *aval = new AvalancheMicroscopic();
  aval->SetSensor(sensor);
  aval->EnableSignalCalculation();
  aval->SetCollisionSteps(100);
  aval->EnableAvalancheSizeLimit(18000);
  aval->EnablePlotting(vDrift);
 for (int i = n_events; i--;)
  {
    clock_t begin_time = clock();

    double z0 = zmax - driftzone;
    std::cout << "Z0: " << z0 << std::endl;
    double x0 = 0.5 * (RndmUniform() - 1);
    double y0 = 0.5 * (RndmUniform() - 1);
    const double t0 = 0.;
    const double dx0 = 0.;
    const double dy0 = 0.;
    const double dz0 = -1;

    track->NewTrack(x0, y0, z0, t0, dx0, dy0, dz0);

    double xc, yc, zc, tc, extra;
    int nc = 0;
    double ec = 0.;
    while (track->GetCluster(xc, yc, zc, tc, nc, ec, extra))
    {

      for (int j = 0; j < nc; j++)
      {

        double xe1, ye1, ze1, te1, e1, dxe, dye, dze;
        track->GetElectron(j, xe1, ye1, ze1, te1, e1, dxe, dye, dze);

        aval->AvalancheElectron(x0, y0, z0, te1, 0., 0., 0., 0.);

        int ne = 0, ni = 0, nf = 0;
        int np = aval->GetNumberOfElectronEndpoints();
        np += ne;

       
        double xe2, ye2, ze2, te2, e2;
        double xe3, ye3, ze3, te3, e3;
        int status;

        for (int j = np; j--;)
        {
          aval->GetElectronEndpoint(j, xe2, ye2, ze2, te2, e2,
                                    xe3, ye3, ze3, te3, e3, status);

          if (ze3 <= zmin + 0.0001)
          {
            nf += 1;
          }
        }
      }
    }
  }
}