v3.2 - Pastebin.com

#include <cstdio>

#include <cstdlib>

#include <cstring>

#include <cstdbool>

#include <Eigen/Sparse>

#include <cmath>

#include <ctime>

#include <random>

#include <vector>

#include <string>

#include <fstream>

#include <sstream>

#include <algorithm>

#include <stdexcept>

#include <iostream>

#include <iomanip>

#include <sys/stat.h>

#include <errno.h>

#include <stdio.h>

#include <unordered_map>

using namespace::std;

using namespace Eigen;

// constants

const double PI = 3.141592653589793; // mathematical constant Pi

const double TWO_PI = 2.0 * PI; // two times Pi

const double E_CHARGE = 1.60217662e-19; // electron charge [C]

const double EV_TO_J = E_CHARGE; // eV <-> Joule conversion factor

const double E_MASS = 9.109e-31; // mass of electron [kg]

const double HE_MASS = 6.67e-27; // mass of He atom [kg]

const double MU_HEHE = HE_MASS / 2.0; // reduced mass of two He atoms [kg]

const double K_BOLTZMANN = 1.38064852e-23; // Boltzmann's constant [J/K]

const double EPSILON0 = 8.85418781e-12; // permittivity of free space [F/m]

// simulation parameters

const int N_G = 501; // number of grid points

const int N_T = 2000; // time steps within an RF period

const double FREQUENCY = 13.56e6; // driving frequency [Hz]

const double VOLTAGE = 350.0; // voltage amplitude [V]

const double L = 0.04; // electrode gap [m]

const double PRESSURE = 13.3322; // gas pressure [Pa] // n*k*T to match Turner's case

const double T_neutral = 300.0; // background gas temperature [K] (also ion temperature)

const double T_electron = 30000.0; // initial electron temperatutre [K]

const double WEIGHT = 41875.0; // weight of superparticles

const double ELECTRODE_AREA = 1.6e-4; // (fictive) electrode area [m^2]

const int N_INIT = 65536; // number of initial electrons and ions

// additional (derived) constants

const double PERIOD = 1.0 / FREQUENCY; // RF period length [s]

const double DT_E = PERIOD / (double)(N_T); // electron time step [s]

const int N_SUB = 5; // ions move only in these cycles (subcycling)

const int N_avg = 1; // number of cycles to average rates --- Artem

int counter = 0; // cycles since update of rates --- Artem

const double DT_I = N_SUB * DT_E; // ion time step [s]

const double DX = L / (double)(N_G - 1); // spatial grid division [m]

const double INV_DX = 1.0 / DX; // inverse of spatial grid size [1/m]

const double GAS_DENSITY = PRESSURE / (K_BOLTZMANN * T_neutral); // background gas density [1/m^3]

const double OMEGA = TWO_PI * FREQUENCY; // angular frequency [rad/s]

const double Diff_coeff = 8.992E-6*pow(T_neutral,1.5)/(PRESSURE/133.3); // diffusion coefficient for metastable states

const double kT_over_m = K_BOLTZMANN*T_electron/E_MASS; // KT/m for calculations

// electron and ion cross sections

const int N_CS = 8; // total number of processes / cross sections

const int E_ELA = 0; // process identifier: electron/elastic

const int E_ION = 1; // process identifier: electron/ionization

const int I_ISO = 2; // process identifier: ion/elastic/isotropic

const int I_BACK = 3; // process identifier: ion/elastic/backscattering

const int E_EXC_1 = 4; // process identifier: electron/excitation1

const int E_EXC_2 = 5; // process identifier: electron/excitation1

const int E_SUPER_1 = 6; // triplet superelastic - Artem

const int E_SUPER_2 = 7; // singlet superelastic - Artem

const double E_EXC_TH_1 = 19.82; // electron impact excitation threshold [eV]

const double E_EXC_TH_2 = 20.61; // electron impact excitation threshold [eV]

const double E_ION_TH = 24.587; // electron impact ionization threshold [eV]

const int CS_RANGES = 1000000; // number of entries in cross section arrays

const double DE_CS = 0.001; // energy division in cross section arrays [eV]

typedef float cross_section[CS_RANGES]; // cross section array

cross_section sigma[N_CS]; // set of cross section arrays

cross_section sigma_tot_e; // total macroscopic cross section of electrons

cross_section sigma_tot_i; // total macroscopic cross section of ions

double nu_avg; // average nu for electrons thorugh 1 cycle

// particle coordinates

const int MAX_N_P = 10000000; // maximum number of particles (electrons / ions)

typedef double particle_vector[MAX_N_P]; // array for particle properties

int N_e = 0; // number of electrons

int N_i = 0; // number of ions

int N_e1 = 0; // number of singlet excited states

int N_e2 = 0; // number of triplet excited states

particle_vector x_e, vx_e, vy_e, vz_e; // coordinates of electrons (one spatial, three velocity components)

particle_vector x_i, vx_i, vy_i, vz_i; // coordinates of ions (one spatial, three velocity components)

typedef double xvector[N_G]; // array for quantities defined at gird points

xvector efield,pot; // electric field and potential

xvector e_density,i_density; // electron and ion densities

xvector cumul_e_density,cumul_i_density; // cumulative densities

//excited states handling

xvector e1_density;

xvector e2_density;

xvector sum_electron_density = {0.0}; xvector avg_electron_density = {0.0};

xvector sum_rate1f = {0.0}; xvector sum_rate1b = {0.0}; xvector sum_rate2f = {0.0}; xvector sum_rate2b = {0.0};

xvector avg_rate1f = {0.0}; xvector avg_rate1b = {0.0}; xvector avg_rate2f = {0.0}; xvector avg_rate2b = {0.0};

// array for Triplet excitation rates!!! Artem

xvector S1 = {0.0};

xvector S2 = {0.0}; // source terms for DRR module -- Artem

typedef unsigned long long int Ullong; // compact name for 64 bit unsigned integer

Ullong N_e_abs_pow = 0; // counter for electrons absorbed at the powered electrode

Ullong N_e_abs_gnd = 0; // counter for electrons absorbed at the grounded electrode

Ullong N_i_abs_pow = 0; // counter for ions absorbed at the powered electrode

Ullong N_i_abs_gnd = 0; // counter for ions absorbed at the grounded electrode

// electron energy probability function

const int N_EEPF = 200000; // number of energy bins in Electron Energy Probability Function (EEPF)

const double DE_EEPF = 0.05; // resolution of EEPF [eV]

typedef double eepf_vector[N_EEPF]; // array for EEPF

eepf_vector eepf = {0.0}; // time integrated EEPF in the center of the plasma

// ION VELOCITY and ENERGY PROBABILITY FUNCTION --- Artem

//spatial resolution

const int N_sp = 100; // number of spatial bins of the discharge regions

const double DX_PF = L/N_sp; // corresponded spatial step for probability functions

//IVDF

const int N_IVDF = 50000;

const double DV_IVDF = 10.0;

std::vector<std::vector<double>>

ivdf_pos(N_sp, std::vector<double>(N_IVDF, 0.0)), // positive part of IVDF

ivdf_neg(N_sp, std::vector<double>(N_IVDF, 0.0)); // negative part of IVDF

//IEPF

const int N_IEPF = 200000; // number of energy bins in Ion Energy Probability Function (EEPF)

const double DE_IEPF = 0.001; // resolution of IEPF [eV]

std::vector<std::vector<double>>

iepf(N_sp, std::vector<double>(N_IEPF, 0.0)); // IEDF array

// ion flux-energy distributions

const int N_IFED = 200; // number of energy bins in Ion Flux-Energy Distributions (IFEDs)

const double DE_IFED = 1.0; // resolution of IFEDs [eV]

typedef int ifed_vector[N_IFED]; // array for IFEDs

ifed_vector ifed_pow = {0}; // IFED at the powered electrode

ifed_vector ifed_gnd = {0}; // IFED at the grounded electrode

double mean_i_energy_pow; // mean ion energy at the powered electrode

double mean_i_energy_gnd; // mean ion energy at the grounded electrode

// spatio-temporal (XT) distributions

const int N_BIN = 10; // number of time steps binned for the XT distributions

const int N_XT = N_T / N_BIN; // number of spatial bins for the XT distributions

typedef double xt_distr[N_G][N_XT]; // array for XT distributions (decimal numbers)

xt_distr pot_xt = {0.0}; // XT distribution of the potential

xt_distr efield_xt = {0.0}; // XT distribution of the electric field

xt_distr ne_xt = {0.0}; // XT distribution of the electron density

xt_distr ni_xt = {0.0}; // XT distribution of the ion density

xt_distr ue_xt = {0.0}; // XT distribution of the mean electron velocity

xt_distr ui_xt = {0.0}; // XT distribution of the mean ion velocity

xt_distr je_xt = {0.0}; // XT distribution of the electron current density

xt_distr ji_xt = {0.0}; // XT distribution of the ion current density

xt_distr powere_xt = {0.0}; // XT distribution of the electron powering (power absorption) rate

xt_distr poweri_xt = {0.0}; // XT distribution of the ion powering (power absorption) rate

xt_distr meanee_xt = {0.0}; // XT distribution of the mean electron energy

xt_distr meanei_xt = {0.0}; // XT distribution of the mean ion energy

xt_distr counter_e_xt = {0.0}; // XT counter for electron properties

xt_distr counter_i_xt = {0.0}; // XT counter for ion properties

xt_distr ioniz_rate_xt = {0.0}; // XT distribution of the ionisation rate

xt_distr e1_xt = {0.0}; // XT distribution of triplet excited states density - Artem

xt_distr e2_xt = {0.0}; // XT distribution of singlet excited states density - Artem

double mean_energy_accu_center = 0; // mean electron energy accumulator in the center of the gap

Ullong mean_energy_counter_center = 0; // mean electron energy counter in the center of the gap

Ullong N_e_coll = 0; // counter for electron collisions

Ullong N_i_coll = 0; // counter for ion collisions

double Time; // total simulated time (from the beginning of the simulation)

int cycle,no_of_cycles,cycles_done; // current cycle and total cycles in the run, cycles completed

int arg1; // used for reading command line arguments

char st0[80]; // used for reading command line arguments

FILE *datafile; // used for saving data

bool measurement_mode; // flag that controls measurements and data saving

const double gamma_i = 0.0; // yielding coefficient for ions bombarding the surface

const double gamma_m = 0.3; // yielding coefficient for metastables bombarding the surface

const double gamma_p = -10.0; // yielding coefficient for photons bombarding the surface

const double gamma_k = 0.5; // recombination coefficient

const double r = 0.0; // probabilty for the electron to reflect at the surface

const double T_SEE = 300; // initial temperature of SEE electrons [K]

const double Ry = 13.6057; // Rydberg constant, eV

