#include <cstdlib>
#include <iostream>

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

#include <Garfield/MediumMagboltz.hh>
#include <Garfield/MediumConductor.hh>
#include <Garfield/GeometrySimple.hh>
#include <Garfield/SolidSphere.hh>
#include <Garfield/SolidWire.hh>
#include <Garfield/ComponentNeBem3d.hh>
#include <Garfield/Sensor.hh>
#include <Garfield/TrackHeed.hh>
#include <Garfield/ViewField.hh>
#include <Garfield/ViewGeometry.hh>
#include <Garfield/ViewDrift.hh>

using namespace Garfield;

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

    // Spherical shell
    constexpr double innerR = 2.0; // [cm]
    constexpr double outerR = innerR + 0.01; // [cm]
    constexpr double sphereV = 0.0; // [V]
    SolidSphere sphere(0.0, 0.0, 0.0, // x0, y0, z0 (centroid)
                       innerR, outerR // inner radius, outer radius
                      );
    sphere.SetBoundaryPotential(sphereV);
    // sphere.SetMeridians(6); // UV sphere approximation

    // Solid wire
    constexpr double wireR = 0.001; // [cm]
    constexpr double wireV = 1e3; // [V]
    SolidWire wire(0.0, 0.0, 0.0, // x0, y0, z0 (centroid)
                   wireR, outerR * 1.2, // radius, half-length
                   0.0, 0.0, 1.0 // dx, dy, dz (orientation)
                   );
    wire.SetBoundaryPotential(wireV);

    MediumMagboltz gas;
    MediumConductor metal;

    GeometrySimple geo;
    geo.SetMedium(&gas);
    geo.AddSolid(&sphere, &metal);
    geo.AddSolid(&wire, &metal);

    ComponentNeBem3d cmp;
    cmp.SetGeometry(&geo);
    cmp.UseLUInversion();
    cmp.SetNumberOfThreads(36);
    cmp.SetTargetElementSize(1);
    cmp.SetMinMaxNumberOfElements(1, 5); 
    cmp.Initialise();

    Sensor sensor(&cmp);
    sensor.SetArea(-outerR, -outerR, -outerR,
                    outerR,  outerR,  outerR);

    ViewDrift driftView;

    TrackHeed track(&sensor);
    track.EnablePlotting(&driftView);
    track.NewTrack(-innerR + 1e-4, 0.0, 0.0, 0.0, 
                    1.732 * 0.5, 0.5, 0.0); // plot a particle track

    // Plotting limits, x-axis
    constexpr double xlim0 = -outerR;
    constexpr double xlim1 =  outerR;
    // Plotting limits, y-axis
    constexpr double ylim0 = -outerR;
    constexpr double ylim1 =  outerR;
    // Plotting limits, z-axis
    constexpr double zlim0 = -outerR;
    constexpr double zlim1 =  outerR;

    // YZ view of the electric potential field
    ViewField fieldViewYZ(&cmp);
    fieldViewYZ.SetArea(ylim0, zlim0,
                        ylim1, zlim1);
    fieldViewYZ.SetPlaneYZ();
    fieldViewYZ.PlotContour("v");
    fieldViewYZ.GetCanvas()->SetTitle("Electric potential contour (Y-Z view)");
    fieldViewYZ.GetCanvas()->Modified();
    fieldViewYZ.GetCanvas()->Update();

    // XY view of the electric potential field
    ViewField fieldViewXY(&cmp);
    fieldViewXY.SetArea(xlim0, ylim0,
                        xlim1, ylim1);
    fieldViewXY.SetPlaneXY();
    fieldViewXY.PlotContour("v");
    fieldViewXY.GetCanvas()->SetTitle("Electric potential contour (X-Y view)");
    fieldViewXY.GetCanvas()->Modified();
    fieldViewXY.GetCanvas()->Update();

    // 3D view of the geometry
    ViewGeometry geoView(&geo);
    geoView.Plot3d();
    driftView.SetCanvas(geoView.GetCanvas());
    driftView.Plot3d();
    geoView.GetCanvas()->SetTitle("Geometry view");
    geoView.GetCanvas()->Modified();
    geoView.GetCanvas()->Update();

    std::cout << "Done." << std::endl;

    app.Run(true);

    return EXIT_SUCCESS;
}