v3.0 - Hagelaar forms

#include <iostream>
#include <random>
#include <fstream>
#include <assert.h>

#include <math.h>
#include <time.h>
#include <iomanip>  // For std::fixed and std::setprecision

#include <algorithm>  // For std::shuffle
#include <numeric>    // For std::iota

//physical constants

#define m_e 9.1093837E-31 // electron mass in kg
#define M_n 6.6464731E-27 // Helium atom mass
#define k_b 1.380649E-23 // Boltzmann constant
#define q 1.602176634E-19 // elementary charge    - eV -> J transfer param
#define Coulomb_log 15.87 // Coulomb logarithm
#define epsilon_0 8.854188E-12 // Vacuum permittivity
#define Coulomb_const pow(q,4)/(pow(4.0*M_PI*epsilon_0,2)) // e^4/(4*pi*eps0)^2
#define thresh0 0.0 //elastic threshold
#define thresh1 19.82 // threshold energy excitation tripet state
#define thresh2 20.61 // threshold energy excitation singlet state
#define tau_singlet 0.0195

// simulation parameters

#define n_e 100000
// #define N_He 1000000 // Helium neutrals number
#define T_n 2.0 // Helium neutral temperature in eV
#define T_e 10.0    // electron Maxwell initial distribution
#define Emin 0.0
#define Emax 3000.0
#define Volume 1.0E-12 // Volume to calculate netral density and collision frequency
#define total_time 1.0E-2 // 500 microsec time to equalibrate the system
#define dopant 1.0E-5 // addition to avoid zero
#define E_reduced 1.0 // constant electrin field along z-axis

#define bin_width 0.05 // keep energy step ~ this to maintain cross-section clarity (Ramsauer minimum etc)
#define N ( (int)((Emax-Emin)/bin_width) + 1) // add 1 to include E_max if needed?

// handling final energy bin

#define bin_width_smooth 0.05 // energy bin for smooth final distribution
#define N_smooth ( (int)((Emax-Emin)/bin_width_smooth) )



double solve_A(double s) { // Netwon method solver

    if (s > 3) {
        return 3*exp(-s);
    }
    if (s < 0.01) {
        return 1.0/s;
    }

    double A0 = 0.01; // initial guess
    double A = A0;  //starting with initial guess
    double tol = 1.0E-7; // accuracy


    for (int i = 0; i < 1000; i++){

        double tanhA = tanh(A);
        // Avoid division by an extremely small tanh(A)
        if (fabs(tanhA) < 1e-12) {
            std::cerr << "tanh(A) too small, returning fallback at iteration " << i << "\n";
            return 1.0E-7;
        }

        double f = 1.0 / tanhA - 1.0 / A - exp(-s);
        if (fabs(f) < tol)
            break;

        double sinhA = sinh(A);
        if (fabs(sinhA) < 1e-12) {
            std::cerr << "sinh(A) too small, returning fallback at iteration " << i << "\n";
            return 1.0E-5;
        }

        double dfdA = -1.0/(sinh(A)*sinh(A)) + 1.0/(A*A);

        // Check if derivative is too close to zero to avoid huge updates
        if (fabs(dfdA) < 1e-12) {
            std::cerr << "dfdA is too small at iteration " << i << ", returning fallback\n";
            if (s < 0.01) {
//                std::cout << "Small s! Huge A!" << "\n";
                return 1.0/s;
            }
            if (s > 3) {
                return 3.0*exp(-s);
            }
        }

        A -= f/dfdA;

        // Early check for numerical issues
        if (std::isnan(A) || std::isinf(A)) {
            std::cerr << "Numerical error detected, invalid A at iteration " << i << "\n";
            return (A > 0) ? 1.0E-5 : -1.0E-5;  // Fallback value based on sign
        }


    }

    return A;
}

struct Electron {

    //velocity components
    double vx = 0.0;
    double vy = 0.0;
    double vz = 0.0;
    //energy in eV
    double energy = 0.0;
    //Collision flag
    bool collided_en = false;
    bool collided_ee = false;

    // initializing Maxwell-Boltzmann distribution with T_e
    void initialize(std::mt19937& gen, std::uniform_real_distribution<double>& dis, std::gamma_distribution<double>& maxwell) {

        double R = dis(gen);

        // velocity angles in spherical coordinates
        double phi = 2*M_PI*dis(gen);
        double cosTheta = 2.0*dis(gen) - 1.0;
        double sinTheta = sqrt(1.0 - cosTheta*cosTheta);


        energy = maxwell(gen); // neutrals energies sampled as Maxwell distribution in eV

        double speed = sqrt(2*energy*q/m_e);

        //velocity components of neutrals in m/s
        vx = speed * sinTheta * cos(phi);
        vy = speed * sinTheta * sin(phi);
        vz = speed * cosTheta;
    }


};

struct CrossSection {
    double energy;
    double sigma;
};

double interpolate (double energy, const std::vector<CrossSection>& CS) {


    if (energy < CS.front().energy) {
//        std::cout << " required energy value lower than range of cross-section data at energy: " << energy << "\n";
        return 0.0;
    }
    if (energy > CS.back().energy) {
//        std::cout << " required energy value higher than range of cross-section data" << "\n";
        return 0.0;
    }

    int step = 0;
        while (step < CS.size() && energy > CS[step].energy) {
            step++;
        }

    double k = (CS[step].sigma - CS[step-1].sigma)/(CS[step].energy - CS[step-1].energy);
    double m = CS[step].sigma - k*CS[step].energy;

    return k*energy + m;
}


struct NeutralParticle {

    double energy = 0.0;
    double vx = 0.0;
    double vy = 0.0;
    double vz = 0.0;

