#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <map>
#include <regex>
#include <set>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>
#include "TApplication.h"
#include "Garfield/MediumGas.hh"
#include "Garfield/MediumMagboltz.hh"
#include "Garfield/ViewMedium.hh"


namespace {

constexpr double kBoltzmann = 1.380649e-23;          // J/K
constexpr double kTorrToPa  = 133.32236842105263;    // Pa / Torr
constexpr double kPi        = 3.14159265358979323846;

// ---------- utilities ----------

std::string Trim(const std::string& s) {
  size_t i = 0, j = s.size();
  while (i < j && std::isspace(static_cast<unsigned char>(s[i]))) ++i;
  while (j > i && std::isspace(static_cast<unsigned char>(s[j - 1]))) --j;
  return s.substr(i, j - i);
}

std::string Lower(std::string s) {
  for (char& c : s) c = static_cast<char>(std::tolower(static_cast<unsigned char>(c)));
  return s;
}

std::string NormalizeLabel(std::string s) {
  s = Lower(Trim(s));
  s = std::regex_replace(s, std::regex(R"(\s+)"), " ");
  return s;
}

struct TransportRow {
  double en_td;
  double muN;
  double alpha_over_N;
  double eta_over_N;
  double dL_N;
  double dT_N;
};

bool ParseDoubleToken(const std::string& token, double& out) {
  char* end = nullptr;
  out = std::strtod(token.c_str(), &end);
  return end != token.c_str() && *end == '\0';
}

std::vector<double> ExtractNumbers(const std::string& line) {
  std::vector<double> values;
  std::stringstream ss(line);
  std::string tok;
  while (ss >> tok) {
    double x = 0.;
    if (ParseDoubleToken(tok, x)) values.push_back(x);
  }
  return values;
}

bool IsMostlyNumericLine(const std::string& line) {
  auto nums = ExtractNumbers(line);
  if (nums.empty()) return false;
  // Require at least one numeric token and no alphabetic chars.
  for (char c : line) {
    if (std::isalpha(static_cast<unsigned char>(c))) return false;
  }
  return true;
}

double GasNumberDensity(double pressure_torr, double temperature_K) {
  const double p_pa = pressure_torr * kTorrToPa;
  return p_pa / (kBoltzmann * temperature_K);   // m^-3
}

// BOLSIG+ reduced field: E/N in Td, with 1 Td = 1e-21 V m^2.
// Garfield++ grid needs E in V/cm.
double ENTdToEfieldVcm(double en_td, double N) {
  return en_td * 1.e-23 * N;
}

// BOLSIG+ mobility*N is in 1/(m V s).
// v = mu * E = (muN / N) * (E/N * N) = muN * (E/N).
double DriftVelocityCmNs(double muN, double en_td) {
  const double en_si = en_td * 1.e-21;     // V m^2
  const double v_m_s = muN * en_si;        // m/s
  return v_m_s * 1.e-7;                    // cm/ns
}

// alpha/N, eta/N in m^2 -> alpha, eta in cm^-1
double CoeffCmInv(double coeff_over_N, double N) {
  return coeff_over_N * N * 1.e-2;
}

// D*N in 1/(m s) -> D in m^2/s
double ReducedDiffusionToSI(double DN, double N) {
  return DN / N;
}

// Garfield++ wants diffusion coefficients in sqrt(cm). The user guide gives
// dl, dt in sqrt(cm). A consistent conversion from a conventional diffusion
// constant D uses sqrt(2D/v). The units match Garfield's transport interface. :contentReference[oaicite:2]{index=2}
double GarfieldDiffusionSqrtCm(double D_m2_s, double v_cm_ns) {
  if (v_cm_ns <= 0.) return 0.;
  const double D_cm2_ns = D_m2_s * 1.e-5;
  const double coeff_cm = 2.0 * D_cm2_ns / v_cm_ns;
  return coeff_cm > 0. ? std::sqrt(coeff_cm) : 0.;
}

// ---------- block parser ----------
// files are block-style:
//   Label line
//   one or more numeric lines
//   next label line
//
// parse into: label -> concatenated numeric array.
std::map<std::string, std::vector<double>> ParseParameterBlocks(const std::string& filename) {
  std::ifstream fin(filename);
  if (!fin) throw std::runtime_error("Cannot open file: " + filename);

  std::map<std::string, std::vector<double>> blocks;
  std::string line, currentLabel;

  while (std::getline(fin, line)) {
    line = Trim(line);
    if (line.empty()) continue;

    auto nums = ExtractNumbers(line);

    std::stringstream ss(line);
    std::string tok;
    int nTokens = 0;
    int nNumeric = 0;

    while (ss >> tok) {
      ++nTokens;
      double x = 0.;
      if (ParseDoubleToken(tok, x)) ++nNumeric;
    }

    // If almost all tokens are numeric, treat as data.
    if (nTokens > 0 && nNumeric == nTokens) {
      if (!currentLabel.empty()) {
        auto& dest = blocks[currentLabel];
        dest.insert(dest.end(), nums.begin(), nums.end());
      }
      continue;
    }

    // Otherwise treat as a label.
    currentLabel = NormalizeLabel(line);
    if (!blocks.count(currentLabel)) {
      blocks[currentLabel] = {};
    }
  }

  return blocks;
}

const std::vector<double>& GetRequiredBlock(
    const std::map<std::string, std::vector<double>>& blocks,
    const std::vector<std::string>& candidateLabels,
    const std::string& fileTag) {

  for (const auto& label : candidateLabels) {
    auto it = blocks.find(NormalizeLabel(label));
    if (it != blocks.end() && !it->second.empty()) return it->second;
  }

  std::ostringstream os;
  os << "Could not find required block in " << fileTag << ". Tried labels:";
  for (const auto& s : candidateLabels) os << "\n  - " << s;

  os << "\n\nAvailable parsed labels in " << fileTag << ":";
  for (const auto& kv : blocks) {
    os << "\n  * [" << kv.first << "]  (" << kv.second.size() << " values)";
  }

  throw std::runtime_error(os.str());
}


void CheckSameSize(const std::vector<std::pair<std::string, size_t>>& sizes) {
  if (sizes.empty()) return;
  const size_t n = sizes.front().second;
  for (const auto& x : sizes) {
    if (x.second != n) {
      std::ostringstream os;
      os << "Array size mismatch:\n";
      for (const auto& y : sizes) os << "  " << y.first << ": " << y.second << "\n";
      throw std::runtime_error(os.str());
    }
  }
}

void Usage(const char* argv0) {
  std::cerr
      << "Usage:\n"
      << "  " << argv0
      << " output_PT.dat output_SST.dat output_gn_PT.dat output.gas pressure_torr temperature_K\n\n"
      << "Example:\n"
      << "  " << argv0
      << " output_PT.dat output_SST.dat output_gn_PT.dat custom.gas 760 293.15\n";
}

} // namespace

