e-n collisions(v10) e-e collisions added

#include <iostream>
#include <random>
#include <fstream>

#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.1E-31 // electron mass in kg
#define M_n 6.6464731E-27 // Helium atom mass
#define k_b 1.38E-23 // Boltzmann constant
#define q 1.602176634E-19 // elementary charge    - eV -> J transfer param
#define Coulomb_log 10 // 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

// simulation parameters

#define n_e 500000
#define N_He 1000000 // Helium neutrals number
#define T_n 10.0 // Helium neutral temperature in eV
#define T_e 50.0    // electron Maxwell initial distribution
#define Emin 0.0
#define Emax 450.0
#define Volume 1.0E-12 // Volume to calculate netral density and collision frequency
#define time 6.0E-4 // 500 microsec time to equalibrate the system

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

// handling final energy bin

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

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;
    }

};




int main() {

    clock_t start = clock();

    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);


    std::vector<int> histo_random(N, 0); // initialize N size zero-vector for random (initial) histogram
    std::vector<int> histo_maxwell(N_smooth, 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<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::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::uniform_int_distribution<int> pair(0, n_e-1);
    std::uniform_int_distribution<int> neutral_pair(0, N_He-1);


    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;
    }

    // --- starts reading cross section datafiles

    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;

    std::vector<CrossSection> inelastic1_CS_temp;

    while (excitation1_cs_dat >> energy >> sigma) {
        inelastic1_CS_temp.push_back({energy, sigma});
    }
    excitation1_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");

    // Initialize all electrons
    for (auto& e : electrons) {
        e.initialize(gen, dis, maxwell_electron);
    }
    // initialize all nenutrals
    for (auto&n : neutrals) {
        n.initialize(gen, dis, maxwell_neutral);
    }
    // precalculate elastic cross-section for each energy bin
    for (int i = 0; i < N; i++){
        elastic_vec[i] = interpolate(bin_width*(i+0.5), elastic_CS_temp);
    }
    // precalculate inelastic cross-section (triplet) for each energy bin
    for (int i = 0; i < N; i++){
        inelastic1_vec[i] = interpolate(bin_width*(i+0.5), inelastic1_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
    }
    // -----initial electrons energy distribution ends------------////

    // -----neutrals Maxwell-Boltzmann distribution starts------------////
    for (int i = 0; i < N_He; i++){
        int bin = (int)( (neutrals[i].energy - Emin)/bin_width );
        if (bin >=0 && bin < histo_neutral.size())
            histo_neutral[bin]++;
    }

    for (int i = 0; i < histo_neutral.size(); i++){
        double bin_center = Emin + (i + 0.5) * bin_width;
        file5 << bin_center << " " << static_cast<double>(histo_neutral[i])/(neutrals.size()*bin_width) << "\n"; // this is real f(E) - normalized distribution
        file6 << bin_center << " " << bin_center*static_cast<double>(histo_neutral[i])/(neutrals.size()*bin_width) << "\n"; // this should be E*f(E)

    }
    // -----neutrals Maxwell-Boltzmann distribution starts------------////

    // -----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++){
        sigma_total = elastic_vec[i] + inelastic1_vec[i];
        nu_max_temp = (N_He/Volume)*sigma_total * sqrt(2.0*(i*bin_width + bin_width/2.0)*q/m_e);
        file7 << i << " " << nu_max_temp << "\n";
        if (nu_max_temp > nu_max)
            nu_max = nu_max_temp;
    }
    // -----calculating nu-max for null-collision method ends ------------////

    //----- 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 << dt << "\n";

    // ---- e-e Coulomb collisions variables starts ---/////
    double nu_ee = 0.0;
    // // calculating initial mean energy
    // double sum = 0.0;
    // for (const auto& e : electrons) sum += e.energy;
    // double mean_energy = sum / n_e;
    // nu_ee = (Coulomb_const*(n_e/Volume)*Coulomb_log)/(sqrt(m_e)*pow(mean_energy*q,1.5));

    //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;