const double Bohr = 0.8797e-20; // atomic consant related to Bohr Radius, m^2

//---------------------------------------------------------------------------//

// C++ Mersenne Twister 19937 generator //

// R01(MTgen) will genarate uniform distribution over [0,1) interval //

// RMB(MTgen) will generate Maxwell-Boltzmann distribution (of gas atoms) //

//---------------------------------------------------------------------------//

std::random_device rd{};

std::mt19937 MTgen(rd());

std::uniform_real_distribution<> R01(0.0, 1.0);

std::normal_distribution<> RMB_n(0.0,sqrt(K_BOLTZMANN * T_neutral / HE_MASS));

std::normal_distribution<> RMB_e(0.0,sqrt(K_BOLTZMANN * T_electron / E_MASS));

std::normal_distribution<> RMB_SEE(0.0,sqrt(K_BOLTZMANN * T_SEE / E_MASS));

enum class ProcKind { ELASTIC, IONIZATION, EXCITATION, SUPERELASTIC };

struct Process {

ProcKind kind;

std::string name; // e.g. "1^1S→2^3S"

int csSlot; // index into sigma (filled later)

int init_Lvl; // index in levels vector

int final_Lvl; // "

double threshold_eV; // 0 for elastic

double A1 = 0.0, A2 = 0.0, A3 = 0.0, A4 = 0.0, A5 = 0.0, A6 =0.0; // coefficients for cross-section calculation [Ralchenko2008]

char type; // type of a process - dipole allowed/forbidden etc..

std::vector<double> cum_rate; // storing cumulative rate of this process, feed to avg_rate() function

std::vector<double> avg_rate; //stoding average rate of this process, feed to DRR

};

struct level {

std::string name;

double th_energy;

int weight;

std::vector<double> density;

};

struct transition {

std::string lvl_init;

std::string lvl_final;

char type;

double thresh_energy;

double A1, A2, A3, A4, A5, A6;

int is_superelastic;

};

