Thermalization_1

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

#define n_e 100000
#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.5
#define m_e 9.1E-31 // electron mass in kg
#define k_b 1.38E-23 // Boltzmann constant
#define q 1.6E-19 // elementary charge
#define N ( (int)((Emax-Emin)/bin_width) ) // add 1 to include E_max if needed?
#define P 0.005
#define steps 100

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 initialize(std::mt19937 gen, std::uniform_real_distribution<double> dis) {
        // velocity angles in spherical coordinates
        double phi = 2*M_PI*dis(gen);
        double sinTheta = dis(gen);
        double cosTheta = sqrt(1.0 - sinTheta*sinTheta);

        energy = Emax*dis(gen);

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

        vx = speed * sinTheta * cos(phi);
        vy = speed * sinTheta * sin(phi);
        vz = speed * cosTheta;

    }

};

int main() {

    double electron_vel[3*n_e];
    double electron_e[n_e];
    int histo_random[N] = {0};
    int histo_maxwell[N] = {0};
    bool collided[n_e] = {0};

    clock_t start = clock();

    for (int i = 0; i < 3*n_e; i++){
        electron_vel[i] = 0.0;
    }

    for (int i = 0; i < n_e; i++){
        electron_e[i] = 0.0;
    }


    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::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution<double> dis(0.0, 1.0);

    std::random_device rd1;
    std::mt19937 gen1(rd1());
    std::uniform_int_distribution<int> pair(0, n_e-1);


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

        double phi = 2*M_PI*dis(gen);
        double sinTheta = dis(gen);
        double cosTheta = sqrt(1.0 - sinTheta*sinTheta);

        double energy = Emax*dis(gen);
        double speed = sqrt(2*energy*k_b/m_e);

        electron_vel[i*3] = speed * sinTheta * cos(phi); // x component
        electron_vel[i*3+1] = speed * sinTheta * sin(phi); // y component
        electron_vel[i*3+2] = speed * cosTheta; // z component

        electron_e[i] = energy;

        file0  << i << " " << electron_vel[i*3] << " " << electron_vel[i*3+1] << " " << electron_vel[i*3+2]  << "\n";

        file1 << i << " " << electron_e[i] << "\n";

        energy = 0.0;
    }

    for (int i = 0; i < n_e; i++){
        int bin = (int)( (electron_e[i] - Emin)/bin_width );
        if (bin >=0 && bin < N)
            histo_random[bin]++;
    }

    for (int i = 0; i < N; i++){
        double bin_start = Emin + i*bin_width;
//        printf("%5.1f - %5.1f\t%d\n", bin_start, bin_start + bin_width, histo_random[i]);
        file3 << i*bin_width << " " <<  static_cast<double>(histo_random[i]) << "\n"; // dividing by n_e to get normalized distribution
    }


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

        std::ostringstream filename;
        filename << "data/distribution_" << std::setw(4) << std::setfill('0') << t << ".dat";

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

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

        for (int i = 0 ; i < n_e; i++){
            if (dis(gen) < P) {
                int partner = -1;
                int attempts = 0;

                while (attempts < n_e){
                    int candidate = pair(gen1);
                    if (candidate != i && !collided[candidate]){
                        partner = candidate;
                        break;
                    }
                    attempts++;
                }

                if (partner != -1) {
                    collided[i] = true;
                    collided[partner] = true;
                }

                // ----   Collision energy redistribution module

                // first particle X Y Z initial velocities
                double V0_x_1 = electron_vel[i*3 + 0];
                double V0_y_1 = electron_vel[i*3 + 1];
                double V0_z_1 = electron_vel[i*3 + 2];

                // second particle X Y Z initial velocities
                double V0_x_2 = electron_vel[partner*3 + 0];
                double V0_y_2 = electron_vel[partner*3 + 1];
                double V0_z_2 = electron_vel[partner*3 + 2];

                // 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)/2.0;
                double V0_rel_y = (V0_y_1 - V0_y_2)/2.0;
                double V0_rel_z = (V0_z_1 - V0_z_2)/2.0;

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

                // calculating final energies of first and second colliding particles

                electron_e[i] = V_1*V_1*m_e/(2.0*k_b);
                electron_e[partner] = V_2*V_2*m_e/(2.0*k_b);

                // --- collision energy redistrubution module ends

            }
        }

        for (int i = 0; i < n_e; i++){
//        filename << i << " " << electron_e[i] << "\n";
        int bin = (int)( (electron_e[i] - Emin)/bin_width );
        if (bin >=0 && bin < N)
            histo_maxwell[bin]++;
        }

        for (int i = 0; i < N; i++){
            double bin_start = Emin + i*bin_width;
            file << i*bin_width << " " <<  static_cast<double>(histo_maxwell[i]) << "\n"; // later need to divide by total partcles number to get normalized distribution
            histo_maxwell[i] = 0;
        }

        file.close();

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

//         file2 << i << " " << electron_e[i] << "\n";

//         int bin = (int)( (electron_e[i] - Emin)/bin_width );
//         if (bin >=0 && bin < N)
//             histo_maxwell[bin]++;
//     }

//     for (int i = 0; i < N; i++){
//         double bin_start = Emin + i*bin_width;
// //        printf("%5.1f - %5.1f\t%d\n", bin_start, bin_start + bin_width, histo_maxwell[i]);
//         file4 << i*bin_width << " " <<  static_cast<double>(histo_maxwell[i]) << "\n"; // dividing by partcles number to get normalized distribution
//     }


    file0.close();
    file1.close();
    file2.close();
    file3.close();
    file4.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;


}