//    int Ne_collided = n_e*0.98;  // in case I want to check smth


    // Generate shuffled list of electron indices
    std::vector<int> electron_indices(n_e);
    std::iota(electron_indices.begin(), electron_indices.end(), 0); // fill with index
    std::shuffle(electron_indices.begin(), electron_indices.end(), gen); // shuffle the indexes
    int reshuffle_interval = 1;
    int print_interval = 100;
    int el_coll_counter = 0; // track all elastic collisions
    int exc1_coll_counter = 0; // track all excitation collisions
    int null_coll_counter = 0; // track null-collisions
    int ee_coll_counter = 0; //track e-e Coulomb collisions

    for (int t = 0; t < steps; t++){

        std::cout << "timestep remains: " << steps - t << "\n";

        // calculating mean energy
        double total_energy = 0.0;
        for (const auto& e : electrons) total_energy += e.energy;
        double mean_energy = total_energy / n_e;
        file8 << t*dt << " " << mean_energy << "\n";

        // calculating e-e Coulomb collisions frequency and number of collised electrons

        double nu_ee = (Coulomb_const*(n_e/Volume)*Coulomb_log)/(sqrt(m_e)*pow(mean_energy*q,1.5));
        int EE_collided_pairs = static_cast<int>(0.5 * nu_ee * n_e * dt); // maximum number of electron pairs collided through e-e collisions

        //reshuffle the indices
        if (t % reshuffle_interval == 0) {
            std::shuffle(electron_indices.begin(), electron_indices.end(), gen);
        }

        // 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


        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 electron_energy = e.energy;
            int bin_energy = static_cast<int>(electron_energy / bin_width);
            double nu_elastic = (N_He/Volume) * elastic_vec[bin_energy] * sqrt(2.0*electron_energy*q/m_e);
            double nu_inelastic1 = (N_He/Volume) * inelastic1_vec[bin_energy] * sqrt(2.0*electron_energy*q/m_e);

            double r = dis(gen);

            if (r < nu_elastic/nu_max) {

                // elastic collision happens

                // ----   Collision energy redistribution module

                // electron particle X Y Z initial velocities and energy
                double V0_x_1 = e.vx;
                double V0_y_1 = e.vy;
                double V0_z_1 = e.vz;

                // neutral particle X Y Z initial velocities

                // int k = neutral_pair(gen);

                // double V0_x_2 = neutrals[k].vx;
                // double V0_y_2 = neutrals[k].vy;
                // double V0_z_2 = neutrals[k].vz;

                // randomize particles each collision
                NeutralParticle tmp_neutral;
                tmp_neutral.initialize(gen, dis, maxwell_neutral);
                double V0_x_2 = tmp_neutral.vx;
                double V0_y_2 = tmp_neutral.vy;
                double V0_z_2 = tmp_neutral.vz;

                // initial relative velocity X Y Z (must be equal to final relative velocity in center-of-mass frame)

                double V0_rel_x = (V0_x_1 - V0_x_2);
                double V0_rel_y = (V0_y_1 - V0_y_2);
                double V0_rel_z = (V0_z_1 - V0_z_2);

                double V0_rel = sqrt(V0_rel_x*V0_rel_x + V0_rel_y*V0_rel_y + V0_rel_z*V0_rel_z);

                // center-of-mass frame initial velocity (magnitude of it must be equal to the counterpart in this frame)

                double V_cm_x = (m_e*V0_x_1 + M_n*V0_x_2)/(m_e + M_n);
                double V_cm_y = (m_e*V0_y_1 + M_n*V0_y_2)/(m_e + M_n);
                double V_cm_z = (m_e*V0_z_1 + M_n*V0_z_2)/(m_e + M_n);

                // 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;

                //calculating final relative velocity with random direction

                double V_rel_x = V0_rel*sin(theta)*cos(phi);
                double V_rel_y = V0_rel*sin(theta)*sin(phi);
                double V_rel_z = V0_rel*cos(theta);

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

                //calculating final velocity of electron

                double V_x_1 = V_cm_x + V_rel_x * (M_n/(m_e + M_n));
                double V_y_1 = V_cm_y + V_rel_y * (M_n/(m_e + M_n));
                double V_z_1 = V_cm_z + V_rel_z * (M_n/(m_e + M_n));

                double V_1 = sqrt(V_x_1*V_x_1 + V_y_1*V_y_1 + V_z_1*V_z_1);

                //updating electron energy and velocities

                e.energy = m_e*V_1*V_1/(2.0*q);
                e.vx = V_x_1;
                e.vy = V_y_1;
                e.vz = V_z_1;

                collision_counter_en++;
                el_coll_counter++;

                e.collided_en = true;
            }

            else if (r < (nu_elastic + nu_inelastic1)/nu_max) {

                //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;

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

                // 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 = (2.0 + E_0 - 2.0*pow((1+E_0), R1))/E_0;
                double sin_khi = sqrt(1.0 - cos_khi*cos_khi);

                double phi = 2.0*M_PI*R2;
                double cos_theta = V0_x/V0;
                double sin_theta = sqrt(1.0 - cos_theta*cos_theta);
                // //calculating final relative velocity with random direction

                //calculating final velocity of electron

                double i_scat = (V0_x/V0)*cos_khi + (1.0 - (V0_x/V0)*(V0_x/V0))*(sin_khi*cos(phi)/sin_theta);
                double j_scat = (V0_y/V0)*cos_khi +  (V0_z/V0)*sin_khi*sin(phi)/sin_theta - (V0_x/V0)*(V0_y/V0)*sin_khi*cos(phi)/sin_theta;
                double k_scat = (V0_z/V0)*cos_khi - (V0_y/V0)*sin_khi*sin(phi)/sin_theta - (V0_x/V0)*(V0_z/V0)*sin_khi*cos(phi)/sin_theta;

                //updating electron energy and velocities

                double delta_E = 19.82; //triplet excitation energy

                if (e.energy < delta_E) {
                    null_coll_counter++;
                    collision_counter_en++;
                    e.collided_en = true;
                    continue;
                }
                else {
                    e.energy = E_0 - delta_E;

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

                    e.vx = speed*i_scat;
                    e.vy = speed*j_scat;
                    e.vz = speed*k_scat;

                    collision_counter_en++;
                    exc1_coll_counter++;

                    e.collided_en = true;
                }
            }

            else {
                collision_counter_en++;
                null_coll_counter++;
                e.collided_en = true;
            }
        }

        // ----- -------now begin e-e collisions ------ /////

        // Reshuffle electron indices for random pairing for e-e collisions
        std::shuffle(electron_indices.begin(), electron_indices.end(), gen);

        int max_pairs = std::min(EE_collided_pairs, n_e/2); // handling the case when # of collided pairs ~ ne/2

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

            int id1 = electron_indices[2 * i];
            int id2 = electron_indices[2 * i + 1];
            if (id1 >= n_e || id2 >= n_e) continue; // Handle edge case

            Electron& e1 = electrons[id1];
            Electron& e2 = electrons[id2];

            if (e1.collided_ee || e2.collided_ee) continue; //handle already collided cases

            double E_initial = e1.energy + e2.energy; // total initial energy of pair to check the energy conservation

            // ----   Collision energy redistribution module

            // first particle X Y Z initial velocities
            double V0_x_1 = e1.vx;
            double V0_y_1 = e1.vy;
            double V0_z_1 = e1.vz;

            // second particle X Y Z initial velocities
            double V0_x_2 = e2.vx;
            double V0_y_2 = e2.vy;
            double V0_z_2 = e2.vz;

            // initial relative velocity X Y Z (must be equal to final relative velocity in center-of-mass frame)

            double V0_rel_x = (V0_x_1 - V0_x_2);
            double V0_rel_y = (V0_y_1 - V0_y_2);
            double V0_rel_z = (V0_z_1 - V0_z_2);

            double V0_rel = sqrt(V0_rel_x*V0_rel_x + V0_rel_y*V0_rel_y + V0_rel_z*V0_rel_z);

            // center-of-mass frame initial velocity (magnitude of it must be equal to the counterpart in this frame)

            double V_cm_x = (V0_x_1 + V0_x_2)/2.0;
            double V_cm_y = (V0_y_1 + V0_y_2)/2.0;
            double V_cm_z = (V0_z_1 + V0_z_2)/2.0;

            double V_cm = sqrt(V_cm_x*V_cm_x + V_cm_y*V_cm_y + V_cm_z*V_cm_z);

            // 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;

            //calculating final relative velocity with random direction

            double V_rel_x = V0_rel*sin(theta)*cos(phi);
            double V_rel_y = V0_rel*sin(theta)*sin(phi);
            double V_rel_z = V0_rel*cos(theta);

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

            //calculating final velocity of first particle

            double V_x_1 = V_cm_x + V_rel_x/2.0;
            double V_y_1 = V_cm_y + V_rel_y/2.0;
            double V_z_1 = V_cm_z + V_rel_z/2.0;

            double V_1 = sqrt(V_x_1*V_x_1 + V_y_1*V_y_1 + V_z_1*V_z_1);

            //calculating final velocity of second particle

            double V_x_2 = V_cm_x - V_rel_x/2.0;
            double V_y_2 = V_cm_y - V_rel_y/2.0;
            double V_z_2 = V_cm_z - V_rel_z/2.0;

            double V_2 = sqrt(V_x_2*V_x_2 + V_y_2*V_y_2 + V_z_2*V_z_2);

            // updating velocities

            e1.vx = V_x_1;
            e1.vy = V_y_1;
            e1.vz = V_z_1;

            e2.vx = V_x_2; // Update velocity components
            e2.vy = V_y_2;
            e2.vz = V_z_2;

            // calculating final energies of first and second colliding particles

            e1.energy = V_1*V_1*m_e/(2.0*q);
            e2.energy = V_2*V_2*m_e/(2.0*q);

            double E_final = e1.energy + e2.energy;

            if(fabs(E_final - E_initial) > 1e-6) {
                std::cerr << "Energy conservation violation: " << E_final - E_initial << " eV\n";
            }

            // --- collision energy redistrubution module ends

            // collision counters handling

            collision_counter_ee++;
            ee_coll_counter++;
            e1.collided_ee = true;
            e2.collided_ee = true;

        }
        //////----------------------e-e coulomb collision ends --------------/////////////////

    }

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

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

        int bin = (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 < histo_maxwell.size(); 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";

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

        //recalculating average frequencies to adjust validation constants later
    for (int i = 0; i < histo_maxwell.size(); i++){
        double bin_center = Emin + (i + 0.5) * bin_width;
        file11 << bin_center << " " << (N_He/Volume)*elastic_vec[i] * sqrt(2.0*bin_center*q/m_e)*static_cast<double>(histo_maxwell[i])/(electrons.size()*bin_width) << "\n";
        file12 << bin_center << " " << (N_He/Volume)*inelastic1_vec[i] * sqrt(2.0*bin_center*q/m_e)*static_cast<double>(histo_maxwell[i])/(electrons.size()*bin_width) << "\n";
    }
    //----- calculating nu-average from our electron 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();

    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 triplet excitation collisions per timestep: " << static_cast<int>(exc1_coll_counter/steps) << "\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 << "Average e-e collisions per timestep: " << static_cast<int>(ee_coll_counter/steps) << "\n";

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


    return 0;

}