API tools faq
paste
Advertisement
phystota

e-n collisions(v7) - elastic_decay_final

phystota
Apr 1st, 2025
629
0
Never
1
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 17.39 KB | None | 0 0
raw download clone embed print report
  1. #include <iostream>
  2. #include <random>
  3. #include <fstream>
  4.  
  5. #include <math.h>
  6. #include <time.h>
  7. #include <iomanip>  // For std::fixed and std::setprecision
  8.  
  9. #include <algorithm>  // For std::shuffle
  10. #include <numeric>    // For std::iota
  11.  
  12.  
  13. #define n_e 100000
  14. #define Emin 0.0
  15. #define Emax 450.0
  16. #define E_mean 100 // mean electron energy, initial distribution
  17. #define bin_width 0.1 // keep energy step ~ this to maintain cross-section clarity (Ramsauer minimum etc)
  18. #define m_e 9.1E-31 // electron mass in kg
  19. #define k_b 1.38E-23 // Boltzmann constant
  20. #define q 1.602176634E-19 // elementary charge    - eV -> J transfer param
  21. #define N ( (int)((Emax-Emin)/bin_width) ) // add 1 to include E_max if needed?
  22. #define T_n 10.0 // Helium neutral temperature in eV
  23. #define T_e 50.0    // electron Maxwell initial distribution
  24. #define M_n 6.6464731E-27 // Helium atom mass
  25. #define N_He 1000000 // Helium neutrals number
  26. #define Volume 1.0E-12 // Volume to calculate netral density and collision frequency
  27. #define time 1.0E-4 // 500 microsec time to equalibrate the system
  28.  
  29. struct Electron {
  30.  
  31.     //velocity components
  32.     double vx = 0.0;
  33.     double vy = 0.0;
  34.     double vz = 0.0;
  35.     //energy in eV
  36.     double energy = 0.0;
  37.     //Collision flag
  38.     bool collided = false;
  39.  
  40.     // initializing Maxwell-Boltzmann distribution with T_e
  41.     void initialize(std::mt19937& gen, std::uniform_real_distribution<double>& dis, std::gamma_distribution<double>& maxwell) {
  42.  
  43.         double R = dis(gen);
  44.  
  45.         // velocity angles in spherical coordinates
  46.         double phi = 2*M_PI*dis(gen);
  47.         double cosTheta = 2.0*dis(gen) - 1.0;
  48.         double sinTheta = sqrt(1.0 - cosTheta*cosTheta);
  49.  
  50.            
  51.         energy = maxwell(gen); // neutrals energies sampled as Maxwell distribution in eV
  52.            
  53.         double speed = sqrt(2*energy*q/m_e);
  54.  
  55.         //velocity components of neutrals in m/s
  56.         vx = speed * sinTheta * cos(phi);
  57.         vy = speed * sinTheta * sin(phi);
  58.         vz = speed * cosTheta;
  59.     }
  60.  
  61.  
  62. };
  63.  
  64.  
  65. struct CrossSection {
  66.     double energy;
  67.     double sigma;
  68. };
  69.  
  70. double interpolate (double energy, const std::vector<CrossSection>& elastic_CS) {
  71.  
  72.  
  73.     if (energy < elastic_CS.front().energy) {
  74.         std::cout << " required energy value lower than range of cross-section data" << "\n";
  75.         return 0.0;
  76.     }
  77.     if (energy > elastic_CS.back().energy) {
  78.         std::cout << " required energy value higher than range of cross-section data" << "\n";
  79.         return 0.0;        
  80.     }
  81.  
  82.     int step = 0;  
  83.         while (step < elastic_CS.size() && energy > elastic_CS[step].energy) {
  84.             step++;
  85.         }
  86.  
  87.     double k = (elastic_CS[step].sigma - elastic_CS[step-1].sigma)/(elastic_CS[step].energy - elastic_CS[step-1].energy);
  88.     double m = elastic_CS[step].sigma - k*elastic_CS[step].energy;
  89.    
  90.     return k*energy + m;
  91. }
  92.  
  93.  
  94. struct NeutralParticle {
  95.  
  96.     double energy = 0.0;
  97.     double vx = 0.0;
  98.     double vy = 0.0;
  99.     double vz = 0.0;
  100.  
  101.     void initialize(std::mt19937& gen, std::uniform_real_distribution<double>& dis, std::gamma_distribution<double>& maxwell) {
  102.  
  103.         double R = dis(gen);
  104.  
  105.         // velocity angles in spherical coordinates
  106.         double phi = 2*M_PI*dis(gen);
  107.         double cosTheta = 2.0*dis(gen) - 1.0;
  108.         double sinTheta = sqrt(1.0 - cosTheta*cosTheta);
  109.  
  110.            
  111.         energy = maxwell(gen); // neutrals energies sampled as Maxwell distribution in eV
  112.            
  113.         double speed = sqrt(2*energy*q/M_n);
  114.  
  115.         //velocity components of neutrals in m/s
  116.         vx = speed * sinTheta * cos(phi);
  117.         vy = speed * sinTheta * sin(phi);
  118.         vz = speed * cosTheta;
  119.     }
  120.    
  121. };
  122.  
  123.  
  124.  
  125.  
  126. int main() {
  127.  
  128.     clock_t start = clock();
  129.  
  130.     std::vector<Electron> electrons(n_e); // better to use vector instead of simple array as it's dynamically allocated (beneficial for ionization)
  131.     std::vector<NeutralParticle> neutrals(N_He);
  132.  
  133.  
  134.     std::vector<int> histo_random(N, 0); // initialize N size zero-vector for random (initial) histogram
  135.     std::vector<int> histo_maxwell(N, 0); // initialize N size zero-vector for maxwellian histogram
  136.     std::vector<int> histo_neutral(N, 0); // initialize N size zero-vector for neutral distribution histogram
  137.  
  138.     std::vector<double> elastic_vec(N, 0); // precompiled elastic cross-section-energy vector
  139.     std::vector<double> inelastic1_vec(N, 0); // precompiled inelastic(triplet excitation) cross-section-energy vector
  140.  
  141.     std::random_device rd;
  142.     std::mt19937 gen(rd());
  143.     std::uniform_real_distribution<double> dis(0.0, 1.0);
  144.     std::gamma_distribution<double> maxwell_neutral(1.5, T_n);
  145.     std::gamma_distribution<double> maxwell_electron(1.5, T_e);
  146.  
  147.     std::uniform_int_distribution<int> pair(0, n_e-1);
  148.     std::uniform_int_distribution<int> neutral_pair(0, N_He-1);    
  149.  
  150.  
  151.     std::ifstream elastic_cs_dat("cross_sections/elastic.dat");
  152.     if (!elastic_cs_dat.is_open()) {
  153.         std::cerr << "Error opening elastic cross-sections file!" << std::endl;
  154.         return 1;
  155.     }    
  156.  
  157.     std::ifstream excitation1_cs_dat("cross_sections/inelastic_triplet.dat");
  158.     if (!elastic_cs_dat.is_open()) {
  159.         std::cerr << "Error opening inelastic triplet cross-sections file!" << std::endl;
  160.         return 1;
  161.     }  
  162.  
  163.     // --- starts reading cross section datafiles
  164.  
  165.     std::vector<CrossSection> elastic_CS_temp;
  166.  
  167.     double energy, sigma;
  168.  
  169.     while (elastic_cs_dat >> energy >> sigma) {
  170.         elastic_CS_temp.push_back({energy, sigma});
  171.     }    
  172.     elastic_cs_dat.close();
  173.  
  174.     energy = 0.0;
  175.     sigma = 0.0;
  176.  
  177.     std::vector<CrossSection> inelastic1_CS_temp;
  178.  
  179.     while (excitation1_cs_dat >> energy >> sigma) {
  180.         inelastic1_CS_temp.push_back({energy, sigma});
  181.     }    
  182.     excitation1_cs_dat.close();    
  183.  
  184.     // --- finish reading cross-section datafiles  
  185.  
  186.     std::ofstream file0("velocities.dat");    
  187.     std::ofstream file1("energies.dat");        
  188.     std::ofstream file2("energies_final.dat");    
  189.     std::ofstream file3("histo_random.dat");    
  190.     file3 << std::fixed << std::setprecision(10);
  191.    
  192.     std::ofstream file4("histo_maxwell.dat");
  193.     file4 << std::fixed << std::setprecision(10);          
  194.    
  195.     std::ofstream file5("neutral_distribution.dat");    
  196.     std::ofstream file6("E*f(E).dat");    
  197.     std::ofstream file7("nu_max.dat");
  198.     std::ofstream file8("electron_mean_energy.dat");
  199.     std::ofstream file9("nu.dat");
  200.     std::ofstream file10("inelastic_CS_check.dat");        
  201.  
  202.     // Initialize all electrons
  203.     for (auto& e : electrons) {
  204.         e.initialize(gen, dis, maxwell_electron);
  205.     }
  206.     // initialize all nenutrals
  207.     for (auto&n : neutrals) {
  208.         n.initialize(gen, dis, maxwell_neutral);
  209.     }
  210.     // precalculate elastic cross-section for each energy bin
  211.     for (int i = 0; i < N; i++){
  212.         elastic_vec[i] = interpolate(bin_width*(i+0.5), elastic_CS_temp);
  213.     }
  214.     // precalculate inelastic cross-section (triplet) for each energy bin
  215.     for (int i = 0; i < N; i++){
  216.         inelastic1_vec[i] = interpolate(bin_width*(i+0.5), inelastic1_CS_temp);
  217.     }
  218.  
  219.     for (int i = 0; i < n_e; i++){
  220.         file1 << i << " " << electrons.at(i).energy << "\n";
  221.         file0 << i << " " << electrons[i].vx << " " << electrons[i].vy << " " << electrons[i].vz << "\n";
  222.     }
  223.  
  224.     // -----initial electrons energy distribution starts------------////
  225.     for (int i = 0; i < n_e; i++){
  226.         int bin = (int)( (electrons[i].energy - Emin)/bin_width );
  227.         if (bin >=0 && bin < histo_random.size())
  228.             histo_random[bin]++;
  229.     }
  230.  
  231.     for (int i = 0; i < histo_random.size(); i++){
  232.         double bin_center = Emin + (i + 0.5) * bin_width;
  233.         file3 << bin_center << " " <<  bin_center*static_cast<double>(histo_random[i])/(electrons.size()*bin_width) << "\n"; // this is electron normalized distribution function
  234.     }
  235.     // -----initial electrons energy distribution ends------------////    
  236.  
  237.     // -----neutrals Maxwell-Boltzmann distribution starts------------////
  238.     for (int i = 0; i < N_He; i++){
  239.         int bin = (int)( (neutrals[i].energy - Emin)/bin_width );
  240.         if (bin >=0 && bin < histo_neutral.size())
  241.             histo_neutral[bin]++;
  242.     }    
  243.  
  244.     for (int i = 0; i < histo_neutral.size(); i++){
  245.         double bin_center = Emin + (i + 0.5) * bin_width;
  246.         file5 << bin_center << " " << static_cast<double>(histo_neutral[i])/(neutrals.size()*bin_width) << "\n"; // this is real f(E) - normalized distribution
  247.         file6 << bin_center << " " << bin_center*static_cast<double>(histo_neutral[i])/(neutrals.size()*bin_width) << "\n"; // this should be E*f(E)
  248.  
  249.     }
  250.     // -----neutrals Maxwell-Boltzmann distribution starts------------////      
  251.  
  252.     // -----calculating nu-max for null-collision method starts ------------////
  253.     double nu_max = 0.0;
  254.     double nu_max_temp = 0.0;
  255.     double sigma_total = 0.0;
  256.    
  257.     for (int i = 0; i < N; i++){
  258.         nu_max_temp = (N_He/Volume)*elastic_vec[i] * sqrt(2.0*(i*bin_width + bin_width/2.0)*q/m_e);
  259.         file7 << i << " " << nu_max_temp << "\n";
  260.         if (nu_max_temp > nu_max)
  261.             nu_max = nu_max_temp;
  262.     }
  263.     // -----calculating nu-max for null-collision method ends ------------////
  264.  
  265.     //----- calculating nu-average from our electron distribution starts---------///
  266.     for (int i = 0; i < N; i++){
  267.         double bin_center = Emin + (i + 0.5) * bin_width;
  268.         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";
  269.     }
  270.     //----- calculating nu-average from our electron distribution ends ---------///    
  271.  
  272.     std::cout << nu_max << "\n";
  273.  
  274.     double dt = 0.1/nu_max;   // minimum should be 0.1/nu_max to get acceptable numerical error range see Vahedi Surrendra 1995
  275.     double steps = static_cast<int>(time/dt);
  276.  
  277. //    std::cout << steps << "\n";
  278.  
  279.     //using  null-collision technique, getting the number of particles colliding each step: P_collision = 1 - exp(-nu_max*dt)
  280.     int Ne_collided = (1.0-exp(-1.0*dt*nu_max))*n_e;
  281. //    int Ne_collided = n_e*0.98;  // in case I want to check smth
  282.  
  283.  
  284.     // Generate shuffled list of electron indices
  285.     std::vector<int> electron_indices(n_e);
  286.     std::iota(electron_indices.begin(), electron_indices.end(), 0); // fill with index
  287.     std::shuffle(electron_indices.begin(), electron_indices.end(), gen); // shuffle the indexes    
  288.     int reshuffle_interval = 1;
  289.     int print_interval = 100;
  290.  
  291.     for (int t = 0; t < steps; t++){
  292.         std::cout << "timestep remains: " << steps - t << "\n";
  293.  
  294.         //reshuffle the indices
  295.         if (t % reshuffle_interval == 0) {
  296.             std::shuffle(electron_indices.begin(), electron_indices.end(), gen);
  297.         }
  298.  
  299.         // setting flags to false each timestep
  300.         for (auto& e : electrons) e.collided = false;
  301.  
  302.         int collision_counter = 0;
  303.  
  304.  
  305.         for (int idx : electron_indices) {
  306.  
  307.             if (collision_counter >= Ne_collided) break; // quit if reached all collisions
  308.  
  309.             Electron& e = electrons[idx];
  310.             if (e.collided) continue;  // Skip already collided electrons
  311.  
  312.             double electron_energy = e.energy;
  313.             int bin_energy = static_cast<int>(electron_energy / bin_width);
  314.             double nu_elastic = (N_He/Volume) * elastic_vec[bin_energy] * sqrt(2.0*electron_energy*q/m_e);
  315.  
  316.             if (dis(gen) < nu_elastic/nu_max) {
  317.  
  318.                 // ----   Collision energy redistribution module
  319.  
  320.                 // electron particle X Y Z initial velocities and energy
  321.                 double V0_x_1 = e.vx;
  322.                 double V0_y_1 = e.vy;
  323.                 double V0_z_1 = e.vz;
  324.  
  325.                 // neutral particle X Y Z initial velocities
  326.  
  327.                 // int k = neutral_pair(gen);
  328.  
  329.                 // double V0_x_2 = neutrals[k].vx;
  330.                 // double V0_y_2 = neutrals[k].vy;
  331.                 // double V0_z_2 = neutrals[k].vz;
  332.  
  333.                 // randomize particles each collision
  334.                 NeutralParticle tmp_neutral;
  335.                 tmp_neutral.initialize(gen, dis, maxwell_neutral);
  336.                 double V0_x_2 = tmp_neutral.vx;
  337.                 double V0_y_2 = tmp_neutral.vy;
  338.                 double V0_z_2 = tmp_neutral.vz;
  339.  
  340.                 // initial relative velocity X Y Z (must be equal to final relative velocity in center-of-mass frame)
  341.  
  342.                 double V0_rel_x = (V0_x_1 - V0_x_2);
  343.                 double V0_rel_y = (V0_y_1 - V0_y_2);
  344.                 double V0_rel_z = (V0_z_1 - V0_z_2);
  345.  
  346.                 double V0_rel = sqrt(V0_rel_x*V0_rel_x + V0_rel_y*V0_rel_y + V0_rel_z*V0_rel_z);
  347.  
  348.                 // center-of-mass frame initial velocity (magnitude of it must be equal to the counterpart in this frame)
  349.  
  350.                 double V_cm_x = (m_e*V0_x_1 + M_n*V0_x_2)/(m_e + M_n);
  351.                 double V_cm_y = (m_e*V0_y_1 + M_n*V0_y_2)/(m_e + M_n);
  352.                 double V_cm_z = (m_e*V0_z_1 + M_n*V0_z_2)/(m_e + M_n);                    
  353.  
  354.                 // generating random variables to calculate random direction of center-of-mass after the collision
  355.  
  356.                 double R1 = dis(gen);
  357.                 double R2 = dis(gen);
  358.  
  359.                 // calculating spherical angles for center-of-mass random direction
  360.                 double theta = acos(1.0- 2.0*R1);
  361.                 double phi = 2*M_PI*R2;
  362.  
  363.                 //calculating final relative velocity with random direction
  364.  
  365.                 double V_rel_x = V0_rel*sin(theta)*cos(phi);
  366.                 double V_rel_y = V0_rel*sin(theta)*sin(phi);
  367.                 double V_rel_z = V0_rel*cos(theta);
  368.  
  369.                 double V_rel = sqrt(V_rel_x*V_rel_x + V_rel_y*V_rel_y + V_rel_z*V_rel_z);
  370.  
  371.                 //calculating final velocity of electron
  372.  
  373.                 double V_x_1 = V_cm_x + V_rel_x * (M_n/(m_e + M_n));
  374.                 double V_y_1 = V_cm_y + V_rel_y * (M_n/(m_e + M_n));
  375.                 double V_z_1 = V_cm_z + V_rel_z * (M_n/(m_e + M_n));
  376.  
  377.                 double V_1 = sqrt(V_x_1*V_x_1 + V_y_1*V_y_1 + V_z_1*V_z_1);
  378.  
  379.                 //updating electron energy and velocities
  380.                
  381.                 e.energy = m_e*V_1*V_1/(2.0*q);
  382.                 e.vx = V_x_1;
  383.                 e.vy = V_y_1;
  384.                 e.vz = V_z_1;
  385.  
  386.                 collision_counter++;
  387.  
  388.                 e.collided = true;
  389.             }            
  390.         }
  391.                 // if (t%print_interval == 0){
  392.                 // // open datafiles to write each time step to see evolution
  393.                 // std::ostringstream filename;
  394.                 // filename << "data/distribution_" << std::setw(4) << std::setfill('0') << t << ".dat";
  395.  
  396.                 // std::ofstream file(filename.str());
  397.                 // if (!file.is_open()){
  398.                 // std::cerr << "Error opening file: " << filename.str() << std::endl;
  399.                 // return 1;
  400.                 // }
  401.                 // // end opening datafiles for each timestep
  402.                
  403.                 // // creating histogram each timestep
  404.                 // for (int i = 0; i < n_e; i++){
  405.                 //     int bin = (int)( (electrons[i].energy - Emin)/bin_width );
  406.                 //     if (bin >=0 && bin < N)
  407.                 //     histo_maxwell[bin]++;
  408.                 // }
  409.  
  410.                 // // writing data each time step
  411.                 // for (int i = 0; i < N; i++){
  412.                 //     double bin_start = Emin + i*bin_width;
  413.                 //     file << i*bin_width << " " <<  static_cast<double>(histo_maxwell[i])/(electrons.size()*bin_width) << "\n"; // later need to divide by total partcles number to get normalized distribution
  414.                 //     histo_maxwell[i] = 0;
  415.                 // }
  416.  
  417.                 // file.close();
  418.                 // // end writing data each timestep
  419.  
  420.                 // calculating mean energy
  421.                 double total_energy = 0.0;
  422.                 for (const auto& e : electrons) total_energy += e.energy;
  423.                 double mean_energy = total_energy / n_e;
  424.                 file8 << t*dt << " " << mean_energy << "\n";                
  425.             }
  426.  
  427.     // ----- final electron energies distribution begins
  428.     for (int i = 0; i < n_e; i++){
  429.  
  430.         file2 << i << " " << electrons[i].energy << "\n";
  431.  
  432.         int bin = (int)( (electrons[i].energy - Emin)/bin_width );
  433.         if (bin >=0 && bin < histo_maxwell.size())
  434.             histo_maxwell[bin]++;
  435.     }
  436.  
  437.     int check = 0;
  438.     for (int i = 0; i < histo_maxwell.size(); i++){
  439.         check += histo_maxwell[i];
  440.         double bin_center = Emin + (i + 0.5) * bin_width;
  441.         file4 << bin_center << " " <<  bin_center*static_cast<double>(histo_maxwell[i])/(electrons.size()*bin_width) << "\n"; // getting f(E)*E
  442.     }
  443.     std::cout << "Total # of electrons in histo: " << check << "\n";
  444.  
  445.     // ----- final electron energies distribution begins    
  446.  
  447.  
  448.     file0.close();
  449.     file1.close();
  450.     file2.close();
  451.     file3.close();
  452.     file4.close();
  453.     file5.close();
  454.     file6.close();
  455.     file7.close();
  456.     file8.close();
  457.     file9.close();
  458.  
  459.     clock_t end = clock();
  460.  
  461.     double elapsed = (double)(end - start) / CLOCKS_PER_SEC;
  462.  
  463.     std::cout << "Ne collided each timesteps:" << Ne_collided << "\n";
  464.     std::cout << "Energies written successfuly\n";
  465.     std::cout << "Elapsed time: %f seconds " << elapsed << "\n";
  466.  
  467.  
  468.     return 0;
  469.  
  470. }
Advertisement
Comments
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!