Coulomb_Nanbu

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

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

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

#define n_e 10000
#define V_0 30000.0     // max velocity using to generate random distribution ---- doesn't work -> produces skew distribution???
#define Emin 0.0
#define Emax 200.0
#define E_average 100.0 // electron sampling energy
#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.99
#define time 5.0E-6
#define dt 1.0E-8
#define steps time/dt
#define Coulomb_log 15.0 // Coulomb logarithm calculated based on Nanmbu paper
#define epsilon_0 8.854188E-12 // Vacuum permittivity
#define Volume 1.0E-14 // volume, m-3

double solve_A(double s) {

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

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

        double f = 1.0/tanh(A) - 1.0/A - exp(-s);

        if (abs(f) < error) break;

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

        A -= f/dfdA;

    }

    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 = 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 * dis(gen) - 1.0;
        double sinTheta = sqrt(1.0 - cosTheta*cosTheta);

        energy = E_average*dis(gen);

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

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

    }

};

int main() {

    std::vector<Electron> electrons(n_e); // better to use vector instead of simple array as it's dynamically allocated (beneficial for ionization)

    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

    clock_t start = clock();

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

    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("mean_energy.dat");



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

    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; i++){
    //     std::cout << i << " " << histo_random.at(i) << "\n"; // using vector.at allows to access elements with boundary checkings
    // }


    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 << bin_start << " " <<  bin_start*static_cast<double>(histo_random[i])/(electrons.size()*bin_width) << "\n"; // E*f(E)
    }

    // for (int i = 0; i < static_cast<int>(4/0.1)+1; i++){
    //     double s = i*0.1;
    //     std::cout << "s = " << s << ", " << "Solution: A = " << solve_A(s) << "\n";
    // }

    double energy_sum = 0.0;

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

        std::cout << "timestep remains: " << steps - t << " " << "deltaE: " << "\n";
        std::ostringstream filename;
        filename << "Coulomb_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++){
            electrons.at(j).collided = 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 && !electrons.at(candidate).collided){
                        partner = candidate;
                        break;
                    }
                    attempts++;
                }

                if (partner != -1) {
                    electrons.at(i).collided = true;
                    electrons.at(partner).collided = true;


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

                // generating random variables to calculate random direction of center-of-mass after the collision

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

                // ----   Collision energy redistribution module

                // --- initial energy of partiecles

                double E0_1 = electrons[i].energy;
                double E0_2 = electrons[partner].energy;

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


                // second particle X Y Z initial velocities
                double V0_x_2 = electrons[partner].vx;
                double V0_y_2 = electrons[partner].vy;
                double V0_z_2 = electrons[partner].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);
                double V0_rel_normal = sqrt(V0_rel_y*V0_rel_y + V0_rel_z*V0_rel_z);


                // calcluating h for equations 20a, 20b (Nanbu1995)

                double eps = 2*M_PI*R1;

                double h_x = V0_rel_normal*cos(eps);
                double h_y = -(V0_rel_y*V0_rel_x*cos(eps) + V0_rel*V0_rel_z*sin(eps))/V0_rel_normal;
                double h_z = -(V0_rel_z*V0_rel_x*cos(eps) - V0_rel*V0_rel_y*sin(eps))/V0_rel_normal;


                //  calculating s (Nanbu1995 eq 19)

                double s = Coulomb_log/(4.0*M_PI) * pow((q*q/(epsilon_0*(m_e/2))),2) * (n_e/Volume) * pow(V0_rel,-3) * dt;
                double A = solve_A(s);




                // calculating cos(khi) (Nanbu1995 eq 17)
                double cos_khi = 0.0;
                double sin_khi = 0.0;

                if (s < 1.0E-2) {// taking care of small s
                    cos_khi = 1.0 + s*log(R1);
                }
                else {
                    cos_khi = (1.0/A)*log(exp(-A) + 2.0*R1*sinh(A));
                    // std::cout << cos_khi << "\n";
                }

                sin_khi = sqrt(1.0 - cos_khi*cos_khi);



                //calculating final velocity of first particle

                double V_x_1 = V0_x_1 - 0.5*(V0_rel_x*(1.0-cos_khi) + h_x*sin_khi);
                double V_y_1 = V0_y_1 - 0.5*(V0_rel_y*(1.0-cos_khi) + h_y*sin_khi);
                double V_z_1 = V0_z_1 - 0.5*(V0_rel_z*(1.0-cos_khi) + h_z*sin_khi);

                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 = V0_x_2 + 0.5*(V0_rel_x*(1.0-cos_khi) + h_x*sin_khi);
                double V_y_2 = V0_y_2 + 0.5*(V0_rel_y*(1.0-cos_khi) + h_y*sin_khi);
                double V_z_2 = V0_z_2 + 0.5*(V0_rel_z*(1.0-cos_khi) + h_z*sin_khi);

                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

                //redistributing energy
                electrons[i].energy = V_1*V_1*m_e/(2.0*q);
                electrons[partner].energy = V_2*V_2*m_e/(2.0*q);

                //redistributing velocity
                electrons[i].vx = V_x_1;
                electrons[i].vy = V_y_1;
                electrons[i].vz = V_z_1;

                electrons[partner].vx = V_x_2;
                electrons[partner].vy = V_y_2;
                electrons[partner].vz = V_z_2;

                // --- collision energy redistrubution module ends

                energy_sum += E0_1 + E0_2 - electrons[i].energy - electrons[partner].energy;

                }

            }
        }

        double total_energy = 0.0;
        for (int i = 0; i < n_e; i++){
            total_energy += electrons[i].energy;
            int bin = (int)( (electrons[i].energy - Emin)/bin_width );
            if (bin >=0 && bin < N)
            histo_maxwell[bin]++;
        }

        file5 << t*dt << " " << total_energy << "\n";

        for (int i = 0; i < N; i++){
            double bin_start = Emin + i*bin_width;
            file << bin_start << " " <<  static_cast<double>(histo_maxwell[i])/(electrons.size()*bin_width) << "\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 << " " << 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;
//        printf("%5.1f - %5.1f\t%d\n", bin_start, bin_start + bin_width, histo_maxwell[i]);
        file4 << bin_start << " " <<  bin_start*static_cast<double>(histo_maxwell[i])/(electrons.size()*bin_width) << "\n"; // E*f(E)
    }


    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";
    std::cout << "Total energy change: " << energy_sum << "\n";


    return 0;


}