    void initialize(std::mt19937& gen, std::uniform_real_distribution<double>& dis, std::gamma_distribution<double>& maxwell) {

        double R = dis(gen);

        // velocity angles in spherical coordinates
        double phi = 2*M_PI*dis(gen);
        double cosTheta = 2.0*dis(gen) - 1.0;
        double sinTheta = sqrt(1.0 - cosTheta*cosTheta);


        energy = maxwell(gen); // neutrals energies sampled as Maxwell distribution in eV

        double speed = sqrt(2*energy*q/M_n);

        //velocity components of neutrals in m/s
        vx = speed * sinTheta * cos(phi);
        vy = speed * sinTheta * sin(phi);
        vz = speed * cosTheta;
    }

};

struct Excited_neutral {

    double energy;
    double vx;
    double vy;
    double vz;

};



int main() {

    clock_t start = clock();

    int N_He = 100000000;

    std::vector<Electron> electrons(n_e); // better to use vector instead of simple array as it's dynamically allocated (beneficial for ionization)
//    std::vector<NeutralParticle> neutrals(N_He); // I don't need a vector of neutrals bcs it's like a backhround in MCC-simulation

    std::vector<int> histo_random(N, 0); // initialize N size zero-vector for random (initial) histogram
    std::vector<int> histo_maxwell(N, 0); // initialize N size zero-vector for maxwellian histogram
    std::vector<int> histo_neutral(N, 0); // initialize N size zero-vector for neutral distribution histogram
    std::vector<int> histo_excited(N, 0); // initialize N size zero-vector for excited He distribution histogram

    std::vector<double> elastic_vec(N, 0); // precompiled elastic cross-section-energy vector
    std::vector<double> inelastic1_vec(N, 0); // precompiled inelastic(triplet excitation) cross-section-energy vector
    std::vector<double> inelastic2_vec(N, 0); // precompiled inelastic(singlet excitation) cross-section-energy vector
    std::vector<double> superelastic1_vec(N, 0); // precompiled superelastic(triplet de-excitation) cross-section-energy vector
    std::vector<double> superelastic2_vec(N, 0); // precompiled superelastic(triplet de-excitation) cross-section-energy vector

    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<double> dis(0.0, 1.0);
    std::gamma_distribution<double> maxwell_neutral(1.5, T_n);
    std::gamma_distribution<double> maxwell_electron(1.5, T_e);

    std::ifstream elastic_cs_dat("cross_sections/elastic.dat");
    if (!elastic_cs_dat.is_open()) {
        std::cerr << "Error opening elastic cross-sections file!" << std::endl;
        return 1;
    }

    std::ifstream excitation1_cs_dat("cross_sections/inelastic_triplet.dat");
    if (!excitation1_cs_dat.is_open()) {
        std::cerr << "Error opening inelastic triplet cross-sections file!" << std::endl;
        return 1;
    }

    std::ifstream excitation2_cs_dat("cross_sections/inelastic_singlet.dat");
    if (!excitation2_cs_dat.is_open()) {
        std::cerr << "Error opening inelastic singlet cross-sections file!" << std::endl;
        return 1;
    }

    // --- starts reading cross section datafiles

//-----------------elastic---------------------------//
    std::vector<CrossSection> elastic_CS_temp;

    double energy, sigma;

    while (elastic_cs_dat >> energy >> sigma) {
        elastic_CS_temp.push_back({energy, sigma});
    }
    elastic_cs_dat.close();

    energy = 0.0;
    sigma = 0.0;
//-----------------triplet excitation---------------------------//
    std::vector<CrossSection> inelastic1_CS_temp;

    while (excitation1_cs_dat >> energy >> sigma) {
        inelastic1_CS_temp.push_back({energy, sigma});
    }
    excitation1_cs_dat.close();
//-----------------singlet excitation---------------------------//
    energy = 0.0;
    sigma = 0.0;

    std::vector<CrossSection> inelastic2_CS_temp;

    while (excitation2_cs_dat >> energy >> sigma) {
        inelastic2_CS_temp.push_back({energy, sigma});
    }
    excitation2_cs_dat.close();

    // --- finish reading cross-section datafiles

    std::ofstream file0("output_files/velocities.dat");
    std::ofstream file1("output_files/energies.dat");
    std::ofstream file2("output_files/energies_final.dat");
    std::ofstream file3("output_files/histo_random.dat");
    file3 << std::fixed << std::setprecision(10);

    std::ofstream file4("output_files/histo_maxwell.dat");
    file4 << std::fixed << std::setprecision(10);

    std::ofstream file5("output_files/neutral_distribution.dat");
    std::ofstream file6("output_files/E*f(E).dat");
    std::ofstream file7("output_files/nu_max.dat");
    std::ofstream file8("output_files/electron_mean_energy.dat");
    std::ofstream file9("output_files/nu_elastic_average_initial.dat");
    std::ofstream file10("output_files/nu_inelastic1_average_initial.dat");
    std::ofstream file11("output_files/nu_elastic_average_final.dat");
    std::ofstream file12("output_files/nu_inelastic1_average_final.dat");
    std::ofstream file13("output_files/log_output.dat");
    std::ofstream file14("output_files/excited_energies.dat");
    std::ofstream file15("output_files/excited_histo.dat");
    std::ofstream file_temp("output_files/collision_rates.dat");
    std::ofstream file16("output_files/energy_gain.dat");

    // Initialize all electrons
    for (auto& e : electrons) {
        e.initialize(gen, dis, maxwell_electron);
    }

    // precalculate cross-sections for each energy bin
    for (int i = 0; i < N; i++){
        elastic_vec[i] = interpolate(bin_width*(i+0.5), elastic_CS_temp); //elastiuc
        inelastic1_vec[i] = interpolate(bin_width*(i+0.5), inelastic1_CS_temp); //triplet excitation
        inelastic2_vec[i] = interpolate(bin_width*(i+0.5), inelastic2_CS_temp); //singlet excitation
    }

    // precalculate superelastic cross-section (triplet -> ground) for each energy bin
    // detailed balance gives: sigma_j_i(E) = (g_i/g_j)*((E+theshold)/E)*sigma_i_j(E+theshold)
    for (int i = 0; i < N; i++){
        double energy = Emin + (i + 0.5) * bin_width;
        int thresh_bin = (int)( (thresh1 - Emin)/bin_width ); // calculating bin for threshold energy
        superelastic1_vec[i] = (1.0/3.0)*((energy + thresh1)/energy)*interpolate(energy + thresh1, inelastic1_CS_temp); // using detailed balance, calculating backward deexcitation cross-section
        superelastic2_vec[i] = (1.0/1.0)*((energy + thresh2)/energy)*interpolate(energy + thresh2, inelastic2_CS_temp);
    }

    for (int i = 0; i < n_e; i++){
        file1 << i << " " << electrons.at(i).energy << "\n";
        file0 << i << " " << electrons[i].vx << " " << electrons[i].vy << " " << electrons[i].vz << "\n";
    }

    // -----initial electrons energy distribution starts------------////
    for (int i = 0; i < n_e; i++){
        int bin = (int)( (electrons[i].energy - Emin)/bin_width );
        if (bin >=0 && bin < histo_random.size())
            histo_random[bin]++;
    }

    for (int i = 0; i < histo_random.size(); i++){
        double bin_center = Emin + (i + 0.5) * bin_width;
        file3 << bin_center << " " <<  static_cast<double>(histo_random[i])/(electrons.size()*bin_width) << "\n"; // this is electron normalized distribution function
    }



    //calculating excited specied population

    /*From Boltzman distribution y_k = g_k*exp(-E_k/kT)/norm, where g_k - stat weight of k-state,
    E_k - threshold energy for k-state, norm is a total partition function or normaliztion factor     */

    double part_ground = 1.0*exp(-0.0/T_n); // partition function for ground state
    double part_triplet = 3.0*exp(-thresh1/T_n); // partition function for triplet excited state
    double part_singlet = 1.0*exp(-thresh2/T_n); // partition function for singlet excited state
    double part_func_total = part_ground + part_triplet + part_singlet; // total partition function
    double N_trpilet = (part_triplet/part_func_total)*N_He; // population of tripet state
    double N_singlet = (part_singlet/part_func_total)*N_He; // population of singlet state

    std::vector<Excited_neutral> exc_1(static_cast<int>(N_trpilet));  // vector to track triplet excited atoms of Helium
    std::vector<Excited_neutral> exc_2(static_cast<int>(N_singlet));  // vector to track singlet excited atoms of Helium

    // adjusting neutrals number:

    N_He -= (N_trpilet + N_singlet);

    std::cout << N_He << "\n";

    // initializing excited species with Maxwellian distribution

    for (auto& exc : exc_1) {
    NeutralParticle tmp_neutral;
    tmp_neutral.initialize(gen, dis, maxwell_neutral);
    exc.energy = tmp_neutral.energy;
    exc.vx = tmp_neutral.vx;
    exc.vy = tmp_neutral.vy;
    exc.vz = tmp_neutral.vz;
    }

    for (auto& exc : exc_2) {
    NeutralParticle tmp_neutral;
    tmp_neutral.initialize(gen, dis, maxwell_neutral);
    exc.energy = tmp_neutral.energy;
    exc.vx = tmp_neutral.vx;
    exc.vy = tmp_neutral.vy;
    exc.vz = tmp_neutral.vz;
    }

    std::cout << "Triplet population initialized: " << exc_1.size() << "\n";
    std::cout << "Singlet population initialized: " << exc_2.size() << "\n";

    // calculating excited specied population finished //


    //----- calculating number to calculate nu-average (both elastic/inelastic )from our electron distribution starts---------///
    // --- calculating nu(E)*f(E) for later external integration, using initial f(E)
    for (int i = 0; i < N; i++){
        double bin_center = Emin + (i + 0.5) * bin_width;
        file9 << bin_center << " " << (N_He/Volume)*elastic_vec[i] * sqrt(2.0*bin_center*q/m_e)*static_cast<double>(histo_random[i])/(electrons.size()*bin_width) << "\n";
        file10 << bin_center << " " << (N_He/Volume)*inelastic1_vec[i] * sqrt(2.0*bin_center*q/m_e)*static_cast<double>(histo_random[i])/(electrons.size()*bin_width) << "\n";
    }
    //----- calculating nu-average from our electron distribution ends ---------///

    // double dt = 0.1/nu_max;   // minimum should be 0.1/nu_max to get acceptable numerical error range see Vahedi Surrendra 1995
    // double steps = static_cast<int>(time/dt);

    // std::cout << steps << "\n";
    // std::cout << dt << "\n";

    // //using  null-collision technique, getting the number of particles colliding each step: P_collision = 1 - exp(-nu_max*dt)
    // int Ne_collided = (1.0-exp(-1.0*dt*nu_max))*n_e;

    // std::cout << "Ne_collided:" << Ne_collided << "\n";

    int print_interval = 10;
    int el_coll_counter = 0; // track all elastic collisions
    int exc1_coll_counter = 0; // track all triplet excitation collisions
    int exc2_coll_counter = 0; // track all singlet excitation collisions
    int null_coll_counter = 0; // track null-collisions
    int ee_coll_counter = 0; //track e-e Coulomb collisions
    int super1_coll_counter = 0; // track superelastic triplet collisions
    int super2_coll_counter = 0; // track superelastic triplet collisions


    double a_z = ((-1.0)*q * E_reduced) / m_e;
    double mass_ratio = 2.0*(m_e/M_n);
    double charge_mass_ratio = 0.5*m_e/q;
    double sqrt_charge_mass = sqrt(2*q/m_e);
    double C1 = -1.0*q*E_reduced;
    double C2 = 0.5*C1*C1/m_e;

    double t_elapsed = 0.0;

    std::cout << C1 << "    " << C2 << "\n";


    // -----calculating nu-max for null-collision method starts ------------////
    double nu_max = 0.0;
    double nu_max_temp = 0.0;
    double sigma_total = 0.0;

    for (int i = 0; i < N; i++){
        // Get initial densities
        double n_ground = N_He / Volume;
        double n_excited1 = exc_1.size() / Volume;
        double n_excited2 = exc_2.size() / Volume;

        double energy = Emin + (i + 0.5) * bin_width;

        // Total collision frequency for this energy bin
        double sigma_total =
            elastic_vec[i] * n_ground +
            inelastic1_vec[i] * n_ground +
            inelastic2_vec[i] * n_ground +
            superelastic1_vec[i] * n_excited1 +
            superelastic2_vec[i] * n_excited2;

        double v = sqrt(2 * energy * q / m_e);
        double nu_temp = sigma_total * v;

        if (nu_temp > nu_max) nu_max = nu_temp;
    }

    std::cout << "initial nu_max: " <<nu_max << "\n";
    // -----calculating nu-max for null-collision method ends ------------////

    double dt = 0.001/nu_max;   // minimum should be 0.1/nu_max to get acceptable numerical error range see Vahedi Surrendra 1995


    while (t_elapsed < total_time) {
        // Handle edge case for final step
        if (t_elapsed + dt > total_time) {
            dt = total_time - t_elapsed;
        }


        //using  null-collision technique, getting the number of particles colliding each step: P_collision = 1 - exp(-nu_max*dt)
        int Ne_collided = (1.0-exp(-1.0*dt*nu_max))*n_e;

        // Generate shuffled list of electron indices
        int reshuffle_interval = 1;
        std::vector<int> electron_indices(n_e);
        std::iota(electron_indices.begin(), electron_indices.end(), 0); // fill with index


        for (int i = 0; i < Ne_collided; ++i) {
            int j = i + std::uniform_int_distribution<int>(0, n_e - i - 1)(gen);
            std::swap(electron_indices[i], electron_indices[j]);
        }

        int exc1_coll_counter_temp = 0;
        int super1_coll_counter_temp = 0;
        int exc2_coll_counter_temp = 0;
        int super2_coll_counter_temp = 0;
        int null_coll_counter_temp = 0;

        double energy_exc = 0.0; // calculating excitation losses each timestep
        double energy_sup = 0.0; // calculating superelastic gains each timestep
        double energy_Efield = 0.0; // calculating field gains/losses each timestep


        // std::cout << "Progress: " << (t_elapsed/total_time)*100 << "\n";

        // setting flags to false each timestep
        for (auto& e : electrons) e.collided_en = false;
        for (auto& e : electrons) e.collided_ee = false;

        int collision_counter_en = 0; // electron-neutral collision counter
        int collision_counter_ee = 0; // e-e collisoin counter

        /// -- electrin field heating along E-Z axis begin--- /// -- first half!!!
        for (int idx : electron_indices) {
            double half_dt = dt/2.0;
            energy_Efield += ( C1*electrons[idx].vz*half_dt + C2*half_dt*half_dt )/q; // dividing by q to get eV
            // Update velocity component due to electric field
            // double a_z = ((-1.0)*q * E_reduced) / m_e; // acceleration in z-direction, m/s^2
            electrons[idx].vz += a_z * (dt*0.5); // only half timestep

            // Recalculate energy from updated velocity
            double vx = electrons[idx].vx;
            double vy = electrons[idx].vy;
            double vz = electrons[idx].vz;
            electrons[idx].energy = (vx*vx + vy*vy + vz*vz)*charge_mass_ratio;
        }
        // -------------------------------------------- filed heating ends ------------------------//


        for (int idx : electron_indices) {

            if (collision_counter_en >= Ne_collided) break; // quit if reached all collisions

            Electron& e = electrons[idx];
            if (e.collided_en) continue;  // Skip already collided electrons

            double dens_neutrals = (N_He/Volume);
            double dens_exc_1 = (exc_1.size()/Volume);
            double dens_exc_2 = (exc_2.size()/Volume);
            double speed = sqrt_charge_mass*sqrt(e.energy);

            int bin_energy = static_cast<int>(e.energy / bin_width);
            double nu_elastic = dens_neutrals * elastic_vec[bin_energy] * speed;
            double nu_inelastic1 = dens_neutrals * inelastic1_vec[bin_energy] * speed;
            double nu_superelastic1 = dens_exc_1 * superelastic1_vec[bin_energy] * speed;
            double nu_inelastic2 = dens_neutrals * inelastic2_vec[bin_energy] * speed;
            double nu_superelastic2 = dens_exc_2 * superelastic2_vec[bin_energy] * speed;

            double r = dis(gen);

            double P0 = nu_elastic/nu_max;
            double P1 = (nu_elastic + nu_inelastic1)/nu_max;
            double P2 = (nu_elastic + nu_inelastic1 + nu_superelastic1)/nu_max;
            double P3 = (nu_elastic + nu_inelastic1 + nu_superelastic1 + nu_inelastic2)/nu_max;
            double P4 = (nu_elastic + nu_inelastic1 + nu_superelastic1 + nu_inelastic2 + nu_superelastic2)/nu_max;

            if (r < P0) {

                // elastic collision happens

                // ----   Collision energy redistribution module

                // electron particle X Y Z initial velocities and energy
                double V0_x = e.vx;
                double V0_y = e.vy;
                double V0_z = e.vz;
                double E_0 = e.energy;

                // neutral that collides with electron

                // randomize particles each collision

                NeutralParticle tmp_neutral;
                tmp_neutral.initialize(gen, dis, maxwell_neutral);
                double V_x_n = tmp_neutral.vx;
                double V_y_n = tmp_neutral.vy;
                double V_z_n = tmp_neutral.vz;
                double E_n = tmp_neutral.energy;


                double V0 = sqrt(V0_x*V0_x + V0_y*V0_y + V0_z*V0_z);

                // initial relative velocity
                double V_rel_x = V0_x - V_x_n;
                double V_rel_y = V0_y - V_y_n;
                double V_rel_z = V0_z - V_z_n;

                double V_rel = sqrt(V_rel_x*V_rel_x + V_rel_y*V_rel_y + V_rel_z*V_rel_z);
                if (std::isnan(V_rel) || std::isinf(V_rel) || fabs(V_rel) < 1e-12) {
                    std::cerr << "Invalid V_rel copmuted: " << V_rel << " at timestep " << t_elapsed << std::endl;
                }
                //initial relative velocity unit vectors components
                double e1 = V_rel_x/V_rel;
                double e2 = V_rel_y/V_rel;
                double e3 = V_rel_z/V_rel;

                double C = -V_rel/(1.0+(m_e/M_n));
                double beta = sqrt(1.0 - thresh0/e.energy);
                // generating random variables to calculate random direction of center-of-mass after the collision

                double R1 = dis(gen);
                double R2 = dis(gen);

                //// calculating spherical angles for center-of-mass random direction
                // double theta = acos(1.0- 2.0*R1);
                // double phi = 2*M_PI*R2;
                double cos_khi = 1.0 - 2.0*R1;
                double sin_khi = sqrt(1.0 - cos_khi*cos_khi);
                if (std::isnan(sin_khi) || std::isinf(sin_khi)) {
                    std::cerr << "Invalid sin_khi copmuted: " << sin_khi << " at timestep " << t_elapsed << std::endl;
                }

                double phi = 2.0*M_PI*R2;

                //calculating velocity change of electron

                double dv_x = C*( (1.0-beta*cos_khi)*e1 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(-e2) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e1*e3)  );
                double dv_y = C*( (1.0-beta*cos_khi)*e2 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(e1) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e2*e3)  );
                double dv_z = C*( (1.0-beta*cos_khi)*e3 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(0.0) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(e1*e1 + e2*e2)  );

                //updating electron energy and velocities

                //updating electron energy and velocities

                e.vx += dv_x;
                e.vy += dv_y;
                e.vz += dv_z;

                e.energy = (e.vx*e.vx + e.vy*e.vy + e.vz*e.vz) * charge_mass_ratio;

                collision_counter_en++;
                el_coll_counter++;

                e.collided_en = true;
            }

            else if (r < P1) {

                //inelastic 1(triplet) collision happens

                // ----   Collision energy redistribution module

                // electron particle X Y Z initial velocities and energy
                double V0_x = e.vx;
                double V0_y = e.vy;
                double V0_z = e.vz;
                double E_0 = e.energy;

                // neutral that collides with electron

                // randomize particles each collision

                NeutralParticle tmp_neutral;
                tmp_neutral.initialize(gen, dis, maxwell_neutral);
                double V_x_n = tmp_neutral.vx;
                double V_y_n = tmp_neutral.vy;
                double V_z_n = tmp_neutral.vz;
                double E_n = tmp_neutral.energy;


                double V0 = sqrt(V0_x*V0_x + V0_y*V0_y + V0_z*V0_z);

                // initial relative velocity
                double V_rel_x = V0_x - V_x_n;
                double V_rel_y = V0_y - V_y_n;
                double V_rel_z = V0_z - V_z_n;

                double V_rel = sqrt(V_rel_x*V_rel_x + V_rel_y*V_rel_y + V_rel_z*V_rel_z);

                if (std::isnan(V_rel) || std::isinf(V_rel) || fabs(V_rel) < 1e-12) {
                    std::cerr << "Invalid V_rel copmuted: " << V_rel << " at timestep " << t_elapsed << std::endl;
                }

                //initial relative velocity unit vectors components
                double e1 = V_rel_x/V_rel;
                double e2 = V_rel_y/V_rel;
                double e3 = V_rel_z/V_rel;

                double C = -V_rel/(1.0+(m_e/M_n));
                double beta = sqrt(abs(1.0 - thresh1/e.energy));
                // generating random variables to calculate random direction of center-of-mass after the collision

                double R1 = dis(gen);
                double R2 = dis(gen);

                //// calculating spherical angles for center-of-mass random direction
                // double theta = acos(1.0- 2.0*R1);
                // double phi = 2*M_PI*R2;
                double cos_khi = 1.0 - 2.0*R1;
                double sin_khi = sqrt(1.0 - cos_khi*cos_khi);
                if (std::isnan(sin_khi) || std::isinf(sin_khi)) {
                    std::cerr << "Invalid sin_khi copmuted: " << sin_khi << " at timestep " << t_elapsed << std::endl;
                }
                double phi = 2.0*M_PI*R2;

                //calculating velocity change of electron

                double dv_x = C*( (1.0-beta*cos_khi)*e1 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(-e2) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e1*e3)  );
                double dv_y = C*( (1.0-beta*cos_khi)*e2 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(e1) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e2*e3)  );
                double dv_z = C*( (1.0-beta*cos_khi)*e3 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(0.0) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(e1*e1 + e2*e2)  );

                //updating electron energy and velocities

                if (e.energy < thresh1) {
                    null_coll_counter++;
                    collision_counter_en++;
                    e.collided_en = true;
                    continue;
                }
                else {
                    e.vx += dv_x;
                    e.vy += dv_y;
                    e.vz += dv_z;

                    e.energy = (e.vx*e.vx + e.vy*e.vy + e.vz*e.vz) * charge_mass_ratio;

                    collision_counter_en++;
                    exc1_coll_counter++;
                    exc1_coll_counter_temp++;

                    e.collided_en = true;

                    // pushing this neutral to an array of excited species exc_1
                    if (N_He > 0) {
                        exc_1.push_back({E_n, V_x_n, V_y_n, V_z_n});
                        N_He--;
                    }
                }
            }

            else if (r < P2) {

                //superelastic 1(triplet -> ground state) collision happens

                if (exc_1.empty()) {
                    null_coll_counter++;
                    collision_counter_en++;
                    e.collided_en = true;
                    continue;
                }

                // ----   Collision energy redistribution module

                // electron particle X Y Z initial velocities and energy
                double V0_x = e.vx;
                double V0_y = e.vy;
                double V0_z = e.vz;
                double E_0 = e.energy;

                // neutral that collides with electron
                // taking particles from dynamic array of excited neutrals

                int index = std::uniform_int_distribution<int>(0, exc_1.size()-1)(gen);
                Excited_neutral& exc = exc_1[index];
                double V_x_n = exc.vx;
                double V_y_n = exc.vy;
                double V_z_n = exc.vz;
                double E_n = exc.energy;

                double V0 = sqrt(V0_x*V0_x + V0_y*V0_y + V0_z*V0_z);

                // initial relative velocity
                double V_rel_x = V0_x - V_x_n;
                double V_rel_y = V0_y - V_y_n;
                double V_rel_z = V0_z - V_z_n;

                double V_rel = sqrt(V_rel_x*V_rel_x + V_rel_y*V_rel_y + V_rel_z*V_rel_z);

                if (std::isnan(V_rel) || std::isinf(V_rel) || fabs(V_rel) < 1e-12) {
                    std::cerr << "Invalid V_rel copmuted: " << V_rel << " at timestep " << t_elapsed << std::endl;
                }

                //initial relative velocity unit vectors components
                double e1 = V_rel_x/V_rel;
                double e2 = V_rel_y/V_rel;
                double e3 = V_rel_z/V_rel;

                double C = -V_rel/(1.0+(m_e/M_n));
                double beta = sqrt(1.0 + thresh1/e.energy);
                // generating random variables to calculate random direction of center-of-mass after the collision

                double R1 = dis(gen);
                double R2 = dis(gen);

                //// calculating spherical angles for center-of-mass random direction
                // double theta = acos(1.0- 2.0*R1);
                // double phi = 2*M_PI*R2;
                double cos_khi = 1.0 - 2.0*R1;
                double sin_khi = sqrt(1.0 - cos_khi*cos_khi);
                if (std::isnan(sin_khi) || std::isinf(sin_khi)) {
                    std::cerr << "Invalid sin_khi copmuted: " << sin_khi << " at timestep " << t_elapsed << std::endl;
                }
                double phi = 2.0*M_PI*R2;

                //calculating velocity change of electron

                double dv_x = C*( (1.0-beta*cos_khi)*e1 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(-e2) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e1*e3)  );
                double dv_y = C*( (1.0-beta*cos_khi)*e2 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(e1) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e2*e3)  );
                double dv_z = C*( (1.0-beta*cos_khi)*e3 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(0.0) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(e1*e1 + e2*e2)  );

                //updating electron energy and velocities

                e.vx += dv_x;
                e.vy += dv_y;
                e.vz += dv_z;

                e.energy = (e.vx*e.vx + e.vy*e.vy + e.vz*e.vz) * charge_mass_ratio;

                //counting collisions, working with flags, popping atom out of the vector
                if (!exc_1.empty()) {
                    std::swap(exc_1[index], exc_1.back());
                    exc_1.pop_back();
                    N_He++;
                }
                collision_counter_en++;
                super1_coll_counter++;
                super1_coll_counter_temp++;
                energy_sup += thresh1;

                e.collided_en = true;
            }

            else if (r < P3) {

                //inelastic 1(singlet) excitation collision happens

                // ----   Collision energy redistribution module

                // electron particle X Y Z initial velocities and energy
                double V0_x = e.vx;
                double V0_y = e.vy;
                double V0_z = e.vz;
                double E_0 = e.energy;

                // neutral that collides with electron

                // randomize particles each collision

                NeutralParticle tmp_neutral;
                tmp_neutral.initialize(gen, dis, maxwell_neutral);
                double V_x_n = tmp_neutral.vx;
                double V_y_n = tmp_neutral.vy;
                double V_z_n = tmp_neutral.vz;
                double E_n = tmp_neutral.energy;


                double V0 = sqrt(V0_x*V0_x + V0_y*V0_y + V0_z*V0_z);

                // initial relative velocity
                double V_rel_x = V0_x - V_x_n;
                double V_rel_y = V0_y - V_y_n;
                double V_rel_z = V0_z - V_z_n;

                double V_rel = sqrt(V_rel_x*V_rel_x + V_rel_y*V_rel_y + V_rel_z*V_rel_z);

                if (std::isnan(V_rel) || std::isinf(V_rel) || fabs(V_rel) < 1e-12) {
                    std::cerr << "Invalid V_rel copmuted: " << V_rel << " at timestep " << t_elapsed << std::endl;
                }

                //initial relative velocity unit vectors components
                double e1 = V_rel_x/V_rel;
                double e2 = V_rel_y/V_rel;
                double e3 = V_rel_z/V_rel;

                double C = -V_rel/(1.0+(m_e/M_n));
                double beta = sqrt(abs(1.0 - thresh2/e.energy));
                // generating random variables to calculate random direction of center-of-mass after the collision

                double R1 = dis(gen);
                double R2 = dis(gen);

                //// calculating spherical angles for center-of-mass random direction
                // double theta = acos(1.0- 2.0*R1);
                // double phi = 2*M_PI*R2;
                double cos_khi = 1.0 - 2.0*R1;
                double sin_khi = sqrt(1.0 - cos_khi*cos_khi);
                if (std::isnan(sin_khi) || std::isinf(sin_khi)) {
                    std::cerr << "Invalid sin_khi copmuted: " << sin_khi << " at timestep " << t_elapsed << std::endl;
                }
                double phi = 2.0*M_PI*R2;

                //calculating velocity change of electron

                double dv_x = C*( (1.0-beta*cos_khi)*e1 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(-e2) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e1*e3)  );
                double dv_y = C*( (1.0-beta*cos_khi)*e2 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(e1) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e2*e3)  );
                double dv_z = C*( (1.0-beta*cos_khi)*e3 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(0.0) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(e1*e1 + e2*e2)  );

                //updating electron energy and velocities

                if (e.energy < thresh1) {
                    null_coll_counter++;
                    collision_counter_en++;
                    e.collided_en = true;
                    continue;
                }
                else {
                    e.vx += dv_x;
                    e.vy += dv_y;
                    e.vz += dv_z;

                    e.energy = (e.vx*e.vx + e.vy*e.vy + e.vz*e.vz) * charge_mass_ratio;

                    collision_counter_en++;
                    exc2_coll_counter++;
                    exc2_coll_counter_temp++;

                    e.collided_en = true;

                    // pushing this neutral to an array of excited species exc_2

                    if (N_He > 0) {
                        exc_2.push_back({E_n, V_x_n, V_y_n, V_z_n});
                        N_He--;
                    }
                }
            }

            else if (r < P4) {

                //supernelastic 1(singlet -> ground state) collision happens

                if (exc_2.empty()) {
                    null_coll_counter++;
                    collision_counter_en++;
                    e.collided_en = true;
                    continue;
                }

                // ----   Collision energy redistribution module

                // electron particle X Y Z initial velocities and energy
                double V0_x = e.vx;
                double V0_y = e.vy;
                double V0_z = e.vz;
                double E_0 = e.energy;

                // neutral that collides with electron
                // taking particles from dynamic array of excited neutrals

                int index = std::uniform_int_distribution<int>(0, exc_2.size()-1)(gen);
                Excited_neutral& exc = exc_2[index];
                double V_x_n = exc.vx;
                double V_y_n = exc.vy;
                double V_z_n = exc.vz;
                double E_n = exc.energy;

                double V0 = sqrt(V0_x*V0_x + V0_y*V0_y + V0_z*V0_z);

                // initial relative velocity
                double V_rel_x = V0_x - V_x_n;
                double V_rel_y = V0_y - V_y_n;
                double V_rel_z = V0_z - V_z_n;

                double V_rel = sqrt(V_rel_x*V_rel_x + V_rel_y*V_rel_y + V_rel_z*V_rel_z);

                if (std::isnan(V_rel) || std::isinf(V_rel) || fabs(V_rel) < 1e-12) {
                    std::cerr << "Invalid V_rel copmuted: " << V_rel << " at timestep " << t_elapsed << std::endl;
                }

                //initial relative velocity unit vectors components
                double e1 = V_rel_x/V_rel;
                double e2 = V_rel_y/V_rel;
                double e3 = V_rel_z/V_rel;

                double C = -V_rel/(1.0+(m_e/M_n));
                double beta = sqrt(1.0 + thresh2/e.energy);
                // generating random variables to calculate random direction of center-of-mass after the collision

                double R1 = dis(gen);
                double R2 = dis(gen);

                //// calculating spherical angles for center-of-mass random direction
                // double theta = acos(1.0- 2.0*R1);
                // double phi = 2*M_PI*R2;
                double cos_khi = 1.0 - 2.0*R1;
                double sin_khi = sqrt(1.0 - cos_khi*cos_khi);
                if (std::isnan(sin_khi) || std::isinf(sin_khi)) {
                    std::cerr << "Invalid sin_khi copmuted: " << sin_khi << " at timestep " << t_elapsed << std::endl;
                }
                double phi = 2.0*M_PI*R2;

                //calculating velocity change of electron

                double dv_x = C*( (1.0-beta*cos_khi)*e1 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(-e2) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e1*e3)  );
                double dv_y = C*( (1.0-beta*cos_khi)*e2 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(e1) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e2*e3)  );
                double dv_z = C*( (1.0-beta*cos_khi)*e3 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(0.0) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(e1*e1 + e2*e2)  );

                //updating electron energy and velocities

                e.vx += dv_x;
                e.vy += dv_y;
                e.vz += dv_z;

                e.energy = (e.vx*e.vx + e.vy*e.vy + e.vz*e.vz) * charge_mass_ratio;
                //counting collisions, working with flags, popping atom out of the vector

                if (!exc_2.empty()) {
                    std::swap(exc_2[index], exc_2.back());
                    exc_2.pop_back();
                    N_He++;
                }

                collision_counter_en++;
                super2_coll_counter++;
                super2_coll_counter_temp++;
                energy_sup += thresh2;

                e.collided_en = true;
            }

            else {
                // null-collision
                collision_counter_en++;
                null_coll_counter++;
                e.collided_en = true;
            }
        }


        /// -- electrin field heating along E-Z axis begin--- /// -- second half!!!
        for (int idx : electron_indices) {
            double half_dt = dt/2.0;
            energy_Efield += ( C1*electrons[idx].vz*half_dt + C2*half_dt*half_dt )/q; //dividing by q to get eV
            // Update velocity component due to electric field
            // double a_z = ((-1.0)*q * E_reduced) / m_e; // acceleration in z-direction, m/s^2
            electrons[idx].vz += a_z * (dt*0.5); // only half timestep

            // Recalculate energy from updated velocity
            double vx = electrons[idx].vx;
            double vy = electrons[idx].vy;
            double vz = electrons[idx].vz;
            electrons[idx].energy = (vx*vx + vy*vy + vz*vz) * charge_mass_ratio;
        }
        // -------------------------------------------- filed heating ends ------------------------////////////////

        int decay_counter = 0;

        // // Iterate backwards to safely remove elements
        // for (int i = exc_2.size() - 1; i >= 0; --i) {
        //     if (dis(gen) < dt / tau_singlet) {
        //         // Swap with last element and pop (like your superelastic code)
        //         std::swap(exc_2[i], exc_2.back());
        //         exc_2.pop_back();
        //         N_He++;
        //         decay_counter++;
        //     }
        // }


        t_elapsed += dt; // Advance time

        // Recalculate nu_max periodically (e.g., every 100 steps)
        static int recalc_counter = 0;
        if (++recalc_counter >= 1000) {

            recalc_counter = 0;

            // Recalculate nu_max with CURRENT densities
            nu_max = 0.0;
            for (int i = 0; i < N; i++) {
                double energy = Emin + (i + 0.5) * bin_width;

                // Get current densities
                double n_ground = N_He / Volume;
                double n_excited1 = exc_1.size() / Volume;
                double n_excited2 = exc_2.size() / Volume;

                // Total collision frequency for this energy bin
                double sigma_total =
                    elastic_vec[i] * n_ground +
                    inelastic1_vec[i] * n_ground +
                    inelastic2_vec[i] * n_ground +
                    superelastic1_vec[i] * n_excited1 +
                    superelastic2_vec[i] * n_excited2;

                double speed = sqrt_charge_mass*sqrt(energy);
                double nu_temp = sigma_total * speed;

                if (nu_temp > nu_max) nu_max = nu_temp;
            }

            // Update dt based on new nu_max
            dt = 0.001 / nu_max;
        }

        // calculating mean energy
        if (static_cast<int>(t_elapsed/dt)%print_interval == 0) {
            double total_energy = 0.0;
            for (const auto& e : electrons) total_energy += e.energy;
            double mean_energy = total_energy / n_e;
            file8 << t_elapsed << " " << mean_energy << "\n";
            file_temp << t_elapsed << " " << exc_1.size() << " " << exc_2.size() << "\n";
            std::cout << "Progress: " << (t_elapsed/total_time)*100 << "%" << " ";
            std::cout << "   nu_max: " << nu_max << "    " << "dt: " << dt << " " << "decay counter: " << decay_counter <<   "\n";
            file16 << t_elapsed << " " << energy_Efield/n_e << " " << energy_sup/n_e << "\n";
        }

    }

    // ----- final electron energies distribution begins
    for (int i = 0; i < n_e; i++){

        file2 << i << " " << electrons[i].energy << "\n";

        int bin = static_cast<int>( (electrons[i].energy - Emin)/bin_width_smooth);
        if (bin >=0 && bin < histo_maxwell.size())
            histo_maxwell[bin]++;
    }

    int check = 0;
    for (int i = 0; i < N_smooth; i++){
        check += histo_maxwell[i];
        double bin_center = Emin + (i + 0.5) * bin_width_smooth;
        file4 << bin_center << " " <<  static_cast<double>(histo_maxwell[i])/(electrons.size()*bin_width_smooth) << "\n"; // getting f(E)
    }

        std::cout << "Total # of electrons in a final histogram: " << check << "\n";
        std::cout << "Final nu max: " << nu_max << "\n";

    // ----- final electron energies distribution ends


    file0.close();
    file1.close();
    file2.close();
    file3.close();
    file4.close();
    file5.close();
    file6.close();
    file7.close();
    file8.close();
    file9.close();
    file10.close();
    file11.close();
    file12.close();
    file13.close();
    file14.close();
    file15.close();
    file_temp.close();
    file16.close();

    clock_t end = clock();

    double elapsed = (double)(end - start) / CLOCKS_PER_SEC;

    // std::cout << "# of steps: " << steps << "\n";
    // std::cout << "# of electrons collided each timesteps:" << Ne_collided << "\n";

    // std::cout << "Average elastic collisions per timestep: " << static_cast<int>(el_coll_counter/steps) << "\n";
    // std::cout << "Average null collisions per timestep: " << static_cast<int>(null_coll_counter/steps) << "\n";
    // std::cout << "\n";

    // std::cout << "triplet:________" << "\n";
    // std::cout << "Average triplet excitation collisions per timestep: " << static_cast<int>(exc1_coll_counter/steps) << "\n";
    // std::cout << "\n";
    // std::cout << "Average superelastic triplet collisions per timestep: " << static_cast<int>(super1_coll_counter/steps) << "\n";
    // std::cout << "\n";

    // std::cout << "singlet:________" << "\n";
    // std::cout << "Average singlet excitation collisions per timestep: " << static_cast<int>(exc2_coll_counter/steps) << "\n";
    // std::cout << "\n";
    // std::cout << "Average superelastic singlet collisions per timestep: " << static_cast<int>(super2_coll_counter/steps) << "\n";
    // std::cout << "\n";

    // std::cout << "Average e-e collisions per timestep: " << static_cast<int>(ee_coll_counter/steps) << "\n";

    std::cout << "Elapsed time: %f seconds " << elapsed << "\n";


    return 0;

}