#include <cstdlib>
#include <iostream>
#include <thread>
#include <cmath>
#include <vector>

#include <Garfield/SolidWire.hh>
#include <Garfield/GeometrySimple.hh>
#include <Garfield/ComponentNeBem3d.hh>
#include <Garfield/MediumConductor.hh>
#include <Garfield/MediumMagboltz.hh>

using namespace Garfield;

int main(int argc, char **argv) {
    MediumConductor metal;
    MediumMagboltz gas;

    GeometrySimple geo;
    geo.SetMedium(&gas);

    size_t nWire = 100;
    if (argc == 2) { nWire = std::max(1, std::stoi(argv[1])); } // got tired of compiling over and over again
    std::cout << "Adding " << nWire << " wires" << std::endl;;

    const double radius = 0.0025;
    const double length = 10.;
    const double x0 = 0.1;

    std::vector<SolidWire> wires;    
    wires.reserve(nWire);
    for (size_t i = 0; i < nWire; ++i) { 
        wires.emplace_back(x0 * i, 0.0, 0.0, radius, length * 0.5, 0.0, 1.0, 0.0);
        SolidWire *wire = &wires[i];
        if (!wire) { 
            std::cerr << "Invalid memory address for wire #" << i + 1 << std::endl; 
            return EXIT_FAILURE; 
        }
        wire->SetBoundaryPotential(0.);
        wire->SetLabel("wire" + std::to_string(wire->GetId())); 
        geo.AddSolid(wire, &metal);
    }

    ComponentNeBem3d cmp;
    cmp.EnableDebugging();
    cmp.SetGeometry(&geo);
    cmp.SetTargetElementSize(0.01);
    cmp.SetMinMaxNumberOfElements(1, 3);
    cmp.UseLUInversion();
    // const int njobs = std::max(1, (int) std::thread::hardware_concurrency());
    // std::cout << "Using " << njobs << " threads" << std::endl;
    // cmp.SetNumberOfThreads(njobs);
    if (!cmp.Initialise()) { return EXIT_FAILURE; }

    std::cout << "Exiting..." << std::endl;
    return EXIT_SUCCESS;
}