#ifndef FIELDS_H
#define FIELDS_H

#include "nufi/parameters.h"
#include "poisson_problem.h"
#include <cmath>
#include <deal.II/base/function.h>
#include <deal.II/base/point.h>

using namespace dealii;

// inline std::vector<int> Indices_of_points(const std::vector<double> &points,
// double x_min, double x_max, double dx, int grid_type=0)
// {
//     // grid type:
//     // 0 => uniform
//     // 1 => non uniform (TODO)
//
//     if (dx <= 0.0) {
//         throw std::invalid_argument("dx must be positive");
//     }
//     if (x_max <= x_min) {
//         throw std::invalid_argument("x_max must be > x_min");
//     }
//
//     std::vector<int> indices;
//     indices.reserve(points.size());
//
//     switch (grid_type) {
//         case 0:
//             {
//             const double L = x_max - x_min;
//             const int N = std::floor(L/dx);
//
//
//             for (double x : points) //GPT loop, to check
//             {
//                 x-= x_min;
//                 x = x - L * std::floor(x/L);
//
//                 int i = static_cast<int>(std::floor(x / dx));
//
//                 // safety: handle rare edge case due to floating precision
//                 if (i == N) i = 0;
//
//                 indices.push_back(i);
//             }
//             }
//         case 1:
//             {
//             throw std::invalid_argument("Case for non uniform grid is not
//             completed");
//             }
//         default:
//             throw std::invalid_argument("Invalid grid_type argument");
//
//     }
//     return indices;
// }

inline double f0(const double x, const double v,
                 const double eps = Parameters::EPS,
                 const double k = Parameters::WAVE_NR) {
  const double prefactor =
      Parameters::F0_FACTOR * (1.0 + eps * std::cos(k * x));
  const double gaussian = v * v * std::exp(-0.5 * v * v);

  return prefactor * gaussian;
}

inline double compute_rho(const double x,
                          const unsigned int Nv = Parameters::NV) {
  const double dv =
      (Parameters::V_DOMAIN_RIGHT - Parameters::V_DOMAIN_LEFT) / Nv;

  double integral = 0.0;

  for (unsigned int i = 0; i < Nv; ++i) {
    const double v = Parameters::V_DOMAIN_LEFT + (i + 0.5) * dv;
    integral += f0(x, v) * dv;
  }

  return 1.0 - integral;
}

inline double eval(double x, const PoissonProblem<1> &poisson,
                   const Vector<double> &solution) noexcept {

  x -= Parameters::X_DOMAIN_LEFT;

  x = x - Parameters::LX * std::floor(x * Parameters::LX_INV);
  return eval_point<1>(poisson.get_mapping(), poisson.get_dof_handler(),
                       solution, Point<1>(x));
}

inline double integral_space_vector(const PoissonProblem<1> &poisson,
                                    const Vector<double> &solution,
                                    double dx = Parameters::PLOT_DX,
                                    size_t Nx = Parameters::PLOT_NX) {
  double integral = 0.0;
  double xmin = Parameters::X_DOMAIN_LEFT;
#pragma omp parallel for reduction(+ : integral)
  for (size_t i = 0; i < Nx; ++i) {
    double x = xmin + i * dx;
    integral += eval(x, poisson, solution);
  }
  return integral * dx;
};

inline double integral_space_vector_squared(const PoissonProblem<1> &poisson,
                                            const Vector<double> &solution,
                                            double dx = Parameters::PLOT_DX,
                                            size_t Nx = Parameters::PLOT_NX) {
  double integral = 0.0;
  double xmin = Parameters::X_DOMAIN_LEFT;
#pragma omp parallel for reduction(+ : integral)
  for (size_t i = 0; i < Nx; ++i) {
    double x = xmin + i * dx;
    double val = eval(x, poisson, solution);
    integral += val * val;
  }
  return integral * dx;
};

// class Gradient {
// public:
//   Gradient(double xmin, double xmax, unsigned int Nx)
//       : xmin_(xmin), xmax_(xmax), Nx_(Nx) {
//     if (xmax_ <= xmin_) {
//       throw std::invalid_argument("xmax must be greater than xmin");
//     }
//   }
//
//   std::vector<double> compute(const std::vector<double> &values) const {
//     size_t n = values.size();
//     if (n < 2) {
//       throw std::invalid_argument("Need at least 2 points");
//     }
//
//     std::vector<double> grad(n);
//
//     double dx = (xmax_ - xmin_) / (n - 1);
//     // periodic boundaries
//     grad[0] = -(values[1] - values[n - 1]) / (2.0 * dx);
//     grad[n - 1] = -(values[0] - values[n - 2]) / (2.0 * dx);
//
//     for (size_t i = 1; i < n - 1; ++i) {
//       grad[i] = -(values[i + 1] - values[i - 1]) / (2.0 * dx);
//     }
//
//     return grad;
//   }
//
// private:
//   double xmin_;
//   double xmax_;
//   [[maybe_unused]] unsigned int Nx_;
// };

template <int dim>
class ChargeDensity : public Function<dim> // only uses f0
{
public:
  ChargeDensity(double eps, double k, unsigned int Nv)
      : Function<dim>(1), eps(eps), k(k), Nv(Nv) {}

  virtual double
  value(const Point<dim> &p,
        [[maybe_unused]] const unsigned int component = 0) const override {
    return compute_rho(p[0], Nv);
  }

private:
  const double eps;
  const double k;
  const unsigned int Nv;
};

#endif