int main(int argc, char* argv[]) {
  if (argc != 7) {
    Usage(argv[0]);
    return 1;
  }

  const std::string filePT   = argv[1];
  const std::string fileSST  = argv[2];
  const std::string fileGNPT = argv[3];
  const std::string outGas   = argv[4];
  const double pressure_torr = std::stod(argv[5]);
  const double temperature_K = std::stod(argv[6]);

  try {
    const auto pt   = ParseParameterBlocks(filePT);
    const auto sst  = ParseParameterBlocks(fileSST);
    const auto gnpt = ParseParameterBlocks(fileGNPT);

    // PT file
    const auto& en_td = GetRequiredBlock(
        pt,
        {
          "Electric field / N (Td)",
          "Electric field / N",
          "Reduced electric field (Td)"
        },
        "output_PT.dat");

    const auto& muN = GetRequiredBlock(
        pt,
        {
          "Mobility *N (1/m/V/s)",
          "Mobility *N"
        },
        "output_PT.dat");

    // SST file
    const auto& alpha_over_N = GetRequiredBlock(
        sst,
        {
          "Townsend ioniz. coef. alpha/N (m2)",
          "Townsend ioniz. coef. alpha/N",
          "Townsend ioniz. coef. alpha/N (m2)",
          "Townsend ionization coefficient alpha/N",
          "Alpha/N (m2)"
        },
        "output_SST.dat");

    const auto& eta_over_N = GetRequiredBlock(
        sst,
        {
          "Townsend attach. coef. eta/N (m2)",
          "Townsend attach. coef. eta/N",
          "Townsend attach. coef. eta/N (m2)",
          "Townsend attachment coefficient eta/N",
          "Eta/N (m2)"
        },
        "output_SST.dat");

    // gn_PT file
    const auto& dL_N = GetRequiredBlock(
        gnpt,
        {
          "Longitudinal diffusion coefficient *N (1/m/s)",
          "Long. diffusion coefficient *N (1/m/s)",
          "Longitud. diffusion coef. *N (1/m/s)",
          "Long. diffusion coefficient *N"
        },
        "output_gnPT.dat");

    const auto& dT_N = GetRequiredBlock(
        gnpt,
        {
          "Transverse diffusion coefficient *N (1/m/s)",
          "Trans. diffusion coefficient *N (1/m/s)",
          "Transverse diffusion coefficient *N",
          "Trans. diffusion coefficient *N",
          "Diffusion coefficient *N (1/m/s)"
        },
        "output_gnPT.dat");
    
    CheckSameSize({
      {"E/N", en_td.size()},
      {"Mobility*N", muN.size()},
      {"alpha/N", alpha_over_N.size()},
      {"eta/N", eta_over_N.size()},
      {"DL*N", dL_N.size()},
      {"DT*N", dT_N.size()}
    });

    const size_t n = en_td.size();
    const double N = GasNumberDensity(pressure_torr, temperature_K);

    std::vector<TransportRow> rows;
    rows.reserve(n);

    for (size_t i = 0; i < n; ++i) {
      rows.push_back({
        en_td[i],
        muN[i],
        alpha_over_N[i],
        eta_over_N[i],
        dL_N[i],
        dT_N[i]
      });
    }

    std::sort(rows.begin(), rows.end(),
              [](const TransportRow& a, const TransportRow& b) {
                return a.en_td < b.en_td;
              });

    std::vector<TransportRow> uniqueRows;
    uniqueRows.reserve(rows.size());

    for (const auto& r : rows) {
      if (uniqueRows.empty() || r.en_td > uniqueRows.back().en_td) {
        uniqueRows.push_back(r);
      } else if (r.en_td == uniqueRows.back().en_td) {
        std::cerr << "Warning: duplicate E/N value removed: "
                  << r.en_td << " Td\n";
      } else {
        std::cerr << "Warning: non-monotonic E/N value skipped: "
                  << r.en_td << " Td\n";
      }
    }

    std::vector<double> efields;
    efields.reserve(uniqueRows.size());
    for (const auto& r : uniqueRows) {
      efields.push_back(ENTdToEfieldVcm(r.en_td, N));
    }

// Garfield related stuff here
	Garfield::MediumMagboltz gas;
	//gas.LoadGasFile("ar_90_ch4_10_test.gas");
	gas.SetComposition("Ar", 90., "CH4", 10.);
    gas.SetPressure(pressure_torr);
    gas.SetTemperature(temperature_K);
    gas.PrintGas();
        
    for (size_t i = 1; i < efields.size(); ++i) {
      if (efields[i] <= efields[i - 1]) {
        std::cerr << "Non-ascending E-field at i=" << i
              << ": " << efields[i - 1] << " then " << efields[i] << "\n";
      }
    }
    gas.SetFieldGrid(efields, std::vector<double>{0.0}, std::vector<double>{0.5 * kPi});

    for (size_t i = 0; i < uniqueRows.size(); ++i) {
      const auto& r = uniqueRows[i];
      const double vE    = DriftVelocityCmNs(r.muN, r.en_td);
      const double vB    = 0.0;
      const double vExB  = 0.0;
      const double DL_SI = ReducedDiffusionToSI(r.dL_N, N);
      const double DT_SI = ReducedDiffusionToSI(r.dT_N, N);
      const double dl    = GarfieldDiffusionSqrtCm(DL_SI, vE);
      const double dt    = GarfieldDiffusionSqrtCm(DT_SI, vE);
      const double alpha = CoeffCmInv(r.alpha_over_N, N);
      const double eta   = CoeffCmInv(r.eta_over_N, N);

      const size_t ib = 0;
      const size_t ia = 0;
      std::cout << "E index: " << i << "\n";
      std::cout << "E field: " << efields[i] << "\n";

      std::cout << "drift: " << vE << "\n";
      std::cout << "alpha: " << alpha << "\n";
      double v = 0;

      if (!gas.SetElectronVelocityE(i, ib, ia, vE)) {
        throw std::runtime_error("SetElectronVelocityE failed at row " + std::to_string(i));
      }
      
      gas.GetElectronVelocityE(i, ib, ia, v);
      std::cout << "test drift: " << v << "\n";

      
      
      if (!gas.SetElectronVelocityB(i, ib, ia, vB)) {
        throw std::runtime_error("SetElectronVelocityB failed at row " + std::to_string(i));
      }


      if (!gas.SetElectronVelocityExB(i, ib, ia, vExB)) {
        throw std::runtime_error("SetElectronVelocityExB failed at row " + std::to_string(i));
      }
      if (!gas.SetElectronLongitudinalDiffusion(i, ib, ia, dl)) {
        throw std::runtime_error("SetElectronLongitudinalDiffusion failed at row " + std::to_string(i));
      }
      if (!gas.SetElectronTransverseDiffusion(i, ib, ia, dt)) {
        throw std::runtime_error("SetElectronTransverseDiffusion failed at row " + std::to_string(i));
      }
      if (!gas.SetElectronTownsend(i, ib, ia, alpha)) {
        throw std::runtime_error("SetElectronTownsend failed at row " + std::to_string(i));
      }

      if (!gas.SetElectronAttachment(i, ib, ia, eta)) {
        throw std::runtime_error("SetElectronAttachment failed at row " + std::to_string(i));
      }
    }
    
    if (!gas.WriteGasFile(outGas)) {
      throw std::runtime_error("WriteGasFile failed for " + outGas);
    }
    gas.PrintGas();
    gas.WriteGasFile(outGas);
    
    double vcheck;
    gas.GetElectronVelocityE(10., 0., 0., vcheck);
    std::cout << "test " << vcheck<<"\n";
    std::cout << "Wrote " << outGas << "\n";
    std::cout << "Rows: " << n << "\n";
    std::cout << "Pressure [Torr]: " << pressure_torr << "\n";
    std::cout << "Temperature [K]: " << temperature_K << "\n";
    std::cout << "Number density [m^-3]: " << std::setprecision(12) << N << "\n";
    std::cout << "E/N range [Td]: " << en_td.front() << " -> " << en_td.back() << "\n";
    std::cout << "E range [V/cm]: " << efields.front() << " -> " << efields.back() << "\n";
    
    Garfield::MediumMagboltz gas1;
	gas1.LoadGasFile("ar_90_ch4_10_test.gas");
    gas1.SetPressure(pressure_torr);
    gas1.SetTemperature(temperature_K);

    
    
    TApplication app("app",&argc,argv);
    Garfield::ViewMedium view;
    TCanvas c1;
    view.SetCanvas(&c1);
    view.SetMedium(&gas);
    view.PlotElectronVelocity('e');
    std::vector<std::string> labels = {"P10 Bolsig+","P10 magboltz"};
    view.SetLabels(labels);
    view.SetMedium(&gas1);
    view.PlotElectronVelocity('e',true);

    app.Run(true);
    
  

  } catch (const std::exception& e) {
    std::cerr << "Error: " << e.what() << "\n";
    return 2;
  }

  return 0;
  
}