particles.hpp

#pragma once
 
#include "Eigen/Core"
#include <random>
 
/**
 * @class Particles
 * @brief Stores and advances a dynamically sized set of 1D superparticles.
 *
 * Handles particle emission, leapfrog pushing under electric fields, charge
 * deposition onto grid nodes, and absorbing boundary flux diagnostics.
 */
class Particles {
 
private:
    /**
     * @brief Return index of the cell containing position x.
     * @throws std::out_of_range if x lies outside the domain.
     */
    int find_cell_index(const Eigen::VectorXd& node_pos, double x) const;
    Eigen::VectorXd get_kinetic_energies() const;
 
public:
 
    /** @brief Particle positions (1D). */
    Eigen::VectorXd m_positions;
 
    /** @brief Previous positions for boundary‐crossing detection. */
    Eigen::VectorXd m_pos_prev;
 
    /** @brief Half-step velocities for leapfrog integration. */
    Eigen::VectorXd m_vel_half;
 
    /** @brief Charge of a single physical particle (C). */
    double m_charge;
 
    /** @brief Mass of a single physical particle (kg). */
    double m_mass;
 
    /** @brief Superparticle statistical weight. */
    double m_weight = 1.0;
 
    // ---------------------------------------------------------------------
    // Constructors
    // ---------------------------------------------------------------------
 
    /** @brief Default constructor, creates an empty particle set. */
    Particles() = default;
 
    /**
     * @brief Construct a particle set with N superparticles.
     * @param N Number of particles.
     * @param q Charge per physical particle.
     * @param m Mass per physical particle.
     * @param w Statistical weight per superparticle.
     */
    Particles(int N, double q, double m, double w);
 
    /**
     * @brief Resize particle arrays to contain N superparticles.
     * @param N New total particle count.
     */
    void resize(int N);
 
    double interp_linear(const Eigen::VectorXd& x_nodes, const Eigen::VectorXd& E_nodes, double x);
    // ---------------------------------------------------------------------
    // Particle creation and pushing
    // ---------------------------------------------------------------------
 
    /**
     * @brief Emit N particles at position xpos with Maxwellian velocity distribution.
     * @param N Number of emitted particles.
     * @param xpos Injection position.
     * @param dt Injection time intervall.
     * @param energy_eV Thermal energy (eV) for Maxwellian sampling.
     * @param rng Random generator for velocity/time sampling.
     */
    void emit(    
        int N, double xpos, double dt, double energy_eV, 
        const Eigen::VectorXd& node_pos, const Eigen::VectorXd& efield, std::mt19937& rng);
 
 
    /**
     * @brief Advance particles using leapfrog (velocity half-step scheme).
     * @param dt Time step.
     * @param efield Cell averaged electric fields.
     * @param node_pos Nodal positions.
     */
    void push_leapfrog(double dt, const Eigen::VectorXd& node_pos, const Eigen::VectorXd& efield);
 
    // ---------------------------------------------------------------------
    // Diagnostics and charge deposition
    // ---------------------------------------------------------------------
 
    /**
     * @brief Compute absorbed particle fluxes on four electrodes.
     *
     * Flux order:
     *   0: filament  
     *   1: control electrode  
     *   2: acceleration grid  
     *   3: ion collector  
     *
     * Randomized absorption at the control electrode and acceleration grid
     * is controlled by probabilities p_ce and p_ag.
     *
     * @param el_pos Positions of the four electrodes.
     * @param rng Random generator for absorption events.
     * @param p_ce Absorption probability at control electrode.
     * @param p_ag Absorption probability at acceleration grid.
     * @return Vector of absorbed charges per boundary.
     */
    Eigen::Vector4d fluxes(const Eigen::Vector4d& el_pos,
                           std::mt19937& rng,
                           double p_ce,
                           double p_ag);
 
    /**
     * @brief Deposit particle charge onto grid nodes (CIC scheme).
     * @param node_pos Node positions.
     * @param dx Cell widths.
     * @return Charge density vector (C/m) at nodes.
     */
    Eigen::VectorXd deposit_charge(const Eigen::VectorXd& node_pos,
                                   const Eigen::VectorXd& dx) const;
 
 
    Eigen::VectorXd energy_distribution(int n_bins, double W_max) const;
};