#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 thresh1 19.82 // threshold energy excitation tripet state
#define thresh2 20.61 // threshold energy excitation singlet state
// 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 5.0 // electron Maxwell initial distribution
#define Emin 0.0
#define Emax 4000.0
#define Volume 1.0E-12 // Volume to calculate netral density and collision frequency
#define time 3.0E-4 // 500 microsec time to equalibrate the system
#define dopant 1.0E-5 // addition to avoid zero
#define E_reduced 0.1 // 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";
// Linear extrapolation using last two points
const auto& last = CS[CS.size() - 1];
const auto& prev = CS[CS.size() - 2];
double slope = (last.sigma - prev.sigma) / (last.energy - prev.energy);
return std::max(0.0, last.sigma + slope * (energy - last.energy));
}
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 = 10000000;
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");
// 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
}
// -----initial electrons energy distribution ends------------////
//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);
// 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 //
int print_interval = 100;
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 total_time = 3.0E-3; // Your original `time` constant
double t_elapsed = 0.0;
// -----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 << nu_max << "\n";
// -----calculating nu-max for null-collision method ends ------------////
double dt = 0.1/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;
}
std::cout << "Progress: " << (t_elapsed/total_time)*100 << "%" << " ";
//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
std::shuffle(electron_indices.begin(), electron_indices.end(), gen); // shuffle the indexes
// Generate shuffled list of triplet excited atoms indices
std::vector<int> excited1_indices(exc_1.size());
std::iota(excited1_indices.begin(), excited1_indices.end(), 0); // fill with index
std::shuffle(excited1_indices.begin(), excited1_indices.end(), gen); // shuffle the indexes
// Generate shuffled list of singlet excited atoms indices
std::vector<int> excited2_indices(exc_2.size());
std::iota(excited2_indices.begin(), excited2_indices.end(), 0); // fill with index
std::shuffle(excited2_indices.begin(), excited2_indices.end(), gen); // shuffle the indexes
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;
// 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_elapsed << " " << mean_energy << "\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--- ///
for (int idx : electron_indices) {
// 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;
// Recalculate energy from updated velocity
double vx = electrons[idx].vx;
double vy = electrons[idx].vy;
double vz = electrons[idx].vz;
electrons[idx].energy = 0.5 * m_e * (vx*vx + vy*vy + vz*vz) / q;
}
// -------------------------------------------- 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 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 nu_superelastic1 = (exc_1.size()/Volume) * superelastic1_vec[bin_energy] * sqrt(2.0*electron_energy*q/m_e);
double nu_inelastic2 = (N_He/Volume) * inelastic2_vec[bin_energy] * sqrt(2.0*electron_energy*q/m_e);
double nu_superelastic2 = (exc_2.size()/Volume) * superelastic2_vec[bin_energy] * sqrt(2.0*electron_energy*q/m_e);
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);
// 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;
if (cos_khi >= 1)
cos_khi = 1.0 - dopant;
if (cos_khi <= -1)
cos_khi = -1.0 + dopant;
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 = 2.0*(m_e/M_n)*(1.0 - cos_khi)*E_0;
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++;
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);
// 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
if (e.energy < thresh1) {
null_coll_counter++;
collision_counter_en++;
e.collided_en = true;
continue;
}
else {
e.energy = E_0 - thresh1;
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++;
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;
double V0 = sqrt(V0_x*V0_x + V0_y*V0_y + V0_z*V0_z);
// 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 = exc.vx;
double V_y = exc.vy;
double V_z = exc.vz;
double E = exc.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 = (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
e.energy = E_0 + thresh1;
double speed = sqrt(2*e.energy*q/m_e);
e.vx = speed*i_scat;
e.vy = speed*j_scat;
e.vz = speed*k_scat;
//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++;
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);
// 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
if (e.energy < thresh2) {
null_coll_counter++;
collision_counter_en++;
e.collided_en = true;
continue;
}
else {
e.energy = E_0 - thresh2;
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++;
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;
double V0 = sqrt(V0_x*V0_x + V0_y*V0_y + V0_z*V0_z);
// 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 = exc.vx;
double V_y = exc.vy;
double V_z = exc.vz;
double E = exc.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 = (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
e.energy = E_0 + thresh2;
double speed = sqrt(2*e.energy*q/m_e);
e.vx = speed*i_scat;
e.vy = speed*j_scat;
e.vz = speed*k_scat;
//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++;
e.collided_en = true;
}
else {
// null-collision
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 = n_e/2; // each electron collides
// 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
// // 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
// // 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;
// // file13 << "V0_x_1: " << V0_x_1 << " " << "V0_y_1: " << V0_y_1 << " " << " V0_z_1: " << V0_z_1 << " ";
// // file13 << "V0_x_2: " << V0_x_2 << " " << "V0_y_2: " << V0_y_2 << " " << " V0_z_2: " << V0_z_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);
// double V0_rel_y = (V0_y_1 - V0_y_2);
// double V0_rel_z = (V0_z_1 - V0_z_2);
// if(std::isnan(V0_x_1) || std::isinf(V0_x_1) || fabs(V0_x_1) < 1e-12 || std::isnan(V0_y_1) || std::isinf(V0_y_1) || fabs(V0_y_1) < 1e-12 || std::isnan(V0_z_1) || std::isinf(V0_z_1) || fabs(V0_z_1) < 1e-12){
// std::cerr << "Invalid V0_rel_x computed: " << V0_rel_x << " at timestep " << t << std::endl;
// std::cerr << "Components of velocities: Vx, Vy, Vz forr the first electron: " << V0_x_1 << " " << V0_y_1 << " " << V0_z_1 << "\n";
// // continue;
// }
// if(std::isnan(V0_x_2) || std::isinf(V0_x_2) || fabs(V0_x_2) < 1e-12 || std::isnan(V0_y_2) || std::isinf(V0_y_2) || fabs(V0_y_2) < 1e-12 || std::isnan(V0_z_2) || std::isinf(V0_z_2) || fabs(V0_z_2) < 1e-12){
// std::cerr << "Invalid V0_rel_x computed: " << V0_rel_x << " at timestep " << t << std::endl;
// std::cerr << "Components of velocities: Vx, Vy, Vz forr the second electron: " << V0_x_2 << " " << V0_y_2 << " " << V0_z_2 << "\n";
// // continue;
// }
// if(std::isnan(V0_rel_y) || std::isinf(V0_rel_y) || fabs(V0_rel_y) < 1e-12){
// std::cerr << "Invalid V0_rel_y computed: " << V0_rel_y << " at timestep " << t << std::endl;
// continue;
// }
// if(std::isnan(V0_rel_z) || std::isinf(V0_rel_z) || fabs(V0_rel_z) < 1e-12){
// std::cerr << "Invalid V0_rel_z computed: " << V0_rel_z << " at timestep " << t << std::endl;
// continue;
// }
// 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);
// // file13 << "V0_rel: " << V0_rel << " " << "V0_rel_normal: " << V0_rel_normal << " ";
// if(std::isnan(V0_rel) || std::isinf(V0_rel) || fabs(V0_rel) < 1e-12){
// std::cerr << "Invalid V0_rel computed: " << V0_rel << " at timestep " << t << std::endl;
// V0_rel = 1.0E-6;
// e1.collided_ee = true;
// e2.collided_ee = true;
// continue;
// }
// if(std::isnan(V0_rel_normal) || std::isinf(V0_rel_normal) || fabs(V0_rel_normal) < 1e-12){
// std::cerr << "Invalid V0_rel_normal computed: " << V0_rel << " at timestep " << t << std::endl;
// continue;
// }
// // calculating spherical angles for center-of-mass random direction
// double theta = acos(1.0- 2.0*R1);
// double phi = 2*M_PI*R2;
// // 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;
// // file13 << "s: " << s << " ";
// if(std::isnan(s) || std::isinf(s) || fabs(s) < 1e-12){
// std::cerr << "Invalid s computed: " << s << " at timestep " << t << std::endl;
// continue;
// }
// double A = solve_A(s);
// if(std::isnan(A) || std::isinf(A) || fabs(A) < 1e-12){
// // std::cerr << "Invalid A computed: " << A << " at timestep " << t << std::endl;
// A = 1.0E-12;
// // continue;
// }
// // calculating cos(khi) (Nanbu1995 eq 17)
// double cos_khi = 0.0;
// double sin_khi = 0.0;
// if (s < 1.0E-2 & R1 != 0.0) {// 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));
// }
// if (cos_khi > 1.0)
// cos_khi = 1.0;
// 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);
// // 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
// ee_coll_counter++;
// e1.collided_ee = true;
// e2.collided_ee = true;
// }
// //////----------------------e-e coulomb collision ends --------------/////////////////
t_elapsed += dt; // Advance time
// Recalculate nu_max periodically (e.g., every 100 steps)
static int recalc_counter = 0;
if (++recalc_counter >= 1) {
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 v = sqrt(2 * energy * q / m_e);
double nu_temp = sigma_total * v;
if (nu_temp > nu_max) nu_max = nu_temp;
}
// Update dt based on new nu_max
dt = 0.1 / nu_max;
}
std::cout << " nu_max: " << nu_max << " " << "dt: " << dt << "\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";
// ----- final electron energies distribution ends
// // ------ excited atoms histogram --------/////
// for (int i = 0; i < exc_1.size(); i++) {
// file14 << i << " " << exc_1[i].energy << "\n";
// int bin = static_cast<int>( (exc_1[i].energy - Emin)/bin_width);
// if (bin >=0 && bin < histo_excited.size())
// histo_excited[bin]++;
// }
// for (int i = 0; i < histo_excited.size(); i++){
// double bin_center = Emin + (i + 0.5) * bin_width;
// file15 << bin_center << " " << static_cast<double>(histo_excited[i])/(electrons.size()*bin_width) << "\n"; // getting f(E)
// }
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();
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;
}