API tools faq
paste
Advertisement
phystota

v3.1 - stable working - fixed populations - converges with Bolsig+

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