e-n collisions(v3)

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

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

#define n_e 50000
#define V_0 30000.0     // max velocity using to generate random distribution ---- doesn't work -> produces skew distribution???
#define Emin 0.0
#define Emax 100.0
#define bin_width 0.01
#define m_e 9.1E-31 // electron mass in kg
#define k_b 1.38E-23 // Boltzmann constant
#define q 1.602176634E-19 // elementary charge    - eV -> J transfer param
#define N ( (int)((Emax-Emin)/bin_width) ) // add 1 to include E_max if needed?
#define P 0.6 // elastic collision constant probability
#define P_null 0.4 // null-collision
#define steps 500
#define T_n 10.0 // Helium neutral temperature in eV
#define M_n 6.6464731E-29 // Helium atom mass
#define N_He 1000000 // Helium neutrals number

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 = false;

    //initialization function // void func(Type0& t) → means the function expects a reference to variable t of type0
    void initialize(std::mt19937& gen, std::uniform_real_distribution<double>& dis) {
        // 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 = Emax*dis(gen);

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

        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>& elastic_CS) {


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

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

    double k = (elastic_CS[step].sigma - elastic_CS[step-1].sigma)/(elastic_CS[step].energy - elastic_CS[step-1].energy);
    double m = elastic_CS[step].sigma - k*elastic_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<double> neutrals(N_He); // vector for neutrals
    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, 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::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<double> dis(0.0, 1.0);
    std::gamma_distribution<double> maxwell(1.5, T_n);

    std::uniform_int_distribution<int> pair(0, n_e-1);


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

    // reading elastic cross section datafile

    std::vector<CrossSection> elastic_CS;

    double energy, sigma;

    while (elastic_cs >> energy >> sigma) {
        elastic_CS.push_back({energy, sigma});
    }

    elastic_cs.close();


    std::ofstream file0("velocities.dat");
    if (!file0.is_open()) {
        std::cerr << "Error opening file!" << std::endl;
        return 1;
    }

    std::ofstream file1("energies.dat");
    if (!file1.is_open()) {
        std::cerr << "Error opening file!" << std::endl;
        return 1;
    }

    std::ofstream file2("energies_final.dat");
    if (!file2.is_open()) {
        std::cerr << "Error opening file!" << std::endl;
        return 1;
    }

    std::ofstream file3("histo_random.dat");
    if (!file3.is_open()) {
        std::cerr << "Error opening file!" << std::endl;
        return 1;
    }
    file3 << std::fixed << std::setprecision(10);

    std::ofstream file4("histo_maxwell.dat");
    if (!file4.is_open()) {
        std::cerr << "Error opening file!" << std::endl;
        return 1;
    }
    file4 << std::fixed << std::setprecision(10);

    std::ofstream file5("neutral_distribution.dat");
    if (!file5.is_open()) {
        std::cerr << "Error opening file!" << std::endl;
        return 1;
    }

    std::ofstream file6("E*f(E).dat");
    if (!file6.is_open()) {
        std::cerr << "Error opening file!" << std::endl;
        return 1;
    }



    // Initialize all electrons
    for (auto& e : electrons) {
        e.initialize(gen, dis);
    }
    // initialize all nenutrals
    for (auto&n : neutrals) {
        n.initialize(gen, dis, maxwell);
    }


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


    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_start = Emin + i*bin_width;
        file3 << i*bin_width << " " <<  static_cast<double>(histo_random[i])/(electrons.size()*bin_width) << "\n"; // dividing by n_e to get normalized distribution
    }


    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_start = Emin + i*bin_width;
        file5 << i*bin_width << " " << static_cast<double>(histo_neutral[i])/(neutrals.size()*bin_width) << "\n"; // this is real f(E) - normalized distribution
        file6 << i*bin_width << " " << (i*bin_width)*static_cast<double>(histo_neutral[i])/(neutrals.size()*bin_width) << "\n"; // this should be E*f(E)

    }

    int Ne_collided = (1.0-P_null)*n_e; // introducing null-collision technique

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

        // setting flags to false each timestep
        for (int i = 0; i < n_e; i++){
            electrons[i].collided = false;
        }

        int collision_counter = 0;

        while (collision_counter < Ne_collided) {

            int i = pair(gen);

            if (dis(gen) < P && !electrons[i].collided) {

                // ----   Collision energy redistribution module

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

                // neutral particle X Y Z initial velocities

                int k = pair(gen);

                double V0_x_2 = neutrals[k].vx;
                double V0_y_2 = neutrals[k].vy;
                double V0_z_2 = neutrals[k].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

                electrons[i].energy = m_e*V_1*V_1/(2.0*q);
                electrons[i].vx = V_x_1;
                electrons[i].vy = V_y_1;
                electrons[i].vz = V_z_1;

                collision_counter += 1;

                electrons[i]. collided = true;
            }
        }
    }

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

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

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

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


    file0.close();
    file1.close();
    file2.close();
    file3.close();
    file4.close();
    file5.close();
    file6.close();

    clock_t end = clock();

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

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


    return 0;

}