phystota

eduPIC_DRR_v1.3 - fixed_super_CS+more

May 26th, 2025
40
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 70.17 KB | None | 0 0
  1. //-------------------------------------------------------------------//
  2. //         eduPIC : educational 1d3v PIC/MCC simulation code         //
  3. //           version 1.0, release date: March 16, 2021               //
  4. //                       :) Share & enjoy :)                         //
  5. //-------------------------------------------------------------------//
  6. // When you use this code, you are required to acknowledge the       //
  7. // authors by citing the paper:                                      //
  8. // Z. Donko, A. Derzsi, M. Vass, B. Horvath, S. Wilczek              //
  9. // B. Hartmann, P. Hartmann:                                         //
  10. // "eduPIC: an introductory particle based  code for radio-frequency //
  11. // plasma simulation"                                                //
  12. // Plasma Sources Science and Technology, vol 30, pp. 095017 (2021)  //
  13. //-------------------------------------------------------------------//
  14. // Disclaimer: The eduPIC (educational Particle-in-Cell/Monte Carlo  //
  15. // Collisions simulation code), Copyright (C) 2021                   //
  16. // Zoltan Donko et al. is free software: you can redistribute it     //
  17. // and/or modify it under the terms of the GNU General Public License//
  18. // as published by the Free Software Foundation, version 3.          //
  19. // This program is distributed in the hope that it will be useful,   //
  20. // but WITHOUT ANY WARRANTY; without even the implied warranty of    //
  21. // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU  //
  22. // General Public License for more details at                        //
  23. // https://www.gnu.org/licenses/gpl-3.0.html.                        //
  24. //-------------------------------------------------------------------//
  25.  
  26. #include <cstdio>
  27. #include <cstdlib>
  28. #include <cstring>
  29. #include <cstdbool>
  30. #include <cmath>
  31. #include <ctime>
  32. #include <random>
  33. #include <vector>
  34. #include <string>
  35. #include <fstream>
  36. #include <sstream>
  37. #include <algorithm>    
  38. #include <stdexcept>
  39. #include <iostream>
  40.  
  41. using namespace::std;
  42.  
  43. // constants
  44.  
  45. const double     PI             = 3.141592653589793;          // mathematical constant Pi
  46. const double     TWO_PI         = 2.0 * PI;                   // two times Pi
  47. const double     E_CHARGE       = 1.60217662e-19;             // electron charge [C]
  48. const double     EV_TO_J        = E_CHARGE;                   // eV <-> Joule conversion factor
  49. const double     E_MASS         = 9.109e-31;                  // mass of electron [kg]
  50. const double     HE_MASS        = 6.67e-27;                   // mass of He atom [kg]
  51. const double     MU_HEHE        = HE_MASS / 2.0;              // reduced mass of two He atoms [kg]
  52. const double     K_BOLTZMANN    = 1.38064852e-23;             // Boltzmann's constant [J/K]
  53. const double     EPSILON0       = 8.85418781e-12;             // permittivity of free space [F/m]
  54.  
  55. // simulation parameters
  56.  
  57. const int        N_G            = 129;                        // number of grid points
  58. const int        N_T            = 400;                        // time steps within an RF period
  59.  
  60. const double     FREQUENCY      = 13.56e6;                    // driving frequency [Hz]
  61. const double     VOLTAGE        = 450.0;                      // voltage amplitude [V]
  62. const double     L              = 0.067;                      // electrode gap [m]
  63. const double     PRESSURE       = 3.99283551984;              // gas pressure [Pa] // n*k*T to match Turner's case
  64. const double     T_neutral      = 300.0;                      // background gas temperature [K] (also ion temperature)
  65. const double     T_electron     = 30000.0;                    // initial electron temperatutre [K]
  66. const double     WEIGHT         = 41875.0;                  // weight of superparticles
  67. const double     ELECTRODE_AREA = 1.6e-4;                     // (fictive) electrode area [m^2]
  68. const int        N_INIT         = 65536;                      // number of initial electrons and ions
  69.  
  70. // additional (derived) constants
  71.  
  72. const double     PERIOD         = 1.0 / FREQUENCY;                           // RF period length [s]
  73. const double     DT_E           = PERIOD / (double)(N_T);                    // electron time step [s]
  74. const int        N_SUB          = 5;                                         // ions move only in these cycles (subcycling)
  75. const int        N_avg          = 1;                          // number of cycles to average rates --- Artem
  76.       int        counter        = 0;                          // cycles since update of rates --- Artem
  77. const double     DT_I           = N_SUB * DT_E;                              // ion time step [s]
  78. const double     DX             = L / (double)(N_G - 1);                     // spatial grid division [m]
  79. const double     INV_DX         = 1.0 / DX;                                  // inverse of spatial grid size [1/m]
  80. const double     GAS_DENSITY    = PRESSURE / (K_BOLTZMANN * T_neutral);      // background gas density [1/m^3]
  81. const double     OMEGA          = TWO_PI * FREQUENCY;                        // angular frequency [rad/s]
  82.  
  83. // electron and ion cross sections
  84.  
  85. const int        N_CS           = 8;                          // total number of processes / cross sections
  86. const int        E_ELA          = 0;                          // process identifier: electron/elastic
  87. const int        E_EXC_1        = 1;                          // process identifier: electron/excitation1
  88. const int        E_EXC_2        = 2;                          // process identifier: electron/excitation1
  89. const int        E_ION          = 3;                          // process identifier: electron/ionization
  90. const int        I_ISO          = 4;                          // process identifier: ion/elastic/isotropic
  91. const int        I_BACK         = 5;                          // process identifier: ion/elastic/backscattering
  92.  
  93. const int        E_SUPER_1      = 6;                          // triplet excitation - Artem
  94. const int        E_SUPER_2      = 7;                          // singlet excitation - Artem
  95.  
  96. const double     E_EXC_TH_1     = 19.82;                      // electron impact excitation threshold [eV]
  97. const double     E_EXC_TH_2     = 20.61;                      // electron impact excitation threshold [eV]
  98. const double     E_ION_TH       = 24.587;                     // electron impact ionization threshold [eV]
  99. const int        CS_RANGES      = 1000000;                    // number of entries in cross section arrays
  100. const double     DE_CS          = 0.001;                      // energy division in cross section arrays [eV]
  101. typedef float    cross_section[CS_RANGES];                    // cross section array
  102. cross_section    sigma[N_CS];                                 // set of cross section arrays
  103. cross_section    sigma_tot_e;                                 // total macroscopic cross section of electrons
  104. cross_section    sigma_tot_i;                                 // total macroscopic cross section of ions
  105.  
  106. double nu_avg;                                                // average nu for electrons thorugh 1 cycle
  107.  
  108. // particle coordinates
  109.  
  110. const int        MAX_N_P = 10000000;                           // maximum number of particles (electrons / ions)
  111. typedef double   particle_vector[MAX_N_P];                    // array for particle properties
  112. int              N_e = 0;                                     // number of electrons
  113. int              N_i = 0;                                     // number of ions
  114. int              N_e1 = 0;                                    // number of singlet excited states
  115. int              N_e2 = 0;                                    // number of triplet excited states
  116. particle_vector  x_e, vx_e, vy_e, vz_e;                       // coordinates of electrons (one spatial, three velocity components)
  117. particle_vector  x_i, vx_i, vy_i, vz_i;                       // coordinates of ions (one spatial, three velocity components)
  118.  
  119.  
  120.  
  121.  
  122. typedef double   xvector[N_G];                                // array for quantities defined at gird points
  123. xvector          efield,pot;                                  // electric field and potential
  124. xvector          e_density,i_density;                         // electron and ion densities
  125. xvector          cumul_e_density,cumul_i_density;             // cumulative densities
  126.  
  127. //excited states handling
  128. xvector e1_density;
  129. xvector e2_density;
  130. xvector sum_electron_density; xvector avg_electron_density;
  131.  
  132. xvector sum_rate1f = {0.0}; xvector sum_rate1b = {0.0}; xvector sum_rate2f = {0.0}; xvector sum_rate2b = {0.0};
  133.  
  134. xvector avg_rate1f = {0.0}; xvector avg_rate1b = {0.0}; xvector avg_rate2f = {0.0}; xvector avg_rate2b = {0.0};
  135.  
  136.                                        // array for Triplet excitation rates!!! Artem
  137. xvector S1 = {0.0};
  138. xvector S2 = {0.0};                                           // source terms for DRR module -- Artem
  139.  
  140. typedef unsigned long long int Ullong;                        // compact name for 64 bit unsigned integer
  141. Ullong       N_e_abs_pow  = 0;                                // counter for electrons absorbed at the powered electrode
  142. Ullong       N_e_abs_gnd  = 0;                                // counter for electrons absorbed at the grounded electrode
  143. Ullong       N_i_abs_pow  = 0;                                // counter for ions absorbed at the powered electrode
  144. Ullong       N_i_abs_gnd  = 0;                                // counter for ions absorbed at the grounded electrode
  145.  
  146. // electron energy probability function
  147.  
  148. const int    N_EEPF  = 2000;                                 // number of energy bins in Electron Energy Probability Function (EEPF)
  149. const double DE_EEPF = 0.05;                                 // resolution of EEPF [eV]
  150. typedef double eepf_vector[N_EEPF];                          // array for EEPF
  151. eepf_vector eepf     = {0.0};                                // time integrated EEPF in the center of the plasma
  152.  
  153. // ion flux-energy distributions
  154.  
  155. const int    N_IFED   = 200;                                 // number of energy bins in Ion Flux-Energy Distributions (IFEDs)
  156. const double DE_IFED  = 1.0;                                 // resolution of IFEDs [eV]
  157. typedef int  ifed_vector[N_IFED];                            // array for IFEDs
  158. ifed_vector  ifed_pow = {0};                                 // IFED at the powered electrode
  159. ifed_vector  ifed_gnd = {0};                                 // IFED at the grounded electrode
  160. double       mean_i_energy_pow;                              // mean ion energy at the powered electrode
  161. double       mean_i_energy_gnd;                              // mean ion energy at the grounded electrode
  162.  
  163. // spatio-temporal (XT) distributions
  164.  
  165. const int N_BIN                     = 20;                    // number of time steps binned for the XT distributions
  166. const int N_XT                      = N_T / N_BIN;           // number of spatial bins for the XT distributions
  167. typedef double xt_distr[N_G][N_XT];                          // array for XT distributions (decimal numbers)
  168. xt_distr pot_xt                     = {0.0};                 // XT distribution of the potential
  169. xt_distr efield_xt                  = {0.0};                 // XT distribution of the electric field
  170. xt_distr ne_xt                      = {0.0};                 // XT distribution of the electron density
  171. xt_distr ni_xt                      = {0.0};                 // XT distribution of the ion density
  172. xt_distr ue_xt                      = {0.0};                 // XT distribution of the mean electron velocity
  173. xt_distr ui_xt                      = {0.0};                 // XT distribution of the mean ion velocity
  174. xt_distr je_xt                      = {0.0};                 // XT distribution of the electron current density
  175. xt_distr ji_xt                      = {0.0};                 // XT distribution of the ion current density
  176. xt_distr powere_xt                  = {0.0};                 // XT distribution of the electron powering (power absorption) rate
  177. xt_distr poweri_xt                  = {0.0};                 // XT distribution of the ion powering (power absorption) rate
  178. xt_distr meanee_xt                  = {0.0};                 // XT distribution of the mean electron energy
  179. xt_distr meanei_xt                  = {0.0};                 // XT distribution of the mean ion energy
  180. xt_distr counter_e_xt               = {0.0};                 // XT counter for electron properties
  181. xt_distr counter_i_xt               = {0.0};                 // XT counter for ion properties
  182. xt_distr ioniz_rate_xt              = {0.0};                 // XT distribution of the ionisation rate
  183.  
  184. xt_distr e1_xt                      = {0.0};                 // XT distribution of triplet excited states density - Artem
  185. xt_distr e2_xt                      = {0.0};                 // XT distribution of singlet excited states density - Artem
  186.  
  187. double   mean_energy_accu_center    = 0;                     // mean electron energy accumulator in the center of the gap
  188. Ullong   mean_energy_counter_center = 0;                     // mean electron energy counter in the center of the gap
  189. Ullong   N_e_coll                   = 0;                     // counter for electron collisions
  190. Ullong   N_i_coll                   = 0;                     // counter for ion collisions
  191. double   Time;                                               // total simulated time (from the beginning of the simulation)
  192. int      cycle,no_of_cycles,cycles_done;                     // current cycle and total cycles in the run, cycles completed
  193. int      arg1;                                               // used for reading command line arguments
  194. char     st0[80];                                            // used for reading command line arguments
  195. FILE     *datafile;                                          // used for saving data
  196. bool     measurement_mode;                                   // flag that controls measurements and data saving
  197.  
  198. //---------------------------------------------------------------------------//
  199. // C++ Mersenne Twister 19937 generator                                      //
  200. // R01(MTgen) will genarate uniform distribution over [0,1) interval         //
  201. // RMB(MTgen) will generate Maxwell-Boltzmann distribution (of gas atoms)    //
  202. //---------------------------------------------------------------------------//
  203.  
  204. std::random_device rd{};
  205. std::mt19937 MTgen(rd());
  206. std::uniform_real_distribution<> R01(0.0, 1.0);
  207. std::normal_distribution<> RMB_n(0.0,sqrt(K_BOLTZMANN * T_neutral / HE_MASS));
  208. std::normal_distribution<> RMB_e(0.0,sqrt(K_BOLTZMANN * T_electron / E_MASS));
  209.  
  210. //----------------------------------------------------------------------------//
  211. //  electron cross sections: A V Phelps & Z Lj Petrovic, PSST 8 R21 (1999)    //
  212. //----------------------------------------------------------------------------//
  213.  
  214. class CSInterpolator {
  215. public:
  216.   // load "filename" which must have two whitespace‐separated columns:
  217.   //    energy (eV)    cross_section (in m^2 or cm^2 as you prefer)
  218.   CSInterpolator(const std::string &filename) {
  219.     std::ifstream in(filename);
  220.     if (!in) throw std::runtime_error("CSInterpolator: cannot open " + filename);
  221.     double E, sigma;
  222.     std::string line;
  223.     while (std::getline(in, line)) {
  224.       std::istringstream iss(line);
  225.       if (iss >> E >> sigma) {
  226.         E_pts_.push_back(E);
  227.         sigma_pts_.push_back(sigma);
  228.       }
  229.     }
  230.     if (E_pts_.size()<2)
  231.       throw std::runtime_error("CSInterpolator: need at least two data points in " + filename);
  232.   }
  233.  
  234.   // return σ(E) by simple linear interpolation (clamped to end‐points)
  235.   double operator()(double E) const {
  236.     auto it = std::lower_bound(E_pts_.begin(), E_pts_.end(), E);
  237.     if (it == E_pts_.begin()) {
  238.         std::cerr << "Warning: E="<<E<<" below data range, clamping to "<< 0.0 <<"\n";
  239.         return 0.0;
  240.     }  
  241.     if (it == E_pts_.end()){
  242.         std::cerr << "Warning: E="<<E<<" above data range, clamping to "<< sigma_pts_.back() <<"\n";
  243.         return sigma_pts_.back();
  244.     }    
  245.     size_t idx = (it - E_pts_.begin());
  246.     double E1 = E_pts_[idx-1], E2 = E_pts_[idx];
  247.     double s1 = sigma_pts_[idx-1], s2 = sigma_pts_[idx];
  248.     // linear interp
  249.     return s1 + (s2 - s1) * (E - E1)/(E2 - E1);
  250.   }
  251.  
  252. private:
  253.   std::vector<double> E_pts_, sigma_pts_;
  254. };
  255.  
  256. void set_electron_cross_sections_ar(void){
  257.     int    i;
  258.     double en,qmel,qexc_1,qexc_2,qion;
  259.    
  260.     printf(">> eduPIC: Setting e- / He cross sections\n");
  261.  
  262.     // load your four datafiles (make sure these names match your files!)
  263.     CSInterpolator cs_ela  ("CS/He_electron_elastic.dat");   // two‐col: E σ_ela
  264.     CSInterpolator cs_exc1 ("CS/He_electron_exc1.dat");      // two‐col: E σ_exc1
  265.     CSInterpolator cs_exc2 ("CS/He_electron_exc2.dat");      // two‐col: E σ_exc2
  266.     CSInterpolator cs_ion  ("CS/He_electron_ionization.dat");// two‐col: E σ_ion
  267.  
  268.     for(int i=0; i<CS_RANGES; i++){
  269.         // your energy grid
  270.         double en = (i==0 ? DE_CS : DE_CS * i);
  271.  
  272.         // interpolate
  273.         sigma[E_ELA][i]   = cs_ela(en);
  274.         sigma[E_EXC_1][i] = cs_exc1(en);
  275.         sigma[E_EXC_2][i] = cs_exc2(en);
  276.         sigma[E_ION][i]   = cs_ion(en);
  277.  
  278.         // Superelastic for triplet (E_SUPER_1)
  279.         double en_plus_1 = en + E_EXC_TH_1;
  280.         int idx1 = en_plus_1 / DE_CS;
  281.         if (idx1 < CS_RANGES && en > 0) {
  282.             sigma[E_SUPER_1][i] = (1.0 / 3.0) * (en_plus_1 / en) * sigma[E_EXC_1][idx1];
  283.         } else {
  284.             sigma[E_SUPER_1][i] = 0.0;
  285.         }        
  286.  
  287.         // Superelastic for singlet (E_SUPER_2)
  288.         double en_plus_2 = en + E_EXC_TH_2;
  289.         int idx2 = en_plus_2 / DE_CS;
  290.         if (idx2 < CS_RANGES && en > 0) {
  291.             sigma[E_SUPER_2][i] = (1.0 / 1.0) * (en_plus_2 / en) * sigma[E_EXC_2][idx2];
  292.         } else {
  293.             sigma[E_SUPER_2][i] = 0.0;
  294.         }        
  295.     }
  296.  
  297.     for(int i=0; i<CS_RANGES; i++){
  298.         // your energy grid
  299.         double en = (i==0 ? DE_CS : DE_CS * i);
  300.  
  301.         // Superelastic for triplet (E_SUPER_1)
  302.         double en_plus_1 = en + E_EXC_TH_1;
  303.         int idx1 = en_plus_1 / DE_CS;
  304.         if (idx1 < CS_RANGES && en > 0) {
  305.             sigma[E_SUPER_1][i] = (1.0 / 3.0) * (en_plus_1 / en) * sigma[E_EXC_1][idx1];
  306.         } else {
  307.             sigma[E_SUPER_1][i] = 0.0;
  308.         }        
  309.  
  310.         // Superelastic for singlet (E_SUPER_2)
  311.         double en_plus_2 = en + E_EXC_TH_2;
  312.         int idx2 = en_plus_2 / DE_CS;
  313.         if (idx2 < CS_RANGES && en > 0) {
  314.             sigma[E_SUPER_2][i] = (1.0 / 1.0) * (en_plus_2 / en) * sigma[E_EXC_2][idx2];
  315.         } else {
  316.             sigma[E_SUPER_2][i] = 0.0;
  317.         }        
  318.     }    
  319.  
  320. }
  321.  
  322. //------------------------------------------------------------------------------//
  323. //  ion cross sections: A. V. Phelps, J. Appl. Phys. 76, 747 (1994)             //
  324. //------------------------------------------------------------------------------//
  325.  
  326. void set_ion_cross_sections_ar(void){
  327.     int    i;
  328.     double e_com,e_lab,qmom,qback,qiso;
  329.    
  330.     printf(">> eduPIC: Setting He+ / He cross sections\n");
  331.     for(i=0; i<CS_RANGES; i++){
  332.         if (i == 0) {e_com = DE_CS;} else {e_com = DE_CS * i;}             // ion energy in the center of mass frame of reference
  333.         e_lab = 2.0 * e_com;                                               // ion energy in the laboratory frame of reference
  334.         qiso  = 7.63 *pow(10,-20) * pow(e_com, -0.5);
  335.         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 );
  336.         sigma[I_ISO][i]  = qiso;             // cross section for He+ / He isotropic part of elastic scattering
  337.         sigma[I_BACK][i] = qback;            // cross section for He+ / He backward elastic scattering
  338.     }
  339. }
  340.  
  341. //----------------------------------------------------------------------//
  342. //  calculation of total cross sections for electrons and ions          //
  343. //----------------------------------------------------------------------//
  344.  
  345. void calc_total_cross_sections(void){
  346.     int i;
  347.    
  348.     for(i=0; i<CS_RANGES; i++){
  349.         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
  350.         sigma_tot_i[i] = (sigma[I_ISO][i] + sigma[I_BACK][i]) * GAS_DENSITY;                    // total macroscopic cross section of ions
  351.     }
  352. }
  353.  
  354. //----------------------------------------------------------------------//
  355. //  test of cross sections for electrons and ions                       //
  356. //----------------------------------------------------------------------//
  357.  
  358. void test_cross_sections(void){
  359.     FILE  * f;
  360.     int   i,j;
  361.    
  362.     f = fopen("cross_sections.dat","w");        // cross sections saved in data file: cross_sections.dat
  363.     for(i=0; i<CS_RANGES; i++){
  364.         fprintf(f,"%12.4f ",i*DE_CS);
  365.         for(j=0; j<N_CS; j++) fprintf(f,"%14e ",sigma[j][i]);
  366.         fprintf(f,"\n");
  367.     }
  368.     fclose(f);
  369. }
  370.  
  371. //---------------------------------------------------------------------//
  372. // find upper limit of collision frequencies                           //
  373. //---------------------------------------------------------------------//
  374.  
  375. double max_electron_coll_freq (void){
  376.     int i;
  377.     double e,v,nu,nu_max;
  378.     nu_max = 0;
  379.     for(i=0; i<CS_RANGES; i++){
  380.         e  = i * DE_CS;
  381.         v  = sqrt(2.0 * e * EV_TO_J / E_MASS);
  382.         nu = v * sigma_tot_e[i];
  383.         if (nu > nu_max) {nu_max = nu;}
  384.     }
  385.     return nu_max;
  386. }
  387.  
  388. double max_ion_coll_freq (void){
  389.     int i;
  390.     double e,g,nu,nu_max;
  391.     nu_max = 0;
  392.     for(i=0; i<CS_RANGES; i++){
  393.         e  = i * DE_CS;
  394.         g  = sqrt(2.0 * e * EV_TO_J / MU_HEHE);
  395.         nu = g * sigma_tot_i[i];
  396.         if (nu > nu_max) nu_max = nu;
  397.     }
  398.     return nu_max;
  399. }
  400.  
  401. //----------------------------------------------------------------------//
  402. // initialization of the simulation by placing a given number of        //
  403. // electrons and ions at random positions between the electrodes        //
  404. //----------------------------------------------------------------------//
  405.  
  406. // initialization of excited states distribtuion, assuming Maxwellian-Boltzmann balance -- Artem
  407. std::pair<double, double> init_excited_distr() {
  408.     double part_ground = 1.0*exp(-0.0/T_neutral); // partition function for ground state
  409.     double part_triplet = 3.0*exp(-E_EXC_TH_1*EV_TO_J/(K_BOLTZMANN*T_neutral)); // partition function for triplet excited state
  410.     double part_singlet = 1.0*exp(-E_EXC_TH_2*EV_TO_J/(K_BOLTZMANN*T_neutral)); // partition function for singlet excited state
  411.     double part_func_total = part_ground + part_triplet + part_singlet; // total partition function
  412.     double n_triplet = ((part_triplet/part_func_total)*GAS_DENSITY); // denisty population of tripet state
  413.     double n_singlet = ((part_singlet/part_func_total)*GAS_DENSITY); // density population of singlet state
  414.  
  415.     return {n_triplet, n_singlet};
  416. }
  417.  
  418. void print_excitation_densities(void) {
  419.     double total_e1 = 0.0, total_e2 = 0.0;
  420.     const double cell_volume = ELECTRODE_AREA * DX;
  421.    
  422.     // Sum densities across all grid cells
  423.     for (int p = 0; p < N_G; p++) {
  424.         total_e1 += e1_density[p];  // triplet state density
  425.         total_e2 += e2_density[p];  // singlet state density
  426.     }
  427.     printf("Triplet SP = %8.2e | Singlet SP = %8.2e\n", total_e1, total_e2);
  428. }
  429.  
  430. void init(int nseed){
  431.     int i;
  432.    
  433.     for (i=0; i<nseed; i++){
  434.         x_e[i]  = L * R01(MTgen);                                                   // initial random position of the electron
  435.         vx_e[i] = RMB_e(MTgen); vy_e[i] = RMB_e(MTgen); vz_e[i] = RMB_e(MTgen);     // initial velocity components of the electron
  436.         x_i[i]  = L * R01(MTgen);                                                   // initial random position of the ion
  437.         vx_i[i] = RMB_n(MTgen); vy_i[i] = RMB_n(MTgen); vz_i[i] = RMB_n(MTgen);     // initial velocity components of the ion
  438.     }
  439.     N_e = nseed;    // initial number of electrons
  440.     N_i = nseed;    // initial number of ions
  441.  
  442.     double n_min = 1e-15 * GAS_DENSITY;
  443.     auto exc = init_excited_distr();
  444.     for (int p = 0; p < N_G; p++) {
  445.         e1_density[p] = exc.first;
  446.         e2_density[p] = exc.second;
  447.     }
  448. }
  449.  
  450. //----------------------------------------------------------------------//
  451. // e / He collision  (cold gas approximation)                           //
  452. //----------------------------------------------------------------------//
  453.  
  454. void collision_electron (double xe, double *vxe, double *vye, double *vze, int eindex){
  455.     const double F1 = E_MASS  / (E_MASS + HE_MASS);
  456.     const double F2 = HE_MASS / (E_MASS + HE_MASS);
  457.     double t0,t1,t2,t3,t4,t5,rnd;
  458.     double g,g2,gx,gy,gz,wx,wy,wz,theta,phi;
  459.     double chi,eta,chi2,eta2,sc,cc,se,ce,st,ct,sp,cp,energy,e_sc,e_ej;
  460.  
  461.     // - Artem
  462.     // Determine cell p where the electron is
  463.     double c0 = xe * INV_DX;
  464.     int p = std::max(0, std::min(N_G - 1, static_cast<int>(c0)));
  465.  
  466.     // Local densities -- Artem
  467.     double e1_dens = e1_density[p];
  468.     double e2_dens = e2_density[p];
  469.     double ground_dens = GAS_DENSITY - e1_dens - e2_dens;    
  470.  
  471.    
  472.     // calculate relative velocity before collision & velocity of the centre of mass
  473.    
  474.     gx = (*vxe);
  475.     gy = (*vye);
  476.     gz = (*vze);
  477.     g  = sqrt(gx * gx + gy * gy + gz * gz);
  478.     wx = F1 * (*vxe);
  479.     wy = F1 * (*vye);
  480.     wz = F1 * (*vze);
  481.    
  482.     // find Euler angles
  483.    
  484.     if (gx == 0) {theta = 0.5 * PI;}
  485.     else {theta = atan2(sqrt(gy * gy + gz * gz),gx);}
  486.     if (gy == 0) {
  487.         if (gz > 0){phi = 0.5 * PI;} else {phi = - 0.5 * PI;}
  488.     } else {phi = atan2(gz, gy);}
  489.     st  = sin(theta);
  490.     ct  = cos(theta);
  491.     sp  = sin(phi);
  492.     cp  = cos(phi);
  493.    
  494.     // choose the type of collision based on the cross sections
  495.     // take into account energy loss in inelastic collisions
  496.     // generate scattering and azimuth angles
  497.     // in case of ionization handle the 'new' electron
  498.    
  499.     t0   =      sigma[E_ELA][eindex] * ground_dens;
  500.     t1   = t0 + sigma[E_EXC_1][eindex] * ground_dens;
  501.     t2   = t1 + sigma[E_EXC_2][eindex] * ground_dens;
  502.     t3   = t2 + sigma[E_ION][eindex] * ground_dens;
  503.     t4   = t3 + sigma[E_SUPER_1][eindex]  * e1_dens;   // account for superelastic triplet - Artem
  504.     t5   = t4 + sigma[E_SUPER_2][eindex] * e2_dens;   // account for superelastic singlet- Artem
  505.  
  506.     rnd  = R01(MTgen);
  507.     if (rnd < (t0/t5)){                              // elastic scattering
  508.         chi = acos(1.0 - 2.0 * R01(MTgen));          // isotropic scattering
  509.         eta = TWO_PI * R01(MTgen);                   // azimuthal angle
  510.  
  511.     } else if (rnd < (t1/t5)){                       // excitation 1 (triplet)
  512.         energy = 0.5 * E_MASS * g * g;               // electron energy
  513.         energy = fabs(energy - E_EXC_TH_1 * EV_TO_J);  // subtract energy loss for excitation
  514.         g   = sqrt(2.0 * energy / E_MASS);           // relative velocity after energy loss
  515.         chi = acos(1.0 - 2.0 * R01(MTgen));          // isotropic scattering
  516.         eta = TWO_PI * R01(MTgen);                   // azimuthal angle
  517.  
  518.         //add new triplet excited He atom density to this grid point, sample velocities from Maxwellian distribution - Artem
  519.         //no need if we calculate rates and densities outside
  520. //        e1_density[p] += WEIGHT / (ELECTRODE_AREA * DX);
  521.  
  522.     } else if (rnd < (t2/t5)){                       // excitation 2 (singlet)
  523.         energy = 0.5 * E_MASS * g * g;               // electron energy
  524.         energy = fabs(energy - E_EXC_TH_2 * EV_TO_J);  // subtract energy loss for excitation
  525.         g   = sqrt(2.0 * energy / E_MASS);           // relative velocity after energy loss
  526.         chi = acos(1.0 - 2.0 * R01(MTgen));          // isotropic scattering
  527.         eta = TWO_PI * R01(MTgen);                   // azimuthal angle    
  528.  
  529.         //add new singet excited He atom, sample velocities from Maxwellian distribution - Artem
  530.         //no need if we calculate rates and densities outside
  531. //        e2_density[p] += WEIGHT / (ELECTRODE_AREA * DX);
  532.  
  533.     } else if (rnd < (t3/t5)) {                                         // ionization
  534.         energy = 0.5 * E_MASS * g * g;               // electron energy
  535.         energy = fabs(energy - E_ION_TH * EV_TO_J);  // subtract energy loss of ionization
  536.         // e_ej  = 10.0 * tan(R01(MTgen) * atan(energy/EV_TO_J / 20.0)) * EV_TO_J; // energy of the ejected electron
  537.         // e_sc = fabs(energy - e_ej);                  // energy of scattered electron after the collision
  538.         e_ej = 0.5*energy;                          // energy of the ejected electron
  539.         e_sc = fabs(energy - e_ej);                  // energy of scattered electron after the collision        
  540.         g    = sqrt(2.0 * e_sc / E_MASS);            // relative velocity of scattered electron
  541.         g2   = sqrt(2.0 * e_ej / E_MASS);            // relative velocity of ejected electron
  542.         // chi  = acos(sqrt(e_sc / energy));            // scattering angle for scattered electron
  543.         // chi2 = acos(sqrt(e_ej / energy));            // scattering angle for ejected electrons
  544.         chi = acos(1.0 - 2.0 * R01(MTgen));          // isotropic scattering for scattered electron (as in Turner's case)
  545.         chi2 = acos(1.0 - 2.0 * R01(MTgen));          // isotropic scattering for ejected electrons (as in Turner's case)
  546.         eta  = TWO_PI * R01(MTgen);                  // azimuthal angle for scattered electron
  547.         eta2 = eta + PI;                             // azimuthal angle for ejected electron
  548.         sc  = sin(chi2);
  549.         cc  = cos(chi2);
  550.         se  = sin(eta2);
  551.         ce  = cos(eta2);
  552.         gx  = g2 * (ct * cc - st * sc * ce);
  553.         gy  = g2 * (st * cp * cc + ct * cp * sc * ce - sp * sc * se);
  554.         gz  = g2 * (st * sp * cc + ct * sp * sc * ce + cp * sc * se);
  555.         x_e[N_e]  = xe;                              // add new electron
  556.         vx_e[N_e] = wx + F2 * gx;
  557.         vy_e[N_e] = wy + F2 * gy;
  558.         vz_e[N_e] = wz + F2 * gz;
  559.         N_e++;
  560.         x_i[N_i]  = xe;                              // add new ion
  561.         vx_i[N_i] = RMB_n(MTgen);                      // velocity is sampled from background thermal distribution
  562.         vy_i[N_i] = RMB_n(MTgen);
  563.         vz_i[N_i] = RMB_n(MTgen);
  564.         N_i++;
  565.     }
  566.      else if (rnd < (t4/t5)) {                      // accounting for superelastic collisions - triplet - Artem
  567.         energy = 0.5 * E_MASS * g * g;               // electron energy
  568.         energy = fabs(energy + E_EXC_TH_1 * EV_TO_J);  // add energy for deexcitation
  569.         g   = sqrt(2.0 * energy / E_MASS);           // relative velocity after energy loss
  570.         chi = acos(1.0 - 2.0 * R01(MTgen));          // isotropic scattering
  571.         eta = TWO_PI * R01(MTgen);                   // azimuthal angle    
  572.  
  573.         //substract  excited He atom density from the grid point - Artem
  574.         //no need if we calculate rates and densities outside
  575. //        e1_density[p] = std::max(e1_density[p] - WEIGHT / (ELECTRODE_AREA * DX), 0.0);
  576.    
  577.     } else {                                         // account for supeelastic collision - singlet - Artem
  578.         energy = 0.5 * E_MASS * g * g;               // electron energy
  579.         energy = fabs(energy + E_EXC_TH_2 * EV_TO_J);  // add energy for deexcitation
  580.         g   = sqrt(2.0 * energy / E_MASS);           // relative velocity after energy loss
  581.         chi = acos(1.0 - 2.0 * R01(MTgen));          // isotropic scattering
  582.         eta = TWO_PI * R01(MTgen);                   // azimuthal angle    
  583.  
  584.         //substract  excited He atom density from the grid point - Artem
  585.         //no need if we calculate rates and densities outside
  586.  //       e2_density[p] = std::max(e2_density[p] - WEIGHT / (ELECTRODE_AREA * DX), 0.0);
  587.     }
  588.  
  589.    
  590.     // scatter the primary electron
  591.    
  592.     sc = sin(chi);
  593.     cc = cos(chi);
  594.     se = sin(eta);
  595.     ce = cos(eta);
  596.    
  597.     // compute new relative velocity:
  598.    
  599.     gx = g * (ct * cc - st * sc * ce);
  600.     gy = g * (st * cp * cc + ct * cp * sc * ce - sp * sc * se);
  601.     gz = g * (st * sp * cc + ct * sp * sc * ce + cp * sc * se);
  602.    
  603.     // post-collision velocity of the colliding electron
  604.    
  605.     (*vxe) = wx + F2 * gx;
  606.     (*vye) = wy + F2 * gy;
  607.     (*vze) = wz + F2 * gz;
  608. }
  609.  
  610. //----------------------------------------------------------------------//
  611. // He+ / He collision                                                   //
  612. //----------------------------------------------------------------------//
  613.  
  614. void collision_ion (double *vx_1, double *vy_1, double *vz_1,
  615.                     double *vx_2, double *vy_2, double *vz_2, int e_index){
  616.     double   g,gx,gy,gz,wx,wy,wz,rnd;
  617.     double   theta,phi,chi,eta,st,ct,sp,cp,sc,cc,se,ce,t1,t2;
  618.    
  619.     // calculate relative velocity before collision
  620.     // random Maxwellian target atom already selected (vx_2,vy_2,vz_2 velocity components of target atom come with the call)
  621.    
  622.     gx = (*vx_1)-(*vx_2);
  623.     gy = (*vy_1)-(*vy_2);
  624.     gz = (*vz_1)-(*vz_2);
  625.     g  = sqrt(gx * gx + gy * gy + gz * gz);
  626.     wx = 0.5 * ((*vx_1) + (*vx_2));
  627.     wy = 0.5 * ((*vy_1) + (*vy_2));
  628.     wz = 0.5 * ((*vz_1) + (*vz_2));
  629.    
  630.     // find Euler angles
  631.    
  632.     if (gx == 0) {theta = 0.5 * PI;} else {theta = atan2(sqrt(gy * gy + gz * gz),gx);}
  633.     if (gy == 0) {
  634.         if (gz > 0){phi = 0.5 * PI;} else {phi = - 0.5 * PI;}
  635.     } else {phi = atan2(gz, gy);}
  636.    
  637.     // determine the type of collision based on cross sections and generate scattering angle
  638.    
  639.     t1  =      sigma[I_ISO][e_index];
  640.     t2  = t1 + sigma[I_BACK][e_index];
  641.     rnd = R01(MTgen);
  642.     if  (rnd < (t1 /t2)){                        // isotropic scattering
  643.         chi = acos(1.0 - 2.0 * R01(MTgen));      // scattering angle
  644.     } else {                                     // backward scattering
  645.         chi = PI;                                // scattering angle
  646.     }
  647.     eta = TWO_PI * R01(MTgen);                   // azimuthal angle
  648.     sc  = sin(chi);
  649.     cc  = cos(chi);
  650.     se  = sin(eta);
  651.     ce  = cos(eta);
  652.     st  = sin(theta);
  653.     ct  = cos(theta);
  654.     sp  = sin(phi);
  655.     cp  = cos(phi);
  656.    
  657.     // compute new relative velocity
  658.    
  659.     gx = g * (ct * cc - st * sc * ce);
  660.     gy = g * (st * cp * cc + ct * cp * sc * ce - sp * sc * se);
  661.     gz = g * (st * sp * cc + ct * sp * sc * ce + cp * sc * se);
  662.    
  663.     // post-collision velocity of the ion
  664.    
  665.     (*vx_1) = wx + 0.5 * gx;
  666.     (*vy_1) = wy + 0.5 * gy;
  667.     (*vz_1) = wz + 0.5 * gz;
  668. }
  669.  
  670. //-----------------------------------------------------------------//
  671. // solve Poisson equation (Thomas algorithm)                       //
  672. //-----------------------------------------------------------------//
  673.  
  674. void solve_Poisson (xvector rho1, double tt){
  675.     const double A =  1.0;
  676.     const double B = -2.0;
  677.     const double C =  1.0;
  678.     const double S = 1.0 / (2.0 * DX);
  679.     const double ALPHA = -DX * DX / EPSILON0;
  680.     xvector      g, w, f;
  681.     int          i;
  682.    
  683.     // apply potential to the electrodes - boundary conditions
  684.    
  685.     pot[0]     = VOLTAGE * cos(OMEGA * tt);         // potential at the powered electrode
  686.     pot[N_G-1] = 0.0;                               // potential at the grounded electrode
  687.    
  688.     // solve Poisson equation
  689.    
  690.     for(i=1; i<=N_G-2; i++) f[i] = ALPHA * rho1[i];
  691.     f[1] -= pot[0];
  692.     f[N_G-2] -= pot[N_G-1];
  693.     w[1] = C/B;
  694.     g[1] = f[1]/B;
  695.     for(i=2; i<=N_G-2; i++){
  696.         w[i] = C / (B - A * w[i-1]);
  697.         g[i] = (f[i] - A * g[i-1]) / (B - A * w[i-1]);
  698.     }
  699.     pot[N_G-2] = g[N_G-2];
  700.     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
  701.    
  702.     // compute electric field
  703.    
  704.     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
  705.     efield[0]     = (pot[0]     - pot[1])     * INV_DX - rho1[0]     * DX / (2.0 * EPSILON0);   // powered electrode
  706.     efield[N_G-1] = (pot[N_G-2] - pot[N_G-1]) * INV_DX + rho1[N_G-1] * DX / (2.0 * EPSILON0);   // grounded electrode
  707. }
  708.  
  709. //---------------------------------------------------------------------//
  710. // simulation of one radiofrequency cycle                              //
  711. //---------------------------------------------------------------------//
  712.  
  713. void accumulate_rates() {
  714.     double v_sqr, velocity, energy, c0_temp;
  715.     int energy_index, p_temp;
  716.     const double Volume = (ELECTRODE_AREA * DX);
  717.  
  718.     for (int k=0; k<N_e; k++){                              
  719.         v_sqr = vx_e[k] * vx_e[k] + vy_e[k] * vy_e[k] + vz_e[k] * vz_e[k];
  720.         velocity = sqrt(v_sqr);
  721.         energy   = 0.5 * E_MASS * v_sqr / EV_TO_J;
  722.         energy_index = min( int(energy / DE_CS + 0.5), CS_RANGES-1);
  723.  
  724.         c0_temp = x_e[k] * INV_DX;
  725.         p_temp = std::max(0, std::min(N_G - 1, static_cast<int>(c0_temp)));
  726.    
  727.         sum_electron_density[p_temp] += WEIGHT / Volume;
  728.  
  729.         sum_rate1f[p_temp] += sigma[E_EXC_1][energy_index] * velocity * WEIGHT;
  730.         sum_rate2f[p_temp] += sigma[E_EXC_2][energy_index] * velocity * WEIGHT;
  731.  
  732.         sum_rate1b[p_temp] += sigma[E_SUPER_1][energy_index] * velocity * WEIGHT;
  733.         sum_rate2b[p_temp] += sigma[E_SUPER_2][energy_index] * velocity * WEIGHT;
  734.     }    
  735. }
  736.  
  737. // averaging the rates each RF cycle
  738. void average_rates() {
  739.     const double inv_NT = 1.0 / (N_avg * N_T); // averaging through N_avg RF periods each contains N_T cycles
  740.     const double inv_Volume = 1.0/(ELECTRODE_AREA * DX);
  741.     for (int p = 0; p < N_G; p++) {
  742.         avg_rate1f[p] = sum_rate1f[p] * inv_NT * inv_Volume;  
  743.         avg_rate1b[p] = sum_rate1b[p] * inv_NT * inv_Volume;
  744.         avg_rate2f[p] = sum_rate2f[p] * inv_NT * inv_Volume;
  745.         avg_rate2b[p] = sum_rate2b[p] * inv_NT * inv_Volume;
  746.         avg_electron_density[p] = sum_electron_density[p] * inv_NT;
  747.     }    
  748.  
  749.     // for (int i = 0; i < N_G; i++)
  750.     // std::cout << avg_rate2b[i] << "\n";
  751. }
  752.  
  753. void solve_steady_state(int p) {
  754.     double dn0, dn1, dn2;
  755.     double gas_dens_local;
  756.     double tol = 1.0E-5;
  757.  
  758.  
  759.     const int max_iter = static_cast<int>(N_avg * PERIOD / DT_E);
  760.  
  761.     gas_dens_local = GAS_DENSITY - e1_density[p] - e2_density[p];
  762.  
  763.     // Store original values for rollback
  764.     const double original_e1 = e1_density[p];
  765.     const double original_e2 = e2_density[p];
  766.     double original_gas = GAS_DENSITY - original_e1 - original_e2;    
  767.  
  768.     for (int j = 0; j < max_iter; j++) {
  769.    
  770.         dn0 = DT_E*(-(avg_rate1f[p]+avg_rate2f[p])*gas_dens_local + avg_rate1b[p]*e1_density[p] + avg_rate2b[p]*e2_density[p]);
  771.         dn1 = DT_E*(avg_rate1f[p]*gas_dens_local - avg_rate1b[p]*e1_density[p]);
  772.         dn2 = DT_E*(avg_rate2f[p]*gas_dens_local - avg_rate2b[p]*e2_density[p]);
  773.  
  774.         e1_density[p] += dn1;
  775.         e1_density[p] = std::max(e1_density[p], 0.0);
  776.        
  777.         e2_density[p] += dn2;
  778.         e2_density[p] = std::max(e2_density[p], 0.0);
  779.  
  780.         if (p == N_G - 1) {
  781.             e1_density[p] = e1_density[p-1];
  782.             e2_density[p] = e2_density[p-1];
  783.         }
  784.  
  785.         gas_dens_local += dn0;
  786.         gas_dens_local = std::max(gas_dens_local, 0.0);  // Prevent negative ground density
  787.  
  788.     }  
  789.  
  790. }
  791.  
  792. void update_excited_dens() {
  793.     for (int p = 0; p < N_G; p++){
  794.         solve_steady_state(p);
  795.     }
  796. }
  797.  
  798.  
  799. void do_one_cycle (void){
  800.     const double DV       = ELECTRODE_AREA * DX;
  801.     const double FACTOR_W = WEIGHT / DV;
  802.     const double FACTOR_E = DT_E / E_MASS * E_CHARGE;
  803.     const double FACTOR_I = DT_I / HE_MASS * E_CHARGE;
  804.     const double MIN_X    = 0.45 * L;                       // min. position for EEPF collection
  805.     const double MAX_X    = 0.55 * L;                       // max. position for EEPF collection
  806.     int      k, t, p, energy_index;
  807.     double   g, g_sqr, gx, gy, gz, vx_a, vy_a, vz_a, e_x, energy, nu, p_coll, v_sqr, velocity;
  808.     double   mean_v, c0, c1, c2, rate;
  809.     bool     out;
  810.     xvector  rho;
  811.     int      t_index;
  812.  
  813.     nu_avg = 0.0;
  814.    
  815.     for (t=0; t<N_T; t++){          // the RF period is divided into N_T equal time intervals (time step DT_E)
  816.         Time += DT_E;               // update of the total simulated time
  817.         t_index = t / N_BIN;        // index for XT distributions
  818.        
  819.         // step 1: compute densities at grid points
  820.        
  821.         for(p=0; p<N_G; p++) e_density[p] = 0;                             // electron density - computed in every time step
  822.         for(k=0; k<N_e; k++){
  823.  
  824.             if      (p < 0)        p = 0;
  825.             else if (p > N_G - 2)  p = N_G - 2;
  826.             c0 = x_e[k] * INV_DX;
  827.             p  = int(c0);
  828.             e_density[p]   += (p + 1 - c0) * FACTOR_W;
  829.             e_density[p+1] += (c0 - p) * FACTOR_W;
  830.         }
  831.         e_density[0]     *= 2.0; // double at the edge bcs working with half-domain (no left/right neighbour)
  832.         e_density[N_G-1] *= 2.0;
  833.         for(p=0; p<N_G; p++) cumul_e_density[p] += e_density[p];
  834.        
  835.         if ((t % N_SUB) == 0) {                                            // ion density - computed in every N_SUB-th time steps (subcycling)
  836.             for(p=0; p<N_G; p++) i_density[p] = 0;
  837.             for(k=0; k<N_i; k++){
  838.                 c0 = x_i[k] * INV_DX;
  839.                 p  = int(c0);
  840.                 i_density[p]   += (p + 1 - c0) * FACTOR_W;  
  841.                 i_density[p+1] += (c0 - p) * FACTOR_W;
  842.             }
  843.             i_density[0]     *= 2.0; // double at the edge bcs working with half-domain (no left/right neighbour)
  844.             i_density[N_G-1] *= 2.0;
  845.         }
  846.         for(p=0; p<N_G; p++) cumul_i_density[p] += i_density[p];
  847.        
  848.         // step 2: solve Poisson equation
  849.        
  850.         for(p=0; p<N_G; p++) rho[p] = E_CHARGE * (i_density[p] - e_density[p]);  // get charge density
  851.         solve_Poisson(rho,Time);                                                 // compute potential and electric field
  852.        
  853.         // steps 3 & 4: move particles according to electric field interpolated to particle positions
  854.        
  855.         for(k=0; k<N_e; k++){                       // move all electrons in every time step
  856.             c0  = x_e[k] * INV_DX;
  857.             p   = int(c0);
  858.             c1  = p + 1.0 - c0;
  859.             c2  = c0 - p;
  860.             e_x = c1 * efield[p] + c2 * efield[p+1];
  861.            
  862.             if (measurement_mode) {
  863.                
  864.                 // measurements: 'x' and 'v' are needed at the same time, i.e. old 'x' and mean 'v'
  865.                
  866.                 mean_v = vx_e[k] - 0.5 * e_x * FACTOR_E;
  867.                 counter_e_xt[p][t_index]   += c1;
  868.                 counter_e_xt[p+1][t_index] += c2;
  869.                 ue_xt[p][t_index]   += c1 * mean_v;
  870.                 ue_xt[p+1][t_index] += c2 * mean_v;
  871.                 v_sqr  = mean_v * mean_v + vy_e[k] * vy_e[k] + vz_e[k] * vz_e[k];
  872.                 energy = 0.5 * E_MASS * v_sqr / EV_TO_J;
  873.                 meanee_xt[p][t_index]   += c1 * energy;
  874.                 meanee_xt[p+1][t_index] += c2 * energy;
  875.                 energy_index = min( int(energy / DE_CS + 0.5), CS_RANGES-1);
  876.                 velocity = sqrt(v_sqr);
  877.                 double local_neut_dens = GAS_DENSITY - e1_density[p] - e2_density[p];
  878.                 rate = sigma[E_ION][energy_index] * velocity * DT_E * local_neut_dens;
  879.                 ioniz_rate_xt[p][t_index]   += c1 * rate;
  880.                 ioniz_rate_xt[p+1][t_index] += c2 * rate;
  881.  
  882.                 // measure EEPF in the center
  883.                
  884.                 if ((MIN_X < x_e[k]) && (x_e[k] < MAX_X)){
  885.                     energy_index = (int)(energy / DE_EEPF);
  886.                     if (energy_index < N_EEPF) {eepf[energy_index] += 1.0;}
  887.                     mean_energy_accu_center += energy;
  888.                     mean_energy_counter_center++;
  889.                 }
  890.             }
  891.            
  892.             // update velocity and position
  893.            
  894.             vx_e[k] -= e_x * FACTOR_E;
  895.             x_e[k]  += vx_e[k] * DT_E;
  896.         }
  897.        
  898.         if ((t % N_SUB) == 0) {                       // move all ions in every N_SUB-th time steps (subcycling)
  899.             for(k=0; k<N_i; k++){
  900.                 c0  = x_i[k] * INV_DX;
  901.                 p   = int(c0);
  902.                 c1  = p + 1 - c0;
  903.                 c2  = c0 - p;
  904.                 e_x = c1 * efield[p] + c2 * efield[p+1];
  905.                
  906.                 if (measurement_mode) {
  907.                    
  908.                     // measurements: 'x' and 'v' are needed at the same time, i.e. old 'x' and mean 'v'
  909.  
  910.                     mean_v = vx_i[k] + 0.5 * e_x * FACTOR_I;
  911.                     counter_i_xt[p][t_index]   += c1;
  912.                     counter_i_xt[p+1][t_index] += c2;
  913.                     ui_xt[p][t_index]   += c1 * mean_v;
  914.                     ui_xt[p+1][t_index] += c2 * mean_v;
  915.                     v_sqr  = mean_v * mean_v + vy_i[k] * vy_i[k] + vz_i[k] * vz_i[k];
  916.                     energy = 0.5 * HE_MASS * v_sqr / EV_TO_J;
  917.                     meanei_xt[p][t_index]   += c1 * energy;
  918.                     meanei_xt[p+1][t_index] += c2 * energy;
  919.                 }
  920.                
  921.                 // update velocity and position and accumulate absorbed energy
  922.                
  923.                 vx_i[k] += e_x * FACTOR_I;
  924.                 x_i[k]  += vx_i[k] * DT_I;
  925.             }
  926.         }
  927.        
  928.         // step 5: check boundaries
  929.        
  930.         k = 0;
  931.         while(k < N_e) {    // check boundaries for all electrons in every time step
  932.             out = false;
  933.             if (x_e[k] < 0) {N_e_abs_pow++; out = true;}    // the electron is out at the powered electrode
  934.             if (x_e[k] > L) {N_e_abs_gnd++; out = true;}    // the electron is out at the grounded electrode
  935.             if (out) {                                      // remove the electron, if out
  936.                 x_e [k] = x_e [N_e-1]; // pushing last element on a vacant place
  937.                 vx_e[k] = vx_e[N_e-1];
  938.                 vy_e[k] = vy_e[N_e-1];
  939.                 vz_e[k] = vz_e[N_e-1];
  940.                 N_e--;
  941.             } else k++;
  942.         }
  943.        
  944.         if ((t % N_SUB) == 0) {   // check boundaries for all ions in every N_SUB-th time steps (subcycling)
  945.             k = 0;
  946.             while(k < N_i) {
  947.                 out = false;
  948.                 if (x_i[k] < 0) {       // the ion is out at the powered electrode
  949.                     N_i_abs_pow++;
  950.                     out    = true;
  951.                     v_sqr  = vx_i[k] * vx_i[k] + vy_i[k] * vy_i[k] + vz_i[k] * vz_i[k];
  952.                     energy = 0.5 * HE_MASS *  v_sqr/ EV_TO_J;
  953.                     energy_index = (int)(energy / DE_IFED);
  954.                     if (energy_index < N_IFED) {ifed_pow[energy_index]++;}       // save IFED at the powered electrode
  955.                 }
  956.                 if (x_i[k] > L) {       // the ion is out at the grounded electrode
  957.                     N_i_abs_gnd++;
  958.                     out    = true;
  959.                     v_sqr  = vx_i[k] * vx_i[k] + vy_i[k] * vy_i[k] + vz_i[k] * vz_i[k];
  960.                     energy = 0.5 * HE_MASS * v_sqr / EV_TO_J;
  961.                     energy_index = (int)(energy / DE_IFED);
  962.                     if (energy_index < N_IFED) {ifed_gnd[energy_index]++;}        // save IFED at the grounded electrode
  963.                 }
  964.                 if (out) {  // delete the ion, if out
  965.                     x_i [k] = x_i [N_i-1];
  966.                     vx_i[k] = vx_i[N_i-1];
  967.                     vy_i[k] = vy_i[N_i-1];
  968.                     vz_i[k] = vz_i[N_i-1];
  969.                     N_i--;
  970.                 } else k++;
  971.             }
  972.         }
  973.        
  974.         // step 6: collisions
  975.  
  976.        
  977.        
  978.         for (k=0; k<N_e; k++){                              // checking for occurrence of a collision for all electrons in every time step
  979.             v_sqr = vx_e[k] * vx_e[k] + vy_e[k] * vy_e[k] + vz_e[k] * vz_e[k];
  980.             velocity = sqrt(v_sqr);
  981.             energy   = 0.5 * E_MASS * v_sqr / EV_TO_J;
  982.             energy_index = min( int(energy / DE_CS + 0.5), CS_RANGES-1);
  983.  
  984.             // Artem  - adding superelastic impact on total collisoinal probability//
  985.  
  986.            
  987.             int e_crdnt = std::max(0, std::min(N_G - 1, static_cast<int>(x_e[k] * INV_DX)));
  988.             double sigma_super_e = sigma[E_SUPER_1][energy_index] * e1_density[e_crdnt] + sigma[E_SUPER_2][energy_index] * e2_density[e_crdnt];
  989.             double ground_dens_local = GAS_DENSITY - e1_density[e_crdnt] - e2_density[e_crdnt];
  990.             double sigma_ground = (sigma[E_ELA][energy_index] + sigma[E_EXC_1][energy_index] +
  991.                                 sigma[E_EXC_2][energy_index] + sigma[E_ION][energy_index]) * ground_dens_local;
  992.             nu = (sigma_ground + sigma_super_e) * velocity;
  993.  
  994.             nu_avg += nu;
  995.            
  996.             p_coll = 1 - exp(- nu * DT_E);                  // collision probability for electrons
  997.             if (R01(MTgen) < p_coll) {                      // electron collision takes place
  998.                 collision_electron(x_e[k], &vx_e[k], &vy_e[k], &vz_e[k], energy_index);
  999.                 N_e_coll++;
  1000.             }
  1001.         }
  1002.        
  1003.         if ((t % N_SUB) == 0) {                             // checking for occurrence of a collision for all ions in every N_SUB-th time steps (subcycling)
  1004.             for (k=0; k<N_i; k++){
  1005.                 vx_a = RMB_n(MTgen);                          // pick velocity components of a random target gas atom
  1006.                 vy_a = RMB_n(MTgen);
  1007.                 vz_a = RMB_n(MTgen);
  1008.                 gx   = vx_i[k] - vx_a;                       // compute the relative velocity of the collision partners
  1009.                 gy   = vy_i[k] - vy_a;
  1010.                 gz   = vz_i[k] - vz_a;
  1011.                 g_sqr = gx * gx + gy * gy + gz * gz;
  1012.                 g = sqrt(g_sqr);
  1013.                 energy = 0.5 * MU_HEHE * g_sqr / EV_TO_J;
  1014.                 energy_index = min( int(energy / DE_CS + 0.5), CS_RANGES-1);
  1015.                 nu = sigma_tot_i[energy_index] * g;
  1016.                 p_coll = 1 - exp(- nu * DT_I);              // collision probability for ions
  1017.                 if (R01(MTgen)< p_coll) {                   // ion collision takes place
  1018.                     collision_ion (&vx_i[k], &vy_i[k], &vz_i[k], &vx_a, &vy_a, &vz_a, energy_index);
  1019.                     N_i_coll++;
  1020.                 }
  1021.             }
  1022.         }
  1023.  
  1024.         //step 7: accumulate rates
  1025.         accumulate_rates();
  1026.  
  1027.         if (measurement_mode) {
  1028.            
  1029.             // collect 'xt' data from the grid
  1030.            
  1031.             for (p=0; p<N_G; p++) {
  1032.                 pot_xt   [p][t_index] += pot[p];
  1033.                 efield_xt[p][t_index] += efield[p];
  1034.                 ne_xt    [p][t_index] += e_density[p];
  1035.                 ni_xt    [p][t_index] += i_density[p];
  1036.                 e1_xt    [p][t_index] += e1_density[p];  // Artem
  1037.                 e2_xt    [p][t_index] += e2_density[p];  // Artem
  1038.             }
  1039.         }
  1040.        
  1041.         if ((t % 1000) == 0) printf(" c = %8d  t = %8d  #e = %8d  #i = %8d\n", cycle,t,N_e,N_i);
  1042.     }
  1043.  
  1044.     counter++;
  1045.  
  1046.     // updating denisites each N_avg cycles: --- Artem
  1047.     if (counter%N_avg == 0) {
  1048.         // compute average rates over a cycle
  1049.         average_rates();
  1050.         // updating densities
  1051.         update_excited_dens();
  1052.         // reset accumulators
  1053.         memset(sum_rate1f, 0, sizeof(sum_rate1f));
  1054.         memset(sum_rate1b, 0, sizeof(sum_rate1b));
  1055.         memset(sum_rate2f, 0, sizeof(sum_rate2f));
  1056.         memset(sum_rate2b, 0, sizeof(sum_rate2b));
  1057.         memset(sum_electron_density, 0, sizeof(sum_electron_density));  
  1058.     }    
  1059.  
  1060.  
  1061.     //calculate nu electron average:
  1062.  
  1063.     nu_avg /= (N_T*N_e);
  1064.  
  1065.     fprintf(datafile,"%8d  %8d  %8d\n",cycle,N_e,N_i);
  1066.     print_excitation_densities();
  1067. }
  1068.  
  1069. //---------------------------------------------------------------------//
  1070. // save particle coordinates                                           //
  1071. //---------------------------------------------------------------------//
  1072.  
  1073. void save_particle_data(){
  1074.     double   d;
  1075.     FILE   * f;
  1076.     char fname[80];
  1077.    
  1078.     strcpy(fname,"picdata.bin");
  1079.     f = fopen(fname,"wb");
  1080.     fwrite(&Time,sizeof(double),1,f);
  1081.     d = (double)(cycles_done);
  1082.     fwrite(&d,sizeof(double),1,f);
  1083.     d = (double)(N_e);
  1084.     fwrite(&d,sizeof(double),1,f);
  1085.     fwrite(x_e, sizeof(double),N_e,f);
  1086.     fwrite(vx_e,sizeof(double),N_e,f);
  1087.     fwrite(vy_e,sizeof(double),N_e,f);
  1088.     fwrite(vz_e,sizeof(double),N_e,f);
  1089.     d = (double)(N_i);
  1090.     fwrite(&d,sizeof(double),1,f);
  1091.     fwrite(x_i, sizeof(double),N_i,f);
  1092.     fwrite(vx_i,sizeof(double),N_i,f);
  1093.     fwrite(vy_i,sizeof(double),N_i,f);
  1094.     fwrite(vz_i,sizeof(double),N_i,f);
  1095.  
  1096.     // saving excited state densities - Artem
  1097.  
  1098.     fwrite(e1_density, sizeof(double), N_G, f);
  1099.     fwrite(e2_density, sizeof(double), N_G, f);    
  1100.  
  1101.     fclose(f);
  1102.     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);
  1103. }
  1104.  
  1105. //---------------------------------------------------------------------//
  1106. // load particle coordinates                                           //
  1107. //---------------------------------------------------------------------//
  1108.  
  1109. void load_particle_data(){
  1110.     double   d;
  1111.     FILE   * f;
  1112.     char fname[80];
  1113.    
  1114.     strcpy(fname,"picdata.bin");    
  1115.     f = fopen(fname,"rb");
  1116.     if (f==NULL) {printf(">> eduPIC: ERROR: No particle data file found, try running initial cycle using argument '0'\n"); exit(0); }
  1117.     fread(&Time,sizeof(double),1,f);
  1118.     fread(&d,sizeof(double),1,f);
  1119.     cycles_done = int(d);
  1120.     fread(&d,sizeof(double),1,f);
  1121.     N_e = int(d);
  1122.     fread(x_e, sizeof(double),N_e,f);
  1123.     fread(vx_e,sizeof(double),N_e,f);
  1124.     fread(vy_e,sizeof(double),N_e,f);
  1125.     fread(vz_e,sizeof(double),N_e,f);
  1126.     fread(&d,sizeof(double),1,f);
  1127.     N_i = int(d);
  1128.     fread(x_i, sizeof(double),N_i,f);
  1129.     fread(vx_i,sizeof(double),N_i,f);
  1130.     fread(vy_i,sizeof(double),N_i,f);
  1131.     fread(vz_i,sizeof(double),N_i,f);
  1132.  
  1133.     // reading excited states densities -- Artem
  1134.  
  1135.     fread(e1_density, sizeof(double), N_G, f);
  1136.     fread(e2_density, sizeof(double), N_G, f);    
  1137.  
  1138.     fclose(f);
  1139.     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);
  1140. }
  1141.  
  1142. //---------------------------------------------------------------------//
  1143. // save density data                                                   //
  1144. //---------------------------------------------------------------------//
  1145.  
  1146. void save_density(void){
  1147.     FILE *f;
  1148.     double c;
  1149.     int m;
  1150.    
  1151.     f = fopen("density.dat","w");
  1152.     c = 1.0 / (double)(no_of_cycles) / (double)(N_T);
  1153.     for(m=0; m<N_G; m++){
  1154.         fprintf(f,"%8.5f  %12e  %12e\n",m * DX, cumul_e_density[m] * c, cumul_i_density[m] * c);
  1155.     }
  1156.     fclose(f);
  1157. }
  1158.  
  1159. // save Final excited states densities - Artem
  1160.  
  1161. void save_final_excited_densities(void) {
  1162.     FILE *f = fopen("excited_densities.dat", "w");
  1163.     for(int p=0; p<N_G; p++) {
  1164.         fprintf(f, "%8.5f %12e %12e\n", p*DX, e1_density[p], e2_density[p]);
  1165.     }
  1166.     fclose(f);
  1167. }
  1168.  
  1169. //---------------------------------------------------------------------//
  1170. // save EEPF data                                                      //
  1171. //---------------------------------------------------------------------//
  1172.  
  1173. void save_eepf(void) {
  1174.     FILE   *f;
  1175.     int    i;
  1176.     double h,energy;
  1177.    
  1178.     h = 0.0;
  1179.     for (i=0; i<N_EEPF; i++) {h += eepf[i];}
  1180.     h *= DE_EEPF;
  1181.     f = fopen("eepf.dat","w");
  1182.     for (i=0; i<N_EEPF; i++) {
  1183.         energy = (i + 0.5) * DE_EEPF;
  1184.         fprintf(f,"%e  %e\n", energy, eepf[i] / h / sqrt(energy));
  1185.     }
  1186.     fclose(f);
  1187. }
  1188.  
  1189. //---------------------------------------------------------------------//
  1190. // save IFED data                                                      //
  1191. //---------------------------------------------------------------------//
  1192.  
  1193. void save_ifed(void) {
  1194.     FILE   *f;
  1195.     int    i;
  1196.     double h_pow,h_gnd,energy;
  1197.    
  1198.     h_pow = 0.0;
  1199.     h_gnd = 0.0;
  1200.     for (i=0; i<N_IFED; i++) {h_pow += ifed_pow[i]; h_gnd += ifed_gnd[i];}
  1201.     h_pow *= DE_IFED;
  1202.     h_gnd *= DE_IFED;
  1203.     mean_i_energy_pow = 0.0;
  1204.     mean_i_energy_gnd = 0.0;
  1205.     f = fopen("ifed.dat","w");
  1206.     for (i=0; i<N_IFED; i++) {
  1207.         energy = (i + 0.5) * DE_IFED;
  1208.         fprintf(f,"%6.2f %10.6f %10.6f\n", energy, (double)(ifed_pow[i])/h_pow, (double)(ifed_gnd[i])/h_gnd);
  1209.         mean_i_energy_pow += energy * (double)(ifed_pow[i]) / h_pow;
  1210.         mean_i_energy_gnd += energy * (double)(ifed_gnd[i]) / h_gnd;
  1211.     }
  1212.     fclose(f);
  1213. }
  1214.  
  1215. //--------------------------------------------------------------------//
  1216. // save XT data                                                       //
  1217. //--------------------------------------------------------------------//
  1218.  
  1219. void save_xt_1(xt_distr distr, char *fname) {
  1220.     FILE   *f;
  1221.     int    i, j;
  1222.    
  1223.     f = fopen(fname,"w");
  1224.     for (i=0; i<N_G; i++){
  1225.         for (j=0; j<N_XT; j++){
  1226.             fprintf(f,"%e  ", distr[i][j]);
  1227.         }
  1228.         fprintf(f,"\n");
  1229.     }
  1230.     fclose(f);
  1231. }
  1232.  
  1233. void norm_all_xt(void){
  1234.     double f1, f2;
  1235.     int    i, j;
  1236.    
  1237.     // normalize all XT data
  1238.    
  1239.     f1 = (double)(N_XT) / (double)(no_of_cycles * N_T);
  1240.     f2 = WEIGHT / (ELECTRODE_AREA * DX) / (no_of_cycles * (PERIOD / (double)(N_XT)));
  1241.    
  1242.     for (i=0; i<N_G; i++){
  1243.         for (j=0; j<N_XT; j++){
  1244.             pot_xt[i][j]    *= f1;
  1245.             efield_xt[i][j] *= f1;
  1246.             ne_xt[i][j]     *= f1;
  1247.             ni_xt[i][j]     *= f1;
  1248.             e1_xt[i][j]     *= f1;   // Artem
  1249.             e2_xt[i][j]     *= f1;   // Artem
  1250.             if (counter_e_xt[i][j] > 0) {
  1251.                 ue_xt[i][j]     =  ue_xt[i][j] / counter_e_xt[i][j];
  1252.                 je_xt[i][j]     = -ue_xt[i][j] * ne_xt[i][j] * E_CHARGE;
  1253.                 meanee_xt[i][j] =  meanee_xt[i][j] / counter_e_xt[i][j];
  1254.                 ioniz_rate_xt[i][j] *= f2;
  1255.              } else {
  1256.                 ue_xt[i][j]         = 0.0;
  1257.                 je_xt[i][j]         = 0.0;
  1258.                 meanee_xt[i][j]     = 0.0;
  1259.                 ioniz_rate_xt[i][j] = 0.0;
  1260.             }
  1261.             if (counter_i_xt[i][j] > 0) {
  1262.                 ui_xt[i][j]     = ui_xt[i][j] / counter_i_xt[i][j];
  1263.                 ji_xt[i][j]     = ui_xt[i][j] * ni_xt[i][j] * E_CHARGE;
  1264.                 meanei_xt[i][j] = meanei_xt[i][j] / counter_i_xt[i][j];
  1265.             } else {
  1266.                 ui_xt[i][j]     = 0.0;
  1267.                 ji_xt[i][j]     = 0.0;
  1268.                 meanei_xt[i][j] = 0.0;
  1269.             }
  1270.             powere_xt[i][j] = je_xt[i][j] * efield_xt[i][j];
  1271.             poweri_xt[i][j] = ji_xt[i][j] * efield_xt[i][j];
  1272.         }
  1273.     }
  1274. }
  1275.  
  1276. void save_all_xt(void){
  1277.     char fname[80];
  1278.    
  1279.     strcpy(fname,"pot_xt.dat");     save_xt_1(pot_xt, fname);
  1280.     strcpy(fname,"efield_xt.dat");  save_xt_1(efield_xt, fname);
  1281.     strcpy(fname,"ne_xt.dat");      save_xt_1(ne_xt, fname);
  1282.     strcpy(fname,"ni_xt.dat");      save_xt_1(ni_xt, fname);
  1283.     strcpy(fname,"je_xt.dat");      save_xt_1(je_xt, fname);
  1284.     strcpy(fname,"ji_xt.dat");      save_xt_1(ji_xt, fname);
  1285.     strcpy(fname,"powere_xt.dat");  save_xt_1(powere_xt, fname);
  1286.     strcpy(fname,"poweri_xt.dat");  save_xt_1(poweri_xt, fname);
  1287.     strcpy(fname,"meanee_xt.dat");  save_xt_1(meanee_xt, fname);
  1288.     strcpy(fname,"meanei_xt.dat");  save_xt_1(meanei_xt, fname);
  1289.     strcpy(fname,"ioniz_xt.dat");   save_xt_1(ioniz_rate_xt, fname);
  1290.     strcpy(fname,"e1_xt.dat");      save_xt_1(e1_xt, fname);
  1291.     strcpy(fname,"e2_xt.dat");      save_xt_1(e2_xt, fname);
  1292. }
  1293.  
  1294. //---------------------------------------------------------------------//
  1295. // simulation report including stability and accuracy conditions       //
  1296. //---------------------------------------------------------------------//
  1297.  
  1298. void check_and_save_info(void){
  1299.     FILE     *f;
  1300.     double   plas_freq, meane, kT, debye_length, density, ecoll_freq, icoll_freq, sim_time, e_max, v_max, power_e, power_i, c;
  1301.     int      i,j;
  1302.     bool     conditions_OK;
  1303.    
  1304.     density    = cumul_e_density[N_G / 2] / (double)(no_of_cycles) / (double)(N_T);  // e density @ center
  1305.     plas_freq  = E_CHARGE * sqrt(density / EPSILON0 / E_MASS);                       // e plasma frequency @ center
  1306.     meane      = mean_energy_accu_center / (double)(mean_energy_counter_center);     // e mean energy @ center
  1307.     kT         = 2.0 * meane * EV_TO_J / 3.0;                                        // k T_e @ center (approximate)
  1308.     sim_time   = (double)(no_of_cycles) / FREQUENCY;                                 // simulated time
  1309.     ecoll_freq = (double)(N_e_coll) / sim_time / (double)(N_e);                      // e collision frequency
  1310.     icoll_freq = (double)(N_i_coll) / sim_time / (double)(N_i);                      // ion collision frequency
  1311.     debye_length = sqrt(EPSILON0 * kT / density) / E_CHARGE;                         // e Debye length @ center
  1312.    
  1313.     f = fopen("info.txt","w");
  1314.     fprintf(f,"########################## eduPIC simulation report ############################\n");
  1315.     fprintf(f,"Simulation parameters:\n");
  1316.     fprintf(f,"Gap distance                          = %12.3e [m]\n",  L);
  1317.     fprintf(f,"# of grid divisions                   = %12d\n",      N_G);
  1318.     fprintf(f,"Frequency                             = %12.3e [Hz]\n", FREQUENCY);
  1319.     fprintf(f,"# of time steps / period              = %12d\n",      N_T);
  1320.     fprintf(f,"# of electron / ion time steps        = %12d\n",      N_SUB);
  1321.     fprintf(f,"Voltage amplitude                     = %12.3e [V]\n",  VOLTAGE);
  1322.     fprintf(f,"Pressure (Ar)                         = %12.3e [Pa]\n", PRESSURE);
  1323.     fprintf(f,"Temperature                           = %12.3e [K]\n",  T_neutral);
  1324.     fprintf(f,"Superparticle weight                  = %12.3e\n",      WEIGHT);
  1325.     fprintf(f,"# of simulation cycles in this run    = %12d\n",      no_of_cycles);
  1326.     fprintf(f,"--------------------------------------------------------------------------------\n");
  1327.     fprintf(f,"Plasma characteristics:\n");
  1328.     fprintf(f,"Electron density @ center             = %12.3e [m^{-3}]\n", density);
  1329.     fprintf(f,"Plasma frequency @ center             = %12.3e [rad/s]\n",  plas_freq);
  1330.     fprintf(f,"Debye length @ center                 = %12.3e [m]\n",      debye_length);
  1331.     fprintf(f,"Electron collision frequency          = %12.3e [1/s]\n",    ecoll_freq);
  1332.     fprintf(f,"Ion collision frequency               = %12.3e [1/s]\n",    icoll_freq);
  1333.     fprintf(f,"--------------------------------------------------------------------------------\n");
  1334.     fprintf(f,"Stability and accuracy conditions:\n");
  1335.     conditions_OK = true;
  1336.     c = plas_freq * DT_E;
  1337.     fprintf(f,"Plasma frequency @ center * DT_E      = %12.3f (OK if less than 0.20)\n", c);
  1338.     if (c > 0.2) {conditions_OK = false;}
  1339.     c = DX / debye_length;
  1340.     fprintf(f,"DX / Debye length @ center            = %12.3f (OK if less than 1.00)\n", c);
  1341.     if (c > 1.0) {conditions_OK = false;}
  1342.     c = max_electron_coll_freq() * DT_E;
  1343.     fprintf(f,"Max. electron coll. frequency * DT_E  = %12.3f (OK if less than 0.05)\n", c);
  1344.     if (c > 0.05) {conditions_OK = false;}
  1345.     c = max_ion_coll_freq() * DT_I;
  1346.     fprintf(f,"Max. ion coll. frequency * DT_I       = %12.3f (OK if less than 0.05)\n", c);
  1347.     if (c > 0.05) {conditions_OK = false;}
  1348.     if (conditions_OK == false){
  1349.         fprintf(f,"--------------------------------------------------------------------------------\n");
  1350.         fprintf(f,"** STABILITY AND ACCURACY CONDITION(S) VIOLATED - REFINE SIMULATION SETTINGS! **\n");
  1351.         fprintf(f,"--------------------------------------------------------------------------------\n");
  1352.         fclose(f);
  1353.         printf(">> eduPIC: ERROR: STABILITY AND ACCURACY CONDITION(S) VIOLATED!\n");
  1354.         printf(">> eduPIC: for details see 'info.txt' and refine simulation settings!\n");
  1355.     }
  1356.     else
  1357.     {
  1358.         // calculate maximum energy for which the Courant-Friedrichs-Levy condition holds:
  1359.        
  1360.         v_max = DX / DT_E;
  1361.         e_max = 0.5 * E_MASS * v_max * v_max / EV_TO_J;
  1362.         fprintf(f,"Max e- energy for CFL condition       = %12.3f [eV]\n", e_max);
  1363.         fprintf(f,"Check EEPF to ensure that CFL is fulfilled for the majority of the electrons!\n");
  1364.         fprintf(f,"--------------------------------------------------------------------------------\n");
  1365.        
  1366.         // saving of the following data is done here as some of the further lines need data
  1367.         // that are computed / normalized in these functions
  1368.        
  1369.         printf(">> eduPIC: saving diagnostics data\n");
  1370.         save_density();
  1371.         save_final_excited_densities();   // Artem
  1372.         save_eepf();
  1373.         save_ifed();
  1374.         norm_all_xt();
  1375.         save_all_xt();
  1376.         fprintf(f,"Particle characteristics at the electrodes:\n");
  1377.         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));
  1378.         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));
  1379.         fprintf(f,"Mean ion energy at powered electrode  = %12.3e [eV]\n", mean_i_energy_pow);
  1380.         fprintf(f,"Mean ion energy at grounded electrode = %12.3e [eV]\n", mean_i_energy_gnd);
  1381.         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));
  1382.         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));
  1383.         fprintf(f,"--------------------------------------------------------------------------------\n");
  1384.        
  1385.         // calculate spatially and temporally averaged power absorption by the electrons and ions
  1386.        
  1387.         power_e = 0.0;
  1388.         power_i = 0.0;
  1389.         for (i=0; i<N_G; i++){
  1390.             for (j=0; j<N_XT; j++){
  1391.                 power_e += powere_xt[i][j];
  1392.                 power_i += poweri_xt[i][j];
  1393.             }
  1394.         }
  1395.         power_e /= (double)(N_XT * N_G);
  1396.         power_i /= (double)(N_XT * N_G);
  1397.         fprintf(f,"Absorbed power calculated as <j*E>:\n");
  1398.         fprintf(f,"Electron power density (average)      = %12.3e [W m^{-3}]\n", power_e);
  1399.         fprintf(f,"Ion power density (average)           = %12.3e [W m^{-3}]\n", power_i);
  1400.         fprintf(f,"Total power density(average)          = %12.3e [W m^{-3}]\n", power_e + power_i);
  1401.         fprintf(f,"--------------------------------------------------------------------------------\n");
  1402.         fclose(f);
  1403.     }
  1404. }
  1405.  
  1406. //------------------------------------------------------------------------------------------//
  1407. // main                                                                                     //
  1408. // command line arguments:                                                                  //
  1409. // [1]: number of cycles (0 for init)                                                       //
  1410. // [2]: "m" turns on data collection and saving                                             //
  1411. //------------------------------------------------------------------------------------------//
  1412.  
  1413. int main (int argc, char *argv[]){
  1414.     printf(">> eduPIC: starting...\n");
  1415.     printf(">> eduPIC: **************************************************************************\n");
  1416.     printf(">> eduPIC: Copyright (C) 2021 Z. Donko et al.\n");
  1417.     printf(">> eduPIC: This program comes with ABSOLUTELY NO WARRANTY\n");
  1418.     printf(">> eduPIC: This is free software, you are welcome to use, modify and redistribute it\n");
  1419.     printf(">> eduPIC: according to the GNU General Public License, https://www.gnu.org/licenses/\n");
  1420.     printf(">> eduPIC: **************************************************************************\n");
  1421.  
  1422.     if (argc == 1) {
  1423.         printf(">> eduPIC: error = need starting_cycle argument\n");
  1424.         return 1;
  1425.     } else {
  1426.         strcpy(st0,argv[1]);
  1427.         arg1 = atol(st0);
  1428.         if (argc > 2) {
  1429.             if (strcmp (argv[2],"m") == 0){
  1430.                 measurement_mode = true;                  // measurements will be done
  1431.             } else {
  1432.                 measurement_mode = false;
  1433.             }
  1434.         }
  1435.     }
  1436.     if (measurement_mode) {
  1437.         printf(">> eduPIC: measurement mode: on\n");
  1438.     } else {
  1439.         printf(">> eduPIC: measurement mode: off\n");
  1440.     }
  1441.     set_electron_cross_sections_ar();
  1442.     set_ion_cross_sections_ar();
  1443.     calc_total_cross_sections();
  1444.  
  1445.     auto exc = init_excited_distr();
  1446.     printf("Excited state population - superparticles!!!\n");
  1447.     printf("%lf %lf\n", exc.first, exc.second);
  1448.  
  1449.     printf("Sanity check \n");
  1450.     fflush(stdout);
  1451.  
  1452.     // for (int i = 0; i < 10; ++i) {
  1453.     //     double E = i * DE_CS;
  1454.     //     printf("σ_super1(%.3f eV) = %e  σ_super2(%.3f eV) = %e\n",
  1455.     //         E, sigma[E_SUPER_1][i],
  1456.     //         E, sigma[E_SUPER_2][i]);
  1457.     // }
  1458.  
  1459.     std::ofstream file0("nu_avg.dat");
  1460.  
  1461.     //test_cross_sections(); return 1;
  1462.     datafile = fopen("conv.dat","a");
  1463.     if (arg1 == 0) {
  1464.         if (FILE *file = fopen("picdata.bin", "r")) { fclose(file);
  1465.             printf(">> eduPIC: Warning: Data from previous calculation are detected.\n");
  1466.             printf("           To start a new simulation from the beginning, please delete all output files before running ./eduPIC 0\n");
  1467.             printf("           To continue the existing calculation, please specify the number of cycles to run, e.g. ./eduPIC 100\n");
  1468.             exit(0);
  1469.         }
  1470.         no_of_cycles = 1;
  1471.         cycle = 1;                                        // init cycle
  1472.         init(N_INIT);                                     // seed initial electrons & ions
  1473.         printf(">> eduPIC: running initializing cycle\n");
  1474.         Time = 0;
  1475.         do_one_cycle();
  1476.         print_excitation_densities();
  1477.         file0 << cycle << " " << nu_avg << "\n";
  1478.         cycles_done = 1;
  1479.     } else {
  1480.         no_of_cycles = arg1;                              // run number of cycles specified in command line
  1481.         load_particle_data();                            // read previous configuration from file
  1482.         printf(">> eduPIC: running %d cycle(s)\n",no_of_cycles);
  1483.         for (cycle=cycles_done+1;cycle<=cycles_done+no_of_cycles;cycle++) {
  1484.             do_one_cycle();
  1485.             file0 << cycle << " " << nu_avg << "\n";
  1486.         }
  1487.         cycles_done += no_of_cycles;
  1488.     }
  1489.     fclose(datafile);
  1490.     save_particle_data();
  1491.     if (measurement_mode) {
  1492.         check_and_save_info();
  1493.     }
  1494.     printf(">> eduPIC: simulation of %d cycle(s) is completed.\n",no_of_cycles);
  1495.     file0.close();
  1496. }
  1497.  
Add Comment
Please, Sign In to add comment