#include <iostream>
#include <fstream>

#include <TApplication.h>

#include "Garfield/SolidBox.hh" 
#include "Garfield/SolidHole.hh" 
#include "Garfield/GeometrySimple.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/MediumConductor.hh"
#include "Garfield/MediumPlastic.hh"
#include "Garfield/ComponentNeBem3d.hh"
#include "Garfield/ViewGeometry.hh"
#include "Garfield/ViewField.hh"
#include "Garfield/Medium.hh"

using namespace Garfield;

int main(int argc, char *argv[]) {
  TApplication app("app", &argc, argv);

  // Define gas medium

    MediumMagboltz gas;
    const double gasT = 293.15; //Temperature of gas in Kelvin
    const double gasP = 1 * AtmosphericPressure;
    gas.SetTemperature(gasT);
    gas.SetPressure(gasP);
    gas.LoadGasFile("Ar97p7_iso2p3_1atm.gas");

  MediumConductor metal;
  MediumPlastic dielectric;
  dielectric.SetDielectricConstant(4.0);

  // Define geometry
  GeometrySimple geo;
  geo.SetMedium(&gas);
  
  double drftgap = 1.0;
  double wiredia = 0.004;
  double holedia = 0.006;
  double wirepitch = wiredia + holedia;
  double andplz = 0.0;
  double andplLZ = 0.0010;
  double thickZ1 = 0.0005;
  double thickZ = 0.0045;
  double wirez = 0.00275;
  double wirez2 = 0.00525;
  double lenLX = wirepitch;
  double lenLY = wirepitch;
  double drftplLZ = 0.0010;
  double drftplz = wirez + wiredia / 2 + drftgap + drftplLZ / 2;
  double drftplV = -250.0;
  double andplV = 345.0;
  double wireV1 = 0.0;

  // Drift plane
  SolidBox drift(0.0, 0.0, drftplz, lenLX / 2, lenLY / 2, drftplLZ / 2);
  drift.SetBoundaryPotential(drftplV);
  geo.AddSolid(&drift, &metal);

  // Hole at 0.00525
  SolidHole hole1(0.0, 0.0, wirez2, holedia / 2, holedia / 2, lenLX / 2, lenLY / 2, thickZ1 / 2);
  hole1 .SetSectors(4);
  hole1.SetBoundaryPotential(wireV1);
  geo.AddSolid(&hole1, &metal);

  // Hole at 0.00275 with dielectric
  SolidHole hole2(0.0, 0.0, wirez, holedia / 2, holedia / 2, lenLX / 2, lenLY / 2, thickZ / 2);
  hole2 .SetSectors(4);
  hole2.SetBoundaryDielectric();
  geo.AddSolid(&hole2, &dielectric);

  // Anode plane
  SolidBox anode(0.0, 0.0, andplz, lenLX / 2, lenLY / 2, andplLZ / 2);
  anode.SetBoundaryPotential(andplV);
  geo.AddSolid(&anode, &metal);

  // NeBEM component setup
  
    double tgtElSize = 10.e-4;
  int minEl = 3, maxEl = 5;
  int xcopy = 10, ycopy = 10, zcopy = 0;
  ComponentNeBem3d nebem;
  nebem.SetGeometry(&geo);
  nebem.SetTargetElementSize(tgtElSize);
  nebem.SetMinMaxNumberOfElements(minEl, maxEl);
  nebem.SetPeriodicityX(lenLX);
  nebem.SetPeriodicityY(lenLY);
  nebem.SetPeriodicCopies(xcopy, ycopy, zcopy);
  nebem.UseLUInversion();
  nebem.SetNumberOfThreads(6);
  nebem.Initialise();
  
  ViewGeometry geomView2d;
  geomView2d.SetGeometry(&geo);
  geomView2d.SetArea(-2*lenLX, -2*lenLY, 0.0, 2*lenLX, 2*lenLY, 0.01);
  geomView2d.SetPlane(0, 1, 0, 0, 0, 0.0);
  geomView2d.Plot2d();
  
  
  // Plot device geometry in 3D
  ViewGeometry geomView3d;
  geomView3d.SetGeometry(&geo);
  geomView3d.Plot(); 

 
    ViewField fieldView4;
    fieldView4.SetComponent(&nebem);
    fieldView4.SetPlane(0, -1, 0, 0, 0, 0); // Plane along y = 0(xz plane)
//    fieldView4.SetArea(-0.05, -0.05, 0, 0.05, 0.05, 0.01); // Visualization area
    fieldView4.SetArea(-0.01, -0.01, 0.0, 0.01, 0.01, 0.01); // Visualization area

    fieldView4.SetNumberOfContours(10);
    fieldView4.PlotContour();
    std::vector<double> x0;
    std::vector<double> y0;
    std::vector<double> z0;

    std::vector<double> x1;
    std::vector<double> y1;
    std::vector<double> z1;

    fieldView4.SetElectricFieldRange(10.,300.);
    
    fieldView4.EqualFluxIntervals(-0.008,0.,0.007,0.008,0.,0.007,x1,y1,z1,10);
    fieldView4.PlotFieldLines(x1,y1,z1,true,false);
 

 
 
 
 
 
 
  
 app.Run(true);
}