#include <TApplication.h>
 #include <TSystem.h>
 
#include "Garfield/ComponentNeBem3d.hh"
#include "Garfield/GeometrySimple.hh"
#include "Garfield/MediumConductor.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/MediumPlastic.hh"
#include "Garfield/SolidBox.hh"
#include "Garfield/SolidTube.hh"
#include "Garfield/SolidWire.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/ViewGeometry.hh"

using namespace Garfield;

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

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

  MediumMagboltz gas;
  gas.LoadGasFile("../ar_93_co2_7.gas");
  gas.LoadIonMobility("IonMobility_Ar+_Ar.txt");
  gas.SetPressure(760.);
  gas.SetTemperature(293.15);

  MediumConductor Cu;
  MediumPlastic glass;
  glass.SetDielectricConstant(10.0);
  
  // Geometry.
  GeometrySimple geo;

  double half_tck = 0.00015; 

  const double cat_v=-1500.;
  const double sense_v=+10.;
  const double field_v=-1700.;
  const double backplane_v=-1700.;

  const double hx = 0.4;
  const double hy = 0.5;
  SolidBox box1(0, 0, 0.5, hx, hy, 0.0);
  box1.SetBoundaryPotential(cat_v);

  SolidBox box2(0, 0, -0.1 - half_tck, hx, hy, 0.1);
  box2.SetBoundaryDielectric();

  SolidBox box3(0, 0, -0.2 - half_tck, hx, hy, 0.0);
  box3.SetBoundaryPotential(backplane_v);

  geo.AddSolid(&box1, &Cu);
  geo.AddSolid(&box2, &glass);
  geo.AddSolid(&box3, &Cu);

  double passo = 0.4;

  // adding sense strips
  double xS = -hx;
  const double rS = 0.0015;
  std::vector<SolidBox*> stripsS(3, nullptr);
  for (int i = 0; i < 3; ++i) {
    if (i == 0) {
      stripsS[i] = new SolidBox(xS + rS/2, 0.0, 0.0, rS/2, hy, half_tck);
    } else if (i == 2) {
      stripsS[i] = new SolidBox(xS - rS/2, 0.0, 0.0, rS/2, hy, half_tck);
    } else {
      stripsS[i] = new SolidBox(xS, 0.0, 0.0, rS, hy, half_tck);
    }
    stripsS[i]->SetBoundaryPotential(sense_v);
    stripsS[i]->SetLabel("readout");
    geo.AddSolid(stripsS[i], &Cu);
    xS += passo;
  }

  // adding field strips
  double xF = -0.5 * hx;
  const double rF = 0.005;
  std::vector<SolidBox*> stripsF(2, nullptr);
  for (int i = 0; i < 2; ++i) {
    stripsF[i] = new SolidBox(xF, 0.0, 0.0, rF, hy, half_tck); 
    stripsF[i]->SetBoundaryPotential(field_v);
    geo.AddSolid(stripsF[i], &Cu);
    xF += passo;
  }

  geo.PrintSolids();
  geo.SetMedium(&gas);

  ComponentNeBem3d nebem;
  nebem.SetGeometry(&geo);
  nebem.SetTargetElementSize(0.2);
  nebem.SetMinMaxNumberOfElements(1, 100); 
  nebem.SetPeriodicityX(2 * hx); 
  nebem.SetPeriodicityY(2 * hy);
  // nebem.SetPeriodicCopies(3, 5, 0);
  nebem.SetPeriodicCopies(1, 1, 0);
  nebem.Initialise();

  // Plot geometry
  ViewGeometry geomView2d2(&geo);
  TCanvas c2("c2", "", 600, 600);
  geomView2d2.SetCanvas(&c2);
  geomView2d2.SetPlaneXY();
  geomView2d2.Plot3d();
  geomView2d2.PlotPanels();

  ///*
  TCanvas c3("c3", "", 600, 600);
  // Plot field contour
  ViewField fieldView(&nebem);
  fieldView.SetCanvas(&c3);
  fieldView.SetArea(-1., -hx, -0.2, 1., hx, 0.5);
  fieldView.SetPlaneXZ();
  fieldView.PlotContour("v");
 
  TCanvas c4("c4", "", 600, 600);
  ViewField fieldView1(&nebem);
  fieldView1.SetCanvas(&c4);
  // Plot field contour
  fieldView1.SetPlaneXY();
  fieldView1.PlotContour("v");
  //*/
  gSystem->ProcessEvents();
  app.Run(true);
}