API tools faq
paste
Advertisement
phystota

e-n collisions(v9) clean

phystota
Apr 2nd, 2025
436
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 20.78 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 6.0E-4 // 500 microsec time to equalibrate the system
  28.  
  29. // handling final energy bin
  30.  
  31. #define bin_width_smooth 1.0 // energy bin for smooth final distribution
  32. #define N_smooth ( (int)((Emax-Emin)/bin_width_smooth) )
  33.  
  34. struct Electron {
  35.  
  36.     //velocity components
  37.     double vx = 0.0;
  38.     double vy = 0.0;
  39.     double vz = 0.0;
  40.     //energy in eV
  41.     double energy = 0.0;
  42.     //Collision flag
  43.     bool collided = false;
  44.  
  45.     // initializing Maxwell-Boltzmann distribution with T_e
  46.     void initialize(std::mt19937& gen, std::uniform_real_distribution<double>& dis, std::gamma_distribution<double>& maxwell) {
  47.  
  48.         double R = dis(gen);
  49.  
  50.         // velocity angles in spherical coordinates
  51.         double phi = 2*M_PI*dis(gen);
  52.         double cosTheta = 2.0*dis(gen) - 1.0;
  53.         double sinTheta = sqrt(1.0 - cosTheta*cosTheta);
  54.  
  55.            
  56.         energy = maxwell(gen); // neutrals energies sampled as Maxwell distribution in eV
  57.            
  58.         double speed = sqrt(2*energy*q/m_e);
  59.  
  60.         //velocity components of neutrals in m/s
  61.         vx = speed * sinTheta * cos(phi);
  62.         vy = speed * sinTheta * sin(phi);
  63.         vz = speed * cosTheta;
  64.     }
  65.  
  66.  
  67. };
  68.  
  69.  
  70. struct CrossSection {
  71.     double energy;
  72.     double sigma;
  73. };
  74.  
  75. double interpolate (double energy, const std::vector<CrossSection>& CS) {
  76.  
  77.  
  78.     if (energy < CS.front().energy) {
  79.         std::cout << " required energy value lower than range of cross-section data" << "\n";
  80.         return 0.0;
  81.     }
  82.     if (energy > CS.back().energy) {
  83.         std::cout << " required energy value higher than range of cross-section data" << "\n";
  84.         return 0.0;        
  85.     }
  86.  
  87.     int step = 0;  
  88.         while (step < CS.size() && energy > CS[step].energy) {
  89.             step++;
  90.         }
  91.  
  92.     double k = (CS[step].sigma - CS[step-1].sigma)/(CS[step].energy - CS[step-1].energy);
  93.     double m = CS[step].sigma - k*CS[step].energy;
  94.    
  95.     return k*energy + m;
  96. }
  97.  
  98.  
  99. struct NeutralParticle {
  100.  
  101.     double energy = 0.0;
  102.     double vx = 0.0;
  103.     double vy = 0.0;
  104.     double vz = 0.0;
  105.  
  106.     void initialize(std::mt19937& gen, std::uniform_real_distribution<double>& dis, std::gamma_distribution<double>& maxwell) {
  107.  
  108.         double R = dis(gen);
  109.  
  110.         // velocity angles in spherical coordinates
  111.         double phi = 2*M_PI*dis(gen);
  112.         double cosTheta = 2.0*dis(gen) - 1.0;
  113.         double sinTheta = sqrt(1.0 - cosTheta*cosTheta);
  114.  
  115.            
  116.         energy = maxwell(gen); // neutrals energies sampled as Maxwell distribution in eV
  117.            
  118.         double speed = sqrt(2*energy*q/M_n);
  119.  
  120.         //velocity components of neutrals in m/s
  121.         vx = speed * sinTheta * cos(phi);
  122.         vy = speed * sinTheta * sin(phi);
  123.         vz = speed * cosTheta;
  124.     }
  125.    
  126. };
  127.  
  128.  
  129.  
  130.  
  131. int main() {
  132.  
  133.     clock_t start = clock();
  134.  
  135.     std::vector<Electron> electrons(n_e); // better to use vector instead of simple array as it's dynamically allocated (beneficial for ionization)
  136.     std::vector<NeutralParticle> neutrals(N_He);
  137.  
  138.  
  139.     std::vector<int> histo_random(N, 0); // initialize N size zero-vector for random (initial) histogram
  140.     std::vector<int> histo_maxwell(N, 0); // initialize N size zero-vector for maxwellian histogram
  141.     std::vector<int> histo_neutral(N, 0); // initialize N size zero-vector for neutral distribution histogram
  142.  
  143.     std::vector<double> elastic_vec(N, 0); // precompiled elastic cross-section-energy vector
  144.     std::vector<double> inelastic1_vec(N, 0); // precompiled inelastic(triplet excitation) cross-section-energy vector
  145.  
  146.     std::random_device rd;
  147.     std::mt19937 gen(rd());
  148.     std::uniform_real_distribution<double> dis(0.0, 1.0);
  149.     std::gamma_distribution<double> maxwell_neutral(1.5, T_n);
  150.     std::gamma_distribution<double> maxwell_electron(1.5, T_e);
  151.  
  152.     std::uniform_int_distribution<int> pair(0, n_e-1);
  153.     std::uniform_int_distribution<int> neutral_pair(0, N_He-1);    
  154.  
  155.  
  156.     std::ifstream elastic_cs_dat("cross_sections/elastic.dat");
  157.     if (!elastic_cs_dat.is_open()) {
  158.         std::cerr << "Error opening elastic cross-sections file!" << std::endl;
  159.         return 1;
  160.     }    
  161.  
  162.     std::ifstream excitation1_cs_dat("cross_sections/inelastic_triplet.dat");
  163.     if (!excitation1_cs_dat.is_open()) {
  164.         std::cerr << "Error opening inelastic triplet cross-sections file!" << std::endl;
  165.         return 1;
  166.     }  
  167.  
  168.     // --- starts reading cross section datafiles
  169.  
  170.     std::vector<CrossSection> elastic_CS_temp;
  171.  
  172.     double energy, sigma;
  173.  
  174.     while (elastic_cs_dat >> energy >> sigma) {
  175.         elastic_CS_temp.push_back({energy, sigma});
  176.     }    
  177.     elastic_cs_dat.close();
  178.  
  179.     energy = 0.0;
  180.     sigma = 0.0;
  181.  
  182.     std::vector<CrossSection> inelastic1_CS_temp;
  183.  
  184.     while (excitation1_cs_dat >> energy >> sigma) {
  185.         inelastic1_CS_temp.push_back({energy, sigma});
  186.     }    
  187.     excitation1_cs_dat.close();    
  188.  
  189.     // --- finish reading cross-section datafiles  
  190.  
  191.     std::ofstream file0("output_files/velocities.dat");    
  192.     std::ofstream file1("output_files/energies.dat");        
  193.     std::ofstream file2("output_files/energies_final.dat");    
  194.     std::ofstream file3("output_files/histo_random.dat");    
  195.     file3 << std::fixed << std::setprecision(10);
  196.    
  197.     std::ofstream file4("output_files/histo_maxwell.dat");
  198.     file4 << std::fixed << std::setprecision(10);          
  199.    
  200.     std::ofstream file5("output_files/neutral_distribution.dat");    
  201.     std::ofstream file6("output_files/E*f(E).dat");    
  202.     std::ofstream file7("output_files/nu_max.dat");
  203.     std::ofstream file8("output_files/electron_mean_energy.dat");
  204.     std::ofstream file9("output_files/nu_elastic_average_initial.dat");
  205.     std::ofstream file10("output_files/nu_inelastic1_average_initial.dat");
  206.     std::ofstream file11("output_files/nu_elastic_average_final.dat");
  207.     std::ofstream file12("output_files/nu_inelastic1_average_final.dat");                    
  208.  
  209.     // Initialize all electrons
  210.     for (auto& e : electrons) {
  211.         e.initialize(gen, dis, maxwell_electron);
  212.     }
  213.     // initialize all nenutrals
  214.     for (auto&n : neutrals) {
  215.         n.initialize(gen, dis, maxwell_neutral);
  216.     }
  217.     // precalculate elastic cross-section for each energy bin
  218.     for (int i = 0; i < N; i++){
  219.         elastic_vec[i] = interpolate(bin_width*(i+0.5), elastic_CS_temp);
  220.     }
  221.     // precalculate inelastic cross-section (triplet) for each energy bin
  222.     for (int i = 0; i < N; i++){
  223.         inelastic1_vec[i] = interpolate(bin_width*(i+0.5), inelastic1_CS_temp);
  224.     }
  225.  
  226.     for (int i = 0; i < n_e; i++){
  227.         file1 << i << " " << electrons.at(i).energy << "\n";
  228.         file0 << i << " " << electrons[i].vx << " " << electrons[i].vy << " " << electrons[i].vz << "\n";
  229.     }
  230.  
  231.     // -----initial electrons energy distribution starts------------////
  232.     for (int i = 0; i < n_e; i++){
  233.         int bin = (int)( (electrons[i].energy - Emin)/bin_width );
  234.         if (bin >=0 && bin < histo_random.size())
  235.             histo_random[bin]++;
  236.     }
  237.  
  238.     for (int i = 0; i < histo_random.size(); i++){
  239.         double bin_center = Emin + (i + 0.5) * bin_width;
  240.         file3 << bin_center << " " <<  static_cast<double>(histo_random[i])/(electrons.size()*bin_width) << "\n"; // this is electron normalized distribution function
  241.     }
  242.     // -----initial electrons energy distribution ends------------////    
  243.  
  244.     // -----neutrals Maxwell-Boltzmann distribution starts------------////
  245.     for (int i = 0; i < N_He; i++){
  246.         int bin = (int)( (neutrals[i].energy - Emin)/bin_width );
  247.         if (bin >=0 && bin < histo_neutral.size())
  248.             histo_neutral[bin]++;
  249.     }    
  250.  
  251.     for (int i = 0; i < histo_neutral.size(); i++){
  252.         double bin_center = Emin + (i + 0.5) * bin_width;
  253.         file5 << bin_center << " " << static_cast<double>(histo_neutral[i])/(neutrals.size()*bin_width) << "\n"; // this is real f(E) - normalized distribution
  254.         file6 << bin_center << " " << bin_center*static_cast<double>(histo_neutral[i])/(neutrals.size()*bin_width) << "\n"; // this should be E*f(E)
  255.  
  256.     }
  257.     // -----neutrals Maxwell-Boltzmann distribution starts------------////      
  258.  
  259.     // -----calculating nu-max for null-collision method starts ------------////
  260.     double nu_max = 0.0;
  261.     double nu_max_temp = 0.0;
  262.     double sigma_total = 0.0;
  263.    
  264.     for (int i = 0; i < N; i++){
  265.         sigma_total = elastic_vec[i] + inelastic1_vec[i];
  266.         nu_max_temp = (N_He/Volume)*sigma_total * sqrt(2.0*(i*bin_width + bin_width/2.0)*q/m_e);
  267.         file7 << i << " " << nu_max_temp << "\n";
  268.         if (nu_max_temp > nu_max)
  269.             nu_max = nu_max_temp;
  270.     }
  271.     // -----calculating nu-max for null-collision method ends ------------////
  272.  
  273.     //----- calculating number to calculate nu-average (both elastic/inelastic )from our electron distribution starts---------///
  274.     // --- calculating nu(E)*f(E) for later external integration, using initial f(E)
  275.     for (int i = 0; i < N; i++){
  276.         double bin_center = Emin + (i + 0.5) * bin_width;
  277.         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";
  278.         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";
  279.     }
  280.     //----- calculating nu-average from our electron distribution ends ---------///    
  281.  
  282.     std::cout << nu_max << "\n";
  283.  
  284.     double dt = 0.1/nu_max;   // minimum should be 0.1/nu_max to get acceptable numerical error range see Vahedi Surrendra 1995
  285.     double steps = static_cast<int>(time/dt);
  286.  
  287. //    std::cout << steps << "\n";
  288.  
  289.     //using  null-collision technique, getting the number of particles colliding each step: P_collision = 1 - exp(-nu_max*dt)
  290.     int Ne_collided = (1.0-exp(-1.0*dt*nu_max))*n_e;
  291. //    int Ne_collided = n_e*0.98;  // in case I want to check smth
  292.  
  293.  
  294.     // Generate shuffled list of electron indices
  295.     std::vector<int> electron_indices(n_e);
  296.     std::iota(electron_indices.begin(), electron_indices.end(), 0); // fill with index
  297.     std::shuffle(electron_indices.begin(), electron_indices.end(), gen); // shuffle the indexes    
  298.     int reshuffle_interval = 1;
  299.     int print_interval = 100;
  300.     int el_coll_counter = 0; // track all elastic collisions
  301.     int exc1_coll_counter = 0; // track all excitation collisions
  302.     int null_coll_counter = 0; // track null-collisions
  303.  
  304.     for (int t = 0; t < steps; t++){
  305.         std::cout << "timestep remains: " << steps - t << "\n";
  306.  
  307.         //reshuffle the indices
  308.         if (t % reshuffle_interval == 0) {
  309.             std::shuffle(electron_indices.begin(), electron_indices.end(), gen);
  310.         }
  311.  
  312.         // setting flags to false each timestep
  313.         for (auto& e : electrons) e.collided = false;
  314.  
  315.         int collision_counter = 0;
  316.  
  317.  
  318.         for (int idx : electron_indices) {
  319.  
  320.             if (collision_counter >= Ne_collided) break; // quit if reached all collisions
  321.  
  322.             Electron& e = electrons[idx];
  323.             if (e.collided) continue;  // Skip already collided electrons
  324.  
  325.             double electron_energy = e.energy;
  326.             int bin_energy = static_cast<int>(electron_energy / bin_width);
  327.             double nu_elastic = (N_He/Volume) * elastic_vec[bin_energy] * sqrt(2.0*electron_energy*q/m_e);
  328.             double nu_inelastic1 = (N_He/Volume) * inelastic1_vec[bin_energy] * sqrt(2.0*electron_energy*q/m_e);
  329.  
  330.             double r = dis(gen);
  331.  
  332.             if (r < nu_elastic/nu_max) {
  333.  
  334.                 // elastic collision happens
  335.  
  336.                 // ----   Collision energy redistribution module
  337.  
  338.                 // electron particle X Y Z initial velocities and energy
  339.                 double V0_x_1 = e.vx;
  340.                 double V0_y_1 = e.vy;
  341.                 double V0_z_1 = e.vz;
  342.  
  343.                 // neutral particle X Y Z initial velocities
  344.  
  345.                 // int k = neutral_pair(gen);
  346.  
  347.                 // double V0_x_2 = neutrals[k].vx;
  348.                 // double V0_y_2 = neutrals[k].vy;
  349.                 // double V0_z_2 = neutrals[k].vz;
  350.  
  351.                 // randomize particles each collision
  352.                 NeutralParticle tmp_neutral;
  353.                 tmp_neutral.initialize(gen, dis, maxwell_neutral);
  354.                 double V0_x_2 = tmp_neutral.vx;
  355.                 double V0_y_2 = tmp_neutral.vy;
  356.                 double V0_z_2 = tmp_neutral.vz;
  357.  
  358.                 // initial relative velocity X Y Z (must be equal to final relative velocity in center-of-mass frame)
  359.  
  360.                 double V0_rel_x = (V0_x_1 - V0_x_2);
  361.                 double V0_rel_y = (V0_y_1 - V0_y_2);
  362.                 double V0_rel_z = (V0_z_1 - V0_z_2);
  363.  
  364.                 double V0_rel = sqrt(V0_rel_x*V0_rel_x + V0_rel_y*V0_rel_y + V0_rel_z*V0_rel_z);
  365.  
  366.                 // center-of-mass frame initial velocity (magnitude of it must be equal to the counterpart in this frame)
  367.  
  368.                 double V_cm_x = (m_e*V0_x_1 + M_n*V0_x_2)/(m_e + M_n);
  369.                 double V_cm_y = (m_e*V0_y_1 + M_n*V0_y_2)/(m_e + M_n);
  370.                 double V_cm_z = (m_e*V0_z_1 + M_n*V0_z_2)/(m_e + M_n);                    
  371.  
  372.                 // generating random variables to calculate random direction of center-of-mass after the collision
  373.  
  374.                 double R1 = dis(gen);
  375.                 double R2 = dis(gen);
  376.  
  377.                 // calculating spherical angles for center-of-mass random direction
  378.                 double theta = acos(1.0- 2.0*R1);
  379.                 double phi = 2*M_PI*R2;
  380.  
  381.                 //calculating final relative velocity with random direction
  382.  
  383.                 double V_rel_x = V0_rel*sin(theta)*cos(phi);
  384.                 double V_rel_y = V0_rel*sin(theta)*sin(phi);
  385.                 double V_rel_z = V0_rel*cos(theta);
  386.  
  387.                 double V_rel = sqrt(V_rel_x*V_rel_x + V_rel_y*V_rel_y + V_rel_z*V_rel_z);
  388.  
  389.                 //calculating final velocity of electron
  390.  
  391.                 double V_x_1 = V_cm_x + V_rel_x * (M_n/(m_e + M_n));
  392.                 double V_y_1 = V_cm_y + V_rel_y * (M_n/(m_e + M_n));
  393.                 double V_z_1 = V_cm_z + V_rel_z * (M_n/(m_e + M_n));
  394.  
  395.                 double V_1 = sqrt(V_x_1*V_x_1 + V_y_1*V_y_1 + V_z_1*V_z_1);
  396.  
  397.                 //updating electron energy and velocities
  398.                
  399.                 e.energy = m_e*V_1*V_1/(2.0*q);
  400.                 e.vx = V_x_1;
  401.                 e.vy = V_y_1;
  402.                 e.vz = V_z_1;
  403.  
  404.                 collision_counter++;
  405.                 el_coll_counter++;
  406.  
  407.                 e.collided = true;
  408.             }        
  409.  
  410.             else if (r < (nu_elastic + nu_inelastic1)/nu_max) {
  411.  
  412.                 //inelastic 1(triplet) collision happens
  413.  
  414.                 // ----   Collision energy redistribution module
  415.  
  416.                 // electron particle X Y Z initial velocities and energy
  417.                 double V0_x = e.vx;
  418.                 double V0_y = e.vy;
  419.                 double V0_z = e.vz;
  420.                 double E_0 = e.energy;
  421.  
  422.                 double V0 = sqrt(V0_x*V0_x + V0_y*V0_y + V0_z*V0_z);
  423.                
  424.                 // generating random variables to calculate random direction of center-of-mass after the collision
  425.  
  426.                 double R1 = dis(gen);
  427.                 double R2 = dis(gen);
  428.  
  429.                 //// calculating spherical angles for center-of-mass random direction
  430.                 // double theta = acos(1.0- 2.0*R1);
  431.                 // double phi = 2*M_PI*R2;
  432.                 double cos_khi = (2.0 + E_0 - 2.0*pow((1+E_0), R1))/E_0;
  433.                 double sin_khi = sqrt(1.0 - cos_khi*cos_khi);
  434.  
  435.                 double phi = 2.0*M_PI*R2;
  436.                 double cos_theta = V0_x/V0;
  437.                 double sin_theta = sqrt(1.0 - cos_theta*cos_theta);
  438.                 // //calculating final relative velocity with random direction
  439.  
  440.                 //calculating final velocity of electron
  441.  
  442.                 double i_scat = (V0_x/V0)*cos_khi + (1.0 - (V0_x/V0)*(V0_x/V0))*(sin_khi*cos(phi)/sin_theta);
  443.                 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;
  444.                 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;
  445.  
  446.                 //updating electron energy and velocities
  447.  
  448.                 double delta_E = 19.82; //triplet excitation energy              
  449.                
  450.                 if (e.energy < delta_E) {
  451.                     null_coll_counter++;
  452.                     collision_counter++;
  453.                     e.collided = true;
  454.                     continue;
  455.                 }
  456.                 else {
  457.                     e.energy = E_0 - delta_E;
  458.  
  459.                     double speed = sqrt(2*e.energy*q/m_e);
  460.  
  461.                     e.vx = speed*i_scat;
  462.                     e.vy = speed*j_scat;
  463.                     e.vz = speed*k_scat;
  464.  
  465.                     collision_counter++;  
  466.                     exc1_coll_counter++;
  467.  
  468.                     e.collided = true;
  469.                 }
  470.             }    
  471.  
  472.             else {
  473.                 collision_counter++;
  474.                 null_coll_counter++;
  475.                 e.collided = true;
  476.             }
  477.         }
  478.                 // calculating mean energy
  479.                 double total_energy = 0.0;
  480.                 for (const auto& e : electrons) total_energy += e.energy;
  481.                 double mean_energy = total_energy / n_e;
  482.                 file8 << t*dt << " " << mean_energy << "\n";                
  483.     }
  484.  
  485.     // ----- final electron energies distribution begins
  486.     for (int i = 0; i < n_e; i++){
  487.  
  488.         file2 << i << " " << electrons[i].energy << "\n";
  489.  
  490.         int bin = (int)( (electrons[i].energy - Emin)/bin_width_smooth );
  491.         if (bin >=0 && bin < histo_maxwell.size())
  492.             histo_maxwell[bin]++;
  493.     }
  494.  
  495.     int check = 0;
  496.     for (int i = 0; i < histo_maxwell.size(); i++){
  497.         check += histo_maxwell[i];
  498.         double bin_center = Emin + (i + 0.5) * bin_width;
  499.         file4 << bin_center << " " <<  static_cast<double>(histo_maxwell[i])/(electrons.size()*bin_width) << "\n"; // getting f(E)
  500.     }
  501.  
  502.     //recalculating average frequencies to adjust validation constants later
  503.     for (int i = 0; i < histo_maxwell.size(); i++){
  504.         double bin_center = Emin + (i + 0.5) * bin_width;
  505.         file11 << bin_center << " " << (N_He/Volume)*elastic_vec[i] * sqrt(2.0*bin_center*q/m_e)*static_cast<double>(histo_maxwell[i])/(electrons.size()*bin_width) << "\n";
  506.         file12 << bin_center << " " << (N_He/Volume)*inelastic1_vec[i] * sqrt(2.0*bin_center*q/m_e)*static_cast<double>(histo_maxwell[i])/(electrons.size()*bin_width) << "\n";
  507.     }
  508.     //----- calculating nu-average from our electron distribution ends ---------///        
  509.  
  510.  
  511.         std::cout << "Total # of electrons in a final histogram: " << check << "\n";
  512.  
  513.     // ----- final electron energies distribution begins    
  514.  
  515.  
  516.     file0.close();
  517.     file1.close();
  518.     file2.close();
  519.     file3.close();
  520.     file4.close();
  521.     file5.close();
  522.     file6.close();
  523.     file7.close();
  524.     file8.close();
  525.     file9.close();
  526.     file10.close();
  527.     file11.close();
  528.     file12.close();
  529.  
  530.     clock_t end = clock();
  531.  
  532.     double elapsed = (double)(end - start) / CLOCKS_PER_SEC;
  533.  
  534.     std::cout << "# of steps: " << steps << "\n";
  535.     std::cout << "# of electrons collided each timesteps:" << Ne_collided << "\n";
  536.  
  537.     std::cout << "Average triplet excitation collsisions per timestep: " << exc1_coll_counter/steps << "\n";
  538.     std::cout << "Average elastic collsisions per timestep: " << el_coll_counter/steps << "\n";
  539.     std::cout << "Average null collsisions per timestep: " << null_coll_counter/steps << "\n";
  540.  
  541.     std::cout << "Elapsed time: %f seconds " << elapsed << "\n";
  542.  
  543.  
  544.     return 0;
  545.  
  546. }
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!