API tools faq
paste
phystota

v3.4 - Clean, full

phystota
Apr 25th, 2025
57
0
Never
1
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 50.51 KB | None | 0 0
raw download clone embed print report
  1. #include <iostream>
  2. #include <random>
  3. #include <fstream>
  4. #include <assert.h>
  5.  
  6.  
  7. #include <omp.h>
  8. #include <chrono>
  9.  
  10. #include <math.h>
  11. #include <time.h>
  12. #include <iomanip>  // For std::fixed and std::setprecision
  13.  
  14. #include <algorithm>  // For std::shuffle
  15. #include <numeric>    // For std::iota
  16.  
  17. #include <random>
  18. #include <limits>
  19.  
  20. inline double fast_uniform(std::mt19937 &g) {
  21.     // g() returns uint32_t in [0, 2^32−1].
  22.     // Divide by max to get a double in [0,1].
  23.     return double(g()) * (1.0 / double(std::numeric_limits<uint32_t>::max()));
  24. }
  25.  
  26. //physical constants
  27.  
  28. #define m_e 9.1093837E-31 // electron mass in kg
  29. #define M_n 6.6464731E-27 // Helium atom mass
  30. #define k_b 1.380649E-23 // Boltzmann constant
  31. #define q 1.602176634E-19 // elementary charge    - eV -> J transfer param
  32. #define Coulomb_log 15.87 // Coulomb logarithm
  33. #define epsilon_0 8.854188E-12 // Vacuum permittivity
  34. #define Coulomb_const pow(q,4)/(pow(4.0*M_PI*epsilon_0,2)) // e^4/(4*pi*eps0)^2
  35. #define thresh0 0.0 //elastic threshold
  36. #define thresh1 19.82 // threshold energy excitation tripet state
  37. #define thresh2 20.61 // threshold energy excitation singlet state
  38. #define thresh3 24.587 // threshold energy ionization HeI
  39. #define tau_singlet 0.0195
  40.  
  41. // simulation parameters
  42.  
  43.  
  44. // #define N_He 1000000 // Helium neutrals number
  45. #define T_n 2.0 // Helium neutral temperature in eV
  46. #define T_e 10.0    // electron Maxwell initial distribution
  47. #define Emin 0.0
  48. #define Emax 3000.0
  49. #define Volume 1.0E-12 // Volume to calculate netral density and collision frequency
  50. #define total_time 8.0E-4 // 500 microsec time to equalibrate the system
  51. #define dopant 1.0E-5 // addition to avoid zeros
  52. #define E_reduced 0.6 // constant electrin field along z-axis
  53.  
  54. #define bin_width 0.05 // keep energy step ~ this to maintain cross-section clarity (Ramsauer minimum etc)
  55. #define N ( (int)((Emax-Emin)/bin_width) + 1) // add 1 to include E_max if needed?
  56.  
  57. // handling final energy bin
  58.  
  59. #define bin_width_smooth 0.5 // energy bin for smooth final distribution
  60. #define N_smooth ( (int)((Emax-Emin)/bin_width_smooth) )
  61.  
  62.  
  63.  
  64. double solve_A(double s) { // Netwon method solver
  65.  
  66.     if (s > 3) {
  67.         return 3*exp(-s);
  68.     }
  69.     if (s < 0.01) {
  70.         return 1.0/s;
  71.     }
  72.    
  73.     double A0 = 0.01; // initial guess
  74.     double A = A0;  //starting with initial guess
  75.     double tol = 1.0E-7; // accuracy
  76.  
  77.              
  78.     for (int i = 0; i < 1000; i++){
  79.  
  80.         double tanhA = tanh(A);
  81.         // Avoid division by an extremely small tanh(A)
  82.         if (fabs(tanhA) < 1e-12) {
  83.             std::cerr << "tanh(A) too small, returning fallback at iteration " << i << "\n";
  84.             return 1.0E-7;
  85.         }        
  86.  
  87.         double f = 1.0 / tanhA - 1.0 / A - exp(-s);
  88.         if (fabs(f) < tol)
  89.             break;
  90.  
  91.         double sinhA = sinh(A);
  92.         if (fabs(sinhA) < 1e-12) {
  93.             std::cerr << "sinh(A) too small, returning fallback at iteration " << i << "\n";
  94.             return 1.0E-5;
  95.         }
  96.  
  97.         double dfdA = -1.0/(sinh(A)*sinh(A)) + 1.0/(A*A);
  98.  
  99.         // Check if derivative is too close to zero to avoid huge updates
  100.         if (fabs(dfdA) < 1e-12) {
  101.             std::cerr << "dfdA is too small at iteration " << i << ", returning fallback\n";
  102.             if (s < 0.01) {
  103. //                std::cout << "Small s! Huge A!" << "\n";
  104.                 return 1.0/s;
  105.             }
  106.             if (s > 3) {
  107.                 return 3.0*exp(-s);
  108.             }
  109.         }        
  110.  
  111.         A -= f/dfdA;
  112.  
  113.         // Early check for numerical issues
  114.         if (std::isnan(A) || std::isinf(A)) {
  115.             std::cerr << "Numerical error detected, invalid A at iteration " << i << "\n";
  116.             return (A > 0) ? 1.0E-5 : -1.0E-5;  // Fallback value based on sign
  117.         }        
  118.  
  119.  
  120.     }
  121.  
  122.     return A;
  123. }
  124.  
  125. struct Electron {
  126.  
  127.     //velocity components
  128.     double vx = 0.0;
  129.     double vy = 0.0;
  130.     double vz = 0.0;
  131.     //energy in eV
  132.     double energy = 0.0;
  133.     //Collision flag
  134.     bool collided_en = false;
  135.     bool collided_ee = false;
  136.  
  137.     // initializing Maxwell-Boltzmann distribution with T_e
  138.     void initialize(std::mt19937& gen, std::uniform_real_distribution<double>& dis, std::gamma_distribution<double>& maxwell) {
  139.  
  140.         double R = dis(gen);
  141.  
  142.         // velocity angles in spherical coordinates
  143.         double phi = 2*M_PI*dis(gen);
  144.         double cosTheta = 2.0*dis(gen) - 1.0;
  145.         double sinTheta = sqrt(1.0 - cosTheta*cosTheta);
  146.  
  147.            
  148.         energy = maxwell(gen); // neutrals energies sampled as Maxwell distribution in eV
  149.            
  150.         double speed = sqrt(2*energy*q/m_e);
  151.  
  152.         //velocity components of neutrals in m/s
  153.         vx = speed * sinTheta * cos(phi);
  154.         vy = speed * sinTheta * sin(phi);
  155.         vz = speed * cosTheta;
  156.     }
  157.  
  158.  
  159. };
  160.  
  161. struct CrossSection {
  162.     double energy;
  163.     double sigma;
  164. };
  165.  
  166. double interpolate (double energy, const std::vector<CrossSection>& CS) {
  167.  
  168.  
  169.     if (energy < CS.front().energy) {
  170. //        std::cout << " required energy value lower than range of cross-section data at energy: " << energy << "\n";
  171.         return 0.0;
  172.     }
  173.     if (energy > CS.back().energy) {
  174. //        std::cout << " required energy value higher than range of cross-section data" << "\n";
  175.         return 0.0;        
  176.     }
  177.  
  178.     int step = 0;  
  179.         while (step < CS.size() && energy > CS[step].energy) {
  180.             step++;
  181.         }
  182.  
  183.     double k = (CS[step].sigma - CS[step-1].sigma)/(CS[step].energy - CS[step-1].energy);
  184.     double m = CS[step].sigma - k*CS[step].energy;
  185.    
  186.     return k*energy + m;
  187. }
  188.  
  189.  
  190. struct NeutralParticle {
  191.  
  192.     double energy = 0.0;
  193.     double vx = 0.0;
  194.     double vy = 0.0;
  195.     double vz = 0.0;
  196.  
  197.     void initialize(std::mt19937& gen, std::uniform_real_distribution<double>& dis, std::gamma_distribution<double>& maxwell) {
  198.  
  199.         double R = dis(gen);
  200.  
  201.         // velocity angles in spherical coordinates
  202.         double phi = 2*M_PI*dis(gen);
  203.         double cosTheta = 2.0*dis(gen) - 1.0;
  204.         double sinTheta = sqrt(1.0 - cosTheta*cosTheta);
  205.  
  206.            
  207.         energy = maxwell(gen); // neutrals energies sampled as Maxwell distribution in eV
  208.            
  209.         double speed = sqrt(2*energy*q/M_n);
  210.  
  211.         //velocity components of neutrals in m/s
  212.         vx = speed * sinTheta * cos(phi);
  213.         vy = speed * sinTheta * sin(phi);
  214.         vz = speed * cosTheta;
  215.     }
  216.    
  217. };
  218.  
  219. struct Excited_neutral {
  220.  
  221.     double energy;
  222.     double vx;
  223.     double vy;
  224.     double vz;
  225.    
  226. };
  227.  
  228.  
  229.  
  230. int main() {
  231.  
  232.     clock_t start = clock();
  233.    
  234.     int N_He = 10000000;
  235.     int n_e = 500000;
  236.  
  237.     std::vector<Electron> electrons(n_e); // better to use vector instead of simple array as it's dynamically allocated (beneficial for ionization)
  238. //    std::vector<NeutralParticle> neutrals(N_He); // I don't need a vector of neutrals bcs it's like a backhround in MCC-simulation
  239.  
  240.     std::vector<int> histo_random(N, 0); // initialize N size zero-vector for random (initial) histogram
  241.     std::vector<int> histo_maxwell(N, 0); // initialize N size zero-vector for maxwellian histogram
  242.     std::vector<int> histo_neutral(N, 0); // initialize N size zero-vector for neutral distribution histogram
  243.     std::vector<int> histo_excited(N, 0); // initialize N size zero-vector for excited He distribution histogram
  244.  
  245.     std::vector<double> elastic_vec(N, 0); // precompiled elastic cross-section-energy vector
  246.     std::vector<double> inelastic1_vec(N, 0); // precompiled inelastic(triplet excitation) cross-section-energy vector
  247.     std::vector<double> inelastic2_vec(N, 0); // precompiled inelastic(singlet excitation) cross-section-energy vector    
  248.     std::vector<double> superelastic1_vec(N, 0); // precompiled superelastic(triplet de-excitation) cross-section-energy vector
  249.     std::vector<double> superelastic2_vec(N, 0); // precompiled superelastic(triplet de-excitation) cross-section-energy vector
  250.     std::vector<double> ionization_vec(N, 0); // precompiled ionization cross-section-energy vector
  251.  
  252.     std::random_device rd;
  253.     std::mt19937 gen(rd());
  254.     std::uniform_real_distribution<double> dis(0.0, 1.0);
  255.     std::gamma_distribution<double> maxwell_neutral(1.5, T_n);
  256.     std::gamma_distribution<double> maxwell_electron(1.5, T_e);
  257.  
  258.     std::ifstream elastic_cs_dat("cross_sections/elastic.dat");
  259.     if (!elastic_cs_dat.is_open()) {
  260.         std::cerr << "Error opening elastic cross-sections file!" << std::endl;
  261.         return 1;
  262.     }    
  263.  
  264.     std::ifstream excitation1_cs_dat("cross_sections/inelastic_triplet.dat");
  265.     if (!excitation1_cs_dat.is_open()) {
  266.         std::cerr << "Error opening inelastic triplet cross-sections file!" << std::endl;
  267.         return 1;
  268.     }
  269.  
  270.     std::ifstream excitation2_cs_dat("cross_sections/inelastic_singlet.dat");
  271.     if (!excitation2_cs_dat.is_open()) {
  272.         std::cerr << "Error opening inelastic singlet cross-sections file!" << std::endl;
  273.         return 1;
  274.     }
  275.  
  276.     std::ifstream ionization_cs_dat("cross_sections/ionization.dat");
  277.     if (!ionization_cs_dat.is_open()) {
  278.         std::cerr << "Error opening inelastic singlet cross-sections file!" << std::endl;
  279.         return 1;
  280.     }
  281.  
  282.     // --- starts reading cross section datafiles
  283.  
  284. //-----------------elastic---------------------------//
  285.     std::vector<CrossSection> elastic_CS_temp;
  286.  
  287.     double energy, sigma;
  288.  
  289.     while (elastic_cs_dat >> energy >> sigma) {
  290.         elastic_CS_temp.push_back({energy, sigma});
  291.     }    
  292.     elastic_cs_dat.close();
  293.  
  294.     energy = 0.0;
  295.     sigma = 0.0;
  296. //-----------------triplet excitation---------------------------//
  297.     std::vector<CrossSection> inelastic1_CS_temp;
  298.  
  299.     while (excitation1_cs_dat >> energy >> sigma) {
  300.         inelastic1_CS_temp.push_back({energy, sigma});
  301.     }    
  302.     excitation1_cs_dat.close();    
  303. //-----------------singlet excitation---------------------------//
  304.     energy = 0.0;
  305.     sigma = 0.0;
  306.  
  307.     std::vector<CrossSection> inelastic2_CS_temp;
  308.  
  309.     while (excitation2_cs_dat >> energy >> sigma) {
  310.         inelastic2_CS_temp.push_back({energy, sigma});
  311.     }    
  312.     excitation2_cs_dat.close();    
  313. //-----------------Ionization---------------------------//
  314.     energy = 0.0;
  315.     sigma = 0.0;
  316.  
  317.     std::vector<CrossSection> ionization_CS_temp;
  318.  
  319.     while (ionization_cs_dat >> energy >> sigma) {
  320.         ionization_CS_temp.push_back({energy, sigma});
  321.     }    
  322.     ionization_cs_dat.close();    
  323.     // --- finish reading cross-section datafiles  
  324.  
  325.     std::ofstream file0("output_files/velocities.dat");    
  326.     std::ofstream file1("output_files/energies.dat");        
  327.     std::ofstream file2("output_files/energies_final.dat");    
  328.     std::ofstream file3("output_files/histo_random.dat");    
  329.     file3 << std::fixed << std::setprecision(10);
  330.    
  331.     std::ofstream file4("output_files/histo_maxwell.dat");
  332.     file4 << std::fixed << std::setprecision(10);          
  333.    
  334.     std::ofstream file5("output_files/neutral_distribution.dat");    
  335.     std::ofstream file6("output_files/E*f(E).dat");    
  336.     std::ofstream file7("output_files/nu_max.dat");
  337.     std::ofstream file8("output_files/electron_mean_energy.dat");
  338.     std::ofstream file9("output_files/nu_elastic_average_initial.dat");
  339.     std::ofstream file10("output_files/nu_inelastic1_average_initial.dat");
  340.     std::ofstream file11("output_files/nu_elastic_average_final.dat");
  341.     std::ofstream file12("output_files/nu_inelastic1_average_final.dat");
  342.     std::ofstream file13("output_files/log_output.dat");  
  343.     std::ofstream file14("output_files/excited_energies.dat");      
  344.     std::ofstream file15("output_files/excited_histo.dat");            
  345.     std::ofstream file_temp("output_files/collision_rates.dat");
  346.     std::ofstream file16("output_files/energy_gain.dat");  
  347.     std::ofstream file17("output_files/ionization_check.dat");
  348.  
  349.     // Initialize all electrons
  350.     for (auto& e : electrons) {
  351.         e.initialize(gen, dis, maxwell_electron);
  352.     }
  353.  
  354.     // precalculate cross-sections for each energy bin
  355.     for (int i = 0; i < N; i++){
  356.         elastic_vec[i] = interpolate(bin_width*(i+0.5), elastic_CS_temp); //elastiuc
  357.         inelastic1_vec[i] = interpolate(bin_width*(i+0.5), inelastic1_CS_temp); //triplet excitation
  358.         inelastic2_vec[i] = interpolate(bin_width*(i+0.5), inelastic2_CS_temp); //singlet excitation
  359.         ionization_vec[i] = interpolate(bin_width*(i+0.5), ionization_CS_temp); //ionization
  360.     }
  361.  
  362.     // precalculate superelastic cross-section (triplet -> ground) for each energy bin
  363.     // detailed balance gives: sigma_j_i(E) = (g_i/g_j)*((E+theshold)/E)*sigma_i_j(E+theshold)
  364.     for (int i = 0; i < N; i++){
  365.         double energy = Emin + (i + 0.5) * bin_width;
  366.         int thresh_bin = (int)( (thresh1 - Emin)/bin_width ); // calculating bin for threshold energy
  367.         superelastic1_vec[i] = (1.0/3.0)*((energy + thresh1)/energy)*interpolate(energy + thresh1, inelastic1_CS_temp); // using detailed balance, calculating backward deexcitation cross-section
  368.         superelastic2_vec[i] = (1.0/1.0)*((energy + thresh2)/energy)*interpolate(energy + thresh2, inelastic2_CS_temp);
  369.     }
  370.  
  371.     for (int i = 0; i < n_e; i++){
  372.         file1 << i << " " << electrons.at(i).energy << "\n";
  373.         file0 << i << " " << electrons[i].vx << " " << electrons[i].vy << " " << electrons[i].vz << "\n";
  374.     }
  375.  
  376.     // -----initial electrons energy distribution starts------------////
  377.     for (int i = 0; i < n_e; i++){
  378.         int bin = (int)( (electrons[i].energy - Emin)/bin_width );
  379.         if (bin >=0 && bin < histo_random.size())
  380.             histo_random[bin]++;
  381.     }
  382.  
  383.     for (int i = 0; i < histo_random.size(); i++){
  384.         double bin_center = Emin + (i + 0.5) * bin_width;
  385.         file3 << bin_center << " " <<  static_cast<double>(histo_random[i])/(electrons.size()*bin_width) << "\n"; // this is electron normalized distribution function
  386.     }
  387.  
  388.     //calculating excited specied population
  389.  
  390.     /*From Boltzman distribution y_k = g_k*exp(-E_k/kT)/norm, where g_k - stat weight of k-state,
  391.     E_k - threshold energy for k-state, norm is a total partition function or normaliztion factor     */
  392.  
  393.     double part_ground = 1.0*exp(-0.0/T_n); // partition function for ground state
  394.     double part_triplet = 3.0*exp(-thresh1/T_n); // partition function for triplet excited state
  395.     double part_singlet = 1.0*exp(-thresh2/T_n); // partition function for singlet excited state
  396.     double part_func_total = part_ground + part_triplet + part_singlet; // total partition function
  397.     double N_trpilet = (part_triplet/part_func_total)*N_He; // population of tripet state
  398.     double N_singlet = (part_singlet/part_func_total)*N_He; // population of singlet state
  399.  
  400.     std::vector<Excited_neutral> exc_1(static_cast<int>(N_trpilet));  // vector to track triplet excited atoms of Helium
  401.     std::vector<Excited_neutral> exc_2(static_cast<int>(N_singlet));  // vector to track singlet excited atoms of Helium    
  402.  
  403.     // adjusting neutrals number:
  404.  
  405.     N_He -= (N_trpilet + N_singlet);
  406.  
  407.     std::cout << N_He << "\n";
  408.  
  409.     // initializing excited species with Maxwellian distribution
  410.  
  411.     for (auto& exc : exc_1) {
  412.     NeutralParticle tmp_neutral;
  413.     tmp_neutral.initialize(gen, dis, maxwell_neutral);
  414.     exc.energy = tmp_neutral.energy;
  415.     exc.vx = tmp_neutral.vx;
  416.     exc.vy = tmp_neutral.vy;
  417.     exc.vz = tmp_neutral.vz;
  418.     }
  419.  
  420.     for (auto& exc : exc_2) {
  421.     NeutralParticle tmp_neutral;
  422.     tmp_neutral.initialize(gen, dis, maxwell_neutral);
  423.     exc.energy = tmp_neutral.energy;
  424.     exc.vx = tmp_neutral.vx;
  425.     exc.vy = tmp_neutral.vy;
  426.     exc.vz = tmp_neutral.vz;
  427.     }
  428.  
  429.     std::cout << "Triplet population initialized: " << exc_1.size() << "\n";
  430.     std::cout << "Singlet population initialized: " << exc_2.size() << "\n";    
  431.  
  432.     // calculating excited specied population finished //
  433.  
  434.  
  435.     int print_interval = 10;
  436.     int el_coll_counter = 0; // track all elastic collisions
  437.     int exc1_coll_counter = 0; // track all triplet excitation collisions
  438.     int exc2_coll_counter = 0; // track all singlet excitation collisions
  439.     int null_coll_counter = 0; // track null-collisions
  440.     int ee_coll_counter = 0; //track e-e Coulomb collisions
  441.     int super1_coll_counter = 0; // track superelastic triplet collisions
  442.     int super2_coll_counter = 0; // track superelastic triplet collisions    
  443.     int ionization_counter = 0; //track ionization collisions
  444.  
  445.  
  446.     double a_z = ((-1.0)*q * E_reduced) / m_e;
  447.     double mass_ratio = 2.0*(m_e/M_n);
  448.     double charge_mass_ratio = 0.5*m_e/q;
  449.     double sqrt_charge_mass = sqrt(2*q/m_e);
  450.     double C1 = -1.0*q*E_reduced;
  451.     double C2 = 0.5*C1*C1/m_e;
  452.  
  453.     double t_elapsed = 0.0;
  454.  
  455.     std::cout << C1 << "    " << C2 << "\n";
  456.  
  457.  
  458.     // -----calculating nu-max for null-collision method starts ------------////
  459.     double nu_max = 0.0;
  460.     double nu_max_temp = 0.0;
  461.     double sigma_total = 0.0;
  462.    
  463.     for (int i = 0; i < N; i++){
  464.         // Get initial densities
  465.         double n_ground = N_He / Volume;
  466.         double n_excited1 = exc_1.size() / Volume;
  467.         double n_excited2 = exc_2.size() / Volume;
  468.  
  469.         double energy = Emin + (i + 0.5) * bin_width;
  470.  
  471.         // Total collision frequency for this energy bin
  472.         double sigma_total =
  473.             elastic_vec[i] * n_ground +
  474.             inelastic1_vec[i] * n_ground +
  475.             inelastic2_vec[i] * n_ground +
  476.             superelastic1_vec[i] * n_excited1 +
  477.             superelastic2_vec[i] * n_excited2 +
  478.             ionization_vec[i] * n_ground;
  479.  
  480.         double v = sqrt(2 * energy * q / m_e);
  481.         double nu_temp = sigma_total * v;
  482.        
  483.         if (nu_temp > nu_max) nu_max = nu_temp;
  484.     }
  485.  
  486.     std::cout << "initial nu_max: " <<nu_max << "\n";
  487.     // -----calculating nu-max for null-collision method ends ------------////    
  488.  
  489.     double dt = 0.1/nu_max;   // minimum should be 0.1/nu_max to get acceptable numerical error range see Vahedi Surrendra 1995
  490.  
  491.     std::vector<Electron> newborn;
  492.     newborn.reserve(n_e / 10);  // heuristic
  493.  
  494.     while (t_elapsed < total_time) {
  495.         // Handle edge case for final step
  496.         if (t_elapsed + dt > total_time) {
  497.             dt = total_time - t_elapsed;
  498.         }    
  499.  
  500.  
  501.         //using  null-collision technique, getting the number of particles colliding each step: P_collision = 1 - exp(-nu_max*dt)
  502.         int Ne_collided = (1.0-exp(-1.0*dt*nu_max))*n_e;  
  503.  
  504.  
  505.         newborn.clear();
  506.  
  507.         // Generate shuffled list of electron indices
  508.         std::vector<int> electron_indices(n_e);
  509.         std::iota(electron_indices.begin(), electron_indices.end(), 0); // fill with index
  510.         int reshuffle_interval = 1;
  511.  
  512.  
  513.  
  514.         for (int i = 0; i < Ne_collided; ++i) {
  515.             int j = i + std::uniform_int_distribution<int>(0, n_e - i - 1)(gen);
  516.             std::swap(electron_indices[i], electron_indices[j]);
  517.         }
  518.  
  519.         int exc1_coll_counter_temp = 0;
  520.         int super1_coll_counter_temp = 0;
  521.         int exc2_coll_counter_temp = 0;
  522.         int super2_coll_counter_temp = 0;
  523.         int null_coll_counter_temp = 0;
  524.         int ionization_counter_temp = 0;
  525.  
  526.         double energy_exc = 0.0; // calculating excitation losses each timestep
  527.         double energy_sup = 0.0; // calculating superelastic gains each timestep
  528.         double energy_Efield = 0.0; // calculating field gains/losses each timestep
  529.  
  530.  
  531.         // std::cout << "Progress: " << (t_elapsed/total_time)*100 << "\n";
  532.  
  533.         // setting flags to false each timestep
  534.         for (auto& e : electrons) e.collided_en = false;
  535.         for (auto& e : electrons) e.collided_ee = false;        
  536.  
  537.         int collision_counter_en = 0; // electron-neutral collision counter
  538.         int collision_counter_ee = 0; // e-e collisoin counter
  539.  
  540.         /// -- electrin field heating along E-Z axis begin--- /// -- first half!!!
  541.         for (int idx : electron_indices) {        
  542.             energy_Efield += ( C1*electrons[idx].vz*dt + C2*dt*dt )/q; // dividing by q to get eV            
  543.             // Update velocity component due to electric field
  544.             // double a_z = ((-1.0)*q * E_reduced) / m_e; // acceleration in z-direction, m/s^2
  545.             electrons[idx].vz += a_z * (dt); // only half timestep
  546.  
  547.             // Recalculate energy from updated velocity
  548.             double vx = electrons[idx].vx;
  549.             double vy = electrons[idx].vy;
  550.             double vz = electrons[idx].vz;
  551.             electrons[idx].energy = (vx*vx + vy*vy + vz*vz)*charge_mass_ratio;
  552.         }
  553.         // -------------------------------------------- filed heating ends ------------------------//  
  554.  
  555.  
  556.         for (int ii = 0; ii < Ne_collided; ++ii) {
  557.             int idx = electron_indices[ii];
  558.             // if (collision_counter_en >= Ne_collided) break; // quit if reached all collisions
  559.  
  560.             Electron& e = electrons[idx];
  561.             if (e.collided_en) continue;  // Skip already collided electrons
  562.  
  563.             double dens_neutrals = (N_He/Volume);
  564.             double dens_exc_1 = (exc_1.size()/Volume);
  565.             double dens_exc_2 = (exc_2.size()/Volume);
  566.             double speed = sqrt_charge_mass*sqrt(e.energy);
  567.  
  568.             int bin_energy = static_cast<int>(e.energy / bin_width);
  569.             double nu_elastic = dens_neutrals * elastic_vec[bin_energy] * speed;
  570.             double nu_inelastic1 = dens_neutrals * inelastic1_vec[bin_energy] * speed;
  571.             double nu_superelastic1 = dens_exc_1 * superelastic1_vec[bin_energy] * speed;
  572.             double nu_inelastic2 = dens_neutrals * inelastic2_vec[bin_energy] * speed;
  573.             double nu_superelastic2 = dens_exc_2 * superelastic2_vec[bin_energy] * speed;
  574.             double nu_ionization = dens_neutrals * ionization_vec[bin_energy] * speed;
  575.  
  576.             double r = dis(gen);
  577.  
  578.             double P0 = nu_elastic/nu_max;
  579.             double P1 = (nu_elastic + nu_inelastic1)/nu_max;
  580.             double P2 = (nu_elastic + nu_inelastic1 + nu_superelastic1)/nu_max;
  581.             double P3 = (nu_elastic + nu_inelastic1 + nu_superelastic1 + nu_inelastic2)/nu_max;
  582.             double P4 = (nu_elastic + nu_inelastic1 + nu_superelastic1 + nu_inelastic2 + nu_superelastic2)/nu_max;
  583.             double P5 = (nu_elastic + nu_inelastic1 + nu_superelastic1 + nu_inelastic2 + nu_superelastic2 + nu_ionization)/nu_max;
  584.            
  585.             if (r < P0) {
  586.  
  587.                 // elastic collision happens
  588.  
  589.                 // ----   Collision energy redistribution module
  590.  
  591.                 // electron particle X Y Z initial velocities and energy
  592.                 double V0_x = e.vx;
  593.                 double V0_y = e.vy;
  594.                 double V0_z = e.vz;
  595.                 double E_0 = e.energy;
  596.                
  597.                 // neutral that collides with electron
  598.  
  599.                 // randomize particles each collision
  600.  
  601.                 NeutralParticle tmp_neutral;
  602.                 tmp_neutral.initialize(gen, dis, maxwell_neutral);
  603.                 double V_x_n = tmp_neutral.vx;
  604.                 double V_y_n = tmp_neutral.vy;
  605.                 double V_z_n = tmp_neutral.vz;
  606.                 double E_n = tmp_neutral.energy;
  607.  
  608.  
  609.                 double V0 = sqrt(V0_x*V0_x + V0_y*V0_y + V0_z*V0_z);
  610.  
  611.                 // initial relative velocity
  612.                 double V_rel_x = V0_x - V_x_n;
  613.                 double V_rel_y = V0_y - V_y_n;
  614.                 double V_rel_z = V0_z - V_z_n;
  615.  
  616.                 double V_rel = sqrt(V_rel_x*V_rel_x + V_rel_y*V_rel_y + V_rel_z*V_rel_z);
  617.                 if (std::isnan(V_rel) || std::isinf(V_rel) || fabs(V_rel) < 1e-12) {
  618.                     std::cerr << "Invalid V_rel copmuted: " << V_rel << " at timestep " << t_elapsed << std::endl;
  619.                 }
  620.                 //initial relative velocity unit vectors components
  621.                 double e1 = V_rel_x/V_rel;
  622.                 double e2 = V_rel_y/V_rel;
  623.                 double e3 = V_rel_z/V_rel;
  624.  
  625.                 double C = -V_rel/(1.0+(m_e/M_n));
  626.                 double beta = sqrt(1.0 - thresh0/e.energy);
  627.                 // generating random variables to calculate random direction of center-of-mass after the collision
  628.  
  629.                 double R1 = dis(gen);
  630.                 double R2 = dis(gen);
  631.  
  632.                 //// calculating spherical angles for center-of-mass random direction
  633.                 // double theta = acos(1.0- 2.0*R1);
  634.                 // double phi = 2*M_PI*R2;
  635.                 double cos_khi = 1.0 - 2.0*R1;
  636.                 double sin_khi = sqrt(1.0 - cos_khi*cos_khi);
  637.                 if (std::isnan(sin_khi) || std::isinf(sin_khi)) {
  638.                     std::cerr << "Invalid sin_khi copmuted: " << sin_khi << " at timestep " << t_elapsed << std::endl;
  639.                 }                
  640.  
  641.                 double phi = 2.0*M_PI*R2;
  642.  
  643.                 //calculating velocity change of electron
  644.  
  645.                 double dv_x = C*( (1.0-beta*cos_khi)*e1 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(-e2) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e1*e3)  );
  646.                 double dv_y = C*( (1.0-beta*cos_khi)*e2 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(e1) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e2*e3)  );
  647.                 double dv_z = C*( (1.0-beta*cos_khi)*e3 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(0.0) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(e1*e1 + e2*e2)  );
  648.  
  649.                 //updating electron energy and velocities  
  650.  
  651.                 //updating electron energy and velocities
  652.                
  653.                 e.vx += dv_x;
  654.                 e.vy += dv_y;
  655.                 e.vz += dv_z;
  656.  
  657.                 e.energy = (e.vx*e.vx + e.vy*e.vy + e.vz*e.vz) * charge_mass_ratio;
  658.  
  659.                 collision_counter_en++;
  660.                 el_coll_counter++;
  661.  
  662.                 e.collided_en = true;
  663.             }      
  664.  
  665.             else if (r < P1) {
  666.                 if (e.energy < thresh1) {
  667.                     null_coll_counter++;
  668.                     collision_counter_en++;
  669.                     e.collided_en = true;
  670.                     continue;
  671.                 }
  672.                 else {                
  673.                     //inelastic 1(triplet) collision happens
  674.  
  675.                     // ----   Collision energy redistribution module
  676.    
  677.                     // electron particle X Y Z initial velocities and energy
  678.                     double V0_x = e.vx;
  679.                     double V0_y = e.vy;
  680.                     double V0_z = e.vz;
  681.                     double E_0 = e.energy;
  682.                    
  683.                     // neutral that collides with electron
  684.  
  685.                     // randomize particles each collision
  686.  
  687.                     NeutralParticle tmp_neutral;
  688.                     tmp_neutral.initialize(gen, dis, maxwell_neutral);
  689.                     double V_x_n = tmp_neutral.vx;
  690.                     double V_y_n = tmp_neutral.vy;
  691.                     double V_z_n = tmp_neutral.vz;
  692.                     double E_n = tmp_neutral.energy;
  693.  
  694.    
  695.                     double V0 = sqrt(V0_x*V0_x + V0_y*V0_y + V0_z*V0_z);
  696.  
  697.                     // initial relative velocity
  698.                     double V_rel_x = V0_x - V_x_n;
  699.                     double V_rel_y = V0_y - V_y_n;
  700.                     double V_rel_z = V0_z - V_z_n;
  701.  
  702.                     double V_rel = sqrt(V_rel_x*V_rel_x + V_rel_y*V_rel_y + V_rel_z*V_rel_z);
  703.                    
  704.                     if (std::isnan(V_rel) || std::isinf(V_rel) || fabs(V_rel) < 1e-12) {
  705.                         std::cerr << "Invalid V_rel copmuted: " << V_rel << " at timestep " << t_elapsed << std::endl;
  706.                     }
  707.  
  708.                     //initial relative velocity unit vectors components
  709.                     double e1 = V_rel_x/V_rel;
  710.                     double e2 = V_rel_y/V_rel;
  711.                     double e3 = V_rel_z/V_rel;
  712.  
  713.                     double C = -V_rel/(1.0+(m_e/M_n));
  714.                     double beta = sqrt(abs(1.0 - thresh1/e.energy));
  715.                     // generating random variables to calculate random direction of center-of-mass after the collision
  716.    
  717.                     double R1 = dis(gen);
  718.                     double R2 = dis(gen);
  719.    
  720.                     //// calculating spherical angles for center-of-mass random direction
  721.                     // double theta = acos(1.0- 2.0*R1);
  722.                     // double phi = 2*M_PI*R2;
  723.                     double cos_khi = 1.0 - 2.0*R1;
  724.                     double sin_khi = sqrt(1.0 - cos_khi*cos_khi);
  725.                     if (std::isnan(sin_khi) || std::isinf(sin_khi)) {
  726.                         std::cerr << "Invalid sin_khi copmuted: " << sin_khi << " at timestep " << t_elapsed << std::endl;
  727.                     }      
  728.                     double phi = 2.0*M_PI*R2;
  729.  
  730.                     //calculating velocity change of electron
  731.  
  732.                     double dv_x = C*( (1.0-beta*cos_khi)*e1 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(-e2) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e1*e3)  );
  733.                     double dv_y = C*( (1.0-beta*cos_khi)*e2 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(e1) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e2*e3)  );
  734.                     double dv_z = C*( (1.0-beta*cos_khi)*e3 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(0.0) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(e1*e1 + e2*e2)  );
  735.  
  736.                     //updating electron energy and velocities        
  737.  
  738.                     e.vx += dv_x;
  739.                     e.vy += dv_y;
  740.                     e.vz += dv_z;
  741.  
  742.                     e.energy = (e.vx*e.vx + e.vy*e.vy + e.vz*e.vz) * charge_mass_ratio;
  743.  
  744.                     collision_counter_en++;  
  745.                     exc1_coll_counter++;
  746.                     exc1_coll_counter_temp++;
  747.  
  748.                     e.collided_en = true;
  749.  
  750.                 }
  751.             }    
  752.  
  753.             else if (r < P2) {
  754.  
  755.                 //superelastic 1(triplet -> ground state) collision happens
  756.  
  757.                 if (exc_1.empty()) {
  758.                     null_coll_counter++;
  759.                     collision_counter_en++;
  760.                     e.collided_en = true;
  761.                     continue;            
  762.                 }
  763.  
  764.                 // ----   Collision energy redistribution module
  765.  
  766.                 // electron particle X Y Z initial velocities and energy
  767.                 double V0_x = e.vx;
  768.                 double V0_y = e.vy;
  769.                 double V0_z = e.vz;
  770.                 double E_0 = e.energy;
  771.  
  772.                 // neutral that collides with electron
  773.                 // taking particles from dynamic array of excited neutrals
  774.  
  775.                 int index = std::uniform_int_distribution<int>(0, exc_1.size()-1)(gen);
  776.                 Excited_neutral& exc = exc_1[index];
  777.                 double V_x_n = exc.vx;
  778.                 double V_y_n = exc.vy;
  779.                 double V_z_n = exc.vz;
  780.                 double E_n = exc.energy;
  781.                
  782.                 double V0 = sqrt(V0_x*V0_x + V0_y*V0_y + V0_z*V0_z);
  783.  
  784.                 // initial relative velocity
  785.                 double V_rel_x = V0_x - V_x_n;
  786.                 double V_rel_y = V0_y - V_y_n;
  787.                 double V_rel_z = V0_z - V_z_n;
  788.  
  789.                 double V_rel = sqrt(V_rel_x*V_rel_x + V_rel_y*V_rel_y + V_rel_z*V_rel_z);
  790.  
  791.                 if (std::isnan(V_rel) || std::isinf(V_rel) || fabs(V_rel) < 1e-12) {
  792.                     std::cerr << "Invalid V_rel copmuted: " << V_rel << " at timestep " << t_elapsed << std::endl;
  793.                 }
  794.  
  795.                 //initial relative velocity unit vectors components
  796.                 double e1 = V_rel_x/V_rel;
  797.                 double e2 = V_rel_y/V_rel;
  798.                 double e3 = V_rel_z/V_rel;
  799.  
  800.                 double C = -V_rel/(1.0+(m_e/M_n));
  801.                 double beta = sqrt(1.0 + thresh1/e.energy);
  802.                 // generating random variables to calculate random direction of center-of-mass after the collision
  803.  
  804.                 double R1 = dis(gen);
  805.                 double R2 = dis(gen);
  806.  
  807.                 //// calculating spherical angles for center-of-mass random direction
  808.                 // double theta = acos(1.0- 2.0*R1);
  809.                 // double phi = 2*M_PI*R2;
  810.                 double cos_khi = 1.0 - 2.0*R1;
  811.                 double sin_khi = sqrt(1.0 - cos_khi*cos_khi);
  812.                 if (std::isnan(sin_khi) || std::isinf(sin_khi)) {
  813.                     std::cerr << "Invalid sin_khi copmuted: " << sin_khi << " at timestep " << t_elapsed << std::endl;
  814.                 }      
  815.                 double phi = 2.0*M_PI*R2;
  816.  
  817.                 //calculating velocity change of electron
  818.  
  819.                 double dv_x = C*( (1.0-beta*cos_khi)*e1 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(-e2) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e1*e3)  );
  820.                 double dv_y = C*( (1.0-beta*cos_khi)*e2 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(e1) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e2*e3)  );
  821.                 double dv_z = C*( (1.0-beta*cos_khi)*e3 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(0.0) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(e1*e1 + e2*e2)  );
  822.  
  823.                 //updating electron energy and velocities        
  824.                
  825.                 e.vx += dv_x;
  826.                 e.vy += dv_y;
  827.                 e.vz += dv_z;
  828.  
  829.                 e.energy = (e.vx*e.vx + e.vy*e.vy + e.vz*e.vz) * charge_mass_ratio;
  830.  
  831.                 // //counting collisions, working with flags, popping atom out of the vector
  832.                 // if (!exc_1.empty()) {
  833.                 //     std::swap(exc_1[index], exc_1.back());
  834.                 //     exc_1.pop_back();
  835.                 //     N_He++;
  836.                 // }
  837.                 collision_counter_en++;  
  838.                 super1_coll_counter++;
  839.                 super1_coll_counter_temp++;
  840.                 energy_sup += thresh1;
  841.  
  842.                 e.collided_en = true;
  843.             }  
  844.  
  845.             else if (r < P3) {
  846.  
  847.                 if (e.energy < thresh1) {
  848.                     null_coll_counter++;
  849.                     collision_counter_en++;
  850.                     e.collided_en = true;
  851.                     continue;
  852.                 }
  853.                 else {
  854.  
  855.                     //inelastic 1(singlet) excitation collision happens
  856.  
  857.                     // ----   Collision energy redistribution module
  858.    
  859.                     // electron particle X Y Z initial velocities and energy
  860.                     double V0_x = e.vx;
  861.                     double V0_y = e.vy;
  862.                     double V0_z = e.vz;
  863.                     double E_0 = e.energy;
  864.                    
  865.                     // neutral that collides with electron
  866.  
  867.                     // randomize particles each collision
  868.  
  869.                     NeutralParticle tmp_neutral;
  870.                     tmp_neutral.initialize(gen, dis, maxwell_neutral);
  871.                     double V_x_n = tmp_neutral.vx;
  872.                     double V_y_n = tmp_neutral.vy;
  873.                     double V_z_n = tmp_neutral.vz;
  874.                     double E_n = tmp_neutral.energy;
  875.  
  876.    
  877.                     double V0 = sqrt(V0_x*V0_x + V0_y*V0_y + V0_z*V0_z);
  878.  
  879.                     // initial relative velocity
  880.                     double V_rel_x = V0_x - V_x_n;
  881.                     double V_rel_y = V0_y - V_y_n;
  882.                     double V_rel_z = V0_z - V_z_n;
  883.  
  884.                     double V_rel = sqrt(V_rel_x*V_rel_x + V_rel_y*V_rel_y + V_rel_z*V_rel_z);
  885.  
  886.                     if (std::isnan(V_rel) || std::isinf(V_rel) || fabs(V_rel) < 1e-12) {
  887.                         std::cerr << "Invalid V_rel copmuted: " << V_rel << " at timestep " << t_elapsed << std::endl;
  888.                     }
  889.  
  890.                     //initial relative velocity unit vectors components
  891.                     double e1 = V_rel_x/V_rel;
  892.                     double e2 = V_rel_y/V_rel;
  893.                     double e3 = V_rel_z/V_rel;
  894.  
  895.                     double C = -V_rel/(1.0+(m_e/M_n));
  896.                     double beta = sqrt(abs(1.0 - thresh2/e.energy));
  897.                     // generating random variables to calculate random direction of center-of-mass after the collision
  898.    
  899.                     double R1 = dis(gen);
  900.                     double R2 = dis(gen);
  901.    
  902.                     //// calculating spherical angles for center-of-mass random direction
  903.                     // double theta = acos(1.0- 2.0*R1);
  904.                     // double phi = 2*M_PI*R2;
  905.                     double cos_khi = 1.0 - 2.0*R1;
  906.                     double sin_khi = sqrt(1.0 - cos_khi*cos_khi);
  907.                     if (std::isnan(sin_khi) || std::isinf(sin_khi)) {
  908.                         std::cerr << "Invalid sin_khi copmuted: " << sin_khi << " at timestep " << t_elapsed << std::endl;
  909.                     }      
  910.                     double phi = 2.0*M_PI*R2;
  911.  
  912.                     //calculating velocity change of electron
  913.  
  914.                     double dv_x = C*( (1.0-beta*cos_khi)*e1 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(-e2) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e1*e3)  );
  915.                     double dv_y = C*( (1.0-beta*cos_khi)*e2 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(e1) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e2*e3)  );
  916.                     double dv_z = C*( (1.0-beta*cos_khi)*e3 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(0.0) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(e1*e1 + e2*e2)  );
  917.  
  918.                     //updating electron energy and velocities        
  919.  
  920.                     e.vx += dv_x;
  921.                     e.vy += dv_y;
  922.                     e.vz += dv_z;
  923.  
  924.                     e.energy = (e.vx*e.vx + e.vy*e.vy + e.vz*e.vz) * charge_mass_ratio;
  925.  
  926.                     collision_counter_en++;  
  927.                     exc2_coll_counter++;
  928.                     exc2_coll_counter_temp++;
  929.  
  930.                     e.collided_en = true;
  931.  
  932.                 }
  933.             }
  934.  
  935.             else if (r < P4) {
  936.  
  937.                 //supernelastic 1(singlet -> ground state) collision happens
  938.  
  939.                 if (exc_2.empty()) {
  940.                     null_coll_counter++;
  941.                     collision_counter_en++;
  942.                     e.collided_en = true;
  943.                     continue;            
  944.                 }
  945.  
  946.                 // ----   Collision energy redistribution module
  947.  
  948.                 // electron particle X Y Z initial velocities and energy
  949.                 double V0_x = e.vx;
  950.                 double V0_y = e.vy;
  951.                 double V0_z = e.vz;
  952.                 double E_0 = e.energy;
  953.  
  954.                 // neutral that collides with electron
  955.                 // taking particles from dynamic array of excited neutrals
  956.  
  957.                 int index = std::uniform_int_distribution<int>(0, exc_2.size()-1)(gen);
  958.                 Excited_neutral& exc = exc_2[index];
  959.                 double V_x_n = exc.vx;
  960.                 double V_y_n = exc.vy;
  961.                 double V_z_n = exc.vz;
  962.                 double E_n = exc.energy;
  963.                
  964.                 double V0 = sqrt(V0_x*V0_x + V0_y*V0_y + V0_z*V0_z);
  965.  
  966.                 // initial relative velocity
  967.                 double V_rel_x = V0_x - V_x_n;
  968.                 double V_rel_y = V0_y - V_y_n;
  969.                 double V_rel_z = V0_z - V_z_n;
  970.  
  971.                 double V_rel = sqrt(V_rel_x*V_rel_x + V_rel_y*V_rel_y + V_rel_z*V_rel_z);
  972.  
  973.                 if (std::isnan(V_rel) || std::isinf(V_rel) || fabs(V_rel) < 1e-12) {
  974.                     std::cerr << "Invalid V_rel copmuted: " << V_rel << " at timestep " << t_elapsed << std::endl;
  975.                 }
  976.  
  977.                 //initial relative velocity unit vectors components
  978.                 double e1 = V_rel_x/V_rel;
  979.                 double e2 = V_rel_y/V_rel;
  980.                 double e3 = V_rel_z/V_rel;
  981.  
  982.                 double C = -V_rel/(1.0+(m_e/M_n));
  983.                 double beta = sqrt(1.0 + thresh2/e.energy);
  984.                 // generating random variables to calculate random direction of center-of-mass after the collision
  985.  
  986.                 double R1 = dis(gen);
  987.                 double R2 = dis(gen);
  988.  
  989.                 //// calculating spherical angles for center-of-mass random direction
  990.                 // double theta = acos(1.0- 2.0*R1);
  991.                 // double phi = 2*M_PI*R2;
  992.                 double cos_khi = 1.0 - 2.0*R1;
  993.                 double sin_khi = sqrt(1.0 - cos_khi*cos_khi);
  994.                 if (std::isnan(sin_khi) || std::isinf(sin_khi)) {
  995.                     std::cerr << "Invalid sin_khi copmuted: " << sin_khi << " at timestep " << t_elapsed << std::endl;
  996.                 }      
  997.                 double phi = 2.0*M_PI*R2;
  998.  
  999.                 //calculating velocity change of electron
  1000.  
  1001.                 double dv_x = C*( (1.0-beta*cos_khi)*e1 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(-e2) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e1*e3)  );
  1002.                 double dv_y = C*( (1.0-beta*cos_khi)*e2 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(e1) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e2*e3)  );
  1003.                 double dv_z = C*( (1.0-beta*cos_khi)*e3 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(0.0) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(e1*e1 + e2*e2)  );
  1004.  
  1005.                 //updating electron energy and velocities        
  1006.                
  1007.                 e.vx += dv_x;
  1008.                 e.vy += dv_y;
  1009.                 e.vz += dv_z;
  1010.  
  1011.                 e.energy = (e.vx*e.vx + e.vy*e.vy + e.vz*e.vz) * charge_mass_ratio;
  1012.                 //counting collisions, working with flags, popping atom out of the vector
  1013.  
  1014.                 // if (!exc_2.empty()) {
  1015.                 //     std::swap(exc_2[index], exc_2.back());
  1016.                 //     exc_2.pop_back();
  1017.                 //     N_He++;
  1018.                 // }
  1019.  
  1020.                 collision_counter_en++;  
  1021.                 super2_coll_counter++;
  1022.                 super2_coll_counter_temp++;
  1023.                 energy_sup += thresh2;
  1024.  
  1025.                 e.collided_en = true;
  1026.             }              
  1027.  
  1028.             else if (r < P5) {
  1029.  
  1030.                 //ionization e-n collision happens
  1031.  
  1032.                 if (e.energy < thresh3) {
  1033.                     null_coll_counter++;
  1034.                     collision_counter_en++;
  1035.                     e.collided_en = true;
  1036.                     continue;
  1037.                 }              
  1038.                 else {
  1039.                     // ----   Collision energy redistribution module
  1040.    
  1041.                     // electron particle X Y Z initial velocities and energy
  1042.                     double V0_x = e.vx;
  1043.                     double V0_y = e.vy;
  1044.                     double V0_z = e.vz;
  1045.                     double E_0 = e.energy;
  1046.                    
  1047.                     // neutral that collides with electron
  1048.  
  1049.                     // randomize particles each collision
  1050.  
  1051.                     NeutralParticle tmp_neutral;
  1052.                     tmp_neutral.initialize(gen, dis, maxwell_neutral);
  1053.                     double V_x_n = tmp_neutral.vx;
  1054.                     double V_y_n = tmp_neutral.vy;
  1055.                     double V_z_n = tmp_neutral.vz;
  1056.                     double E_n = tmp_neutral.energy;
  1057.  
  1058.    
  1059.                     double V0 = sqrt(V0_x*V0_x + V0_y*V0_y + V0_z*V0_z);
  1060.  
  1061.                     // initial relative velocity
  1062.                     double V_rel_x = V0_x - V_x_n;
  1063.                     double V_rel_y = V0_y - V_y_n;
  1064.                     double V_rel_z = V0_z - V_z_n;
  1065.  
  1066.                     double V_rel = sqrt(V_rel_x*V_rel_x + V_rel_y*V_rel_y + V_rel_z*V_rel_z);
  1067.                    
  1068.                     if (std::isnan(V_rel) || std::isinf(V_rel) || fabs(V_rel) < 1e-12) {
  1069.                         std::cerr << "Invalid V_rel copmuted: " << V_rel << " at timestep " << t_elapsed << std::endl;
  1070.                     }
  1071.  
  1072.                     //initial relative velocity unit vectors components
  1073.                     double e1 = V_rel_x/V_rel;
  1074.                     double e2 = V_rel_y/V_rel;
  1075.                     double e3 = V_rel_z/V_rel;
  1076.  
  1077.                     double C = -V_rel/(1.0+(m_e/M_n));
  1078.                     double beta = sqrt(abs(1.0 - thresh3/e.energy));
  1079.                     // generating random variables to calculate random direction of center-of-mass after the collision
  1080.    
  1081.                     double R1 = dis(gen);
  1082.                     double R2 = dis(gen);
  1083.    
  1084.                     //// calculating spherical angles for center-of-mass random direction
  1085.                     // double theta = acos(1.0- 2.0*R1);
  1086.                     // double phi = 2*M_PI*R2;
  1087.                     double cos_khi = 1.0 - 2.0*R1;
  1088.                     double sin_khi = sqrt(1.0 - cos_khi*cos_khi);
  1089.                     if (std::isnan(sin_khi) || std::isinf(sin_khi)) {
  1090.                         std::cerr << "Invalid sin_khi copmuted: " << sin_khi << " at timestep " << t_elapsed << std::endl;
  1091.                     }      
  1092.                     double phi = 2.0*M_PI*R2;
  1093.  
  1094.                     //calculating velocity change of electron
  1095.  
  1096.                     double dv_x = C*( (1.0-beta*cos_khi)*e1 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(-e2) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e1*e3)  );
  1097.                     double dv_y = C*( (1.0-beta*cos_khi)*e2 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(e1) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(-e2*e3)  );
  1098.                     double dv_z = C*( (1.0-beta*cos_khi)*e3 + (beta*sin_khi*sin(phi)/sqrt(e1*e1+e2*e2))*(0.0) + (beta*sin_khi*cos(phi)/sqrt(e1*e1+e2*e2))*(e1*e1 + e2*e2)  );
  1099.  
  1100.                     //updating electron energy and velocities        
  1101.  
  1102.                     //new electron emerges, and two electrons shares the energy and obtain an equal directions (no physical background for this assumption,
  1103.                     // should think what to do with it - maybe regenerate velocities)
  1104.                    
  1105.                     //scatterd electron velocity rescaled by a factor of 2
  1106.                     e.vx = (e.vx + dv_x)/sqrt(2.0);
  1107.                     e.vy = (e.vy + dv_y)/sqrt(2.0);
  1108.                     e.vz = (e.vz + dv_z)/sqrt(2.0);
  1109.                     e.energy = (e.vx*e.vx + e.vy*e.vy + e.vz*e.vz) * charge_mass_ratio;
  1110.                     e.collided_en = true;
  1111.  
  1112.                     //ejected electron handling
  1113.  
  1114.                     // velocity angles in spherical coordinates
  1115.                     double phi_ej = 2*M_PI*fast_uniform(gen);
  1116.                     double cosTheta_ej = 2.0*fast_uniform(gen) - 1.0;
  1117.                     double sinTheta_ej = sqrt(1.0 - cosTheta_ej*cosTheta_ej);
  1118.  
  1119.                        
  1120.                     double energy_ej = e.energy; // ejected electron energy is the same
  1121.                        
  1122.                     double speed = sqrt(2*energy_ej*q/m_e);
  1123.                     double vx_ej = speed * sinTheta_ej * cos(phi_ej);
  1124.                     double vy_ej = speed * sinTheta_ej * sin(phi_ej);
  1125.                     double vz_ej = speed * cosTheta_ej;
  1126.                    
  1127.                     Electron ejected_e{vx_ej, vy_ej, vz_ej, energy_ej, false, false};
  1128.                     // electrons.push_back(ejected_e);
  1129.                    
  1130.                     newborn.push_back(ejected_e);
  1131.  
  1132.  
  1133.                     collision_counter_en++;  
  1134.                     ionization_counter++;
  1135.                     ionization_counter_temp++;
  1136.  
  1137.                 }
  1138.             }
  1139.  
  1140.             else {
  1141.                 // null-collision
  1142.                 collision_counter_en++;
  1143.                 null_coll_counter++;
  1144.                 e.collided_en = true;
  1145.             }
  1146.         }
  1147.  
  1148.         electrons.insert(electrons.end(), newborn.begin(), newborn.end());
  1149.         n_e += newborn.size();
  1150.  
  1151.         t_elapsed += dt; // Advance time
  1152.  
  1153.         // Recalculate nu_max periodically (e.g., every 100 steps)
  1154.         static int recalc_counter = 0;
  1155.         if (++recalc_counter >= 1000) {
  1156.            
  1157.             recalc_counter = 0;
  1158.  
  1159.             // Recalculate nu_max with CURRENT densities
  1160.             nu_max = 0.0;
  1161.             for (int i = 0; i < N; i++) {
  1162.                 double energy = Emin + (i + 0.5) * bin_width;
  1163.                
  1164.                 // Get current densities
  1165.                 double n_ground = N_He / Volume;
  1166.                 double n_excited1 = exc_1.size() / Volume;
  1167.                 double n_excited2 = exc_2.size() / Volume;
  1168.                
  1169.                 // Total collision frequency for this energy bin
  1170.                 double sigma_total =
  1171.                     elastic_vec[i] * n_ground +
  1172.                     inelastic1_vec[i] * n_ground +
  1173.                     inelastic2_vec[i] * n_ground +
  1174.                     superelastic1_vec[i] * n_excited1 +
  1175.                     superelastic2_vec[i] * n_excited2 +
  1176.                     ionization_vec[i] * n_ground;
  1177.  
  1178.                 double speed = sqrt_charge_mass*sqrt(energy);
  1179.                 double nu_temp = sigma_total * speed;
  1180.                
  1181.                 if (nu_temp > nu_max) nu_max = nu_temp;
  1182.             }
  1183.  
  1184.             // Update dt based on new nu_max
  1185.             dt = 0.1 / nu_max;        
  1186.         }  
  1187.  
  1188.         // calculating mean energy
  1189.         if (static_cast<int>(t_elapsed/dt)%print_interval == 0) {
  1190.             double total_energy = 0.0;
  1191.             for (const auto& e : electrons) total_energy += e.energy;
  1192.             double mean_energy = total_energy / n_e;
  1193.             file8 << t_elapsed << " " << mean_energy << "\n";            
  1194.             // file_temp << t_elapsed << " " << exc_1.size() << " " << exc_2.size() << "\n";
  1195.             std::cout << "Progress: " << (t_elapsed/total_time)*100 << "%" << "\n";
  1196.             // file16 << t_elapsed << " " << energy_Efield/n_e << " " << energy_sup/n_e << "\n";
  1197.             file17 << t_elapsed << " " << n_e << "\n";
  1198.         }        
  1199.  
  1200.     }
  1201.  
  1202.     // ----- final electron energies distribution begins
  1203.    
  1204.     for (int i = 0; i < n_e; i++){
  1205.  
  1206.         file2 << i << " " << electrons[i].energy << "\n";
  1207.  
  1208.         int bin = static_cast<int>( (electrons[i].energy - Emin)/bin_width_smooth);
  1209.         if (bin >=0 && bin < histo_maxwell.size())
  1210.             histo_maxwell[bin]++;
  1211.     }
  1212.  
  1213.     int check = 0;
  1214.     for (int i = 0; i < N_smooth; i++){
  1215.         check += histo_maxwell[i];
  1216.         double bin_center = Emin + (i + 0.5) * bin_width_smooth;
  1217.         file4 << bin_center << " " <<  static_cast<double>(histo_maxwell[i])/(sqrt(bin_center)*electrons.size()*bin_width_smooth) << "\n"; // getting f(E)
  1218.     }
  1219.  
  1220.         std::cout << "Total # of electrons in a final histogram: " << check << "\n";
  1221.         std::cout << "Final nu max: " << nu_max << "\n";
  1222.  
  1223.     // ----- final electron energies distribution ends
  1224.  
  1225.  
  1226.     file0.close(); file1.close(); file2.close(); file3.close(); file4.close(); file5.close();
  1227.     file6.close(); file7.close(); file8.close(); file9.close(); file10.close(); file11.close();
  1228.     file12.close(); file13.close(); file14.close(); file15.close(); file_temp.close(); file16.close(); file17.close();
  1229.  
  1230.     clock_t end = clock();
  1231.  
  1232.     double elapsed = (double)(end - start) / CLOCKS_PER_SEC;
  1233.  
  1234.     std::cout << "Elapsed time: " << elapsed << " seconds" << "\n";
  1235.  
  1236.     return 0;
  1237.  
  1238. }
Comments
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!