phystota

1D PIC/MCC + DRR - v 1.2. -> 4levels added_reabsorption

Jun 9th, 2025
24
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 100.65 KB | None | 0 0
  1. #include <cstdio>
  2. #include <cstdlib>
  3. #include <cstring>
  4. #include <cstdbool>
  5. #include <cmath>
  6. #include <ctime>
  7. #include <random>
  8. #include <vector>
  9. #include <string>
  10. #include <fstream>
  11. #include <sstream>
  12. #include <algorithm>    
  13. #include <stdexcept>
  14. #include <iostream>
  15. #include <iomanip>
  16.  
  17. using namespace::std;
  18.  
  19. // constants
  20.  
  21. const double     PI             = 3.141592653589793;          // mathematical constant Pi
  22. const double     TWO_PI         = 2.0 * PI;                   // two times Pi
  23. const double     E_CHARGE       = 1.60217662e-19;             // electron charge [C]
  24. const double     EV_TO_J        = E_CHARGE;                   // eV <-> Joule conversion factor
  25. const double     E_MASS         = 9.109e-31;                  // mass of electron [kg]
  26. const double     HE_MASS        = 6.67e-27;                   // mass of He atom [kg]
  27. const double     MU_HEHE        = HE_MASS / 2.0;              // reduced mass of two He atoms [kg]
  28. const double     K_BOLTZMANN    = 1.38064852e-23;             // Boltzmann's constant [J/K]
  29. const double     EPSILON0       = 8.85418781e-12;             // permittivity of free space [F/m]
  30. const double     C_light        = 299792458.0;                 // speed of light [m/s]
  31.  
  32.  
  33. // simulation parameters
  34.  
  35. const int        N_G            = 501;                        // number of grid points
  36. const int        N_T            = 2000;                        // time steps within an RF period
  37.  
  38.  
  39. const double     FREQUENCY      = 13.56e6;                    // driving frequency [Hz]
  40. const double     VOLTAGE        = 350.0;                      // voltage amplitude [V]
  41. const double     L              = 0.04;                      // electrode gap [m]
  42. const double     PRESSURE       = 13.3322;                     // gas pressure [Pa] // n*k*T to match Turner's case
  43. const double     T_neutral      = 300.0;                      // background gas temperature [K] (also ion temperature)
  44. const double     T_electron     = 30000.0;                    // initial electron temperatutre [K]
  45. const double     WEIGHT         = 41875.0;                  // weight of superparticles
  46. const double     ELECTRODE_AREA = 1.6e-4;                     // (fictive) electrode area [m^2]
  47. const int        N_INIT         = 65536;                      // number of initial electrons and ions
  48.  
  49. // additional (derived) constants
  50.  
  51. const double     PERIOD         = 1.0 / FREQUENCY;                              // RF period length [s]
  52. const double     DT_E           = PERIOD / (double)(N_T);                       // electron time step [s]
  53. const int        N_SUB          = 10;                                            // ions move only in these cycles (subcycling)
  54. const int        N_avg          = 5;                                            // number of cycles to average rates --- Artem
  55.       int        counter        = 0;                                            // cycles since update of rates --- Artem
  56. const double     DT_I           = N_SUB * DT_E;                                 // ion time step [s]
  57. const double     DX             = L / (double)(N_G - 1);                        // spatial grid division [m]
  58. const double     INV_DX         = 1.0 / DX;                                     // inverse of spatial grid size [1/m]
  59. const double     GAS_DENSITY    = PRESSURE / (K_BOLTZMANN * T_neutral);         // background gas density [1/m^3]
  60. const double     OMEGA          = TWO_PI * FREQUENCY;                           // angular frequency [rad/s]
  61.  
  62. const double     Diff_coeff     = 8.992E-6*pow(T_neutral,1.5)/(PRESSURE/133.3); // diffusion coefficient for metastable states
  63.  
  64. // electron and ion cross sections
  65.  
  66. const int        N_CS           = 20;                          // total number of processes / cross sections
  67. const int        E_ELA          = 0;                          // process identifier: electron/elastic
  68. const int        E_EXC_1        = 1;                          // process identifier: electron/excitation1
  69. const int        E_EXC_2        = 2;                          // process identifier: electron/excitation1
  70. const int        E_ION          = 3;                          // process identifier: electron/ionization
  71. const int        I_ISO          = 4;                          // process identifier: ion/elastic/isotropic
  72. const int        I_BACK         = 5;                          // process identifier: ion/elastic/backscattering
  73.  
  74. const int        E_SUPER_1      = 6;                          // 1^1S -> 2^3S
  75. const int        E_SUPER_2      = 7;                          // 1^1S -> 2^1S
  76.  
  77. const int        E_EXC_3        = 8;                          // 1^1S -> 2^1P
  78. const int        E_EXC_4        = 9;                          // 1^1S -> 2^3P
  79. const int        E_EXC_12       = 10;                         // 2^1S -> 2^1P
  80. const int        E_EXC_13       = 11;                         // 2^3S -> 2^3P
  81. const int        E_EXC_17       = 12;                         // 2^3S -> 2^1S
  82. const int        E_EXC_19       = 13;                         // 2^1S -> 2^3P
  83.  
  84. const int        E_SUPER_3        = 14;                       // 1^1S -> 2^1P
  85. const int        E_SUPER_4        = 15;                       // 1^1S -> 2^3P
  86. const int        E_SUPER_12       = 16;                       // 2^1S -> 2^1P
  87. const int        E_SUPER_13       = 17;                       // 2^3S -> 2^3P
  88. const int        E_SUPER_17       = 18;                       // 2^3S -> 2^1S
  89. const int        E_SUPER_19       = 19;                       // 2^1S -> 2^3P
  90.  
  91.  
  92. const double     E_EXC_TH_1     = 19.82;                      // 1^1S -> 2^3S threshold [eV]
  93. const double     E_EXC_TH_2     = 20.61;                      // 1^1S -> 2^1S threshold [eV]
  94. const double     E_EXC_TH_3     = 21.218;                     // 1^1S -> 2^1P threshold [eV]
  95. const double     E_EXC_TH_4     = 20.96;                      // 1^1S -> 2^3P threshold [eV]
  96. const double     E_EXC_TH_12    = 0.602;                      // 2^1S -> 2^1P threshold [eV]
  97. const double     E_EXC_TH_13    = 1.1444;                     // 2^3S -> 2^3P threshold [eV]
  98. const double     E_EXC_TH_17    = 0.796;                      // 2^3S -> 2^1S threshold [eV]
  99. const double     E_EXC_TH_19    = 0.3483;                     // 2^1S -> 2^3P threshold [eV]
  100.  
  101. const double     E_ION_TH       = 24.587;                     // electron impact ionization threshold [eV]
  102.  
  103. const int        CS_RANGES      = 1000000;                    // number of entries in cross section arrays
  104. const double     DE_CS          = 0.001;                      // energy division in cross section arrays [eV]
  105. typedef float    cross_section[CS_RANGES];                    // cross section array
  106. cross_section    sigma[N_CS];                                 // set of cross section arrays
  107. cross_section    sigma_tot_e;                                 // total macroscopic cross section of electrons
  108. cross_section    sigma_tot_i;                                 // total macroscopic cross section of ions
  109.  
  110. double nu_avg;                                                // average nu for electrons thorugh 1 cycle
  111.  
  112. // particle coordinates
  113.  
  114. const int        MAX_N_P = 10000000;                           // maximum number of particles (electrons / ions)
  115. typedef double   particle_vector[MAX_N_P];                    // array for particle properties
  116. int              N_e = 0;                                     // number of electrons
  117. int              N_i = 0;                                     // number of ions
  118. int              N_e1 = 0;                                    // number of singlet excited states
  119. int              N_e2 = 0;                                    // number of triplet excited states
  120. particle_vector  x_e, vx_e, vy_e, vz_e;                       // coordinates of electrons (one spatial, three velocity components)
  121. particle_vector  x_i, vx_i, vy_i, vz_i;                       // coordinates of ions (one spatial, three velocity components)
  122.  
  123.  
  124.  
  125.  
  126. typedef double   xvector[N_G];                                // array for quantities defined at gird points
  127. xvector          efield,pot;                                  // electric field and potential
  128. xvector          e_density,i_density;                         // electron and ion densities
  129. xvector          cumul_e_density,cumul_i_density;             // cumulative densities
  130.  
  131. //excited states handling
  132. xvector e1_density = {0.0};                                   // 2^3S state                                              
  133. xvector e2_density = {0.0};                                   // 2^1S state
  134. xvector e3_density = {0.0};                                   // 2^1P state
  135. xvector e4_density = {0.0};                                   // 2^3P state
  136. xvector sum_electron_density = {0.0}; xvector avg_electron_density = {0.0};
  137.  
  138. xvector sum_rate1f  = {0.0}; xvector sum_rate1b  = {0.0}; xvector sum_rate2f  = {0.0}; xvector sum_rate2b  = {0.0};
  139. xvector sum_rate3f  = {0.0}; xvector sum_rate3b  = {0.0}; xvector sum_rate4f  = {0.0}; xvector sum_rate4b  = {0.0};
  140. xvector sum_rate12f = {0.0}; xvector sum_rate12b = {0.0}; xvector sum_rate13f = {0.0}; xvector sum_rate13b = {0.0};
  141. xvector sum_rate17f = {0.0}; xvector sum_rate17b = {0.0}; xvector sum_rate19f = {0.0}; xvector sum_rate19b = {0.0};
  142.  
  143. xvector avg_rate1f  = {0.0}; xvector avg_rate1b  = {0.0}; xvector avg_rate2f  = {0.0}; xvector avg_rate2b  = {0.0};
  144. xvector avg_rate3f  = {0.0}; xvector avg_rate3b  = {0.0}; xvector avg_rate4f  = {0.0}; xvector avg_rate4b  = {0.0};
  145. xvector avg_rate12f = {0.0}; xvector avg_rate12b = {0.0}; xvector avg_rate13f = {0.0}; xvector avg_rate13b = {0.0};
  146. xvector avg_rate17f = {0.0}; xvector avg_rate17b = {0.0}; xvector avg_rate19f = {0.0}; xvector avg_rate19b = {0.0};
  147.  
  148.  
  149. typedef unsigned long long int Ullong;                        // compact name for 64 bit unsigned integer
  150. Ullong       N_e_abs_pow  = 0;                                // counter for electrons absorbed at the powered electrode
  151. Ullong       N_e_abs_gnd  = 0;                                // counter for electrons absorbed at the grounded electrode
  152. Ullong       N_i_abs_pow  = 0;                                // counter for ions absorbed at the powered electrode
  153. Ullong       N_i_abs_gnd  = 0;                                // counter for ions absorbed at the grounded electrode
  154.  
  155. // electron energy probability function
  156.  
  157. const int    N_EEPF  = 2000;                                 // number of energy bins in Electron Energy Probability Function (EEPF)
  158. const double DE_EEPF = 0.05;                                 // resolution of EEPF [eV]
  159. typedef double eepf_vector[N_EEPF];                          // array for EEPF
  160. eepf_vector eepf     = {0.0};                                // time integrated EEPF in the center of the plasma
  161.  
  162. // ion flux-energy distributions
  163.  
  164. const int    N_IFED   = 200;                                 // number of energy bins in Ion Flux-Energy Distributions (IFEDs)
  165. const double DE_IFED  = 1.0;                                 // resolution of IFEDs [eV]
  166. typedef int  ifed_vector[N_IFED];                            // array for IFEDs
  167. ifed_vector  ifed_pow = {0};                                 // IFED at the powered electrode
  168. ifed_vector  ifed_gnd = {0};                                 // IFED at the grounded electrode
  169. double       mean_i_energy_pow;                              // mean ion energy at the powered electrode
  170. double       mean_i_energy_gnd;                              // mean ion energy at the grounded electrode
  171.  
  172. // spatio-temporal (XT) distributions
  173.  
  174. const int N_BIN                     = 20;                    // number of time steps binned for the XT distributions
  175. const int N_XT                      = N_T / N_BIN;           // number of spatial bins for the XT distributions
  176. typedef double xt_distr[N_G][N_XT];                          // array for XT distributions (decimal numbers)
  177. xt_distr pot_xt                     = {0.0};                 // XT distribution of the potential
  178. xt_distr efield_xt                  = {0.0};                 // XT distribution of the electric field
  179. xt_distr ne_xt                      = {0.0};                 // XT distribution of the electron density
  180. xt_distr ni_xt                      = {0.0};                 // XT distribution of the ion density
  181. xt_distr ue_xt                      = {0.0};                 // XT distribution of the mean electron velocity
  182. xt_distr ui_xt                      = {0.0};                 // XT distribution of the mean ion velocity
  183. xt_distr je_xt                      = {0.0};                 // XT distribution of the electron current density
  184. xt_distr ji_xt                      = {0.0};                 // XT distribution of the ion current density
  185. xt_distr powere_xt                  = {0.0};                 // XT distribution of the electron powering (power absorption) rate
  186. xt_distr poweri_xt                  = {0.0};                 // XT distribution of the ion powering (power absorption) rate
  187. xt_distr meanee_xt                  = {0.0};                 // XT distribution of the mean electron energy
  188. xt_distr meanei_xt                  = {0.0};                 // XT distribution of the mean ion energy
  189. xt_distr counter_e_xt               = {0.0};                 // XT counter for electron properties
  190. xt_distr counter_i_xt               = {0.0};                 // XT counter for ion properties
  191. xt_distr ioniz_rate_xt              = {0.0};                 // XT distribution of the ionisation rate
  192.  
  193. xt_distr e1_xt                      = {0.0};                 // XT distribution of 2^3S excited states density - Artem
  194. xt_distr e2_xt                      = {0.0};                 // XT distribution of 2^1S excited states density - Artem
  195. xt_distr e3_xt                      = {0.0};                 // XT distribution of 2^1P excited states density - Artem
  196. xt_distr e4_xt                      = {0.0};                 // XT distribution of 2^3P excited states density - Artem
  197.  
  198. double   mean_energy_accu_center    = 0;                     // mean electron energy accumulator in the center of the gap
  199. Ullong   mean_energy_counter_center = 0;                     // mean electron energy counter in the center of the gap
  200. Ullong   N_e_coll                   = 0;                     // counter for electron collisions
  201. Ullong   N_i_coll                   = 0;                     // counter for ion collisions
  202. double   Time;                                               // total simulated time (from the beginning of the simulation)
  203. int      cycle,no_of_cycles,cycles_done;                     // current cycle and total cycles in the run, cycles completed
  204. int      arg1;                                               // used for reading command line arguments
  205. char     st0[80];                                            // used for reading command line arguments
  206. FILE     *datafile;                                          // used for saving data
  207. bool     measurement_mode;                                   // flag that controls measurements and data saving
  208.  
  209. const double gamma_i                = 0.3;                   // yielding coefficient for ions bombarding the surface
  210. const double gamma_m                = 0.3;                   // yielding coefficient for metastables bombarding the surface - not using now
  211. const double gamma_p                = -10.0;                   // yielding coefficient for photons bombarding the surface - not using now
  212. const double gamma_k                = 0.5;                     // recombination coefficient  
  213. const double r                      = 0.7;                   // probabilty for the electron to reflect at the surface
  214. const double T_SEE                  = 300;                   // initial temperature of SEE electrons [K]
  215.  
  216. const double Ry                     = 13.6057;                // Rydberg constant, eV
  217. const double Bohr                   = 0.8797e-20;             // atomic consant related to Bohr Radius, m^2
  218.  
  219.  
  220. //---------------------------------------------------------------------------//
  221. // C++ Mersenne Twister 19937 generator                                      //
  222. // R01(MTgen) will genarate uniform distribution over [0,1) interval         //
  223. // RMB(MTgen) will generate Maxwell-Boltzmann distribution (of gas atoms)    //
  224. //---------------------------------------------------------------------------//
  225.  
  226. std::random_device rd{};
  227. std::mt19937 MTgen(rd());
  228. std::uniform_real_distribution<> R01(0.0, 1.0);
  229. std::normal_distribution<> RMB_n(0.0,sqrt(K_BOLTZMANN * T_neutral / HE_MASS));
  230. std::normal_distribution<> RMB_e(0.0,sqrt(K_BOLTZMANN * T_electron / E_MASS));
  231. std::normal_distribution<> RMB_SEE(0.0,sqrt(K_BOLTZMANN * T_SEE / E_MASS));
  232.  
  233. double series(double tau){
  234.     int k = 20;             // accuracy of sum calculation, number of series expansion parts
  235.     double sum = 0.0;
  236.     for (int i  = 1; i < k+1; i++){
  237.         double factor = std::tgamma((i+2) + 1);
  238.         // for (int j = 0; j < i+2; j++){
  239.         //     factor *= i + 2 - j;
  240.         // }
  241.        
  242.         sum += pow(-1.0, i) * pow(tau, i+1) / (i*factor*sqrt(i+2));
  243.     }
  244.     return sum;
  245. }
  246.  
  247. double compute_escape_factor(double tau){
  248.     if (tau <= 0) return 1;                 //optically thin plasma
  249.    
  250.     if (tau > 2.5)    return 1/(2.0*tau*sqrt(M_PI))*(sqrt(log(2.0*tau)) + 0.25/sqrt(log(2.0*tau)) + 0.14);
  251.     else              return 1.0 - 0.8293*tau + 0.7071*tau*log(2.0*tau) + series(tau);    
  252. }
  253.  
  254.  
  255. class CSInterpolator {
  256. public:
  257.   // load "filename" which must have two whitespace‐separated columns:
  258.   //    energy (eV)    cross_section (in m^2 or cm^2 as you prefer)
  259.   CSInterpolator(const std::string &filename) {
  260.     std::ifstream in(filename);
  261.     if (!in) throw std::runtime_error("CSInterpolator: cannot open " + filename);
  262.     double E, sigma;
  263.     std::string line;
  264.     while (std::getline(in, line)) {
  265.       std::istringstream iss(line);
  266.       if (iss >> E >> sigma) {
  267.         E_pts_.push_back(E);
  268.         sigma_pts_.push_back(sigma);
  269.       }
  270.     }
  271.     if (E_pts_.size()<2)
  272.       throw std::runtime_error("CSInterpolator: need at least two data points in " + filename);
  273.   }
  274.  
  275.   // return σ(E) by simple linear interpolation (clamped to end‐points)
  276.   double operator()(double E) const {
  277.     auto it = std::lower_bound(E_pts_.begin(), E_pts_.end(), E);
  278.     if (it == E_pts_.begin()) {
  279. //        std::cerr << "Warning: E="<<E<<" below data range, clamping to "<< 0.0 <<"\n";
  280.         return 0.0;
  281.     }  
  282.     if (it == E_pts_.end()){
  283. //        std::cerr << "Warning: E="<<E<<" above data range, clamping to "<< sigma_pts_.back() <<"\n";
  284.         return sigma_pts_.back();
  285.     }    
  286.     size_t idx = (it - E_pts_.begin());
  287.     double E1 = E_pts_[idx-1], E2 = E_pts_[idx];
  288.     double s1 = sigma_pts_[idx-1], s2 = sigma_pts_[idx];
  289.     // linear interp
  290.     return s1 + (s2 - s1) * (E - E1)/(E2 - E1);
  291.   }
  292.  
  293. private:
  294.   std::vector<double> E_pts_, sigma_pts_;
  295. };
  296.  
  297. // Cross-Section function for excitation cross-sections.
  298. // here for type 'd' - dipole-forbidden, 's' - spin-forbidden, 'a' - dipole-allowed
  299. // g - stat weight of initial level (g = 1 for ground state). See paper: RalchenkoJanev2008
  300. double forbidden_CS(double A1, double A2, double A3, double A4, double A5, double A6, double E, double E_th, int g, char type){
  301.     double x = E/E_th;
  302.     double strength, S1, S2, CS;
  303.  
  304.     if (type == 'd'){
  305.         S1 = (A1 + A2/x + A3/(x*x) + A4/(x*x*x));
  306.         S2 = (x*x/(x*x + A5));
  307.     }
  308.     if (type == 's'){
  309.         S1 = (A1 + A2/x + A3/(x*x) + A4/(x*x*x));
  310.         S2 = (1.0/(x*x + A5));
  311.     }
  312.     if (type == 'a'){
  313.         S1 = (A1*log(x) + A2 + A3/x + A4/(x*x) + A5/(x*x*x));
  314.         S2 = (x + 1.0)/(x + A6);
  315.     }
  316.    
  317.     strength = S1*S2;
  318.     CS = Bohr * Ry/(g*E) * strength;
  319.     return CS;
  320. }
  321.  
  322. void set_electron_cross_sections_ar(void){
  323.     int    i;
  324.     double A1_1, A1_2, A1_3, A1_4, A1_5, A1_6, A2_1, A2_2, A2_3, A2_4, A2_5, A2_6;
  325.     double A3_1, A3_2, A3_3, A3_4, A3_5, A3_6, A4_1, A4_2, A4_3, A4_4, A4_5, A4_6;
  326.     double A12_1, A12_2, A12_3, A12_4, A12_5, A12_6, A13_1, A13_2, A13_3, A13_4, A13_5, A13_6;
  327.     double A17_1, A17_2, A17_3, A17_4, A17_5, A17_6, A19_1, A19_2, A19_3, A19_4, A19_5, A19_6;
  328.     int g1, g2, g3, g4, g12, g13, g17, g19;                                                         // stat weights of initial level      
  329.  
  330.     A1_1  = +6.888e-01; A1_2  = +1.975e-01; A1_3  = +7.232e+00; A1_4  = -4.839e+00; A1_5  = +5.003e+01; A1_6  = +0.000e+00;     // 1¹S -> 2^3S  -- s-f
  331.     A2_1  = +1.888e-01; A2_2  = -5.754e-01; A2_3  = +3.439e+00; A2_4  = -2.088e+00; A2_5  = +2.544e+01; A2_6  = +0.000e+00;     // 1¹S -> 2^1S  -- d-f
  332.     A3_1  = +7.087e-01; A3_2  = -9.347e-02; A3_3  = -1.598e+00; A3_4  = +2.986e+00; A3_5  = -1.293e+00; A3_6  = +3.086e-01;     // 1¹S -> 2¹P   -- d-a
  333.     A4_1  = +2.823e-01; A4_2  = +2.048e+00; A4_3  = +5.287e+00; A4_4  = -7.363e+00; A4_5  = +2.728e+01; A4_6 =  +0.000e+00;     // 1¹S -> 2³P   -- s-f
  334.     A12_1 = +3.404e+01; A12_2 = +7.267e+01; A12_3 = +1.710e+02; A12_4 = -7.033e+02; A12_5 = +4.704e+02; A12_6 = +1.194e+01;     // 2^1S -> 2^1P -- d-a
  335.     A13_1 = +7.696e+01; A13_2 = +1.250e+02; A13_3 = +4.938e+01; A13_4 = -4.778e+02; A13_5 = +3.189e+02; A13_6 = +8.157e+00;     // 2^3S -> 2^3P -- d-a
  336.     A17_1 = +5.475e+01; A17_2 = +3.483e+05; A17_3 = +2.345e+05; A17_4 = -5.090e+05; A17_5 = +5.424e+04; A17_6 = +0.000e+00;     // 2^3S -> 2^1S -- s-f
  337.     A19_1 = +5.983e+02; A19_2 = -5.310e+02; A19_3 = +3.348e+02; A19_4 = -2.412e+02; A19_5 = +2.239e+02; A19_6 = +0.000e+00;     // 2^1S -> 2^3P -- s-f
  338.  
  339.     g1 = 1; g2 = 1; g3 = 1; g4 = 1; g12 = 1; g13 = 3; g17 = 3; g19 = 1;
  340.  
  341.     printf(">> eduPIC: Setting e- / He cross sections\n");
  342.  
  343.     // load your four datafiles (make sure these names match your files!)
  344.     CSInterpolator cs_ela  ("CS/He_electron_elastic.dat");   // two‐col: E σ_ela
  345.     CSInterpolator cs_exc1 ("CS/He_electron_exc1.dat");      // two‐col: E σ_exc1
  346.     CSInterpolator cs_exc2 ("CS/He_electron_exc2.dat");      // two‐col: E σ_exc2
  347.     CSInterpolator cs_ion  ("CS/He_electron_ionization.dat");// two‐col: E σ_ion
  348.  
  349.     for(int i=0; i<CS_RANGES; i++){
  350.         // your energy grid
  351.         double en = (i==0 ? DE_CS : DE_CS * i);
  352.  
  353.         // interpolate
  354.         sigma[E_ELA][i]     = cs_ela(en);
  355.         sigma[E_ION][i]     = cs_ion(en);
  356.         // sigma[E_EXC_1][i]   = cs_exc1(en);
  357.         // sigma[E_EXC_2][i]   = cs_exc2(en);
  358.         if (en < E_EXC_TH_1){
  359.             sigma[E_EXC_1][i] = 0.0;
  360.         }
  361.         else {
  362.             sigma[E_EXC_1][i] = forbidden_CS(A1_1, A1_2, A1_3, A1_4, A1_5, A1_6, en, E_EXC_TH_1, g1, 's');
  363.         }
  364.  
  365.         if (en < E_EXC_TH_2){
  366.             sigma[E_EXC_2][i] = 0.0;
  367.         }
  368.         else {
  369.             sigma[E_EXC_2][i] = forbidden_CS(A2_1, A2_2, A2_3, A2_4, A2_5, A2_6, en, E_EXC_TH_2, g2, 'd');
  370.         }
  371.  
  372.         if (en < E_EXC_TH_3){
  373.             sigma[E_EXC_3][i] = 0.0;
  374.         }
  375.         else {
  376.             sigma[E_EXC_3][i] = forbidden_CS(A3_1, A3_2, A3_3, A3_4, A3_5, A3_6, en, E_EXC_TH_3, g3, 'a');
  377.         }  
  378.  
  379.         if (en < E_EXC_TH_4){
  380.             sigma[E_EXC_4][i] = 0.0;
  381.         }
  382.         else {
  383.             sigma[E_EXC_4][i] = forbidden_CS(A4_1, A4_2, A4_3, A4_4, A4_5, A4_6, en, E_EXC_TH_4, g4, 's');
  384.         }        
  385.  
  386.         if (en < E_EXC_TH_12){
  387.             sigma[E_EXC_12][i] = 0.0;
  388.         }
  389.         else {
  390.             sigma[E_EXC_12][i] = forbidden_CS(A12_1, A12_2, A12_3, A12_4, A12_5, A12_6, en, E_EXC_TH_12, g12, 'a');
  391.         }    
  392.  
  393.         if (en < E_EXC_TH_13){
  394.             sigma[E_EXC_13][i] = 0.0;
  395.         }
  396.         else {
  397.             sigma[E_EXC_13][i] = forbidden_CS(A13_1, A13_2, A13_3, A13_4, A13_5, A13_6, en, E_EXC_TH_13, g13, 'a');
  398.         }      
  399.  
  400.         if (en < E_EXC_TH_17){
  401.             sigma[E_EXC_17][i] = 0.0;
  402.         }
  403.         else {
  404.             sigma[E_EXC_17][i] = forbidden_CS(A17_1, A17_2, A17_3, A17_4, A17_5, A17_6, en, E_EXC_TH_17, g17, 's');
  405.         }                    
  406.  
  407.         if (en < E_EXC_TH_19){
  408.             sigma[E_EXC_19][i] = 0.0;
  409.         }
  410.         else {
  411.             sigma[E_EXC_19][i] = forbidden_CS(A19_1, A19_2, A19_3, A19_4, A19_5, A19_6, en, E_EXC_TH_19, g19, 's');
  412.         }                        
  413.     }
  414.  
  415.     for(int i=0; i<CS_RANGES; i++){
  416.         // your energy grid
  417.         double en = (i==0 ? DE_CS : DE_CS * i);
  418.  
  419.         // Superelastic for 1^1S -> 2^3S (E_SUPER_1)
  420.         double en_plus_1 = en + E_EXC_TH_1;
  421.         int idx1 = en_plus_1 / DE_CS;
  422.         if (idx1 < CS_RANGES && en > 0) {
  423.             sigma[E_SUPER_1][i] = (1.0 / 3.0) * (en_plus_1 / en) * sigma[E_EXC_1][idx1];
  424.         } else {
  425.             sigma[E_SUPER_1][i] = 0.0;
  426.         }        
  427.  
  428.         // Superelastic for 1^1S -> 2^1S (E_SUPER_2)
  429.         double en_plus_2 = en + E_EXC_TH_2;
  430.         int idx2 = en_plus_2 / DE_CS;
  431.         if (idx2 < CS_RANGES && en > 0) {
  432.             sigma[E_SUPER_2][i] = (1.0 / 1.0) * (en_plus_2 / en) * sigma[E_EXC_2][idx2];
  433.         } else {
  434.             sigma[E_SUPER_2][i] = 0.0;
  435.         }        
  436.  
  437.         // Superelastic for 1^1S -> 2^1P (E_SUPER_3)
  438.         double en_plus_3 = en + E_EXC_TH_3;
  439.         int idx3 = en_plus_3 / DE_CS;
  440.         if (idx3 < CS_RANGES && en > 0) {
  441.             sigma[E_SUPER_3][i] = (1.0 / 3.0) * (en_plus_3 / en) * sigma[E_EXC_3][idx3];
  442.         } else {
  443.             sigma[E_SUPER_3][i] = 0.0;
  444.         }
  445.  
  446.         // Superelastic for 1^1S -> 2^3P (E_SUPER_4)
  447.         double en_plus_4 = en + E_EXC_TH_4;
  448.         int idx4 = en_plus_4 / DE_CS;
  449.         if (idx4 < CS_RANGES && en > 0) {
  450.             sigma[E_SUPER_4][i] = (1.0 / 9.0) * (en_plus_4 / en) * sigma[E_EXC_4][idx4];
  451.         } else {
  452.             sigma[E_SUPER_4][i] = 0.0;
  453.         }      
  454.  
  455.         // Superelastic for 2^1S -> 2^1P (E_SUPER_12)
  456.         double en_plus_12 = en + E_EXC_TH_12;
  457.         int idx12 = en_plus_12 / DE_CS;
  458.         if (idx12 < CS_RANGES && en > 0) {
  459.             sigma[E_SUPER_12][i] = (1.0 / 3.0) * (en_plus_12 / en) * sigma[E_EXC_12][idx12];
  460.         } else {
  461.             sigma[E_SUPER_12][i] = 0.0;
  462.         }        
  463.  
  464.         // Superelastic for 2^3S -> 2^3P (E_SUPER_13)
  465.         double en_plus_13 = en + E_EXC_TH_13;
  466.         int idx13 = en_plus_13 / DE_CS;
  467.         if (idx13 < CS_RANGES && en > 0) {
  468.             sigma[E_SUPER_13][i] = (3.0 / 9.0) * (en_plus_13 / en) * sigma[E_EXC_13][idx13];
  469.         } else {
  470.             sigma[E_SUPER_13][i] = 0.0;
  471.         }    
  472.  
  473.         // Superelastic for 2^3S -> 2^1S (E_SUPER_17)
  474.         double en_plus_17 = en + E_EXC_TH_17;
  475.         int idx17 = en_plus_17 / DE_CS;
  476.         if (idx17 < CS_RANGES && en > 0) {
  477.             sigma[E_SUPER_17][i] = (3.0 / 1.0) * (en_plus_17 / en) * sigma[E_EXC_17][idx17];
  478.         } else {
  479.             sigma[E_SUPER_17][i] = 0.0;
  480.         }                          
  481.  
  482.         // Superelastic for 2^1S -> 2^3P (E_SUPER_19)
  483.         double en_plus_19 = en + E_EXC_TH_19;
  484.         int idx19 = en_plus_19 / DE_CS;
  485.         if (idx19 < CS_RANGES && en > 0) {
  486.             sigma[E_SUPER_19][i] = (1.0 / 9.0) * (en_plus_19 / en) * sigma[E_EXC_19][idx19];
  487.         } else {
  488.             sigma[E_SUPER_19][i] = 0.0;
  489.         }            
  490.     }    
  491. }
  492.  
  493. //------------------------------------------------------------------------------//
  494. //  ion cross sections: A. V. Phelps, J. Appl. Phys. 76, 747 (1994)             //
  495. //------------------------------------------------------------------------------//
  496.  
  497. void set_ion_cross_sections_ar(void){
  498.     int    i;
  499.     double e_com,e_lab,qmom,qback,qiso;
  500.    
  501.     printf(">> eduPIC: Setting He+ / He cross sections\n");
  502.     for(i=0; i<CS_RANGES; i++){
  503.         if (i == 0) {e_com = DE_CS;} else {e_com = DE_CS * i;}             // ion energy in the center of mass frame of reference
  504.         e_lab = 2.0 * e_com;                                               // ion energy in the laboratory frame of reference
  505.         qiso  = 7.63 *pow(10,-20) * pow(e_com, -0.5);
  506.         qback = 1.0 * pow(10,-19) * pow( (e_com/1000.0), -0.15 ) * pow( (1.0 + e_com/1000.0), -0.25 ) * pow( (1.0 + 5.0/e_com), -0.15 );
  507.         sigma[I_ISO][i]  = qiso;             // cross section for He+ / He isotropic part of elastic scattering
  508.         sigma[I_BACK][i] = qback;            // cross section for He+ / He backward elastic scattering
  509.     }
  510. }
  511.  
  512. //----------------------------------------------------------------------//
  513. //  calculation of total cross sections for electrons and ions          //
  514. //----------------------------------------------------------------------//
  515.  
  516. void calc_total_cross_sections(void){
  517.     int i;
  518.    
  519.     for(i=0; i<CS_RANGES; i++){
  520.         sigma_tot_e[i] = (sigma[E_ELA][i] + sigma[E_EXC_1][i] + sigma[E_EXC_2][i] + sigma[E_ION][i]) * GAS_DENSITY;   // total macroscopic cross section of electrons
  521.         sigma_tot_i[i] = (sigma[I_ISO][i] + sigma[I_BACK][i]) * GAS_DENSITY;                    // total macroscopic cross section of ions
  522.     }
  523. }
  524.  
  525. //----------------------------------------------------------------------//
  526. //  test of cross sections for electrons and ions                       //
  527. //----------------------------------------------------------------------//
  528.  
  529. void test_cross_sections(void){
  530.     FILE  * f;
  531.     int   i,j;
  532.    
  533.     f = fopen("cross_sections.dat","w");        // cross sections saved in data file: cross_sections.dat
  534.     for(i=0; i<CS_RANGES; i++){
  535.         fprintf(f,"%12.4f ",i*DE_CS);
  536.         for(j=0; j<N_CS; j++) fprintf(f,"%14e ",sigma[j][i]);
  537.         fprintf(f,"\n");
  538.     }
  539.     fclose(f);
  540. }
  541.  
  542. //---------------------------------------------------------------------//
  543. // find upper limit of collision frequencies                           //
  544. //---------------------------------------------------------------------//
  545.  
  546. double max_electron_coll_freq (void){
  547.     int i;
  548.     double e,v,nu,nu_max;
  549.     nu_max = 0;
  550.     for(i=0; i<CS_RANGES; i++){
  551.         e  = i * DE_CS;
  552.         v  = sqrt(2.0 * e * EV_TO_J / E_MASS);
  553.         nu = v * sigma_tot_e[i];
  554.         if (nu > nu_max) {nu_max = nu;}
  555.     }
  556.     return nu_max;
  557. }
  558.  
  559. double max_ion_coll_freq (void){
  560.     int i;
  561.     double e,g,nu,nu_max;
  562.     nu_max = 0;
  563.     for(i=0; i<CS_RANGES; i++){
  564.         e  = i * DE_CS;
  565.         g  = sqrt(2.0 * e * EV_TO_J / MU_HEHE);
  566.         nu = g * sigma_tot_i[i];
  567.         if (nu > nu_max) nu_max = nu;
  568.     }
  569.     return nu_max;
  570. }
  571.  
  572. //----------------------------------------------------------------------//
  573. // initialization of the simulation by placing a given number of        //
  574. // electrons and ions at random positions between the electrodes        //
  575. //----------------------------------------------------------------------//
  576.  
  577. // initialization of excited states distribtuion, assuming Maxwellian-Boltzmann balance -- Artem
  578. std::pair<double, double> init_excited_distr() {
  579.     double part_ground = 1.0*exp(-0.0/T_neutral); // partition function for ground state
  580.     double part_triplet = 3.0*exp(-E_EXC_TH_1*EV_TO_J/(K_BOLTZMANN*T_neutral)); // partition function for triplet excited state
  581.     double part_singlet = 1.0*exp(-E_EXC_TH_2*EV_TO_J/(K_BOLTZMANN*T_neutral)); // partition function for singlet excited state
  582.     double part_func_total = part_ground + part_triplet + part_singlet; // total partition function
  583.     double n_triplet = ((part_triplet/part_func_total)*GAS_DENSITY); // denisty population of tripet state
  584.     double n_singlet = ((part_singlet/part_func_total)*GAS_DENSITY); // density population of singlet state
  585.  
  586.     return {n_triplet, n_singlet};
  587. }
  588.  
  589. void print_excitation_densities(void) {
  590.     double total_e1 = 0.0, total_e2 = 0.0, total_e3 = 0.0, total_e4 = 0.0;
  591.     const double cell_volume = ELECTRODE_AREA * DX;
  592.    
  593.     // Sum densities across all grid cells
  594.     for (int p = 0; p < N_G; p++) {
  595.         total_e1 += e1_density[p];  // 2^3S state density
  596.         total_e2 += e2_density[p];  // 2^1S state density
  597.         total_e3 += e3_density[p];  // 2^1P state density
  598.         total_e4 += e4_density[p];  // 2^3P state density        
  599.     }
  600.     std::cout << "average densities of excited states: " <<"\n";
  601.     printf("2^3S = %12.6e | 2^1S = %12.6e\n", total_e1/N_G, total_e2/N_G);
  602.     printf("2^1P = %12.6e | 2^3P = %12.6e\n", total_e3/N_G, total_e4/N_G);    
  603. }
  604.  
  605. void init(int nseed){
  606.     int i;
  607.    
  608.     for (i=0; i<nseed; i++){
  609.         x_e[i]  = L * R01(MTgen);                                                   // initial random position of the electron
  610.         vx_e[i] = RMB_e(MTgen); vy_e[i] = RMB_e(MTgen); vz_e[i] = RMB_e(MTgen);     // initial velocity components of the electron
  611.         x_i[i]  = L * R01(MTgen);                                                   // initial random position of the ion
  612.         vx_i[i] = RMB_n(MTgen); vy_i[i] = RMB_n(MTgen); vz_i[i] = RMB_n(MTgen);     // initial velocity components of the ion
  613.     }
  614.     N_e = nseed;    // initial number of electrons
  615.     N_i = nseed;    // initial number of ions
  616.  
  617.     // auto exc = init_excited_distr();
  618.     // for (int p = 0; p < N_G; p++) {
  619.     //     e1_density[p] = exc.first;
  620.     //     e2_density[p] = exc.second;
  621.     // }
  622.  
  623.     double n0 = 2.7e+9; double a0 = -1.3e+8; double b0 = 5.0e+15;
  624.     for (int p = 0; p < N_G; p++){
  625.         // e1_density[p] = abs(n0 * (a0 *(pow(L,4) - pow((p*DX - 0.5*L),2)) +
  626.         //                     b0 * (pow(L,8) -  pow((p*DX - 0.5*L),8)) ));
  627.         // e2_density[p] = 0.5 * abs(  n0 * (a0 *(pow(L,4) - pow((p*DX - 0.5*L),2)) +
  628.         //                     b0 * (pow(L,8) -  pow((p*DX - 0.5*L),8)) ));    
  629.         e1_density[p] = 1.0e+15;
  630.         e2_density[p] = 1.0e+15;
  631.         e3_density[p] = 1.0e+12;         // initial non-zero safe density for 2^1P
  632.         e4_density[p] = 1.0e+12;         // initial non-zero safe density for 2^3P                                            
  633.     }
  634. }
  635.  
  636. //----------------------------------------------------------------------//
  637. // e / He collision  (cold gas approximation)                           //
  638. //----------------------------------------------------------------------//
  639.  
  640. void collision_electron (double xe, double *vxe, double *vye, double *vze, int eindex){
  641.     const double F1 = E_MASS  / (E_MASS + HE_MASS);
  642.     const double F2 = HE_MASS / (E_MASS + HE_MASS);
  643.     double t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,t11,t12,t13,t14,t15,t16,t17,rnd;
  644.     double g,g2,gx,gy,gz,wx,wy,wz,theta,phi;
  645.     double chi,eta,chi2,eta2,sc,cc,se,ce,st,ct,sp,cp,energy,e_sc,e_ej;
  646.  
  647.     // - Artem
  648.     // Determine cell p where the electron is
  649.     double c0 = xe * INV_DX;
  650.     int p = std::max(0, std::min(N_G - 1, static_cast<int>(c0)));
  651.  
  652.     // Local densities -- Artem
  653.     double e1_dens = e1_density[p]; // 2^3S density
  654.     double e2_dens = e2_density[p];                                                 // 2^1S density
  655.     double e3_dens = e3_density[p];                                                 // 2^1P density
  656.     double e4_dens = e4_density[p];                                                 // 2^3P density
  657.     double ground_dens = GAS_DENSITY - e1_dens - e2_dens - e3_dens - e4_dens;       // should I extract it? let's see     Artem
  658.  
  659.    
  660.     // calculate relative velocity before collision & velocity of the centre of mass
  661.    
  662.     gx = (*vxe);
  663.     gy = (*vye);
  664.     gz = (*vze);
  665.     g  = sqrt(gx * gx + gy * gy + gz * gz);
  666.     wx = F1 * (*vxe);
  667.     wy = F1 * (*vye);
  668.     wz = F1 * (*vze);
  669.    
  670.     // find Euler angles
  671.    
  672.     if (gx == 0) {theta = 0.5 * PI;}
  673.     else {theta = atan2(sqrt(gy * gy + gz * gz),gx);}
  674.     if (gy == 0) {
  675.         if (gz > 0){phi = 0.5 * PI;} else {phi = - 0.5 * PI;}
  676.     } else {phi = atan2(gz, gy);}
  677.     st  = sin(theta);
  678.     ct  = cos(theta);
  679.     sp  = sin(phi);
  680.     cp  = cos(phi);
  681.    
  682.     // choose the type of collision based on the cross sections
  683.     // take into account energy loss in inelastic collisions
  684.     // generate scattering and azimuth angles
  685.     // in case of ionization handle the 'new' electron
  686.    
  687.     t0   =       sigma[E_ELA     ][eindex] * ground_dens;
  688.     t1   = t0  + sigma[E_EXC_1   ][eindex] * ground_dens;     // 1^1S -> 2^3S
  689.     t2   = t1  + sigma[E_EXC_2   ][eindex] * ground_dens;     // 1^1S -> 2^1S
  690.     t3   = t2  + sigma[E_ION     ][eindex] * ground_dens;
  691.     t4   = t3  + sigma[E_SUPER_1 ][eindex] * e1_dens;         // 1^1S -> 2^3S superelastic
  692.     t5   = t4  + sigma[E_SUPER_2 ][eindex] * e2_dens;         // 1^1S -> 2^1S superelastic
  693.     t6   = t5  + sigma[E_EXC_3   ][eindex] * ground_dens;     // 1^1S -> 2^1P excitation
  694.     t7   = t6  + sigma[E_SUPER_3 ][eindex] * e3_dens;         // 1^1S -> 2^1P superelastic
  695.     t8   = t7  + sigma[E_EXC_4   ][eindex] * ground_dens;     // 1^1S -> 2^3P excitation
  696.     t9   = t8  + sigma[E_SUPER_4 ][eindex] * e4_dens;         // 1^1S -> 2^3P superelastic
  697.     t10  = t9  + sigma[E_EXC_12  ][eindex] * e2_dens;         // 2^1S -> 2^1P excitation
  698.     t11  = t10 + sigma[E_SUPER_12][eindex] * e3_dens;         // 2^1S -> 2^1P superelastic        
  699.     t12  = t11 + sigma[E_EXC_13  ][eindex] * e1_dens;         // 2^3S -> 2^3P excitation
  700.     t13  = t12 + sigma[E_SUPER_13][eindex] * e4_dens;         // 2^3S -> 2^3P superelastic  
  701.     t14  = t13 + sigma[E_EXC_17  ][eindex] * e1_dens;         // 2^3S -> 2^1S excitation
  702.     t15  = t14 + sigma[E_SUPER_17][eindex] * e2_dens;         // 2^3S -> 2^1S superelastic  
  703.     t16  = t15 + sigma[E_EXC_19  ][eindex] * e2_dens;         // 2^1S -> 2^3P excitation
  704.     t17  = t16 + sigma[E_SUPER_19][eindex] * e4_dens;         // 2^1S -> 2^3P superelastic      
  705.  
  706.     rnd  = R01(MTgen);
  707.     if (rnd < (t0/t17)){                              // elastic scattering
  708.         chi = acos(1.0 - 2.0 * R01(MTgen));          // isotropic scattering
  709.         eta = TWO_PI * R01(MTgen);                   // azimuthal angle
  710.     } else if (rnd < (t1/t17)){                       // excitation 1 (triplet)
  711.         energy = 0.5 * E_MASS * g * g;               // electron energy
  712.         energy = fabs(energy - E_EXC_TH_1 * EV_TO_J);  // subtract energy loss for excitation
  713.         g   = sqrt(2.0 * energy / E_MASS);           // relative velocity after energy loss
  714.         chi = acos(1.0 - 2.0 * R01(MTgen));          // isotropic scattering
  715.         eta = TWO_PI * R01(MTgen);                   // azimuthal angle
  716.     } else if (rnd < (t2/t17)){                       // excitation 2 (singlet)
  717.         energy = 0.5 * E_MASS * g * g;               // electron energy
  718.         energy = fabs(energy - E_EXC_TH_2 * EV_TO_J);  // subtract energy loss for excitation
  719.         g   = sqrt(2.0 * energy / E_MASS);           // relative velocity after energy loss
  720.         chi = acos(1.0 - 2.0 * R01(MTgen));          // isotropic scattering
  721.         eta = TWO_PI * R01(MTgen);                   // azimuthal angle    
  722.     } else if (rnd < (t3/t17)) {                     // ionization
  723.         energy = 0.5 * E_MASS * g * g;               // electron energy
  724.         energy = fabs(energy - E_ION_TH * EV_TO_J);  // subtract energy loss of ionization
  725.         // e_ej  = 10.0 * tan(R01(MTgen) * atan(energy/EV_TO_J / 20.0)) * EV_TO_J; // energy of the ejected electron
  726.         // e_sc = fabs(energy - e_ej);                  // energy of scattered electron after the collision
  727.         e_ej = 0.5*energy;                          // energy of the ejected electron
  728.         e_sc = fabs(energy - e_ej);                  // energy of scattered electron after the collision        
  729.         g    = sqrt(2.0 * e_sc / E_MASS);            // relative velocity of scattered electron
  730.         g2   = sqrt(2.0 * e_ej / E_MASS);            // relative velocity of ejected electron
  731.         chi = acos(1.0 - 2.0 * R01(MTgen));          // isotropic scattering for scattered electron (as in Turner's case)
  732.         chi2 = acos(1.0 - 2.0 * R01(MTgen));          // isotropic scattering for ejected electrons (as in Turner's case)
  733.         eta  = TWO_PI * R01(MTgen);                  // azimuthal angle for scattered electron
  734.         eta2 = eta + PI;                             // azimuthal angle for ejected electron
  735.         sc  = sin(chi2);
  736.         cc  = cos(chi2);
  737.         se  = sin(eta2);
  738.         ce  = cos(eta2);
  739.         gx  = g2 * (ct * cc - st * sc * ce);
  740.         gy  = g2 * (st * cp * cc + ct * cp * sc * ce - sp * sc * se);
  741.         gz  = g2 * (st * sp * cc + ct * sp * sc * ce + cp * sc * se);
  742.         x_e[N_e]  = xe;                              // add new electron
  743.         vx_e[N_e] = wx + F2 * gx;
  744.         vy_e[N_e] = wy + F2 * gy;
  745.         vz_e[N_e] = wz + F2 * gz;
  746.         N_e++;
  747.         x_i[N_i]  = xe;                                // add new ion
  748.         vx_i[N_i] = RMB_n(MTgen);                      // velocity is sampled from background thermal distribution
  749.         vy_i[N_i] = RMB_n(MTgen);
  750.         vz_i[N_i] = RMB_n(MTgen);
  751.         N_i++;
  752.     }
  753.      else if (rnd < (t4/t17)) {                         // accounting for superelastic collisions - 1^1S -> 2^3S - Artem
  754.         energy = 0.5 * E_MASS * g * g;                  
  755.         energy = fabs(energy + E_EXC_TH_1 * EV_TO_J);  
  756.         g   = sqrt(2.0 * energy / E_MASS);              
  757.         chi = acos(1.0 - 2.0 * R01(MTgen));            
  758.         eta = TWO_PI * R01(MTgen);                          
  759.     } else if (rnd < (t5/t17)){                         // account for superelastic collision - 1^1S -> 2^1S - Artem
  760.         energy = 0.5 * E_MASS * g * g;                  
  761.         energy = fabs(energy + E_EXC_TH_2 * EV_TO_J);  
  762.         g   = sqrt(2.0 * energy / E_MASS);              
  763.         chi = acos(1.0 - 2.0 * R01(MTgen));            
  764.         eta = TWO_PI * R01(MTgen);                          
  765.     } else if (rnd < (t6/t17)) {                    // excitation 1^1S -> 2^1P
  766.         energy = 0.5 * E_MASS * g * g;              
  767.         energy = fabs(energy - E_EXC_TH_3 * EV_TO_J);  
  768.         g   = sqrt(2.0 * energy / E_MASS);          
  769.         chi = acos(1.0 - 2.0 * R01(MTgen));          
  770.         eta = TWO_PI * R01(MTgen);                  
  771.     } else if (rnd < (t7/t17)) {                    // superelastic 1^1S -> 2^1P
  772.         energy = 0.5 * E_MASS * g * g;                  
  773.         energy = fabs(energy + E_EXC_TH_3 * EV_TO_J);  
  774.         g   = sqrt(2.0 * energy / E_MASS);              
  775.         chi = acos(1.0 - 2.0 * R01(MTgen));            
  776.         eta = TWO_PI * R01(MTgen);        
  777.     } else if (rnd < (t8/t17)) {                    // excitation 1^1S -> 2^3P
  778.         energy = 0.5 * E_MASS * g * g;              
  779.         energy = fabs(energy - E_EXC_TH_4 * EV_TO_J);  
  780.         g   = sqrt(2.0 * energy / E_MASS);          
  781.         chi = acos(1.0 - 2.0 * R01(MTgen));          
  782.         eta = TWO_PI * R01(MTgen);                  
  783.     } else if (rnd < (t9/t17)) {                    // superelastic 1^1S -> 2^3P
  784.         energy = 0.5 * E_MASS * g * g;                  
  785.         energy = fabs(energy + E_EXC_TH_4 * EV_TO_J);  
  786.         g   = sqrt(2.0 * energy / E_MASS);              
  787.         chi = acos(1.0 - 2.0 * R01(MTgen));            
  788.         eta = TWO_PI * R01(MTgen);        
  789.     } else if (rnd < (t10/t17)) {                    // excitation 2^1S -> 2^1P
  790.         energy = 0.5 * E_MASS * g * g;              
  791.         energy = fabs(energy - E_EXC_TH_12 * EV_TO_J);  
  792.         g   = sqrt(2.0 * energy / E_MASS);          
  793.         chi = acos(1.0 - 2.0 * R01(MTgen));          
  794.         eta = TWO_PI * R01(MTgen);                  
  795.     } else if (rnd < (t11/t17)) {                    // superelastic 2^1S -> 2^1P
  796.         energy = 0.5 * E_MASS * g * g;                  
  797.         energy = fabs(energy + E_EXC_TH_12 * EV_TO_J);  
  798.         g   = sqrt(2.0 * energy / E_MASS);              
  799.         chi = acos(1.0 - 2.0 * R01(MTgen));            
  800.         eta = TWO_PI * R01(MTgen);        
  801.     } else if (rnd < (t12/t17)) {                    // excitation 2^3S -> 2^3P
  802.         energy = 0.5 * E_MASS * g * g;              
  803.         energy = fabs(energy - E_EXC_TH_13 * EV_TO_J);  
  804.         g   = sqrt(2.0 * energy / E_MASS);          
  805.         chi = acos(1.0 - 2.0 * R01(MTgen));          
  806.         eta = TWO_PI * R01(MTgen);                  
  807.     } else if (rnd < (t13/t17)) {                    // superelastic 2^3S -> 2^3P
  808.         energy = 0.5 * E_MASS * g * g;                  
  809.         energy = fabs(energy + E_EXC_TH_13 * EV_TO_J);  
  810.         g   = sqrt(2.0 * energy / E_MASS);              
  811.         chi = acos(1.0 - 2.0 * R01(MTgen));            
  812.         eta = TWO_PI * R01(MTgen);        
  813.     } else if (rnd < (t14/t17)) {                    // excitation 2^3S -> 2^1S
  814.         energy = 0.5 * E_MASS * g * g;              
  815.         energy = fabs(energy - E_EXC_TH_17 * EV_TO_J);  
  816.         g   = sqrt(2.0 * energy / E_MASS);          
  817.         chi = acos(1.0 - 2.0 * R01(MTgen));          
  818.         eta = TWO_PI * R01(MTgen);                  
  819.     } else if (rnd < (t15/t17)) {                    // superelastic 2^3S -> 2^1S
  820.         energy = 0.5 * E_MASS * g * g;                  
  821.         energy = fabs(energy + E_EXC_TH_17 * EV_TO_J);  
  822.         g   = sqrt(2.0 * energy / E_MASS);              
  823.         chi = acos(1.0 - 2.0 * R01(MTgen));            
  824.         eta = TWO_PI * R01(MTgen);        
  825.     } else if (rnd < (t16/t17)) {                    // excitation 2^1S -> 2^3P
  826.         energy = 0.5 * E_MASS * g * g;              
  827.         energy = fabs(energy - E_EXC_TH_19 * EV_TO_J);  
  828.         g   = sqrt(2.0 * energy / E_MASS);          
  829.         chi = acos(1.0 - 2.0 * R01(MTgen));          
  830.         eta = TWO_PI * R01(MTgen);                  
  831.     } else {                                         // superelastic 2^1S -> 2^3P
  832.         energy = 0.5 * E_MASS * g * g;                  
  833.         energy = fabs(energy + E_EXC_TH_19 * EV_TO_J);  
  834.         g   = sqrt(2.0 * energy / E_MASS);              
  835.         chi = acos(1.0 - 2.0 * R01(MTgen));            
  836.         eta = TWO_PI * R01(MTgen);        
  837.     }
  838.  
  839.     // scatter the primary electron
  840.    
  841.     sc = sin(chi);
  842.     cc = cos(chi);
  843.     se = sin(eta);
  844.     ce = cos(eta);
  845.    
  846.     // compute new relative velocity:
  847.    
  848.     gx = g * (ct * cc - st * sc * ce);
  849.     gy = g * (st * cp * cc + ct * cp * sc * ce - sp * sc * se);
  850.     gz = g * (st * sp * cc + ct * sp * sc * ce + cp * sc * se);
  851.    
  852.     // post-collision velocity of the colliding electron
  853.    
  854.     (*vxe) = wx + F2 * gx;
  855.     (*vye) = wy + F2 * gy;
  856.     (*vze) = wz + F2 * gz;
  857. }
  858.  
  859. //----------------------------------------------------------------------//
  860. // He+ / He collision                                                   //
  861. //----------------------------------------------------------------------//
  862.  
  863. void collision_ion (double *vx_1, double *vy_1, double *vz_1,
  864.                     double *vx_2, double *vy_2, double *vz_2, int e_index){
  865.     double   g,gx,gy,gz,wx,wy,wz,rnd;
  866.     double   theta,phi,chi,eta,st,ct,sp,cp,sc,cc,se,ce,t1,t2;
  867.    
  868.     // calculate relative velocity before collision
  869.     // random Maxwellian target atom already selected (vx_2,vy_2,vz_2 velocity components of target atom come with the call)
  870.    
  871.     gx = (*vx_1)-(*vx_2);
  872.     gy = (*vy_1)-(*vy_2);
  873.     gz = (*vz_1)-(*vz_2);
  874.     g  = sqrt(gx * gx + gy * gy + gz * gz);
  875.     wx = 0.5 * ((*vx_1) + (*vx_2));
  876.     wy = 0.5 * ((*vy_1) + (*vy_2));
  877.     wz = 0.5 * ((*vz_1) + (*vz_2));
  878.    
  879.     // find Euler angles
  880.    
  881.     if (gx == 0) {theta = 0.5 * PI;} else {theta = atan2(sqrt(gy * gy + gz * gz),gx);}
  882.     if (gy == 0) {
  883.         if (gz > 0){phi = 0.5 * PI;} else {phi = - 0.5 * PI;}
  884.     } else {phi = atan2(gz, gy);}
  885.    
  886.     // determine the type of collision based on cross sections and generate scattering angle
  887.    
  888.     t1  =      sigma[I_ISO][e_index];
  889.     t2  = t1 + sigma[I_BACK][e_index];
  890.     rnd = R01(MTgen);
  891.     if  (rnd < (t1 /t2)){                        // isotropic scattering
  892.         chi = acos(1.0 - 2.0 * R01(MTgen));      // scattering angle
  893.     } else {                                     // backward scattering
  894.         chi = PI;                                // scattering angle
  895.     }
  896.     eta = TWO_PI * R01(MTgen);                   // azimuthal angle
  897.     sc  = sin(chi);
  898.     cc  = cos(chi);
  899.     se  = sin(eta);
  900.     ce  = cos(eta);
  901.     st  = sin(theta);
  902.     ct  = cos(theta);
  903.     sp  = sin(phi);
  904.     cp  = cos(phi);
  905.    
  906.     // compute new relative velocity
  907.    
  908.     gx = g * (ct * cc - st * sc * ce);
  909.     gy = g * (st * cp * cc + ct * cp * sc * ce - sp * sc * se);
  910.     gz = g * (st * sp * cc + ct * sp * sc * ce + cp * sc * se);
  911.    
  912.     // post-collision velocity of the ion
  913.    
  914.     (*vx_1) = wx + 0.5 * gx;
  915.     (*vy_1) = wy + 0.5 * gy;
  916.     (*vz_1) = wz + 0.5 * gz;
  917. }
  918.  
  919. //-----------------------------------------------------------------//
  920. // solve Poisson equation (Thomas algorithm)                       //
  921. //-----------------------------------------------------------------//
  922.  
  923. void solve_Poisson (xvector rho1, double tt){
  924.     const double A =  1.0;
  925.     const double B = -2.0;
  926.     const double C =  1.0;
  927.     const double S = 1.0 / (2.0 * DX);
  928.     const double ALPHA = -DX * DX / EPSILON0;
  929.     xvector      g, w, f;
  930.     int          i;
  931.    
  932.     // apply potential to the electrodes - boundary conditions
  933.    
  934.     pot[0]     = VOLTAGE * cos(OMEGA * tt);         // potential at the powered electrode
  935.     pot[N_G-1] = 0.0;                               // potential at the grounded electrode
  936.    
  937.     // solve Poisson equation
  938.    
  939.     for(i=1; i<=N_G-2; i++) f[i] = ALPHA * rho1[i];
  940.     f[1] -= pot[0];
  941.     f[N_G-2] -= pot[N_G-1];
  942.     w[1] = C/B;
  943.     g[1] = f[1]/B;
  944.     for(i=2; i<=N_G-2; i++){
  945.         w[i] = C / (B - A * w[i-1]);
  946.         g[i] = (f[i] - A * g[i-1]) / (B - A * w[i-1]);
  947.     }
  948.     pot[N_G-2] = g[N_G-2];
  949.     for (i=N_G-3; i>0; i--) pot[i] = g[i] - w[i] * pot[i+1];            // potential at the grid points between the electrodes
  950.    
  951.     // compute electric field
  952.    
  953.     for(i=1; i<=N_G-2; i++) efield[i] = (pot[i-1] - pot[i+1]) * S;      // electric field at the grid points between the electrodes
  954.     efield[0]     = (pot[0]     - pot[1])     * INV_DX - rho1[0]     * DX / (2.0 * EPSILON0);   // powered electrode
  955.     efield[N_G-1] = (pot[N_G-2] - pot[N_G-1]) * INV_DX + rho1[N_G-1] * DX / (2.0 * EPSILON0);   // grounded electrode
  956. }
  957.  
  958. //---------------------------------------------------------------------//
  959. // simulation of one radiofrequency cycle                              //
  960. //---------------------------------------------------------------------//
  961.  
  962. void accumulate_rates() {
  963.     double v_sqr, velocity, energy, c0_temp;
  964.     int energy_index, p_temp;
  965.     const double Volume = (ELECTRODE_AREA * DX);
  966.  
  967.     for (int k=0; k<N_e; k++){                              
  968.         v_sqr = vx_e[k] * vx_e[k] + vy_e[k] * vy_e[k] + vz_e[k] * vz_e[k];
  969.         velocity = sqrt(v_sqr);
  970.         energy   = 0.5 * E_MASS * v_sqr / EV_TO_J;
  971.         energy_index = min( int(energy / DE_CS + 0.5), CS_RANGES-1);
  972.  
  973.         c0_temp = x_e[k] * INV_DX;
  974.         p_temp = std::max(0, std::min(N_G-1, static_cast<int>(c0_temp)));
  975.    
  976.         sum_electron_density[p_temp] += WEIGHT / Volume;
  977.  
  978.         // 1^1S -> 2^3S
  979.         sum_rate1f[p_temp] += sigma[E_EXC_1][energy_index] * velocity * WEIGHT;
  980.         sum_rate1b[p_temp] += sigma[E_SUPER_1][energy_index] * velocity * WEIGHT;
  981.         // 1^1S -> 2^1S
  982.         sum_rate2f[p_temp] += sigma[E_EXC_2][energy_index] * velocity * WEIGHT;
  983.         sum_rate2b[p_temp] += sigma[E_SUPER_2][energy_index] * velocity * WEIGHT;
  984.         // 1^1S -> 2^1P
  985.         sum_rate3f[p_temp] += sigma[E_EXC_3][energy_index] * velocity * WEIGHT;
  986.         sum_rate3b[p_temp] += sigma[E_SUPER_3][energy_index] * velocity * WEIGHT;        
  987.         // 1^1S -> 2^3P
  988.         sum_rate4f[p_temp] += sigma[E_EXC_4][energy_index] * velocity * WEIGHT;
  989.         sum_rate4b[p_temp] += sigma[E_SUPER_4][energy_index] * velocity * WEIGHT;    
  990.         // 2^1S -> 2^1P
  991.         sum_rate12f[p_temp] += sigma[E_EXC_12][energy_index] * velocity * WEIGHT;
  992.         sum_rate12b[p_temp] += sigma[E_SUPER_12][energy_index] * velocity * WEIGHT;      
  993.         // 2^3S -> 2^3P
  994.         sum_rate13f[p_temp] += sigma[E_EXC_13][energy_index] * velocity * WEIGHT;
  995.         sum_rate13b[p_temp] += sigma[E_SUPER_13][energy_index] * velocity * WEIGHT;  
  996.         // 2^3S -> 2^1S
  997.         sum_rate17f[p_temp] += sigma[E_EXC_17][energy_index] * velocity * WEIGHT;
  998.         sum_rate17b[p_temp] += sigma[E_SUPER_17][energy_index] * velocity * WEIGHT;      
  999.         // 2^1S -> 2^3P
  1000.         sum_rate19f[p_temp] += sigma[E_EXC_19][energy_index] * velocity * WEIGHT;
  1001.         sum_rate19b[p_temp] += sigma[E_SUPER_19][energy_index] * velocity * WEIGHT;                    
  1002.     }    
  1003. }
  1004.  
  1005. // averaging the rates each RF cycle
  1006. void average_rates() {
  1007.     const double inv_NT = 1.0 / (N_avg * N_T); // averaging through N_avg RF periods each contains N_T cycles
  1008.     const double inv_Volume = 1.0/(ELECTRODE_AREA * DX);
  1009.     for (int p = 0; p < N_G; p++){
  1010.         avg_electron_density[p] = sum_electron_density[p] * inv_NT;
  1011.     }
  1012.     for (int p = 0; p < N_G; p++) {
  1013.         avg_rate1f[p] = sum_rate1f[p] * inv_NT * inv_Volume;  
  1014.         avg_rate1b[p] = sum_rate1b[p] * inv_NT * inv_Volume;
  1015.  
  1016.         avg_rate2f[p] = sum_rate2f[p] * inv_NT * inv_Volume;
  1017.         avg_rate2b[p] = sum_rate2b[p] * inv_NT * inv_Volume;
  1018.  
  1019.         avg_rate3f[p] = sum_rate3f[p] * inv_NT * inv_Volume;
  1020.         avg_rate3b[p] = sum_rate3b[p] * inv_NT * inv_Volume;
  1021.  
  1022.         avg_rate4f[p] = sum_rate4f[p] * inv_NT * inv_Volume;
  1023.         avg_rate4b[p] = sum_rate4b[p] * inv_NT * inv_Volume;
  1024.  
  1025.         avg_rate12f[p] = sum_rate12f[p] * inv_NT * inv_Volume;  
  1026.         avg_rate12b[p] = sum_rate12b[p] * inv_NT * inv_Volume;
  1027.  
  1028.         avg_rate13f[p] = sum_rate13f[p] * inv_NT * inv_Volume;
  1029.         avg_rate13b[p] = sum_rate13b[p] * inv_NT * inv_Volume;
  1030.  
  1031.         avg_rate17f[p] = sum_rate17f[p] * inv_NT * inv_Volume;
  1032.         avg_rate17b[p] = sum_rate17b[p] * inv_NT * inv_Volume;
  1033.  
  1034.         avg_rate19f[p] = sum_rate19f[p] * inv_NT * inv_Volume;
  1035.         avg_rate19b[p] = sum_rate19b[p] * inv_NT * inv_Volume;
  1036.     }    
  1037. }
  1038.  
  1039. void update_excited_dens() {
  1040.  
  1041.     xvector e1_density_bw, e2_density_bw, e3_density_bw, e4_density_bw, gas_dens_local;
  1042.  
  1043.     double S0, S1, S2, S3, S4;                                                   // source terms to account for reaction
  1044.     double diff1, diff2;                                                         // diffusion terms only for metastable states
  1045.  
  1046.     const double v_mean = sqrt(8.0*K_BOLTZMANN*T_neutral/(M_PI*HE_MASS));       // mean spead of He excited states, assuming T~Tg
  1047.     double k20f, k21f, k26f;
  1048.     k20f = 1.799e+09; k21f = 1.976e+06; k26f = 1.022e+07;                       // radiative transitions
  1049.     // k20f = 0.000e+00; k21f = 0.000e+00; k26f = 0.000e+00;                     // radiative transitions
  1050.  
  1051.     bool converged = false;
  1052.  
  1053.     double accuracy = 1.0E-6;                                                   // demanded accuracy
  1054.     double alpha = 0.01;
  1055.     double dt = alpha*DX*DX/Diff_coeff;
  1056.     int max_iter = 1E8;                
  1057.  
  1058.     // constants to calculate reabsorption:
  1059.     double C1 = sqrt(2.0*K_BOLTZMANN*T_neutral/HE_MASS);
  1060.     double C2 = E_CHARGE*E_CHARGE/(4.0*EPSILON0*E_MASS*C_light);
  1061.     // spectral line 2^1P -> 1^1 S
  1062.     double lambda20, d_nu20, phi20, f_lu20, alpha20, tau20, absorb20;               // spectral line constants for for spectral line 2^1P -> 1^1S    
  1063.     lambda20  = 584.334357e-10;                                                     // wavelength of transition
  1064.     d_nu20   =  1.0/lambda20 * C1;                                                  // Doppler broadening delta_nu
  1065.     phi20    = 1.0/(d_nu20*sqrt(M_PI));                                             // Doppler broadening phi
  1066.     f_lu20   = 2.7616e-01;                                                          // oscillator strength                          
  1067.  
  1068.     double lambda21, d_nu21, phi21, f_lu21, alpha21, tau21, absorb21;            // spectral line constants for for spectral line 2^1P -> 2^1S    
  1069.     lambda21  = 20587.0e-10;
  1070.     d_nu21   =  1.0/lambda21 * C1;
  1071.     phi21    = 1.0/(d_nu21*sqrt(M_PI));
  1072.     f_lu21   = 3.764e-01;    
  1073.  
  1074.     double lambda26, d_nu26, phi26, f_lu26, alpha26, tau26, absorb26;            // spectral line constants for for spectral line 2^3P -> 2^3S    
  1075.     lambda26  = 10830.0e-10;
  1076.     d_nu26   =  1.0/lambda26 * C1;
  1077.     phi26    = 1.0/(d_nu26*sqrt(M_PI));
  1078.     f_lu26   = 5.3907e-01;  
  1079.  
  1080.     // creating a helping array for backward step and also works with initial condition
  1081.     std::memcpy(e1_density_bw, e1_density, sizeof(e1_density));      
  1082.     std::memcpy(e2_density_bw, e2_density, sizeof(e2_density));
  1083.     std::memcpy(e3_density_bw, e3_density, sizeof(e3_density));      
  1084.     std::memcpy(e4_density_bw, e4_density, sizeof(e4_density));
  1085.  
  1086.     for (int i = 0; i < N_G; i++){
  1087.         gas_dens_local[i] = GAS_DENSITY;
  1088.     }
  1089.  
  1090.    double rcmb_coeff = gamma_k*v_mean/(2.0*(2.0-gamma_k));
  1091.  
  1092.     for (int idx = 0; idx < max_iter; idx++) {
  1093.  
  1094.         // swapping the arrays
  1095.         std::swap(e1_density_bw, e1_density);
  1096.         std::swap(e2_density_bw, e2_density);            
  1097.         std::swap(e3_density_bw, e3_density);
  1098.         std::swap(e4_density_bw, e4_density);  
  1099.  
  1100.         for (int i = 1; i < N_G-1; i++){
  1101.             // calculate source terms of diffusion equation, S0 for neutral atoms of He
  1102.             // S1 for triplet excited state, S2 for singlet excited state
  1103.             // triplet and singlet are diffusing, neutral atoms are not
  1104.  
  1105.             // S0 = dt*(-(avg_rate1f[i]+avg_rate2f[i] + avg_rate3f[i] + avg_rate4f[i])*gas_dens_local[i]
  1106.             //         + avg_rate1b[i]*e1_density_bw[i] + avg_rate2b[i]*e2_density_bw[i]
  1107.             //         + avg_rate3b[i]*e3_density_bw[i] + avg_rate4b[i]*e4_density_bw[i] + k20f*e3_density_bw[i]);
  1108.  
  1109.             // calculating reabsorption coefficient for 2^1P -> 1^1S
  1110.             alpha20 = C2 * f_lu20 * gas_dens_local[i] * (1.0 - (e3_density_bw[i]*1)/(gas_dens_local[i]*3)) * phi20;
  1111.             tau20 = alpha20*L/2.0;
  1112.             absorb20 = compute_escape_factor(tau20);
  1113.  
  1114.             // calculating reabsorption coefficient for 2^1P -> 2^1S            
  1115.             alpha21 = C2 * f_lu21 * e2_density_bw[i] * (1.0 - (e3_density_bw[i]*1)/(e2_density_bw[i]*3)) * phi21;
  1116.             tau21 = alpha21*L/2.0;
  1117.             absorb21 = compute_escape_factor(tau21);
  1118.  
  1119.             // calculating reabsorption coefficient for 2^3P -> 2^3S            
  1120.             alpha26 = C2 * f_lu26 * e1_density_bw[i] * (1.0 - (e4_density_bw[i]*3)/(e1_density_bw[i]*9)) * phi26;
  1121.             tau26 = alpha26*L/2.0;
  1122.             absorb26 = compute_escape_factor(tau26);
  1123.  
  1124.             S1 = dt*(   avg_rate1f[i]*gas_dens_local[i] + avg_rate13b[i]*e4_density_bw[i] + avg_rate17b[i]*e2_density_bw[i]
  1125.                     - (avg_rate1b[i] + avg_rate13f[i] + avg_rate17f[i])*e1_density_bw[i] + k26f*absorb26*e4_density_bw[i]     )   ;
  1126.  
  1127.             S2 = dt*(avg_rate2f[i]*gas_dens_local[i] + avg_rate12b[i]*e3_density_bw[i]
  1128.                     + avg_rate19b[i]*e4_density_bw[i] + avg_rate17f[i]*e1_density_bw[i]
  1129.                     - (avg_rate2b[i] + avg_rate12f[i] + avg_rate19f[i] + avg_rate17b[i])*e2_density_bw[i]  + k21f*absorb21*e3_density_bw[i]   );
  1130.  
  1131.             S3 = dt*(avg_rate3f[i]*gas_dens_local[i] + avg_rate12f[i]*e2_density_bw[i]
  1132.                     - (avg_rate3b[i] + avg_rate12b[i])*e3_density_bw[i] - (k20f*absorb20+k21f*absorb21)*e3_density_bw[i] );
  1133.  
  1134.             S4 = dt*(avg_rate4f[i]*gas_dens_local[i] + avg_rate13f[i]*e1_density_bw[i] + avg_rate19f[i]*e2_density_bw[i]
  1135.                     - (avg_rate4b[i] + avg_rate13b[i] + avg_rate19b[i])*e4_density_bw[i] - k26f*absorb26*e4_density_bw[i]);        
  1136.  
  1137.             diff1 = e1_density_bw[i]*(1.0 - 2.0*alpha) + alpha*(e1_density_bw[i-1] + e1_density_bw[i+1]);
  1138.             diff2 = e2_density_bw[i]*(1.0 - 2.0*alpha) + alpha*(e2_density_bw[i-1] + e2_density_bw[i+1]);
  1139.  
  1140.  
  1141.             e1_density[i]       = diff1 + S1;
  1142.             e2_density[i]       = diff2 + S2;
  1143.             e3_density[i]       = e3_density_bw[i] + S3;
  1144.             e4_density[i]       = e4_density_bw[i] + S4;            
  1145. //            gas_dens_local[i]  +=         S0;
  1146.  
  1147.             // Prevent negative densities
  1148.             e1_density[i]     = std::max(e1_density[i], 0.0);
  1149.             e2_density[i]     = std::max(e2_density[i], 0.0);
  1150.             e3_density[i]     = std::max(e3_density[i], 0.0);
  1151.             e4_density[i]     = std::max(e4_density[i], 0.0);            
  1152.         }
  1153.  
  1154.         // boundary conditions: now it's zero so no excited states density at the electrodes, yield is zero
  1155.         e1_density[0] = e1_density[1]/(1.0 + rcmb_coeff*DX/Diff_coeff);
  1156.         e1_density[N_G-1] = e1_density[N_G-2]/(1.0 + rcmb_coeff*DX/Diff_coeff);
  1157.         e2_density[0] = e2_density[1]/(1.0 + rcmb_coeff*DX/Diff_coeff);
  1158.         e2_density[N_G-1] = e2_density[N_G-2]/(1.0 + rcmb_coeff*DX/Diff_coeff);  
  1159.  
  1160.         e3_density[0] = e3_density_bw[0] + dt*(avg_rate3f[0]*gas_dens_local[0] + avg_rate12f[0]*e2_density_bw[0]
  1161.                     - (avg_rate3b[0] + avg_rate12b[0])*e3_density_bw[0] - (k20f+k21f)*e3_density_bw[0] );  
  1162.         e3_density[N_G-1] = e3_density_bw[N_G-1] + dt*(avg_rate3f[N_G-1]*gas_dens_local[N_G-1] + avg_rate12f[N_G-1]*e2_density_bw[N_G-1]
  1163.                     - (avg_rate3b[N_G-1] + avg_rate12b[N_G-1])*e3_density_bw[N_G-1] - (k20f+k21f)*e3_density_bw[N_G-1] );
  1164.         e4_density[0] = e4_density_bw[0] + dt*(avg_rate4f[0]*gas_dens_local[0] + avg_rate13f[0]*e1_density_bw[0] + avg_rate19f[0]*e2_density_bw[0]
  1165.                     - (avg_rate4b[0] + avg_rate13b[0] + avg_rate19b[0])*e4_density_bw[0] - k26f*e4_density_bw[0]);  
  1166.         e4_density[N_G-1] = e4_density_bw[N_G-1] + dt*(avg_rate4f[N_G-1]*gas_dens_local[N_G-1] + avg_rate13f[N_G-1]*e1_density_bw[N_G-1] + avg_rate19f[N_G-1]*e2_density_bw[N_G-1]
  1167.                     - (avg_rate4b[N_G-1] + avg_rate13b[N_G-1] + avg_rate19b[N_G-1])*e4_density_bw[N_G-1] - k26f*e4_density_bw[N_G-1]);
  1168.  
  1169.         //checking convergence        
  1170.         double rel1 = 0.0, rel2 = 0.0, rel3 = 0.0, rel4 = 0.0;
  1171.         for (int i = 1; i < N_G-1; ++i) {
  1172.             double abs1 = fabs(e1_density[i] - e1_density_bw[i]);
  1173.             double abs2 = fabs(e2_density[i] - e2_density_bw[i]);
  1174.             double abs3 = fabs(e3_density[i] - e3_density_bw[i]);
  1175.             double abs4 = fabs(e4_density[i] - e4_density_bw[i]);
  1176.             double base1 = std::max(e1_density[i], 1e-12);
  1177.             double base2 = std::max(e2_density[i], 1e-12);
  1178.             double base3 = std::max(e3_density[i], 1e-12);
  1179.             double base4 = std::max(e4_density[i], 1e-12);            
  1180.             rel1 = std::max(rel1, abs1 / e1_density_bw[i]);
  1181.             rel2 = std::max(rel2, abs2 / e2_density_bw[i]);
  1182.             rel3 = std::max(rel3, abs3 / e3_density_bw[i]);
  1183.             rel4 = std::max(rel4, abs4 / e4_density_bw[i]);            
  1184.         }
  1185.         if (rel1 < accuracy && rel2 < accuracy && rel3 < accuracy && rel4 < accuracy) {
  1186.             converged = true;
  1187.             std::cout << "Steady-state reached after " << idx << " iterations.\n";
  1188.             std::cout << "Relative changes are: " << rel1 << " " << rel2 << " " << rel3 << " " << rel4 << "\n";
  1189.             break;
  1190.         }
  1191.     }
  1192.     if (!converged) {
  1193.         std::cerr << "Steady-state not reached after " << max_iter << " iterations.\n";
  1194.     }    
  1195. }
  1196.  
  1197.  
  1198. void do_one_cycle (void){
  1199.     const double DV       = ELECTRODE_AREA * DX;
  1200.     const double FACTOR_W = WEIGHT / DV;
  1201.     const double FACTOR_E = DT_E / E_MASS * E_CHARGE;
  1202.     const double FACTOR_I = DT_I / HE_MASS * E_CHARGE;
  1203.     const double MIN_X    = 0.45 * L;                       // min. position for EEPF collection
  1204.     const double MAX_X    = 0.55 * L;                       // max. position for EEPF collection
  1205.     int      k, t, p, energy_index;
  1206.     double   g, g_sqr, gx, gy, gz, vx_a, vy_a, vz_a, e_x, energy, nu, p_coll, v_sqr, velocity;
  1207.     double   mean_v, c0, c1, c2, rate;
  1208.     bool     out;
  1209.     xvector  rho;
  1210.     int      t_index;
  1211.  
  1212.     nu_avg = 0.0;
  1213.    
  1214.     for (t=0; t<N_T; t++){          // the RF period is divided into N_T equal time intervals (time step DT_E)
  1215.         Time += DT_E;               // update of the total simulated time
  1216.         t_index = t / N_BIN;        // index for XT distributions
  1217.        
  1218.         // step 1: compute densities at grid points
  1219.        
  1220.         for(p=0; p<N_G; p++) e_density[p] = 0;                             // electron density - computed in every time step
  1221.         for(k=0; k<N_e; k++){
  1222.  
  1223.             if      (p < 0)        p = 0;
  1224.             else if (p > N_G - 2)  p = N_G - 2;
  1225.             c0 = x_e[k] * INV_DX;
  1226.             p  = int(c0);
  1227.             e_density[p]   += (p + 1 - c0) * FACTOR_W;
  1228.             e_density[p+1] += (c0 - p) * FACTOR_W;
  1229.         }
  1230.         e_density[0]     *= 2.0; // double at the edge bcs working with half-domain (no left/right neighbour)
  1231.         e_density[N_G-1] *= 2.0;
  1232.         for(p=0; p<N_G; p++) cumul_e_density[p] += e_density[p];
  1233.        
  1234.         if ((t % N_SUB) == 0) {                                            // ion density - computed in every N_SUB-th time steps (subcycling)
  1235.             for(p=0; p<N_G; p++) i_density[p] = 0;
  1236.             for(k=0; k<N_i; k++){
  1237.                 c0 = x_i[k] * INV_DX;
  1238.                 p  = int(c0);
  1239.                 i_density[p]   += (p + 1 - c0) * FACTOR_W;  
  1240.                 i_density[p+1] += (c0 - p) * FACTOR_W;
  1241.             }
  1242.             i_density[0]     *= 2.0; // double at the edge bcs working with half-domain (no left/right neighbour)
  1243.             i_density[N_G-1] *= 2.0;
  1244.         }
  1245.         for(p=0; p<N_G; p++) cumul_i_density[p] += i_density[p];
  1246.        
  1247.         // step 2: solve Poisson equation
  1248.        
  1249.         for(p=0; p<N_G; p++) rho[p] = E_CHARGE * (i_density[p] - e_density[p]);  // get charge density
  1250.         solve_Poisson(rho,Time);                                                 // compute potential and electric field
  1251.        
  1252.         // steps 3 & 4: move particles according to electric field interpolated to particle positions
  1253.        
  1254.         for(k=0; k<N_e; k++){                       // move all electrons in every time step
  1255.             c0  = x_e[k] * INV_DX;
  1256.             p   = int(c0);
  1257.             c1  = p + 1.0 - c0;
  1258.             c2  = c0 - p;
  1259.             e_x = c1 * efield[p] + c2 * efield[p+1];
  1260.            
  1261.             if (measurement_mode) {
  1262.                
  1263.                 // measurements: 'x' and 'v' are needed at the same time, i.e. old 'x' and mean 'v'
  1264.                
  1265.                 mean_v = vx_e[k] - 0.5 * e_x * FACTOR_E;
  1266.                 counter_e_xt[p][t_index]   += c1;
  1267.                 counter_e_xt[p+1][t_index] += c2;
  1268.                 ue_xt[p][t_index]   += c1 * mean_v;
  1269.                 ue_xt[p+1][t_index] += c2 * mean_v;
  1270.                 v_sqr  = mean_v * mean_v + vy_e[k] * vy_e[k] + vz_e[k] * vz_e[k];
  1271.                 energy = 0.5 * E_MASS * v_sqr / EV_TO_J;
  1272.                 meanee_xt[p][t_index]   += c1 * energy;
  1273.                 meanee_xt[p+1][t_index] += c2 * energy;
  1274.                 energy_index = min( int(energy / DE_CS + 0.5), CS_RANGES-1);
  1275.                 velocity = sqrt(v_sqr);
  1276.                 double local_neut_dens = GAS_DENSITY - e1_density[p] - e2_density[p];
  1277.                 rate = sigma[E_ION][energy_index] * velocity * DT_E * local_neut_dens;
  1278.                 ioniz_rate_xt[p][t_index]   += c1 * rate;
  1279.                 ioniz_rate_xt[p+1][t_index] += c2 * rate;
  1280.  
  1281.                 // measure EEPF in the center
  1282.                
  1283.                 if ((MIN_X < x_e[k]) && (x_e[k] < MAX_X)){
  1284.                     energy_index = (int)(energy / DE_EEPF);
  1285.                     if (energy_index < N_EEPF) {eepf[energy_index] += 1.0;}
  1286.                     mean_energy_accu_center += energy;
  1287.                     mean_energy_counter_center++;
  1288.                 }
  1289.             }
  1290.            
  1291.             // update velocity and position
  1292.            
  1293.             vx_e[k] -= e_x * FACTOR_E;
  1294.             x_e[k]  += vx_e[k] * DT_E;
  1295.         }
  1296.        
  1297.         if ((t % N_SUB) == 0) {                       // move all ions in every N_SUB-th time steps (subcycling)
  1298.             for(k=0; k<N_i; k++){
  1299.                 c0  = x_i[k] * INV_DX;
  1300.                 p   = int(c0);
  1301.                 c1  = p + 1 - c0;
  1302.                 c2  = c0 - p;
  1303.                 e_x = c1 * efield[p] + c2 * efield[p+1];
  1304.                
  1305.                 if (measurement_mode) {
  1306.                    
  1307.                     // measurements: 'x' and 'v' are needed at the same time, i.e. old 'x' and mean 'v'
  1308.  
  1309.                     mean_v = vx_i[k] + 0.5 * e_x * FACTOR_I;
  1310.                     counter_i_xt[p][t_index]   += c1;
  1311.                     counter_i_xt[p+1][t_index] += c2;
  1312.                     ui_xt[p][t_index]   += c1 * mean_v;
  1313.                     ui_xt[p+1][t_index] += c2 * mean_v;
  1314.                     v_sqr  = mean_v * mean_v + vy_i[k] * vy_i[k] + vz_i[k] * vz_i[k];
  1315.                     energy = 0.5 * HE_MASS * v_sqr / EV_TO_J;
  1316.                     meanei_xt[p][t_index]   += c1 * energy;
  1317.                     meanei_xt[p+1][t_index] += c2 * energy;
  1318.                 }
  1319.                
  1320.                 // update velocity and position and accumulate absorbed energy
  1321.                
  1322.                 vx_i[k] += e_x * FACTOR_I;
  1323.                 x_i[k]  += vx_i[k] * DT_I;
  1324.             }
  1325.         }
  1326.        
  1327.         // step 5: check boundaries -- changed by Artem
  1328.        
  1329.         k = 0;
  1330.         while(k < N_e) {    // check boundaries for all electrons in every time step
  1331.             out = false;
  1332.             if ((x_e[k] < 0) || (x_e[k] > L)) {out = true;}                   // the electron is out at the electrodes
  1333.             if (out) {
  1334.                 if (R01(MTgen)< r) {                                          // reflecting elastically with r probability
  1335.                     x_e[k] < 0 ? x_e[k] *= -1.0 : x_e[k] = 2.0*L - x_e[k];    // if left/right border -> put on a left/right border
  1336.                     vx_e[k] *= -1.0;                                          // elastic reflection                          
  1337.                     }
  1338.                 else {                                                        // remove the electron
  1339.                     x_e[k] < 0 ? N_e_abs_pow++ : N_e_abs_gnd++;               // absorption calculation at the powered/grounded (left/right) electrode
  1340.                     x_e [k] = x_e [N_e-1];                                    // pushing last element on a vacant place
  1341.                     vx_e[k] = vx_e[N_e-1];
  1342.                     vy_e[k] = vy_e[N_e-1];
  1343.                     vz_e[k] = vz_e[N_e-1];
  1344.                     N_e--;
  1345.                 }
  1346.             } else k++;
  1347.         }
  1348.        
  1349.         if ((t % N_SUB) == 0) {   // check boundaries for all ions in every N_SUB-th time steps (subcycling)
  1350.             k = 0;
  1351.             while(k < N_i) {
  1352.                 out = false;
  1353.                 if (x_i[k] < 0) {       // the ion is out at the powered electrode
  1354.                     N_i_abs_pow++;
  1355.                     out    = true;
  1356.                     v_sqr  = vx_i[k] * vx_i[k] + vy_i[k] * vy_i[k] + vz_i[k] * vz_i[k];
  1357.                     energy = 0.5 * HE_MASS *  v_sqr/ EV_TO_J;
  1358.                     energy_index = (int)(energy / DE_IFED);
  1359.                     if (energy_index < N_IFED) {ifed_pow[energy_index]++;}       // save IFED at the powered electrode
  1360.  
  1361.                     if (R01(MTgen)< gamma_i) {
  1362.                         if (N_e < MAX_N_P){                                                     // adding a new electron on the left electrode
  1363.                             N_e++;                                                              // with energy sampled from Maxwellian distribution
  1364.                             vx_e[N_e-1] = RMB_e(MTgen); vy_e[N_e-1] = RMB_e(MTgen); vz_e[N_e-1] = RMB_e(MTgen);
  1365.                             x_e[N_e-1] = 0.0;
  1366.                         }
  1367.                         else{
  1368.                             std::cerr << "Impossible to emit a new SE at the powered electrode" <<"\n";
  1369.                         }
  1370.                     }
  1371.                 }
  1372.                 if (x_i[k] > L) {       // the ion is out at the grounded electrode
  1373.                     N_i_abs_gnd++;
  1374.                     out    = true;
  1375.                     v_sqr  = vx_i[k] * vx_i[k] + vy_i[k] * vy_i[k] + vz_i[k] * vz_i[k];
  1376.                     energy = 0.5 * HE_MASS * v_sqr / EV_TO_J;
  1377.                     energy_index = (int)(energy / DE_IFED);
  1378.                     if (energy_index < N_IFED) {ifed_gnd[energy_index]++;}        // save IFED at the grounded electrode
  1379.  
  1380.                     if (R01(MTgen)< gamma_i) {
  1381.                         if (N_e < MAX_N_P){                                                     // adding a new electron on the left electrode
  1382.                             N_e++;                                                              // with energy sampled from Maxwellian distribution
  1383.                             vx_e[N_e-1] = RMB_e(MTgen); vy_e[N_e-1] = RMB_e(MTgen); vz_e[N_e-1] = RMB_e(MTgen);
  1384.                             x_e[N_e-1] = L;
  1385.                         }
  1386.                         else{
  1387.                             std::cerr << "Impossible to emit a new SE at the grounded electrode" <<"\n";
  1388.                         }
  1389.                     }                  
  1390.                 }
  1391.  
  1392.                 if (out) {  // delete the ion, if out
  1393.                     x_i [k] = x_i [N_i-1];
  1394.                     vx_i[k] = vx_i[N_i-1];
  1395.                     vy_i[k] = vy_i[N_i-1];
  1396.                     vz_i[k] = vz_i[N_i-1];
  1397.                     N_i--;
  1398.                 } else k++;
  1399.             }
  1400.         }
  1401.        
  1402.         // step 6: collisions
  1403.  
  1404.        
  1405.        
  1406.         for (k=0; k<N_e; k++){                              // checking for occurrence of a collision for all electrons in every time step
  1407.             v_sqr = vx_e[k] * vx_e[k] + vy_e[k] * vy_e[k] + vz_e[k] * vz_e[k];
  1408.             velocity = sqrt(v_sqr);
  1409.             energy   = 0.5 * E_MASS * v_sqr / EV_TO_J;
  1410.             energy_index = min( int(energy / DE_CS + 0.5), CS_RANGES-1);
  1411. //            double ground_dens_local = GAS_DENSITY; // let the ground density be fixed in each cell according to our DRR solver
  1412.             // Artem  - adding superelastic impact on total collisoinal probability//
  1413.  
  1414.            
  1415.             int e_crdnt = std::max(0, std::min(N_G-1, static_cast<int>(x_e[k] * INV_DX)));
  1416.             double sigma_super_e = sigma[E_SUPER_1][energy_index] * e1_density[e_crdnt] + sigma[E_SUPER_2][energy_index] * e2_density[e_crdnt]
  1417.             + sigma[E_SUPER_3][energy_index] * e3_density[e_crdnt] + sigma[E_SUPER_4][energy_index] * e4_density[e_crdnt]
  1418.             + sigma[E_SUPER_12][energy_index] * e3_density[e_crdnt] + sigma[E_SUPER_13][energy_index] * e4_density[e_crdnt]
  1419.             + sigma[E_SUPER_17][energy_index] * e2_density[e_crdnt] + sigma[E_SUPER_19][energy_index] * e4_density[e_crdnt];
  1420.             double ground_dens_local = GAS_DENSITY - e1_density[e_crdnt] - e2_density[e_crdnt] - e3_density[e_crdnt] - e4_density[e_crdnt];
  1421.                                      
  1422.             double sigma_ground         = (sigma[E_ELA][energy_index] + sigma[E_EXC_1][energy_index]
  1423.                                             + sigma[E_EXC_2][energy_index] + sigma[E_ION][energy_index]
  1424.                                             + sigma[E_EXC_3][energy_index] + sigma[E_EXC_4][energy_index]) * ground_dens_local;
  1425.             double sigma_stepwise       = sigma[E_EXC_12][energy_index] * e2_density[e_crdnt] +  sigma[E_EXC_13][energy_index] * e1_density[e_crdnt]
  1426.                                             + sigma[E_EXC_17][energy_index] * e1_density[e_crdnt] + sigma[E_EXC_19][energy_index] * e2_density[e_crdnt];
  1427.  
  1428.             nu = (sigma_ground + sigma_super_e + sigma_stepwise) * velocity;
  1429.  
  1430.             nu_avg += nu;
  1431.            
  1432.             p_coll = 1 - exp(- nu * DT_E);                  // collision probability for electrons
  1433.             if (R01(MTgen) < p_coll) {                      // electron collision takes place
  1434.                 collision_electron(x_e[k], &vx_e[k], &vy_e[k], &vz_e[k], energy_index);
  1435.                 N_e_coll++;
  1436.             }
  1437.         }
  1438.        
  1439.         if ((t % N_SUB) == 0) {                             // checking for occurrence of a collision for all ions in every N_SUB-th time steps (subcycling)
  1440.             for (k=0; k<N_i; k++){
  1441.                 vx_a = RMB_n(MTgen);                          // pick velocity components of a random target gas atom
  1442.                 vy_a = RMB_n(MTgen);
  1443.                 vz_a = RMB_n(MTgen);
  1444.                 gx   = vx_i[k] - vx_a;                       // compute the relative velocity of the collision partners
  1445.                 gy   = vy_i[k] - vy_a;
  1446.                 gz   = vz_i[k] - vz_a;
  1447.                 g_sqr = gx * gx + gy * gy + gz * gz;
  1448.                 g = sqrt(g_sqr);
  1449.                 energy = 0.5 * MU_HEHE * g_sqr / EV_TO_J;
  1450.                 energy_index = min( int(energy / DE_CS + 0.5), CS_RANGES-1);
  1451.                 nu = sigma_tot_i[energy_index] * g;
  1452.                 p_coll = 1 - exp(- nu * DT_I);              // collision probability for ions
  1453.                 if (R01(MTgen)< p_coll) {                   // ion collision takes place
  1454.                     collision_ion (&vx_i[k], &vy_i[k], &vz_i[k], &vx_a, &vy_a, &vz_a, energy_index);
  1455.                     N_i_coll++;
  1456.                 }
  1457.             }
  1458.         }
  1459.  
  1460.         //step 7: accumulate rates
  1461.         accumulate_rates();
  1462.  
  1463.         if (measurement_mode) {
  1464.            
  1465.             // collect 'xt' data from the grid
  1466.            
  1467.             for (p=0; p<N_G; p++) {
  1468.                 pot_xt   [p][t_index] += pot[p];
  1469.                 efield_xt[p][t_index] += efield[p];
  1470.                 ne_xt    [p][t_index] += e_density[p];
  1471.                 ni_xt    [p][t_index] += i_density[p];
  1472.                 e1_xt    [p][t_index] += e1_density[p];  // Artem
  1473.                 e2_xt    [p][t_index] += e2_density[p];  // Artem
  1474.             }
  1475.         }
  1476.        
  1477.         if ((t % 1000) == 0) printf(" c = %8d  t = %8d  #e = %8d  #i = %8d\n", cycle,t,N_e,N_i);
  1478.     }
  1479.  
  1480.     counter++;
  1481.  
  1482.     // updating denisites each N_avg cycles: --- Artem
  1483.     if (counter%N_avg == 0) {
  1484.         // compute average rates over a cycle
  1485.         average_rates();
  1486.         // updating densities
  1487.         update_excited_dens();
  1488.         // reset accumulators
  1489.         memset(sum_rate1f,   0, sizeof sum_rate1f);
  1490.         memset(sum_rate1b,   0, sizeof sum_rate1b);
  1491.         memset(sum_rate2f,   0, sizeof sum_rate2f);
  1492.         memset(sum_rate2b,   0, sizeof sum_rate2b);
  1493.         memset(sum_rate3f,   0, sizeof sum_rate3f);
  1494.         memset(sum_rate3b,   0, sizeof sum_rate3b);
  1495.         memset(sum_rate4f,   0, sizeof sum_rate4f);
  1496.         memset(sum_rate4b,   0, sizeof sum_rate4b);
  1497.         memset(sum_rate12f,  0, sizeof sum_rate12f);
  1498.         memset(sum_rate12b,  0, sizeof sum_rate12b);
  1499.         memset(sum_rate13f,  0, sizeof sum_rate13f);
  1500.         memset(sum_rate13b,  0, sizeof sum_rate13b);
  1501.         memset(sum_rate17f,  0, sizeof sum_rate17f);
  1502.         memset(sum_rate17b,  0, sizeof sum_rate17b);
  1503.         memset(sum_rate19f,  0, sizeof sum_rate19f);
  1504.         memset(sum_rate19b,  0, sizeof sum_rate19b);
  1505.         memset(sum_electron_density, 0, sizeof sum_electron_density);
  1506.     }    
  1507.  
  1508.  
  1509.     //calculate nu electron average:
  1510.  
  1511.     nu_avg /= (N_T*N_e);
  1512.  
  1513.     fprintf(datafile,"%8d  %8d  %8d\n",cycle,N_e,N_i);
  1514.     print_excitation_densities();
  1515. }
  1516.  
  1517. //---------------------------------------------------------------------//
  1518. // save particle coordinates                                           //
  1519. //---------------------------------------------------------------------//
  1520.  
  1521. void save_particle_data(){
  1522.     double   d;
  1523.     FILE   * f;
  1524.     char fname[80];
  1525.    
  1526.     strcpy(fname,"picdata.bin");
  1527.     f = fopen(fname,"wb");
  1528.     fwrite(&Time,sizeof(double),1,f);
  1529.     d = (double)(cycles_done);
  1530.     fwrite(&d,sizeof(double),1,f);
  1531.     d = (double)(N_e);
  1532.     fwrite(&d,sizeof(double),1,f);
  1533.     fwrite(x_e, sizeof(double),N_e,f);
  1534.     fwrite(vx_e,sizeof(double),N_e,f);
  1535.     fwrite(vy_e,sizeof(double),N_e,f);
  1536.     fwrite(vz_e,sizeof(double),N_e,f);
  1537.     d = (double)(N_i);
  1538.     fwrite(&d,sizeof(double),1,f);
  1539.     fwrite(x_i, sizeof(double),N_i,f);
  1540.     fwrite(vx_i,sizeof(double),N_i,f);
  1541.     fwrite(vy_i,sizeof(double),N_i,f);
  1542.     fwrite(vz_i,sizeof(double),N_i,f);
  1543.  
  1544.     // saving excited state densities and rates - Artem
  1545.  
  1546.     fwrite(e1_density,  sizeof(double), N_G, f); fwrite(e2_density,  sizeof(double), N_G, f);
  1547.     fwrite(e3_density,  sizeof(double), N_G, f); fwrite(e4_density,  sizeof(double), N_G, f);
  1548.  
  1549.     fwrite(sum_rate1f,  sizeof(double), N_G, f); fwrite(sum_rate1b,  sizeof(double), N_G, f);
  1550.     fwrite(sum_rate2f,  sizeof(double), N_G, f); fwrite(sum_rate2b,  sizeof(double), N_G, f);  
  1551.     fwrite(sum_rate3f,  sizeof(double), N_G, f); fwrite(sum_rate3b,  sizeof(double), N_G, f);
  1552.     fwrite(sum_rate4f,  sizeof(double), N_G, f); fwrite(sum_rate4b,  sizeof(double), N_G, f);      
  1553.     fwrite(sum_rate12f, sizeof(double), N_G, f); fwrite(sum_rate12b, sizeof(double), N_G, f);
  1554.     fwrite(sum_rate13f, sizeof(double), N_G, f); fwrite(sum_rate13b, sizeof(double), N_G, f);  
  1555.     fwrite(sum_rate17f, sizeof(double), N_G, f); fwrite(sum_rate17b, sizeof(double), N_G, f);
  1556.     fwrite(sum_rate19f, sizeof(double), N_G, f); fwrite(sum_rate19b, sizeof(double), N_G, f);
  1557.  
  1558.     fwrite(avg_rate1f,  sizeof(double), N_G, f); fwrite(avg_rate1b,  sizeof(double), N_G, f);
  1559.     fwrite(avg_rate2f,  sizeof(double), N_G, f); fwrite(avg_rate2b,  sizeof(double), N_G, f);    
  1560.     fwrite(avg_rate3f,  sizeof(double), N_G, f); fwrite(avg_rate3b,  sizeof(double), N_G, f);
  1561.     fwrite(avg_rate4f,  sizeof(double), N_G, f); fwrite(avg_rate4b,  sizeof(double), N_G, f);
  1562.     fwrite(avg_rate12f, sizeof(double), N_G, f); fwrite(avg_rate12b, sizeof(double), N_G, f);
  1563.     fwrite(avg_rate13f, sizeof(double), N_G, f); fwrite(avg_rate13b, sizeof(double), N_G, f);
  1564.     fwrite(avg_rate17f, sizeof(double), N_G, f); fwrite(avg_rate17b, sizeof(double), N_G, f);
  1565.     fwrite(avg_rate19f, sizeof(double), N_G, f); fwrite(avg_rate19b, sizeof(double), N_G, f);
  1566.  
  1567.     fwrite(sum_electron_density, sizeof(double), N_G, f);
  1568.     fwrite(avg_electron_density, sizeof(double), N_G, f);
  1569.  
  1570.     fclose(f);
  1571.     printf(">> eduPIC: data saved : %d electrons %d ions, excited states densities , %d cycles completed, time is %e [s]\n",N_e,N_i,cycles_done,Time);
  1572. }
  1573.  
  1574. //---------------------------------------------------------------------//
  1575. // load particle coordinates                                           //
  1576. //---------------------------------------------------------------------//
  1577.  
  1578. void load_particle_data(){
  1579.     double   d;
  1580.     FILE   * f;
  1581.     char fname[80];
  1582.    
  1583.     strcpy(fname,"picdata.bin");    
  1584.     f = fopen(fname,"rb");
  1585.     if (f==NULL) {printf(">> eduPIC: ERROR: No particle data file found, try running initial cycle using argument '0'\n"); exit(0); }
  1586.     fread(&Time,sizeof(double),1,f);
  1587.     fread(&d,sizeof(double),1,f);
  1588.     cycles_done = int(d);
  1589.     fread(&d,sizeof(double),1,f);
  1590.     N_e = int(d);
  1591.     fread(x_e, sizeof(double),N_e,f);
  1592.     fread(vx_e,sizeof(double),N_e,f);
  1593.     fread(vy_e,sizeof(double),N_e,f);
  1594.     fread(vz_e,sizeof(double),N_e,f);
  1595.     fread(&d,sizeof(double),1,f);
  1596.     N_i = int(d);
  1597.     fread(x_i, sizeof(double),N_i,f);
  1598.     fread(vx_i,sizeof(double),N_i,f);
  1599.     fread(vy_i,sizeof(double),N_i,f);
  1600.     fread(vz_i,sizeof(double),N_i,f);
  1601.  
  1602.     // reading excited states densities and rates info -- Artem
  1603.  
  1604.     fread(e1_density, sizeof(double), N_G, f);
  1605.     fread(e2_density, sizeof(double), N_G, f);    
  1606.     fread(e3_density, sizeof(double), N_G, f);
  1607.     fread(e4_density, sizeof(double), N_G, f);
  1608.  
  1609.     fread(sum_rate1f,  sizeof(double), N_G, f); fread(sum_rate1b,  sizeof(double), N_G, f);
  1610.     fread(sum_rate2f,  sizeof(double), N_G, f); fread(sum_rate2b,  sizeof(double), N_G, f);  
  1611.     fread(sum_rate3f,  sizeof(double), N_G, f); fread(sum_rate3b,  sizeof(double), N_G, f);
  1612.     fread(sum_rate4f,  sizeof(double), N_G, f); fread(sum_rate4b,  sizeof(double), N_G, f);      
  1613.     fread(sum_rate12f, sizeof(double), N_G, f); fread(sum_rate12b, sizeof(double), N_G, f);
  1614.     fread(sum_rate13f, sizeof(double), N_G, f); fread(sum_rate13b, sizeof(double), N_G, f);  
  1615.     fread(sum_rate17f, sizeof(double), N_G, f); fread(sum_rate17b, sizeof(double), N_G, f);
  1616.     fread(sum_rate19f, sizeof(double), N_G, f); fread(sum_rate19b, sizeof(double), N_G, f);
  1617.  
  1618.     fread(avg_rate1f,  sizeof(double), N_G, f); fread(avg_rate1b,  sizeof(double), N_G, f);
  1619.     fread(avg_rate2f,  sizeof(double), N_G, f); fread(avg_rate2b,  sizeof(double), N_G, f);    
  1620.     fread(avg_rate3f,  sizeof(double), N_G, f); fread(avg_rate3b,  sizeof(double), N_G, f);
  1621.     fread(avg_rate4f,  sizeof(double), N_G, f); fread(avg_rate4b,  sizeof(double), N_G, f);
  1622.     fread(avg_rate12f, sizeof(double), N_G, f); fread(avg_rate12b, sizeof(double), N_G, f);
  1623.     fread(avg_rate13f, sizeof(double), N_G, f); fread(avg_rate13b, sizeof(double), N_G, f);
  1624.     fread(avg_rate17f, sizeof(double), N_G, f); fread(avg_rate17b, sizeof(double), N_G, f);
  1625.     fread(avg_rate19f, sizeof(double), N_G, f); fread(avg_rate19b, sizeof(double), N_G, f);    
  1626.  
  1627.     fread(sum_electron_density, sizeof(double), N_G, f);
  1628.     fread(avg_electron_density, sizeof(double), N_G, f);
  1629.  
  1630.     fclose(f);
  1631.     printf(">> eduPIC: data loaded : %d electrons %d ions, excited states densities, %d cycles completed before, time is %e [s]\n",N_e,N_i,cycles_done,Time);
  1632. }
  1633.  
  1634. //---------------------------------------------------------------------//
  1635. // save density data                                                   //
  1636. //---------------------------------------------------------------------//
  1637.  
  1638. void save_density(void){
  1639.     FILE *f;
  1640.     double c;
  1641.     int m;
  1642.    
  1643.     f = fopen("density.dat","w");
  1644.     c = 1.0 / (double)(no_of_cycles) / (double)(N_T);
  1645.     for(m=0; m<N_G; m++){
  1646.         fprintf(f,"%8.5f  %12e  %12e\n",m * DX, cumul_e_density[m] * c, cumul_i_density[m] * c);
  1647.     }
  1648.     fclose(f);
  1649. }
  1650.  
  1651. // save Final excited states densities - Artem
  1652.  
  1653. void save_final_excited_densities(void) {
  1654.     FILE *f = fopen("excited_densities.dat", "w");
  1655.     for(int p=0; p<N_G; p++) {
  1656.         fprintf(f, "%8.5f %12e %12e %12e %12e\n", p*DX, e1_density[p], e2_density[p], e3_density[p], e4_density[p]);
  1657.     }
  1658.     fclose(f);
  1659. }
  1660.  
  1661. //---------------------------------------------------------------------//
  1662. // save EEPF data                                                      //
  1663. //---------------------------------------------------------------------//
  1664.  
  1665. void save_eepf(void) {
  1666.     FILE   *f;
  1667.     int    i;
  1668.     double h,energy;
  1669.    
  1670.     h = 0.0;
  1671.     for (i=0; i<N_EEPF; i++) {h += eepf[i];}
  1672.     h *= DE_EEPF;
  1673.     f = fopen("eepf.dat","w");
  1674.     for (i=0; i<N_EEPF; i++) {
  1675.         energy = (i + 0.5) * DE_EEPF;
  1676.         fprintf(f,"%e  %e\n", energy, eepf[i] / h / sqrt(energy));
  1677.     }
  1678.     fclose(f);
  1679. }
  1680.  
  1681. //---------------------------------------------------------------------//
  1682. // save IFED data                                                      //
  1683. //---------------------------------------------------------------------//
  1684.  
  1685. void save_ifed(void) {
  1686.     FILE   *f;
  1687.     int    i;
  1688.     double h_pow,h_gnd,energy;
  1689.    
  1690.     h_pow = 0.0;
  1691.     h_gnd = 0.0;
  1692.     for (i=0; i<N_IFED; i++) {h_pow += ifed_pow[i]; h_gnd += ifed_gnd[i];}
  1693.     h_pow *= DE_IFED;
  1694.     h_gnd *= DE_IFED;
  1695.     mean_i_energy_pow = 0.0;
  1696.     mean_i_energy_gnd = 0.0;
  1697.     f = fopen("ifed.dat","w");
  1698.     for (i=0; i<N_IFED; i++) {
  1699.         energy = (i + 0.5) * DE_IFED;
  1700.         fprintf(f,"%6.2f %10.6f %10.6f\n", energy, (double)(ifed_pow[i])/h_pow, (double)(ifed_gnd[i])/h_gnd);
  1701.         mean_i_energy_pow += energy * (double)(ifed_pow[i]) / h_pow;
  1702.         mean_i_energy_gnd += energy * (double)(ifed_gnd[i]) / h_gnd;
  1703.     }
  1704.     fclose(f);
  1705. }
  1706.  
  1707. //--------------------------------------------------------------------//
  1708. // save XT data                                                       //
  1709. //--------------------------------------------------------------------//
  1710.  
  1711. void save_xt_1(xt_distr distr, char *fname) {
  1712.     FILE   *f;
  1713.     int    i, j;
  1714.    
  1715.     f = fopen(fname,"w");
  1716.     for (i=0; i<N_G; i++){
  1717.         for (j=0; j<N_XT; j++){
  1718.             fprintf(f,"%e  ", distr[i][j]);
  1719.         }
  1720.         fprintf(f,"\n");
  1721.     }
  1722.     fclose(f);
  1723. }
  1724.  
  1725. void norm_all_xt(void){
  1726.     double f1, f2;
  1727.     int    i, j;
  1728.    
  1729.     // normalize all XT data
  1730.    
  1731.     f1 = (double)(N_XT) / (double)(no_of_cycles * N_T);
  1732.     f2 = WEIGHT / (ELECTRODE_AREA * DX) / (no_of_cycles * (PERIOD / (double)(N_XT)));
  1733.    
  1734.     for (i=0; i<N_G; i++){
  1735.         for (j=0; j<N_XT; j++){
  1736.             pot_xt[i][j]    *= f1;
  1737.             efield_xt[i][j] *= f1;
  1738.             ne_xt[i][j]     *= f1;
  1739.             ni_xt[i][j]     *= f1;
  1740.             e1_xt[i][j]     *= f1;   // Artem
  1741.             e2_xt[i][j]     *= f1;   // Artem
  1742.             if (counter_e_xt[i][j] > 0) {
  1743.                 ue_xt[i][j]     =  ue_xt[i][j] / counter_e_xt[i][j];
  1744.                 je_xt[i][j]     = -ue_xt[i][j] * ne_xt[i][j] * E_CHARGE;
  1745.                 meanee_xt[i][j] =  meanee_xt[i][j] / counter_e_xt[i][j];
  1746.                 ioniz_rate_xt[i][j] *= f2;
  1747.              } else {
  1748.                 ue_xt[i][j]         = 0.0;
  1749.                 je_xt[i][j]         = 0.0;
  1750.                 meanee_xt[i][j]     = 0.0;
  1751.                 ioniz_rate_xt[i][j] = 0.0;
  1752.             }
  1753.             if (counter_i_xt[i][j] > 0) {
  1754.                 ui_xt[i][j]     = ui_xt[i][j] / counter_i_xt[i][j];
  1755.                 ji_xt[i][j]     = ui_xt[i][j] * ni_xt[i][j] * E_CHARGE;
  1756.                 meanei_xt[i][j] = meanei_xt[i][j] / counter_i_xt[i][j];
  1757.             } else {
  1758.                 ui_xt[i][j]     = 0.0;
  1759.                 ji_xt[i][j]     = 0.0;
  1760.                 meanei_xt[i][j] = 0.0;
  1761.             }
  1762.             powere_xt[i][j] = je_xt[i][j] * efield_xt[i][j];
  1763.             poweri_xt[i][j] = ji_xt[i][j] * efield_xt[i][j];
  1764.         }
  1765.     }
  1766. }
  1767.  
  1768. void save_all_xt(void){
  1769.     char fname[80];
  1770.    
  1771.     strcpy(fname,"pot_xt.dat");     save_xt_1(pot_xt, fname);
  1772.     strcpy(fname,"efield_xt.dat");  save_xt_1(efield_xt, fname);
  1773.     strcpy(fname,"ne_xt.dat");      save_xt_1(ne_xt, fname);
  1774.     strcpy(fname,"ni_xt.dat");      save_xt_1(ni_xt, fname);
  1775.     strcpy(fname,"je_xt.dat");      save_xt_1(je_xt, fname);
  1776.     strcpy(fname,"ji_xt.dat");      save_xt_1(ji_xt, fname);
  1777.     strcpy(fname,"powere_xt.dat");  save_xt_1(powere_xt, fname);
  1778.     strcpy(fname,"poweri_xt.dat");  save_xt_1(poweri_xt, fname);
  1779.     strcpy(fname,"meanee_xt.dat");  save_xt_1(meanee_xt, fname);
  1780.     strcpy(fname,"meanei_xt.dat");  save_xt_1(meanei_xt, fname);
  1781.     strcpy(fname,"ioniz_xt.dat");   save_xt_1(ioniz_rate_xt, fname);
  1782.     strcpy(fname,"e1_xt.dat");      save_xt_1(e1_xt, fname);
  1783.     strcpy(fname,"e2_xt.dat");      save_xt_1(e2_xt, fname);
  1784.     strcpy(fname,"e3_xt.dat");      save_xt_1(e3_xt, fname);
  1785.     strcpy(fname,"e4_xt.dat");      save_xt_1(e4_xt, fname);    
  1786.     //add
  1787. }
  1788.  
  1789. //---------------------------------------------------------------------//
  1790. // simulation report including stability and accuracy conditions       //
  1791. //---------------------------------------------------------------------//
  1792.  
  1793. void check_and_save_info(void){
  1794.     FILE     *f;
  1795.     double   plas_freq, meane, kT, debye_length, density, ecoll_freq, icoll_freq, sim_time, e_max, v_max, power_e, power_i, c;
  1796.     int      i,j;
  1797.     bool     conditions_OK;
  1798.    
  1799.     density    = cumul_e_density[N_G / 2] / (double)(no_of_cycles) / (double)(N_T);  // e density @ center
  1800.     plas_freq  = E_CHARGE * sqrt(density / EPSILON0 / E_MASS);                       // e plasma frequency @ center
  1801.     meane      = mean_energy_accu_center / (double)(mean_energy_counter_center);     // e mean energy @ center
  1802.     kT         = 2.0 * meane * EV_TO_J / 3.0;                                        // k T_e @ center (approximate)
  1803.     sim_time   = (double)(no_of_cycles) / FREQUENCY;                                 // simulated time
  1804.     ecoll_freq = (double)(N_e_coll) / sim_time / (double)(N_e);                      // e collision frequency
  1805.     icoll_freq = (double)(N_i_coll) / sim_time / (double)(N_i);                      // ion collision frequency
  1806.     debye_length = sqrt(EPSILON0 * kT / density) / E_CHARGE;                         // e Debye length @ center
  1807.    
  1808.     f = fopen("info.txt","w");
  1809.     fprintf(f,"########################## eduPIC simulation report ############################\n");
  1810.     fprintf(f,"Simulation parameters:\n");
  1811.     fprintf(f,"Gap distance                          = %12.3e [m]\n",  L);
  1812.     fprintf(f,"# of grid divisions                   = %12d\n",      N_G);
  1813.     fprintf(f,"Frequency                             = %12.3e [Hz]\n", FREQUENCY);
  1814.     fprintf(f,"# of time steps / period              = %12d\n",      N_T);
  1815.     fprintf(f,"# of electron / ion time steps        = %12d\n",      N_SUB);
  1816.     fprintf(f,"Voltage amplitude                     = %12.3e [V]\n",  VOLTAGE);
  1817.     fprintf(f,"Pressure (Ar)                         = %12.3e [Pa]\n", PRESSURE);
  1818.     fprintf(f,"Temperature                           = %12.3e [K]\n",  T_neutral);
  1819.     fprintf(f,"Superparticle weight                  = %12.3e\n",      WEIGHT);
  1820.     fprintf(f,"# of simulation cycles in this run    = %12d\n",      no_of_cycles);
  1821.     fprintf(f,"--------------------------------------------------------------------------------\n");
  1822.     fprintf(f,"Plasma characteristics:\n");
  1823.     fprintf(f,"Electron density @ center             = %12.3e [m^{-3}]\n", density);
  1824.     fprintf(f,"Plasma frequency @ center             = %12.3e [rad/s]\n",  plas_freq);
  1825.     fprintf(f,"Debye length @ center                 = %12.3e [m]\n",      debye_length);
  1826.     fprintf(f,"Electron collision frequency          = %12.3e [1/s]\n",    ecoll_freq);
  1827.     fprintf(f,"Ion collision frequency               = %12.3e [1/s]\n",    icoll_freq);
  1828.     fprintf(f,"--------------------------------------------------------------------------------\n");
  1829.     fprintf(f,"Stability and accuracy conditions:\n");
  1830.     conditions_OK = true;
  1831.     c = plas_freq * DT_E;
  1832.     fprintf(f,"Plasma frequency @ center * DT_E      = %12.3f (OK if less than 0.20)\n", c);
  1833.     if (c > 0.2) {conditions_OK = false;}
  1834.     c = DX / debye_length;
  1835.     fprintf(f,"DX / Debye length @ center            = %12.3f (OK if less than 1.00)\n", c);
  1836.     if (c > 1.0) {conditions_OK = false;}
  1837.     c = max_electron_coll_freq() * DT_E;
  1838.     fprintf(f,"Max. electron coll. frequency * DT_E  = %12.3f (OK if less than 0.05)\n", c);
  1839.     if (c > 0.05) {conditions_OK = false;}
  1840.     c = max_ion_coll_freq() * DT_I;
  1841.     fprintf(f,"Max. ion coll. frequency * DT_I       = %12.3f (OK if less than 0.05)\n", c);
  1842.     if (c > 0.05) {conditions_OK = false;}
  1843.     if (conditions_OK == false){
  1844.         fprintf(f,"--------------------------------------------------------------------------------\n");
  1845.         fprintf(f,"** STABILITY AND ACCURACY CONDITION(S) VIOLATED - REFINE SIMULATION SETTINGS! **\n");
  1846.         fprintf(f,"--------------------------------------------------------------------------------\n");
  1847.         fclose(f);
  1848.         printf(">> eduPIC: ERROR: STABILITY AND ACCURACY CONDITION(S) VIOLATED!\n");
  1849.         printf(">> eduPIC: for details see 'info.txt' and refine simulation settings!\n");
  1850.     }
  1851.     else
  1852.     {
  1853.         // calculate maximum energy for which the Courant-Friedrichs-Levy condition holds:
  1854.        
  1855.         v_max = DX / DT_E;
  1856.         e_max = 0.5 * E_MASS * v_max * v_max / EV_TO_J;
  1857.         fprintf(f,"Max e- energy for CFL condition       = %12.3f [eV]\n", e_max);
  1858.         fprintf(f,"Check EEPF to ensure that CFL is fulfilled for the majority of the electrons!\n");
  1859.         fprintf(f,"--------------------------------------------------------------------------------\n");
  1860.        
  1861.         // saving of the following data is done here as some of the further lines need data
  1862.         // that are computed / normalized in these functions
  1863.        
  1864.         printf(">> eduPIC: saving diagnostics data\n");
  1865.         save_density();
  1866.         save_final_excited_densities();   // Artem
  1867.         save_eepf();
  1868.         save_ifed();
  1869.         norm_all_xt();
  1870.         save_all_xt();
  1871.         fprintf(f,"Particle characteristics at the electrodes:\n");
  1872.         fprintf(f,"Ion flux at powered electrode         = %12.3e [m^{-2} s^{-1}]\n", N_i_abs_pow * WEIGHT / ELECTRODE_AREA / (no_of_cycles * PERIOD));
  1873.         fprintf(f,"Ion flux at grounded electrode        = %12.3e [m^{-2} s^{-1}]\n", N_i_abs_gnd * WEIGHT / ELECTRODE_AREA / (no_of_cycles * PERIOD));
  1874.         fprintf(f,"Mean ion energy at powered electrode  = %12.3e [eV]\n", mean_i_energy_pow);
  1875.         fprintf(f,"Mean ion energy at grounded electrode = %12.3e [eV]\n", mean_i_energy_gnd);
  1876.         fprintf(f,"Electron flux at powered electrode    = %12.3e [m^{-2} s^{-1}]\n", N_e_abs_pow * WEIGHT / ELECTRODE_AREA / (no_of_cycles * PERIOD));
  1877.         fprintf(f,"Electron flux at grounded electrode   = %12.3e [m^{-2} s^{-1}]\n", N_e_abs_gnd * WEIGHT / ELECTRODE_AREA / (no_of_cycles * PERIOD));
  1878.         fprintf(f,"--------------------------------------------------------------------------------\n");
  1879.        
  1880.         // calculate spatially and temporally averaged power absorption by the electrons and ions
  1881.        
  1882.         power_e = 0.0;
  1883.         power_i = 0.0;
  1884.         for (i=0; i<N_G; i++){
  1885.             for (j=0; j<N_XT; j++){
  1886.                 power_e += powere_xt[i][j];
  1887.                 power_i += poweri_xt[i][j];
  1888.             }
  1889.         }
  1890.         power_e /= (double)(N_XT * N_G);
  1891.         power_i /= (double)(N_XT * N_G);
  1892.         fprintf(f,"Absorbed power calculated as <j*E>:\n");
  1893.         fprintf(f,"Electron power density (average)      = %12.3e [W m^{-3}]\n", power_e);
  1894.         fprintf(f,"Ion power density (average)           = %12.3e [W m^{-3}]\n", power_i);
  1895.         fprintf(f,"Total power density(average)          = %12.3e [W m^{-3}]\n", power_e + power_i);
  1896.         fprintf(f,"--------------------------------------------------------------------------------\n");
  1897.         fclose(f);
  1898.     }
  1899. }
  1900.  
  1901. //------------------------------------------------------------------------------------------//
  1902. // main                                                                                     //
  1903. // command line arguments:                                                                  //
  1904. // [1]: number of cycles (0 for init)                                                       //
  1905. // [2]: "m" turns on data collection and saving                                             //
  1906. //------------------------------------------------------------------------------------------//
  1907.  
  1908. int main (int argc, char *argv[]){
  1909.  
  1910.     if (argc == 1) {
  1911.         printf(">> eduPIC: error = need starting_cycle argument\n");
  1912.         return 1;
  1913.     } else {
  1914.         strcpy(st0,argv[1]);
  1915.         arg1 = atol(st0);
  1916.         if (argc > 2) {
  1917.             if (strcmp (argv[2],"m") == 0){
  1918.                 measurement_mode = true;                  // measurements will be done
  1919.             } else {
  1920.                 measurement_mode = false;
  1921.             }
  1922.         }
  1923.     }
  1924.     if (measurement_mode) {
  1925.         printf(">> eduPIC: measurement mode: on\n");
  1926.     } else {
  1927.         printf(">> eduPIC: measurement mode: off\n");
  1928.     }
  1929.     set_electron_cross_sections_ar();
  1930.     set_ion_cross_sections_ar();
  1931.     calc_total_cross_sections();
  1932.        
  1933.     auto exc = init_excited_distr();
  1934.     printf("Excited state population initialized!\n");
  1935.  
  1936.  
  1937.     std::ofstream file0("rates.dat"); file0 << std::scientific << std::setprecision(6);
  1938.     std::ofstream file1("excited_densities.dat"); file1 << std::scientific << std::setprecision(6);
  1939.  
  1940.     //test_cross_sections(); return 1;
  1941.     datafile = fopen("conv.dat","a");
  1942.     if (arg1 == 0) {
  1943.         if (FILE *file = fopen("picdata.bin", "r")) { fclose(file);
  1944.             printf(">> eduPIC: Warning: Data from previous calculation are detected.\n");
  1945.             printf("           To start a new simulation from the beginning, please delete all output files before running ./eduPIC 0\n");
  1946.             printf("           To continue the existing calculation, please specify the number of cycles to run, e.g. ./eduPIC 100\n");
  1947.             exit(0);
  1948.         }
  1949.         no_of_cycles = 1;
  1950.         cycle = 1;                                        // init cycle
  1951.         init(N_INIT);                                     // seed initial electrons & ions
  1952.         printf(">> eduPIC: running initializing cycle\n");
  1953.         Time = 0;
  1954.         do_one_cycle();
  1955.         // std::cout << "CHECKING SUM RATES THROUGH CYCLES ---- " << sum_rate1f[250] << "\n";
  1956.         print_excitation_densities();
  1957.         for (int i = 0; i < N_G; i++){
  1958.             double x = i*DX;
  1959.             file0 << avg_rate1f[i] << " " << avg_rate1b[i] << " " << avg_rate2f[i] << " " << avg_rate2b[i] << "\n";
  1960.             file1 << x << " " << e1_density[i] << " " << e2_density[i]
  1961.             << " " << e3_density[i] << " " << e4_density[i]<< "\n";
  1962.         }
  1963.         cycles_done = 1;
  1964.     } else {
  1965.         no_of_cycles = arg1;                              // run number of cycles specified in command line
  1966.         load_particle_data();                            // read previous configuration from file
  1967.         counter = cycles_done % N_avg;
  1968.         printf(">> eduPIC: running %d cycle(s)\n",no_of_cycles);
  1969.         for (cycle=cycles_done+1;cycle<=cycles_done+no_of_cycles;cycle++) {
  1970.             do_one_cycle();
  1971.             // std::cout << "CHECKING SUM RATES THROUGH CYCLES ---- " << sum_rate1f[250] << "|| counter: " << counter << "\n";
  1972.         }
  1973.         cycles_done += no_of_cycles;
  1974.  
  1975.     }
  1976.     fclose(datafile);
  1977.     save_particle_data();
  1978.     if (measurement_mode) {
  1979.         check_and_save_info();
  1980.     }
  1981.     printf(">> eduPIC: simulation of %d cycle(s) is completed.\n",no_of_cycles);
  1982.     file0.close();
  1983.     file1.close();
  1984. }
  1985.  
Add Comment
Please, Sign In to add comment