API tools faq
paste
Advertisement
phystota

superelastic+field_heating(v3)

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