int find_lvl(int init_lvl, int fin_lvl, const std::vector<Process>& eProc){

int check = -1;

string name;

for (auto const& proc : eProc){

if ((proc.kind == ProcKind::EXCITATION) && (proc.init_Lvl == fin_lvl) && (proc.final_Lvl == init_lvl)){

return proc.csSlot;

check += 1;

}

if (check != 0) {

for (auto const& proc : eProc){

if ((proc.kind == ProcKind::SUPERELASTIC) && (proc.init_Lvl == init_lvl) && (proc.final_Lvl == fin_lvl)){

name = proc.name;

}

std::cout << "Can't find an excitation pair for superelastic process: " << name << "\n";

return -1;

}

class CSInterpolator {

public:

// load "filename" which must have two whitespace‐separated columns:

// energy (eV) cross_section (in m^2 or cm^2 as you prefer)

CSInterpolator(const std::string &filename) {

std::ifstream in(filename);

if (!in) throw std::runtime_error("CSInterpolator: cannot open " + filename);

double E, sigma;

std::string line;

while (std::getline(in, line)) {

std::istringstream iss(line);

if (iss >> E >> sigma) {

E_pts_.push_back(E);

sigma_pts_.push_back(sigma);

}

if (E_pts_.size()<2)

throw std::runtime_error("CSInterpolator: need at least two data points in " + filename);

}

// return σ(E) by simple linear interpolation (clamped to end‐points)

double operator()(double E) const {

auto it = std::lower_bound(E_pts_.begin(), E_pts_.end(), E);

if (it == E_pts_.begin()) {

// std::cerr << "Warning: E="<<E<<" below data range, clamping to "<< 0.0 <<"\n";

return 0.0;

}

if (it == E_pts_.end()){

// std::cerr << "Warning: E="<<E<<" above data range, clamping to "<< sigma_pts_.back() <<"\n";

return sigma_pts_.back();

}

size_t idx = (it - E_pts_.begin());

double E1 = E_pts_[idx-1], E2 = E_pts_[idx];

double s1 = sigma_pts_[idx-1], s2 = sigma_pts_[idx];

// linear interp

return s1 + (s2 - s1) * (E - E1)/(E2 - E1);

}

private:

std::vector<double> E_pts_, sigma_pts_;

};

// Cross-Section function for excitation cross-sections.

// here for type 'd' - dipole-forbidden, 's' - spin-forbidden, 'a' - dipole-allowed

// g - stat weight of initial level (g = 1 for ground state). See paper: RalchenkoJanev2008

double forbidden_CS(double A1, double A2, double A3, double A4, double A5, double A6, double E, double E_th, int g, char type){

double x = E/E_th;

double strength, S1, S2, CS;

if (type == 'd'){

S1 = (A1 + A2/x + A3/(x*x) + A4/(x*x*x));

S2 = (x*x/(x*x + A5));

}

if (type == 's'){

S1 = (A1 + A2/x + A3/(x*x) + A4/(x*x*x));

S2 = (1.0/(x*x + A5));

}

if (type == 'a'){

S1 = (A1*log(x) + A2 + A3/x + A4/(x*x) + A5/(x*x*x));

S2 = (x + 1.0)/(x + A6);

}

strength = S1*S2;

CS = Bohr * Ry/(g*E) * strength;

return CS;

}

// N_CS - number of cross-sections

void set_electron_cross_sections_ar(const std::vector<Process>& eProc, const std::vector<level>& levels, std::vector<std::vector<float>>& sigma_temp){

double A1, A2, A3, A4, A5, A6;

double threshold;

char type;

int g;

printf(">> eduPIC: Setting e- / He cross sections\n");

// load your four datafiles (make sure these names match your files!)

CSInterpolator cs_ela ("CS/He_electron_elastic.dat"); // two‐col: E σ_ela

CSInterpolator cs_ion ("CS/He_electron_ionization.dat");// two‐col: E σ_ion

int P = eProc.size();

for(int i=0; i<CS_RANGES; i++){

double en = (i==0 ? DE_CS : DE_CS * i); // energy grid

for (int s = 0; s < P; s++){

switch(eProc[s].kind){

case ProcKind::ELASTIC:{sigma_temp[s][i] = cs_ela(en); break;}

case ProcKind::IONIZATION: {sigma_temp[s][i] = cs_ion(en); break;}

case ProcKind::EXCITATION: {

double threshold = eProc[s].threshold_eV;

if (en < threshold) { // <–– guard reinstated

sigma_temp[s][i] = 0.0f;

break;

}

int g = levels[eProc[s].init_Lvl].weight; // g_init

sigma_temp[s][i] = forbidden_CS(

eProc[s].A1, eProc[s].A2, eProc[s].A3,

eProc[s].A4, eProc[s].A5, eProc[s].A6,

en, threshold, g, eProc[s].type);

if (sigma_temp[s][i] < 0.0f) // numerical safety

sigma_temp[s][i] = 0.0f;

break;

} break;

default: ;

}

for (int i = 0; i<CS_RANGES; i++){

for (int s = 0; s < P; s++){

if (eProc[s].kind != ProcKind::SUPERELASTIC) continue;

double en = (i==0 ? DE_CS : DE_CS * i);

threshold = eProc[s].threshold_eV;

double enPlus = en + threshold;

int idx = enPlus / DE_CS;

int lvl_init = eProc[s].init_Lvl;

int lvl_fin = eProc[s].final_Lvl;

int excitation_pair = find_lvl(lvl_init, lvl_fin, eProc);

sigma_temp[s][i] = (idx<CS_RANGES && en>0)

? (static_cast<double>(levels[eProc[s].final_Lvl].weight) / static_cast<double>(levels[eProc[s].init_Lvl].weight)) *

(enPlus/en) * sigma_temp[excitation_pair][idx] : 0.0;

}

//------------------------------------------------------------------------------//

// ion cross sections: A. V. Phelps, J. Appl. Phys. 76, 747 (1994) //

//------------------------------------------------------------------------------//

void set_ion_cross_sections_ar(void){

int i;

double e_com,e_lab,qmom,qback,qiso;

printf(">> eduPIC: Setting He+ / He cross sections\n");

for(i=0; i<CS_RANGES; i++){

if (i == 0) {e_com = DE_CS;} else {e_com = DE_CS * i;} // ion energy in the center of mass frame of reference

e_lab = 2.0 * e_com; // ion energy in the laboratory frame of reference

qiso = 7.63 *pow(10,-20) * pow(e_com, -0.5);

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 );

sigma[I_ISO][i] = qiso; // cross section for He+ / He isotropic part of elastic scattering

sigma[I_BACK][i] = qback; // cross section for He+ / He backward elastic scattering

}

//----------------------------------------------------------------------//

// calculation of total cross sections for electrons and ions //

//----------------------------------------------------------------------//

void calc_total_cross_sections(void){

int i;

for(i=0; i<CS_RANGES; i++){

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

sigma_tot_i[i] = (sigma[I_ISO][i] + sigma[I_BACK][i]) * GAS_DENSITY; // total macroscopic cross section of ions

}

//----------------------------------------------------------------------//

// test of cross sections for electrons and ions //

//----------------------------------------------------------------------//

void test_cross_sections(void){

FILE * f;

int i,j;

f = fopen("cross_sections.dat","w"); // cross sections saved in data file: cross_sections.dat

for(i=0; i<CS_RANGES; i++){

fprintf(f,"%12.4f ",i*DE_CS);

for(j=0; j<N_CS; j++) fprintf(f,"%14e ",sigma[j][i]);

fprintf(f,"\n");

}

fclose(f);

}

//---------------------------------------------------------------------//

// find upper limit of collision frequencies //

//---------------------------------------------------------------------//

double max_electron_coll_freq (void){

int i;

double e,v,nu,nu_max;

nu_max = 0;

for(i=0; i<CS_RANGES; i++){

e = i * DE_CS;

v = sqrt(2.0 * e * EV_TO_J / E_MASS);

nu = v * sigma_tot_e[i];

if (nu > nu_max) {nu_max = nu;}

}

return nu_max;

}

double max_ion_coll_freq (void){

int i;

double e,g,nu,nu_max;

nu_max = 0;

for(i=0; i<CS_RANGES; i++){

e = i * DE_CS;

g = sqrt(2.0 * e * EV_TO_J / MU_HEHE);

nu = g * sigma_tot_i[i];

if (nu > nu_max) nu_max = nu;

}

return nu_max;

}

//----------------------------------------------------------------------//

// initialization of the simulation by placing a given number of //

// electrons and ions at random positions between the electrodes //

//----------------------------------------------------------------------//

void print_excitation_densities(std::vector<level>& levels) {

double total_e1 = 0.0, total_e2 = 0.0;

const double cell_volume = ELECTRODE_AREA * DX;

// Sum densities across all grid cells

for (int p = 0; p < N_G; p++) {

total_e1 += levels[1].density[p]; // triplet state density

total_e2 += levels[2].density[p]; // singlet state density

}

printf("average Triplet SP = %8.2e | average Singlet SP = %8.2e\n", total_e1/N_G, total_e2/N_G);

}

void init(int nseed, std::vector<level>& levels){

int i;

for (i=0; i<nseed; i++){

x_e[i] = L * R01(MTgen); // initial random position of the electron

vx_e[i] = RMB_e(MTgen); vy_e[i] = RMB_e(MTgen); vz_e[i] = RMB_e(MTgen); // initial velocity components of the electron

x_i[i] = L * R01(MTgen); // initial random position of the ion

vx_i[i] = RMB_n(MTgen); vy_i[i] = RMB_n(MTgen); vz_i[i] = RMB_n(MTgen); // initial velocity components of the ion

}

N_e = nseed; // initial number of electrons

N_i = nseed; // initial number of ions

double n0 = 2.7e+9; double a0 = -1.3e+8; double b0 = 5.0e+15;

for (int p= 0; p < N_G; p++){

for (auto& lvl:levels){

if (lvl.name == "1^1S") lvl.density[p] = GAS_DENSITY;

else if (lvl.name == "2^3S"){

lvl.density[p] = abs(n0 * (a0 *(pow(L,4) - pow((p*DX - 0.5*L),2)) +

b0 * (pow(L,8) - pow((p*DX - 0.5*L),8)) ));

}

else if (lvl.name == "2^1S") {

lvl.density[p] = 0.5 * abs( n0 * (a0 *(pow(L,4) - pow((p*DX - 0.5*L),2)) +

b0 * (pow(L,8) - pow((p*DX - 0.5*L),8)) ));

}

//----------------------------------------------------------------------//

// e / He collision (cold gas approximation) //

//----------------------------------------------------------------------//

void collision_electron (double xe, double *vxe, double *vye, double *vze, int eindex,

const std::vector<level>& levels, const std::vector<Process>& eProc, const std::vector<std::vector<float>>& sigma_temp){

const double F1 = E_MASS / (E_MASS + HE_MASS);

const double F2 = HE_MASS / (E_MASS + HE_MASS);

double t0,t1,t2,t3,t4,t5;

double g,g2,gx,gy,gz,wx,wy,wz,theta,phi;

double chi,eta,chi2,eta2,sc,cc,se,ce,st,ct,sp,cp,energy,threshold,e_sc,e_ej;

double probability_total = 0.0;

double probability_cumulative = 0.0;

size_t sel = 0;

double rnd = R01(MTgen);

// - Artem

// Determine cell p where the electron is

double c0 = xe * INV_DX;

int pos = std::max(0, std::min(N_G - 1, static_cast<int>(c0)));

for (int i = 0; i < eProc.size(); i++){

const Process& ch = eProc[i];

double dens = (ch.kind == ProcKind::SUPERELASTIC) ? levels[ch.init_Lvl].density[pos] : GAS_DENSITY;

probability_total += sigma_temp[i][eindex] * dens;

}

for (int i = 0; i < eProc.size(); i++) {

const Process& ch = eProc[i];

double dens = (ch.kind == ProcKind::SUPERELASTIC) ? levels[ch.init_Lvl].density[pos] : GAS_DENSITY;

probability_cumulative += sigma_temp[i][eindex] * dens;

if (rnd < probability_cumulative/probability_total) {sel = i; break;} // if cumulative probability/total is larger than a random number, stop and save id of a Process

}

const Process& ch = eProc[sel]; // now we've chosen the process, performing kineitcs:

//------------------------------------------------------------------------------------------------/////

// calculate relative velocity before collision & velocity of the centre of mass

gx = (*vxe);

gy = (*vye);

gz = (*vze);

g = sqrt(gx * gx + gy * gy + gz * gz);

wx = F1 * (*vxe);

wy = F1 * (*vye);

wz = F1 * (*vze);

// find Euler angles

if (gx == 0) {theta = 0.5 * PI;}

else {theta = atan2(sqrt(gy * gy + gz * gz),gx);}

if (gy == 0) {

if (gz > 0){phi = 0.5 * PI;} else {phi = - 0.5 * PI;}

} else {phi = atan2(gz, gy);}

st = sin(theta);

ct = cos(theta);

sp = sin(phi);

cp = cos(phi);

//----------------------------------------------------------------------------------------///

switch (ch.kind) {

case ProcKind::ELASTIC:{

chi = acos(1.0 - 2.0 * R01(MTgen)); // isotropic scattering

eta = TWO_PI * R01(MTgen); // azimuthal angle

break;

}

case ProcKind::EXCITATION:{

energy = 0.5 * E_MASS * g * g; // electron energy

threshold = ch.threshold_eV; // threshold energy

energy = fabs(energy - threshold * EV_TO_J); // subtract energy loss for excitation

g = sqrt(2.0 * energy / E_MASS); // relative velocity after energy loss

chi = acos(1.0 - 2.0 * R01(MTgen)); // isotropic scattering

eta = TWO_PI * R01(MTgen); // azimuthal angle

break;

}

case ProcKind::IONIZATION:{

energy = 0.5 * E_MASS * g * g; // electron energy

threshold = ch.threshold_eV; // threshold energy

energy = fabs(energy - threshold * EV_TO_J); // subtract energy loss of ionization

e_ej = 0.5*energy; // energy of the ejected electron

e_sc = fabs(energy - e_ej); // energy of scattered electron after the collision

g = sqrt(2.0 * e_sc / E_MASS); // relative velocity of scattered electron

g2 = sqrt(2.0 * e_ej / E_MASS); // relative velocity of ejected electron

chi = acos(1.0 - 2.0 * R01(MTgen)); // isotropic scattering for scattered electron (as in Turner's case)

chi2 = acos(1.0 - 2.0 * R01(MTgen)); // isotropic scattering for ejected electrons (as in Turner's case)

eta = TWO_PI * R01(MTgen); // azimuthal angle for scattered electron

eta2 = eta + PI; // azimuthal angle for ejected electron

sc = sin(chi2);

cc = cos(chi2);

se = sin(eta2);

ce = cos(eta2);

gx = g2 * (ct * cc - st * sc * ce);

gy = g2 * (st * cp * cc + ct * cp * sc * ce - sp * sc * se);

gz = g2 * (st * sp * cc + ct * sp * sc * ce + cp * sc * se);

x_e[N_e] = xe; // add new electron

vx_e[N_e] = wx + F2 * gx;

vy_e[N_e] = wy + F2 * gy;

vz_e[N_e] = wz + F2 * gz;

N_e++;

x_i[N_i] = xe; // add new ion

vx_i[N_i] = RMB_n(MTgen); // velocity is sampled from background thermal distribution

vy_i[N_i] = RMB_n(MTgen);

vz_i[N_i] = RMB_n(MTgen);

N_i++;

break;

}

case ProcKind::SUPERELASTIC:{

energy = 0.5 * E_MASS * g * g; // electron energy

threshold = ch.threshold_eV; // threshold energy

energy = fabs(energy + threshold * EV_TO_J); // add energy gain after superelastic collision

g = sqrt(2.0 * energy / E_MASS); // relative velocity after energy loss

chi = acos(1.0 - 2.0 * R01(MTgen)); // isotropic scattering

eta = TWO_PI * R01(MTgen); // azimuthal angle

break;

}

// scatter the primary electron

sc = sin(chi);

cc = cos(chi);

se = sin(eta);

ce = cos(eta);

// compute new relative velocity:

gx = g * (ct * cc - st * sc * ce);

gy = g * (st * cp * cc + ct * cp * sc * ce - sp * sc * se);

gz = g * (st * sp * cc + ct * sp * sc * ce + cp * sc * se);

// post-collision velocity of the colliding electron

(*vxe) = wx + F2 * gx;

(*vye) = wy + F2 * gy;

(*vze) = wz + F2 * gz;

}

//----------------------------------------------------------------------//

// He+ / He collision //

//----------------------------------------------------------------------//

void collision_ion (double *vx_1, double *vy_1, double *vz_1,

double *vx_2, double *vy_2, double *vz_2, int e_index){

double g,gx,gy,gz,wx,wy,wz,rnd;

double theta,phi,chi,eta,st,ct,sp,cp,sc,cc,se,ce,t1,t2;

// calculate relative velocity before collision

// random Maxwellian target atom already selected (vx_2,vy_2,vz_2 velocity components of target atom come with the call)

gx = (*vx_1)-(*vx_2);

gy = (*vy_1)-(*vy_2);

gz = (*vz_1)-(*vz_2);

g = sqrt(gx * gx + gy * gy + gz * gz);

wx = 0.5 * ((*vx_1) + (*vx_2));

wy = 0.5 * ((*vy_1) + (*vy_2));

wz = 0.5 * ((*vz_1) + (*vz_2));

// find Euler angles

if (gx == 0) {theta = 0.5 * PI;} else {theta = atan2(sqrt(gy * gy + gz * gz),gx);}

if (gy == 0) {

if (gz > 0){phi = 0.5 * PI;} else {phi = - 0.5 * PI;}

} else {phi = atan2(gz, gy);}

// determine the type of collision based on cross sections and generate scattering angle

t1 = sigma[I_ISO][e_index];

t2 = t1 + sigma[I_BACK][e_index];

rnd = R01(MTgen);

if (rnd < (t1 /t2)){ // isotropic scattering

chi = acos(1.0 - 2.0 * R01(MTgen)); // scattering angle

} else { // backward scattering

chi = PI; // scattering angle

}

eta = TWO_PI * R01(MTgen); // azimuthal angle

sc = sin(chi);

cc = cos(chi);

se = sin(eta);

ce = cos(eta);

st = sin(theta);

ct = cos(theta);

sp = sin(phi);

cp = cos(phi);

// compute new relative velocity

gx = g * (ct * cc - st * sc * ce);

gy = g * (st * cp * cc + ct * cp * sc * ce - sp * sc * se);

gz = g * (st * sp * cc + ct * sp * sc * ce + cp * sc * se);

// post-collision velocity of the ion

(*vx_1) = wx + 0.5 * gx;

(*vy_1) = wy + 0.5 * gy;

(*vz_1) = wz + 0.5 * gz;

}

//-----------------------------------------------------------------//

// solve Poisson equation (Thomas algorithm) //

//-----------------------------------------------------------------//

void solve_Poisson (xvector rho1, double tt){

const double A = 1.0;

const double B = -2.0;

const double C = 1.0;

const double S = 1.0 / (2.0 * DX);

const double ALPHA = -DX * DX / EPSILON0;

xvector g, w, f;

int i;

// apply potential to the electrodes - boundary conditions

pot[0] = VOLTAGE * cos(OMEGA * tt); // potential at the powered electrode

pot[N_G-1] = 0.0; // potential at the grounded electrode

// solve Poisson equation

for(i=1; i<=N_G-2; i++) f[i] = ALPHA * rho1[i];

f[1] -= pot[0];

f[N_G-2] -= pot[N_G-1];

w[1] = C/B;

g[1] = f[1]/B;

for(i=2; i<=N_G-2; i++){

w[i] = C / (B - A * w[i-1]);

g[i] = (f[i] - A * g[i-1]) / (B - A * w[i-1]);

}

pot[N_G-2] = g[N_G-2];

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

// compute electric field

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

efield[0] = (pot[0] - pot[1]) * INV_DX - rho1[0] * DX / (2.0 * EPSILON0); // powered electrode

efield[N_G-1] = (pot[N_G-2] - pot[N_G-1]) * INV_DX + rho1[N_G-1] * DX / (2.0 * EPSILON0); // grounded electrode

}

//---------------------------------------------------------------------//

// simulation of one radiofrequency cycle //

//---------------------------------------------------------------------//

void accumulate_rates(std::vector<Process>& eProc, const std::vector<std::vector<float>>& sigma_temp) {

double v_sqr, velocity, energy, c0_temp;

int energy_index, p_temp;

const double Volume = (ELECTRODE_AREA * DX);

for (int k=0; k<N_e; k++){

v_sqr = vx_e[k] * vx_e[k] + vy_e[k] * vy_e[k] + vz_e[k] * vz_e[k];

velocity = sqrt(v_sqr);

energy = 0.5 * E_MASS * v_sqr / EV_TO_J;

energy_index = min( int(energy / DE_CS + 0.5), CS_RANGES-1);

c0_temp = x_e[k] * INV_DX;

p_temp = std::max(0, std::min(N_G-1, static_cast<int>(c0_temp)));

for (int i = 0; i < eProc.size(); i++){

Process& ch = eProc[i];

if (ch.kind == ProcKind::ELASTIC || ch.kind == ProcKind::IONIZATION) continue;

ch.cum_rate[p_temp] += sigma_temp[ch.csSlot][energy_index] * velocity * WEIGHT;

}

// averaging the rates each RF cycle

void average_rates(std::vector<Process>& eProc) {

const double inv_NT = 1.0 / (N_avg * N_T); // averaging through N_avg RF periods each contains N_T cycles

const double inv_Volume = 1.0/(ELECTRODE_AREA * DX);

for (int p = 0; p < N_G; p++) {

for (int i = 0; i < eProc.size(); i++){

Process& ch = eProc[i];

if (ch.kind == ProcKind::ELASTIC || ch.kind == ProcKind::IONIZATION) continue;

ch.avg_rate[p] = ch.cum_rate[p] * (inv_NT * inv_Volume);

}

void update_excited_dens(const std::vector<Process>& eProc, std::vector<level>& levels) {

using SpMat = Eigen::SparseMatrix<double>;

SpMat A(3*N_G,3*N_G); A.reserve(Eigen::VectorXi::Constant(3*N_G, 9));

SpMat B(3*N_G,3*N_G); B.reserve(Eigen::VectorXi::Constant(3*N_G, 9));

Eigen::VectorXd sol_old(3*N_G), sol_new(3*N_G); // solution for 1^1S, 2^3S, 2^1S

Eigen::VectorXd RHS(3*N_G); // RHS for our new equation

double v_mean = sqrt(8.0*K_BOLTZMANN*T_neutral/(M_PI*HE_MASS));

double alpha = 1.0;

double DT = alpha*2.0*DX*DX/Diff_coeff;

double Rk = (DX*gamma_k*v_mean)/((2.0-gamma_k)*Diff_coeff);

// temporal assignment of rates:

for (int i = 0; i < N_G; i++){

avg_rate1f[i] = eProc[2].avg_rate[i];

avg_rate2f[i] = eProc[3].avg_rate[i];

avg_rate1b[i] = eProc[4].avg_rate[i];

avg_rate2b[i] = eProc[5].avg_rate[i];

}

for (int i = 0; i < N_G; i++){

double Ai = -alpha;

double Ci = -alpha;

double B0i = 1.0 + 2.0*alpha + 0.5*DT*(avg_rate1f[i] + avg_rate2f[i]); //1^1S

double B1i = 1.0 + 2.0*alpha + 0.5*DT*(avg_rate1b[i]); //2^3S

double B2i = 1.0 + 2.0*alpha + 0.5*DT*(avg_rate2b[i]); //2^1S

double D00i = -0.5*DT*(avg_rate1b[i]); //1^1S <- 2^3S

double D01i = -0.5*DT*(avg_rate2b[i]); //1^1S <- 2^1S

double D1i = -0.5*DT*(avg_rate1f[i]); //2^3S

double D2i = -0.5*DT*(avg_rate2f[i]); //2^1S

if (i==0){

// A.insert(i,i) = B0i - Rk*Ci; A.insert(i, i+1) = Ai + Ci;

A.insert(i,i) = 1.0; A.insert(i, i+1) = 0.0;

A.insert(i+N_G,i+N_G) = B1i - Rk*Ci; A.insert(i+N_G, i+N_G+1) = Ai + Ci;

A.insert(i+2*N_G,i+2*N_G) = B2i - Rk*Ci; A.insert(i+2*N_G, i+2*N_G+1) = Ai + Ci;

// A.insert(i, i+N_G) = D00i; A.insert(i, i+2*N_G) = D01i;

A.insert(i, i+N_G) = 0.0; A.insert(i, i+2*N_G) = 0.0;

A.insert(i+N_G, i) = D1i; A.insert(i+N_G, i+2*N_G) = 0.0;

A.insert(i+2*N_G, i) = D2i; A.insert(i+2*N_G, i + N_G) = 0.0;

}

else if (i == N_G-1){

// A.insert(i,i-1) = Ai + Ci; A.insert(i,i) = B0i - Rk*Ai;

A.insert(i,i-1) = 0.0; A.insert(i,i) = 1.0;

A.insert(i+N_G,i+N_G-1) = Ai + Ci; A.insert(i+N_G,i+N_G) = B1i - Rk*Ai;

A.insert(i+2*N_G,i+2*N_G-1) = Ai + Ci; A.insert(i+2*N_G,i+2*N_G) = B2i - Rk*Ai;

A.insert(i, i+N_G) = D00i; A.insert(i, i+2*N_G) = D01i;

A.insert(i+N_G, i) = D1i; A.insert(i+N_G, i+2*N_G) = 0.0;

A.insert(i+2*N_G, i) = D2i; A.insert(i+2*N_G, i + N_G) = 0.0;

}

else {

// A.insert(i, i-1) = Ci; A.insert(i, i) = B0i; A.insert(i, i+1) = Ai;

A.insert(i, i-1) = 0.0; A.insert(i, i) = 1.0; A.insert(i, i+1) = 0.0;

A.insert(i+N_G, i+N_G-1) = Ci; A.insert(i+N_G, i+N_G) = B1i; A.insert(i+N_G, i+N_G+1) = Ai;

A.insert(i+2*N_G, i+2*N_G-1) = Ci; A.insert(i+2*N_G, i+2*N_G) = B2i; A.insert(i+2*N_G, i+2*N_G+1) = Ai;

A.insert(i, i+N_G) = D00i; A.insert(i, i+2*N_G) = D01i;

A.insert(i+N_G, i) = D1i; A.insert(i+N_G, i+2*N_G) = 0.0;

A.insert(i+2*N_G, i) = D2i; A.insert(i+2*N_G, i+N_G) = 0.0;

}

for (int i = 0; i < N_G; i++){

double Ai = alpha;

double Ci = alpha;

double B0i = 1.0 - 2.0*alpha - 0.5*(avg_rate1f[i] + avg_rate2f[i]); //1^1S

double B1i = 1.0 - 2.0*alpha - 0.5*(avg_rate1b[i]); //2^3S

double B2i = 1.0 - 2.0*alpha - 0.5*(avg_rate2b[i]); //2^1S

double D00i = 0.5*DT*(avg_rate1b[i]); //1^1S <- 2^3S

double D01i = 0.5*DT*(avg_rate2b[i]); //1^1S <- 2^1S

double D1i = 0.5*DT*(avg_rate1f[i]); //2^3S

double D2i = 0.5*DT*(avg_rate2f[i]); //2^1S

if (i==0){

B.insert(i,i) = B0i - Rk*Ci; B.insert(i, i+1) = Ai + Ci;

B.insert(i+N_G,i+N_G) = B1i - Rk*Ci; B.insert(i+N_G, i+N_G+1) = Ai + Ci;

B.insert(i+2*N_G,i+2*N_G) = B2i - Rk*Ci; B.insert(i+2*N_G, i+2*N_G+1) = Ai + Ci;

B.insert(i, i+N_G) = D00i; B.insert(i, i+2*N_G) = D01i;

B.insert(i+N_G, i) = D1i; B.insert(i+N_G, i+2*N_G) = 0.0;

B.insert(i+2*N_G, i) = D2i; B.insert(i+2*N_G, i + N_G) = 0.0;

}

else if (i == N_G-1){

B.insert(i,i-1) = Ai + Ci; B.insert(i,i) = B0i - Rk*Ai;

B.insert(i+N_G,i+N_G-1) = Ai + Ci; B.insert(i+N_G,i+N_G) = B1i - Rk*Ai;

B.insert(i+2*N_G,i+2*N_G-1) = Ai + Ci; B.insert(i+2*N_G,i+2*N_G) = B2i - Rk*Ai;

B.insert(i, i+N_G) = D00i; B.insert(i, i+2*N_G) = D01i;

B.insert(i+N_G, i) = D1i; B.insert(i+N_G, i+2*N_G) = 0.0;

B.insert(i+2*N_G, i) = D2i; B.insert(i+2*N_G, i + N_G) = 0.0;

}

else {

B.insert(i, i-1) = Ci; B.insert(i, i) = B0i; B.insert(i, i+1) = Ai;

B.insert(i+N_G, i+N_G-1) = Ci; B.insert(i+N_G, i+N_G) = B1i; B.insert(i+N_G, i+N_G+1) = Ai;

B.insert(i+2*N_G, i+2*N_G-1) = Ci; B.insert(i+2*N_G, i+2*N_G) = B2i; B.insert(i+2*N_G, i+2*N_G+1) = Ai;

B.insert(i, i+N_G) = D00i; B.insert(i, i+2*N_G) = D01i;

B.insert(i+N_G, i) = D1i; B.insert(i+N_G, i+2*N_G) = 0.0;

B.insert(i+2*N_G, i) = D2i; B.insert(i+2*N_G, i+N_G) = 0.0;

}

A.makeCompressed(); B.makeCompressed();

Eigen::SparseLU<SpMat> solver;

solver.analyzePattern(A);

solver.factorize(A);

if(solver.info() != Eigen::Success)

throw std::runtime_error("LU factorization failed");

//initial conditions:

for (int i = 0; i < N_G; i++){

sol_old(i) = GAS_DENSITY;

sol_old(i+N_G) = levels[1].density[i];

sol_old(i+2*N_G) = levels[2].density[i];

}

bool converged = false;

double accuracy = 1.0E-10;

std::cout << "Timestep chosen: " << DT_E << "\n";

//main loop:

for (int n = 0; n < 1000000; n++){

RHS = B * sol_old;

for (int p = 0; p < N_G; ++p) {

RHS[p] = GAS_DENSITY;

}

sol_new = solver.solve(RHS);

if(solver.info() != Eigen::Success)

throw std::runtime_error("Solve failed at step " + std::to_string(n));

//checking convergence

double rel = 0.0;

for (int i = 0; i < 3*N_G; ++i) {

double abs = fabs(sol_new(i) - sol_old(i));

rel = std::max(rel, abs/sol_old(i));

}

if ((rel < accuracy)) {

converged = true;

std::cout << "Steady-state reached after " << n << " iterations.\n";

std::cout << "Relative changes are: " << rel << "\n";

break;

}

sol_old.swap(sol_new);

}

for (int i = 0; i < N_G; i++){

levels[1].density[i] = sol_new(i+N_G); levels[2].density[i] = sol_new(i+2*N_G);

}

void do_one_cycle (std::vector<level>& levels, std::vector<Process>& eProc, const std::vector<std::vector<float>>& sigma_temp){

const double DV = ELECTRODE_AREA * DX;

const double FACTOR_W = WEIGHT / DV;

const double FACTOR_E = DT_E / E_MASS * E_CHARGE;

const double FACTOR_I = DT_I / HE_MASS * E_CHARGE;

const double MIN_X_EEPF = 0.45 * L; // min. position for EEPF collection

const double MAX_X_EEPF = 0.55 * L; // max. position for EEPF collection

int k, t, p, energy_index, velocity_index;

double g, g_sqr, gx, gy, gz, vx_a, vy_a, vz_a, e_x, energy, nu, p_coll, v_sqr, velocity;

double mean_v, c0, c1, c2, rate;

bool out;

xvector rho;

int t_index;

nu_avg = 0.0;

for (t=0; t<N_T; t++){ // the RF period is divided into N_T equal time intervals (time step DT_E)

Time += DT_E; // update of the total simulated time

t_index = t / N_BIN; // index for XT distributions

// step 1: compute densities at grid points

for(p=0; p<N_G; p++) e_density[p] = 0; // electron density - computed in every time step

for(k=0; k<N_e; k++){

if (p < 0) p = 0;

else if (p > N_G - 2) p = N_G - 2;

c0 = x_e[k] * INV_DX;

p = int(c0);

e_density[p] += (p + 1 - c0) * FACTOR_W;

e_density[p+1] += (c0 - p) * FACTOR_W;

}

e_density[0] *= 2.0; // double at the edge bcs working with half-domain (no left/right neighbour)

e_density[N_G-1] *= 2.0;

for(p=0; p<N_G; p++) cumul_e_density[p] += e_density[p];

if ((t % N_SUB) == 0) { // ion density - computed in every N_SUB-th time steps (subcycling)

for(p=0; p<N_G; p++) i_density[p] = 0;

for(k=0; k<N_i; k++){

c0 = x_i[k] * INV_DX;

p = int(c0);

i_density[p] += (p + 1 - c0) * FACTOR_W;

i_density[p+1] += (c0 - p) * FACTOR_W;

}

i_density[0] *= 2.0; // double at the edge bcs working with half-domain (no left/right neighbour)

i_density[N_G-1] *= 2.0;

}

for(p=0; p<N_G; p++) cumul_i_density[p] += i_density[p];

// step 2: solve Poisson equation

for(p=0; p<N_G; p++) rho[p] = E_CHARGE * (i_density[p] - e_density[p]); // get charge density

solve_Poisson(rho,Time); // compute potential and electric field

// steps 3 & 4: move particles according to electric field interpolated to particle positions

for(k=0; k<N_e; k++){ // move all electrons in every time step

c0 = x_e[k] * INV_DX;

p = int(c0);

c1 = p + 1.0 - c0;

c2 = c0 - p;

e_x = c1 * efield[p] + c2 * efield[p+1];

if (measurement_mode) {

// measurements: 'x' and 'v' are needed at the same time, i.e. old 'x' and mean 'v'

mean_v = vx_e[k] - 0.5 * e_x * FACTOR_E;

counter_e_xt[p][t_index] += c1;

counter_e_xt[p+1][t_index] += c2;

ue_xt[p][t_index] += c1 * mean_v;

ue_xt[p+1][t_index] += c2 * mean_v;

v_sqr = mean_v * mean_v + vy_e[k] * vy_e[k] + vz_e[k] * vz_e[k];

energy = 0.5 * E_MASS * v_sqr / EV_TO_J;

meanee_xt[p][t_index] += c1 * energy;

meanee_xt[p+1][t_index] += c2 * energy;

energy_index = min( int(energy / DE_CS + 0.5), CS_RANGES-1);

velocity = sqrt(v_sqr);

rate = sigma[E_ION][energy_index] * velocity * DT_E * GAS_DENSITY;

ioniz_rate_xt[p][t_index] += c1 * rate;

ioniz_rate_xt[p+1][t_index] += c2 * rate;

// measure EEPF in the center

if ((MIN_X_EEPF < x_e[k]) && (x_e[k] < MAX_X_EEPF)){

energy_index = (int)(energy / DE_EEPF);

if (energy_index < N_EEPF) {eepf[energy_index] += 1.0;}

mean_energy_accu_center += energy;

mean_energy_counter_center++;

}

// update velocity and position

vx_e[k] -= e_x * FACTOR_E;

x_e[k] += vx_e[k] * DT_E;

}

if ((t % N_SUB) == 0) { // move all ions in every N_SUB-th time steps (subcycling)

for(k=0; k<N_i; k++){

c0 = x_i[k] * INV_DX;

p = int(c0);

c1 = p + 1 - c0;

c2 = c0 - p;

e_x = c1 * efield[p] + c2 * efield[p+1];

if (measurement_mode) {

// measurements: 'x' and 'v' are needed at the same time, i.e. old 'x' and mean 'v'

mean_v = vx_i[k] + 0.5 * e_x * FACTOR_I;

counter_i_xt[p][t_index] += c1;

counter_i_xt[p+1][t_index] += c2;

ui_xt[p][t_index] += c1 * mean_v;

ui_xt[p+1][t_index] += c2 * mean_v;

v_sqr = mean_v * mean_v + vy_i[k] * vy_i[k] + vz_i[k] * vz_i[k];

energy = 0.5 * HE_MASS * v_sqr / EV_TO_J;

meanei_xt[p][t_index] += c1 * energy;

meanei_xt[p+1][t_index] += c2 * energy;

velocity = sqrt(v_sqr); // actually, it's a speed: a velocity includes direction. A

int b = int(x_i[k]/DX_PF);

int v_idx = int(vx_i[k] / DV_IVDF);

int e_idx = int(energy / DE_IEPF);

if (b < 0 || b >= N_sp) continue;

if (abs(v_idx) > N_IVDF) std::cerr << "Make IVDF larger!" << "\n";

else if (v_idx >= 0 && v_idx < N_IVDF) ivdf_pos[b][v_idx] += 1.0;

else if (v_idx < 0 && -v_idx < N_IVDF) ivdf_neg[b][-v_idx] += 1.0;

if (e_idx < N_IEPF) iepf[b][e_idx] += 1.0;

else std::cerr << "Make IEDF larger!" << "\n";

}

// update velocity and position and accumulate absorbed energy

vx_i[k] += e_x * FACTOR_I;

x_i[k] += vx_i[k] * DT_I;

}

// step 5: check boundaries -- changed by Artem

k = 0;

while(k < N_e) { // check boundaries for all electrons in every time step

out = false;

if ((x_e[k] < 0) || (x_e[k] > L)) {out = true;} // the electron is out at the electrodes

if (out) {

if (R01(MTgen)< r) { // reflecting elastically with r probability

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

vx_e[k] *= -1.0; // elastic reflection

}

else { // remove the electron

x_e[k] < 0 ? N_e_abs_pow++ : N_e_abs_gnd++; // absorption calculation at the powered/grounded (left/right) electrode

x_e [k] = x_e [N_e-1]; // pushing last element on a vacant place

vx_e[k] = vx_e[N_e-1];

vy_e[k] = vy_e[N_e-1];

vz_e[k] = vz_e[N_e-1];

N_e--;

}

} else k++;

}

if ((t % N_SUB) == 0) { // check boundaries for all ions in every N_SUB-th time steps (subcycling)

k = 0;

while(k < N_i) {

out = false;

if (x_i[k] < 0) { // the ion is out at the powered electrode

N_i_abs_pow++;

out = true;

v_sqr = vx_i[k] * vx_i[k] + vy_i[k] * vy_i[k] + vz_i[k] * vz_i[k];

energy = 0.5 * HE_MASS * v_sqr/ EV_TO_J;

energy_index = (int)(energy / DE_IFED);

if (energy_index < N_IFED) {ifed_pow[energy_index]++;} // save IFED at the powered electrode

if (R01(MTgen)< gamma_i) {

if (N_e < MAX_N_P){ // adding a new electron on the left electrode

N_e++;

double u = R01(MTgen); // with energy sampled from Maxwellian distribution

vx_e[N_e-1] = sqrt(-2.0*kT_over_m*log(u)); // half-maxwellian in x direction

vy_e[N_e-1] = RMB_e(MTgen); vz_e[N_e-1] = RMB_e(MTgen);

x_e[N_e-1] = 0.0;

if (std::isnan(vx_e[N_e-1])) {

std::cerr << "Error: SEE electron has Vx NAN!";

vx_e[N_e-1] = RMB_e(MTgen);

}

else{

std::cerr << "Impossible to emit a new SE at the powered electrode" <<"\n";

}

if (x_i[k] > L) { // the ion is out at the grounded electrode

N_i_abs_gnd++;

out = true;

v_sqr = vx_i[k] * vx_i[k] + vy_i[k] * vy_i[k] + vz_i[k] * vz_i[k];

energy = 0.5 * HE_MASS * v_sqr / EV_TO_J;

energy_index = (int)(energy / DE_IFED);

if (energy_index < N_IFED) {ifed_gnd[energy_index]++;} // save IFED at the grounded electrode

if (R01(MTgen)< gamma_i) {

if (N_e < MAX_N_P){ // adding a new electron on the left electrode

N_e++; // with energy sampled from Maxwellian distribution

double u = R01(MTgen);

vx_e[N_e-1] = sqrt(-2.0*kT_over_m*log(u)); // half-maxwellian in x direction

vy_e[N_e-1] = RMB_e(MTgen); vz_e[N_e-1] = RMB_e(MTgen);

x_e[N_e-1] = L;

}

else{

std::cerr << "Impossible to emit a new SE at the grounded electrode" <<"\n";

}

if (out) { // delete the ion, if out

x_i [k] = x_i [N_i-1];

vx_i[k] = vx_i[N_i-1];

vy_i[k] = vy_i[N_i-1];

vz_i[k] = vz_i[N_i-1];

N_i--;

} else k++;

}

// step 6: collisions

for (k=0; k<N_e; k++){ // checking for occurrence of a collision for all electrons in every time step

v_sqr = vx_e[k] * vx_e[k] + vy_e[k] * vy_e[k] + vz_e[k] * vz_e[k];

velocity = sqrt(v_sqr);

energy = 0.5 * E_MASS * v_sqr / EV_TO_J;

energy_index = min( int(energy / DE_CS + 0.5), CS_RANGES-1);

int e_crdnt = std::max(0, std::min(N_G-1, static_cast<int>(x_e[k] * INV_DX)));

double frequency_total = 0.0; //collision frequency sigma*density*velocity for this particular electron

double probability_total = 0.0; //collision probability sigma*density for this particular electron

// Artem - adding superelastic impact on total collisoinal probability//

for (int i = 0; i < eProc.size(); i++){

const Process& ch = eProc[i];

double dens = (ch.kind == ProcKind::SUPERELASTIC) ? levels[ch.init_Lvl].density[e_crdnt] : GAS_DENSITY;

probability_total += sigma_temp[i][energy_index] * dens;

}

frequency_total = probability_total * velocity;

p_coll = 1.0 - exp(- frequency_total * DT_E); // collision probability for electrons

if (R01(MTgen) < p_coll) { // electron collision takes place

collision_electron(x_e[k], &vx_e[k], &vy_e[k], &vz_e[k], energy_index, levels, eProc, sigma_temp);

N_e_coll++;

}

if ((t % N_SUB) == 0) { // checking for occurrence of a collision for all ions in every N_SUB-th time steps (subcycling)

for (k=0; k<N_i; k++){

vx_a = RMB_n(MTgen); // pick velocity components of a random target gas atom

vy_a = RMB_n(MTgen);

vz_a = RMB_n(MTgen);

gx = vx_i[k] - vx_a; // compute the relative velocity of the collision partners

gy = vy_i[k] - vy_a;

gz = vz_i[k] - vz_a;

g_sqr = gx * gx + gy * gy + gz * gz;

g = sqrt(g_sqr);

energy = 0.5 * MU_HEHE * g_sqr / EV_TO_J;

energy_index = min( int(energy / DE_CS + 0.5), CS_RANGES-1);

nu = sigma_tot_i[energy_index] * g;

p_coll = 1 - exp(- nu * DT_I); // collision probability for ions

if (R01(MTgen)< p_coll) { // ion collision takes place

collision_ion (&vx_i[k], &vy_i[k], &vz_i[k], &vx_a, &vy_a, &vz_a, energy_index);

N_i_coll++;

}

//step 7: accumulate rates

accumulate_rates(eProc, sigma_temp);

if (measurement_mode) {

// collect 'xt' data from the grid

for (p=0; p<N_G; p++) {

pot_xt [p][t_index] += pot[p];

efield_xt[p][t_index] += efield[p];

ne_xt [p][t_index] += e_density[p];

ni_xt [p][t_index] += i_density[p];

e1_xt [p][t_index] += levels[1].density[p]; // Artem

e2_xt [p][t_index] += levels[2].density[p]; // Artem

}

if ((t % 1000) == 0) printf(" c = %8d t = %8d #e = %8d #i = %8d\n", cycle,t,N_e,N_i);

}

counter++;

// updating denisites each N_avg cycles: --- Artem

if (counter%N_avg == 0) {

// compute average rates over a cycle

average_rates(eProc);

// updating densities

update_excited_dens(eProc, levels);

// reset accumulators

for (int i = 0; i < eProc.size(); i++){

Process& ch = eProc[i];

if (ch.kind == ProcKind::ELASTIC || ch.kind == ProcKind::IONIZATION) continue;

std::fill(ch.cum_rate.begin(), ch.cum_rate.end(), 0.0);

std::fill(ch.avg_rate.begin(), ch.avg_rate.end(), 0.0);

}

fprintf(datafile,"%8d %8d %8d\n",cycle,N_e,N_i);

print_excitation_densities(levels);

}

//---------------------------------------------------------------------//

// save particle coordinates //

//---------------------------------------------------------------------//

void save_particle_data(std::vector<level>& levels, std::vector<Process>& eProc){

double d;

FILE * f;

char fname[80];

strcpy(fname,"picdata.bin");

f = fopen(fname,"wb");

fwrite(&Time,sizeof(double),1,f);

d = (double)(cycles_done);

fwrite(&d,sizeof(double),1,f);

d = (double)(N_e);

fwrite(&d,sizeof(double),1,f);

fwrite(x_e, sizeof(double),N_e,f);

fwrite(vx_e,sizeof(double),N_e,f);

fwrite(vy_e,sizeof(double),N_e,f);

fwrite(vz_e,sizeof(double),N_e,f);

d = (double)(N_i);

fwrite(&d,sizeof(double),1,f);

fwrite(x_i, sizeof(double),N_i,f);

fwrite(vx_i,sizeof(double),N_i,f);

fwrite(vy_i,sizeof(double),N_i,f);

fwrite(vz_i,sizeof(double),N_i,f);

// saving excited state densities and rates - Artem

for (auto &lvl : levels)

fwrite(lvl.density.data(), sizeof(double), N_G, f);

for (auto &proc : eProc){

if (proc.kind == ProcKind::ELASTIC || proc.kind == ProcKind::IONIZATION) continue;

fwrite(proc.cum_rate.data(), sizeof(double), N_G, f);

fwrite(proc.avg_rate.data(), sizeof(double), N_G, f);

}

fclose(f);

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);

}

//---------------------------------------------------------------------//

// load particle coordinates //

//---------------------------------------------------------------------//

void load_particle_data(std::vector<level>& levels, std::vector<Process>& eProc){

double d;

FILE * f;

char fname[80];

strcpy(fname,"picdata.bin");

f = fopen(fname,"rb");

if (f==NULL) {printf(">> eduPIC: ERROR: No particle data file found, try running initial cycle using argument '0'\n"); exit(0); }

fread(&Time,sizeof(double),1,f);

fread(&d,sizeof(double),1,f);

cycles_done = int(d);

fread(&d,sizeof(double),1,f);

N_e = int(d);

fread(x_e, sizeof(double),N_e,f);

fread(vx_e,sizeof(double),N_e,f);

fread(vy_e,sizeof(double),N_e,f);

fread(vz_e,sizeof(double),N_e,f);

fread(&d,sizeof(double),1,f);

N_i = int(d);

fread(x_i, sizeof(double),N_i,f);

fread(vx_i,sizeof(double),N_i,f);

fread(vy_i,sizeof(double),N_i,f);

fread(vz_i,sizeof(double),N_i,f);

// reading excited states densities and rates info -- Artem

for (auto &lvl : levels)

fread(lvl.density.data(), sizeof(double), N_G, f);

for (auto &proc : eProc){

if (proc.kind == ProcKind::ELASTIC || proc.kind == ProcKind::IONIZATION) continue;

fread(proc.cum_rate.data(), sizeof(double), N_G, f);

fread(proc.avg_rate.data(), sizeof(double), N_G, f);

}

fclose(f);

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);

}

//---------------------------------------------------------------------//

// save density data //

//---------------------------------------------------------------------//

void save_density(void){

FILE *f;

double c;

int m;

f = fopen("density.dat","w");

c = 1.0 / (double)(no_of_cycles) / (double)(N_T);

for(m=0; m<N_G; m++){

fprintf(f,"%8.5f %12e %12e\n",m * DX, cumul_e_density[m] * c, cumul_i_density[m] * c);

}

fclose(f);

}

// save Final excited states densities - Artem

void save_final_excited_densities(std::vector<level>& levels) {

FILE *f = fopen("excited_densities.dat", "w");

for(int p=0; p<N_G; p++) {

fprintf(f, "%8.5f %12e %12e\n", p*DX, levels[1].density[p], levels[2].density[p]);

}

fclose(f);

}

//---------------------------------------------------------------------//

// save EEPF data //

//---------------------------------------------------------------------//

void save_eepf(void) {

FILE *f;

int i;

double h,energy;

h = 0.0;

for (i=0; i<N_EEPF; i++) {h += eepf[i];}

h *= DE_EEPF;

f = fopen("eepf.dat","w");

for (i=0; i<N_EEPF; i++) {

energy = (i + 0.5) * DE_EEPF;

fprintf(f,"%e %e\n", energy, eepf[i] / h / sqrt(energy));

}

fclose(f);

}

//---------------------------------------------------------------------//

// save IVDF data -- Artem //

//---------------------------------------------------------------------//

void save_ivdf(void) {

if (mkdir("IVDF", 0755) != 0 && errno != EEXIST) {

perror("mkdir IVDF");

}

for (int b = 0; b < N_sp; b++){

FILE *f;

int i;

double h,velocity;

h = 0.0;

for (i=0; i<N_IVDF; i++){ h += (ivdf_pos[b][i] + ivdf_neg[b][i]);}

h *= DV_IVDF; // this stems from: f(V) = N(V)/(bin_width*sum_{V}(N(V)))

// build filename: “IVDF/ivdf_0.00–0.01cm.dat”, etc.

double x0 = b * DX_PF,

x1 = x0 + DX_PF;

char fname[128];

snprintf(fname, sizeof(fname),

"IVDF/ivdf_%.4f-%.4fm.dat", x0, x1);

f = fopen(fname,"w");

//printing out negative part:

for (i=0; i<N_IVDF; i++) {

velocity = (-1.0)*(i + 0.5) * DV_IVDF;

fprintf(f,"%e %e\n", velocity, ivdf_neg[b][i] / h );

}

//printing out positive part

for (i=0; i<N_IVDF; i++) {

velocity = (i + 0.5) * DV_IVDF;

fprintf(f,"%e %e\n", velocity, ivdf_pos[b][i] / h );

}

fclose(f);

}

//---------------------------------------------------------------------//

// save IEPF data -- Artem //

//---------------------------------------------------------------------//

void save_iepf(void) {

if (mkdir("IEPF", 0755) != 0 && errno != EEXIST) { perror("mkdir IEPF");}

for (int b = 0; b < N_sp; b++){

FILE *f;

int i;

double h,energy;

h = 0.0;

for (i=0; i<N_IEPF; i++){ h += iepf[b][i];}

h *= DE_IEPF;

if (h < 1.0E-12) continue;

// build filename: “IEPF/iepf_0.00–0.01cm.dat”, etc.

double x0 = b * DX_PF,

x1 = x0 + DX_PF;

char fname[128];

snprintf(fname, sizeof(fname), "IEPF/iepf_%.4f-%.4fm.dat", x0, x1);

//printing out part

f = fopen(fname,"w");

for (i=0; i<N_IEPF; i++) {

energy = (i + 0.5) * DE_IEPF;

fprintf(f,"%e %e\n", energy, iepf[b][i] / h );

}

fclose(f);

}

//---------------------------------------------------------------------//

// save IFED data //

//---------------------------------------------------------------------//

void save_ifed(void) {

FILE *f;

int i;

double h_pow,h_gnd,energy;

h_pow = 0.0;

h_gnd = 0.0;

for (i=0; i<N_IFED; i++) {h_pow += ifed_pow[i]; h_gnd += ifed_gnd[i];}

h_pow *= DE_IFED;

h_gnd *= DE_IFED;

mean_i_energy_pow = 0.0;

mean_i_energy_gnd = 0.0;

f = fopen("ifed.dat","w");

for (i=0; i<N_IFED; i++) {

energy = (i + 0.5) * DE_IFED;

fprintf(f,"%6.2f %10.6f %10.6f\n", energy, (double)(ifed_pow[i])/h_pow, (double)(ifed_gnd[i])/h_gnd);

mean_i_energy_pow += energy * (double)(ifed_pow[i]) / h_pow;

mean_i_energy_gnd += energy * (double)(ifed_gnd[i]) / h_gnd;

}

fclose(f);

}

//--------------------------------------------------------------------//

// save XT data //

//--------------------------------------------------------------------//

void save_xt_1(xt_distr distr, char *fname) {

FILE *f;

int i, j;

f = fopen(fname,"w");

for (i=0; i<N_G; i++){

for (j=0; j<N_XT; j++){

fprintf(f,"%e ", distr[i][j]);

}

fprintf(f,"\n");

}

fclose(f);

}

void norm_all_xt(void){

double f1, f2;

int i, j;

// normalize all XT data

f1 = (double)(N_XT) / (double)(no_of_cycles * N_T);

f2 = WEIGHT / (ELECTRODE_AREA * DX) / (no_of_cycles * (PERIOD / (double)(N_XT)));

for (i=0; i<N_G; i++){

for (j=0; j<N_XT; j++){

pot_xt[i][j] *= f1;

efield_xt[i][j] *= f1;

ne_xt[i][j] *= f1;

ni_xt[i][j] *= f1;

e1_xt[i][j] *= f1; // Artem

e2_xt[i][j] *= f1; // Artem

if (counter_e_xt[i][j] > 0) {

ue_xt[i][j] = ue_xt[i][j] / counter_e_xt[i][j];

je_xt[i][j] = -ue_xt[i][j] * ne_xt[i][j] * E_CHARGE;

meanee_xt[i][j] = meanee_xt[i][j] / counter_e_xt[i][j];

ioniz_rate_xt[i][j] *= f2;

} else {

ue_xt[i][j] = 0.0;

je_xt[i][j] = 0.0;

meanee_xt[i][j] = 0.0;

ioniz_rate_xt[i][j] = 0.0;

}

if (counter_i_xt[i][j] > 0) {

ui_xt[i][j] = ui_xt[i][j] / counter_i_xt[i][j];

ji_xt[i][j] = ui_xt[i][j] * ni_xt[i][j] * E_CHARGE;

meanei_xt[i][j] = meanei_xt[i][j] / counter_i_xt[i][j];

} else {

ui_xt[i][j] = 0.0;

ji_xt[i][j] = 0.0;

meanei_xt[i][j] = 0.0;

}

powere_xt[i][j] = je_xt[i][j] * efield_xt[i][j];

poweri_xt[i][j] = ji_xt[i][j] * efield_xt[i][j];

}

void save_all_xt(void){

char fname[128];

if (mkdir("XT_dat", 0755) != 0 && errno != EEXIST) {

perror("mkdir XT_dat");

exit(1);

}

// POT

// make sure the folder exists (ignore EEXIST)

if (mkdir("XT_dat/XT_pot", 0755) != 0 && errno != EEXIST) {

perror("mkdir XT_pot");

}

// build the numbered filename

snprintf(fname, sizeof(fname),

"XT_dat/XT_pot/pot_xt_%04d.dat",

cycles_done);

save_xt_1(pot_xt, fname);

// NE

if (mkdir("XT_dat/XT_ne", 0755) != 0 && errno != EEXIST) {

perror("mkdir XT_ne");

}

snprintf(fname, sizeof(fname),

"XT_dat/XT_ne/ne_xt_%04d.dat",

cycles_done);

save_xt_1(ne_xt, fname);

// NI

if (mkdir("XT_dat/XT_ni", 0755) != 0 && errno != EEXIST) {

perror("mkdir XT_ni");

}

snprintf(fname, sizeof(fname),

"XT_dat/XT_ni/ni_xt_%04d.dat",

cycles_done);

save_xt_1(ni_xt, fname);

// e1

if (mkdir("XT_dat/XT_e1", 0755) != 0 && errno != EEXIST) {

perror("mkdir XT_e1");

}

snprintf(fname, sizeof(fname),

"XT_dat/XT_e1/e1_xt_%04d.dat",

cycles_done);

save_xt_1(e1_xt, fname);

// e2

if (mkdir("XT_dat/XT_e2", 0755) != 0 && errno != EEXIST) {

perror("mkdir XT_e2");

}

snprintf(fname, sizeof(fname),

"XT_dat/XT_e2/e2_xt_%04d.dat",

cycles_done);

save_xt_1(e2_xt, fname);

// E-Field

if (mkdir("XT_dat/XT_efield", 0755) != 0 && errno != EEXIST) {

perror("mkdir XT_efield");

}

snprintf(fname, sizeof(fname),

"XT_dat/XT_efield/XT_efield_%04d.dat",

cycles_done);

save_xt_1(efield_xt, fname);

// Ioniz rate

if (mkdir("XT_dat/XT_ioniz_rate", 0755) != 0 && errno != EEXIST) {

perror("mkdir XT_ioniz_rate");

}

snprintf(fname, sizeof(fname),

"XT_dat/XT_ioniz_rate/XT_ioniz_rate_%04d.dat",

cycles_done);

save_xt_1(ioniz_rate_xt, fname);

// j-e

if (mkdir("XT_dat/XT_je", 0755) != 0 && errno != EEXIST) {

perror("mkdir XT_je");

}

snprintf(fname, sizeof(fname),

"XT_dat/XT_je/XT_je_%04d.dat",

cycles_done);

save_xt_1(je_xt, fname);

// j-i

if (mkdir("XT_dat/XT_ji", 0755) != 0 && errno != EEXIST) {

perror("mkdir XT_ji");

}

snprintf(fname, sizeof(fname),

"XT_dat/XT_ji/XT_ji_%04d.dat",

cycles_done);

save_xt_1(ji_xt, fname);

// mean electron energy

if (mkdir("XT_dat/XT_meanee", 0755) != 0 && errno != EEXIST) {

perror("mkdir XT_meanee");

}

snprintf(fname, sizeof(fname),

"XT_dat/XT_meanee/XT_meanee_%04d.dat",

cycles_done);

save_xt_1(meanee_xt, fname);

// mean ion energy

if (mkdir("XT_dat/XT_meanei", 0755) != 0 && errno != EEXIST) {

perror("mkdir XT_meanei");

}

snprintf(fname, sizeof(fname),

"XT_dat/XT_meanei/XT_meanei_%04d.dat",

cycles_done);

save_xt_1(meanei_xt, fname);

// power e

if (mkdir("XT_dat/XT_powere", 0755) != 0 && errno != EEXIST) {

perror("mkdir XT_powere");

}

snprintf(fname, sizeof(fname),

"XT_dat/XT_powere/XT_powere_%04d.dat",

cycles_done);

save_xt_1(powere_xt, fname);

// power i

if (mkdir("XT_dat/XT_poweri", 0755) != 0 && errno != EEXIST) {

perror("mkdir XT_poweri");

}

snprintf(fname, sizeof(fname),

"XT_dat/XT_poweri/XT_poweri_%04d.dat",

cycles_done);

save_xt_1(poweri_xt, fname);

}

//---------------------------------------------------------------------//

// simulation report including stability and accuracy conditions //

//---------------------------------------------------------------------//

void check_and_save_info(std::vector<level>& levels){

FILE *f;

double plas_freq, meane, kT, debye_length, density, ecoll_freq, icoll_freq, sim_time, e_max, v_max, power_e, power_i, c;

int i,j;

bool conditions_OK;

density = cumul_e_density[N_G / 2] / (double)(no_of_cycles) / (double)(N_T); // e density @ center

plas_freq = E_CHARGE * sqrt(density / EPSILON0 / E_MASS); // e plasma frequency @ center

meane = mean_energy_accu_center / (double)(mean_energy_counter_center); // e mean energy @ center

kT = 2.0 * meane * EV_TO_J / 3.0; // k T_e @ center (approximate)

sim_time = (double)(no_of_cycles) / FREQUENCY; // simulated time

ecoll_freq = (double)(N_e_coll) / sim_time / (double)(N_e); // e collision frequency

icoll_freq = (double)(N_i_coll) / sim_time / (double)(N_i); // ion collision frequency

debye_length = sqrt(EPSILON0 * kT / density) / E_CHARGE; // e Debye length @ center

f = fopen("info.txt","w");

fprintf(f,"########################## eduPIC simulation report ############################\n");

fprintf(f,"Simulation parameters:\n");

fprintf(f,"Gap distance = %12.3e [m]\n", L);

fprintf(f,"# of grid divisions = %12d\n", N_G);

fprintf(f,"Frequency = %12.3e [Hz]\n", FREQUENCY);

fprintf(f,"# of time steps / period = %12d\n", N_T);

fprintf(f,"# of electron / ion time steps = %12d\n", N_SUB);

fprintf(f,"Voltage amplitude = %12.3e [V]\n", VOLTAGE);

fprintf(f,"Pressure (Ar) = %12.3e [Pa]\n", PRESSURE);

fprintf(f,"Temperature = %12.3e [K]\n", T_neutral);

fprintf(f,"Superparticle weight = %12.3e\n", WEIGHT);

fprintf(f,"# of simulation cycles in this run = %12d\n", no_of_cycles);

fprintf(f,"--------------------------------------------------------------------------------\n");

fprintf(f,"Plasma characteristics:\n");

fprintf(f,"Electron density @ center = %12.3e [m^{-3}]\n", density);

fprintf(f,"Plasma frequency @ center = %12.3e [rad/s]\n", plas_freq);

fprintf(f,"Debye length @ center = %12.3e [m]\n", debye_length);

fprintf(f,"Electron collision frequency = %12.3e [1/s]\n", ecoll_freq);

fprintf(f,"Ion collision frequency = %12.3e [1/s]\n", icoll_freq);

fprintf(f,"--------------------------------------------------------------------------------\n");

fprintf(f,"Stability and accuracy conditions:\n");

conditions_OK = true;

c = plas_freq * DT_E;

fprintf(f,"Plasma frequency @ center * DT_E = %12.3f (OK if less than 0.20)\n", c);

if (c > 0.2) {conditions_OK = false;}

c = DX / debye_length;

fprintf(f,"DX / Debye length @ center = %12.3f (OK if less than 1.00)\n", c);

if (c > 1.0) {conditions_OK = false;}

c = max_electron_coll_freq() * DT_E;

fprintf(f,"Max. electron coll. frequency * DT_E = %12.3f (OK if less than 0.05)\n", c);

if (c > 0.05) {conditions_OK = false;}

c = max_ion_coll_freq() * DT_I;

fprintf(f,"Max. ion coll. frequency * DT_I = %12.3f (OK if less than 0.05)\n", c);

if (c > 0.05) {conditions_OK = false;}

if (conditions_OK == false){

fprintf(f,"--------------------------------------------------------------------------------\n");

fprintf(f,"** STABILITY AND ACCURACY CONDITION(S) VIOLATED - REFINE SIMULATION SETTINGS! **\n");

fprintf(f,"--------------------------------------------------------------------------------\n");

fclose(f);

printf(">> eduPIC: ERROR: STABILITY AND ACCURACY CONDITION(S) VIOLATED!\n");

printf(">> eduPIC: for details see 'info.txt' and refine simulation settings!\n");

}

else

{

// calculate maximum energy for which the Courant-Friedrichs-Levy condition holds:

v_max = DX / DT_E;

e_max = 0.5 * E_MASS * v_max * v_max / EV_TO_J;

fprintf(f,"Max e- energy for CFL condition = %12.3f [eV]\n", e_max);

fprintf(f,"Check EEPF to ensure that CFL is fulfilled for the majority of the electrons!\n");

fprintf(f,"--------------------------------------------------------------------------------\n");

// saving of the following data is done here as some of the further lines need data

// that are computed / normalized in these functions

printf(">> eduPIC: saving diagnostics data\n");

save_density();

save_final_excited_densities(levels); // Artem

save_eepf();

save_ivdf();

save_iepf();

save_ifed();

norm_all_xt();

save_all_xt();

fprintf(f,"Particle characteristics at the electrodes:\n");

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));

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));

fprintf(f,"Mean ion energy at powered electrode = %12.3e [eV]\n", mean_i_energy_pow);

fprintf(f,"Mean ion energy at grounded electrode = %12.3e [eV]\n", mean_i_energy_gnd);

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));

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));

fprintf(f,"--------------------------------------------------------------------------------\n");

// calculate spatially and temporally averaged power absorption by the electrons and ions

power_e = 0.0;

power_i = 0.0;

for (i=0; i<N_G; i++){

for (j=0; j<N_XT; j++){

power_e += powere_xt[i][j];

power_i += poweri_xt[i][j];

}

power_e /= (double)(N_XT * N_G);

power_i /= (double)(N_XT * N_G);

fprintf(f,"Absorbed power calculated as <j*E>:\n");

fprintf(f,"Electron power density (average) = %12.3e [W m^{-3}]\n", power_e);

fprintf(f,"Ion power density (average) = %12.3e [W m^{-3}]\n", power_i);

fprintf(f,"Total power density(average) = %12.3e [W m^{-3}]\n", power_e + power_i);

fprintf(f,"--------------------------------------------------------------------------------\n");

fclose(f);

}

//------------------------------------------------------------------------------------------//

// main //

// command line arguments: //

// [1]: number of cycles (0 for init) //

// [2]: "m" turns on data collection and saving //

//------------------------------------------------------------------------------------------//

int main (int argc, char *argv[]){

if (argc == 1) {

printf(">> eduPIC: error = need starting_cycle argument\n");

return 1;

} else {

strcpy(st0,argv[1]);

arg1 = atol(st0);

if (argc > 2) {

if (strcmp (argv[2],"m") == 0){

measurement_mode = true; // measurements will be done

} else {

measurement_mode = false;

}

if (measurement_mode) {

printf(">> eduPIC: measurement mode: on\n");

} else {

printf(">> eduPIC: measurement mode: off\n");

}

//---------------------reading excitation levels information--------------------//

std::ifstream levels_data("data/levels.dat");

if (!levels_data) {

std::cerr << "Failed to open file\n";

return 1;

}

std::string header1;

std::getline(levels_data, header1);

std::vector<level> levels;

level temp1;

while (levels_data >> temp1.name >> temp1.th_energy >> temp1.weight) {

levels.push_back(temp1);

}

// ----------------------------------------------------------------------------//

std::unordered_map<std::string,int> lvlIndex; // a map that connects levels names and their indices

for (int i=0;i<levels.size(); ++i) lvlIndex[levels[i].name] = i;

for (auto &lvl : levels) {

lvl.density.resize(N_G);

}

// --------- loading transitions information ----------------------------------//

std::ifstream transitions_data("data/transitions.dat");

if (!transitions_data) {

std::cerr << "Failed to open file\n";

return 1;

}

std::string header2;

std::getline(transitions_data, header2);

std::vector<transition> transitions;

transition temp2;

while (transitions_data >> temp2.lvl_init >> temp2.lvl_final >> temp2.type >> temp2.thresh_energy >> temp2.A1 >> temp2.A2

>> temp2.A3 >> temp2.A4 >> temp2.A5 >> temp2.A6 >> temp2.is_superelastic) {

transitions.push_back(temp2);

}

// -----------------------------------------------------------------------------//

std::vector<Process> eProc;

eProc.push_back({ProcKind::ELASTIC, "1^1S -> 1^1S", static_cast<int>(eProc.size()), lvlIndex["1^1S"], lvlIndex["1^1S"], 0.0});

eProc.push_back({ProcKind::IONIZATION, "1^1S -> He+" , static_cast<int>(eProc.size()), lvlIndex["1^1S"], lvlIndex["He+"], E_ION_TH});

for (auto const& transition : transitions) {

Process p;

p.kind = ((transition.is_superelastic != 0) ? ProcKind::SUPERELASTIC : ProcKind::EXCITATION);

p.name = transition.lvl_init + "→" + transition.lvl_final;

p.csSlot = static_cast<int>(eProc.size());

p.init_Lvl = lvlIndex[transition.lvl_init];

p.final_Lvl = lvlIndex[transition.lvl_final];

p.threshold_eV = transition.thresh_energy;

p.A1 = transition.A1;

p.A2 = transition.A2;

p.A3 = transition.A3;

p.A4 = transition.A4;

p.A5 = transition.A5;

p.A6 = transition.A6;

p.type = transition.type;

eProc.push_back(p);

}

for (auto &pr : eProc) {

pr.cum_rate.resize(N_G);

pr.avg_rate.resize(N_G);

}

std::cout << eProc[4].A1 << " " << eProc[4].A2 << " " << eProc[4].A3 << " " << eProc[4].A4 << " " << eProc[4].A5 << " " << eProc[4].A6 << " "

<< eProc[4].csSlot << " " << eProc[4].init_Lvl << " " << eProc[4].final_Lvl << " " << eProc[4].name << " "

<< eProc[4].threshold_eV << " " << eProc[4].type << "\n";

std::cout << eProc[5].A1 << " " << eProc[5].A2 << " " << eProc[5].A3 << " " << eProc[5].A4 << " " << eProc[5].A5 << " " << eProc[5].A6 << " "

<< eProc[5].csSlot << " " << eProc[5].init_Lvl << " " << eProc[5].final_Lvl << " " << eProc[5].name << " "

<< eProc[5].threshold_eV << " " << eProc[5].type << "\n";

const int P = static_cast<int>(eProc.size());

std::vector<std::vector<float>> sigma_temp(P, std::vector<float>(CS_RANGES,0.0f));

set_electron_cross_sections_ar(eProc, levels, sigma_temp);

set_ion_cross_sections_ar();

calc_total_cross_sections();

std::ofstream file0("rates.dat"); file0 << std::scientific << std::setprecision(6);

std::ofstream file1("excited_densities.dat"); file1 << std::scientific << std::setprecision(6);

std::ofstream file2("ions_info.dat"); file2 << std::scientific << std::setprecision(6);

std::ofstream CS_file("CS_check.dat");

for (int i = 0; i < CS_RANGES; i++){

CS_file << i*DE_CS ;

for (int p = 0; p < P; p++){

CS_file << " " << sigma_temp[p][i];

}

CS_file << "\n";

}

CS_file.close();

//test_cross_sections(); return 1;

datafile = fopen("conv.dat","a");

if (arg1 == 0) {

if (FILE *file = fopen("picdata.bin", "r")) { fclose(file);

printf(">> eduPIC: Warning: Data from previous calculation are detected.\n");

printf(" To start a new simulation from the beginning, please delete all output files before running ./eduPIC 0\n");

printf(" To continue the existing calculation, please specify the number of cycles to run, e.g. ./eduPIC 100\n");

exit(0);

}

no_of_cycles = 1;

cycle = 1; // init cycle

init(N_INIT, levels); // seed initial electrons & ions

printf(">> eduPIC: running initializing cycle\n");

Time = 0;

do_one_cycle(levels, eProc, sigma_temp);

print_excitation_densities(levels);

for (int i = 0; i < N_G; i++){

double x = i*DX;

file0 << avg_rate1f[i] << " " << avg_rate1b[i] << " " << avg_rate2f[i] << " " << avg_rate2b[i] << "\n";

}

cycles_done = 1;

} else {

no_of_cycles = arg1; // run number of cycles specified in command line

load_particle_data(levels, eProc); // read previous configuration from file

counter = cycles_done % N_avg;

printf(">> eduPIC: running %d cycle(s)\n",no_of_cycles);

for (cycle=cycles_done+1;cycle<=cycles_done+no_of_cycles;cycle++) {

do_one_cycle(levels, eProc, sigma_temp);

}

cycles_done += no_of_cycles;

for (int i = 0; i < N_i; i++){

file2 << x_i[i] << " " << vx_i[i] << " " << vy_i[i] << " " << vz_i[i] << "\n";

}

fclose(datafile);

save_particle_data(levels, eProc);

if (measurement_mode) {

check_and_save_info(levels);

}

printf(">> eduPIC: simulation of %d cycle(s) is completed.\n",no_of_cycles);

file0.close();

file1.close();

file2.close();

}