eduPIC_DRR_v1.4 - stable

#include <cstdio>

#include <cstdlib>

#include <cstring>

#include <cstdbool>

#include <cmath>

#include <ctime>

#include <random>

#include <vector>

#include <string>

#include <fstream>

#include <sstream>

#include <algorithm>

#include <stdexcept>

#include <iostream>

#include <iomanip>

using namespace::std;

// 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 = 10; // number of cycles to average rates --- Artem

int counter = -1; // 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

// 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_EXC_1 = 1; // process identifier: electron/excitation1

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

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

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

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

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

const int E_SUPER_2 = 7; // singlet excitation - 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; xvector avg_electron_density;

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 = 2000; // 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 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 = 20; // 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.3; // 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.7; // probabilty for the electron to reflect at the surface

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

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

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

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

// electron cross sections: A V Phelps & Z Lj Petrovic, PSST 8 R21 (1999) //

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

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

};

void set_electron_cross_sections_ar(void){

int i;

double en,qmel,qexc_1,qexc_2,qion;

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_exc1 ("CS/He_electron_exc1.dat"); // two‐col: E σ_exc1

CSInterpolator cs_exc2 ("CS/He_electron_exc2.dat"); // two‐col: E σ_exc2

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

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

// your energy grid

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

// interpolate

sigma[E_ELA][i] = cs_ela(en);

sigma[E_EXC_1][i] = cs_exc1(en);

sigma[E_EXC_2][i] = cs_exc2(en);

sigma[E_ION][i] = cs_ion(en);

// Superelastic for triplet (E_SUPER_1)

double en_plus_1 = en + E_EXC_TH_1;

int idx1 = en_plus_1 / DE_CS;

if (idx1 < CS_RANGES && en > 0) {

sigma[E_SUPER_1][i] = (1.0 / 3.0) * (en_plus_1 / en) * sigma[E_EXC_1][idx1];

} else {

sigma[E_SUPER_1][i] = 0.0;

}

// Superelastic for singlet (E_SUPER_2)

double en_plus_2 = en + E_EXC_TH_2;

int idx2 = en_plus_2 / DE_CS;

if (idx2 < CS_RANGES && en > 0) {

sigma[E_SUPER_2][i] = (1.0 / 1.0) * (en_plus_2 / en) * sigma[E_EXC_2][idx2];

} else {

sigma[E_SUPER_2][i] = 0.0;

}

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

// // your energy grid

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

// // Superelastic for triplet (E_SUPER_1)

// double en_plus_1 = en + E_EXC_TH_1;

// int idx1 = en_plus_1 / DE_CS;

// if (idx1 < CS_RANGES && en > 0) {

// sigma[E_SUPER_1][i] = (1.0 / 3.0) * (en_plus_1 / en) * sigma[E_EXC_1][idx1];

// } else {

// sigma[E_SUPER_1][i] = 0.0;

// }

// // Superelastic for singlet (E_SUPER_2)

// double en_plus_2 = en + E_EXC_TH_2;

// int idx2 = en_plus_2 / DE_CS;

// if (idx2 < CS_RANGES && en > 0) {

// sigma[E_SUPER_2][i] = (1.0 / 1.0) * (en_plus_2 / en) * sigma[E_EXC_2][idx2];

// } else {

// sigma[E_SUPER_2][i] = 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 //

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

// initialization of excited states distribtuion, assuming Maxwellian-Boltzmann balance -- Artem

std::pair<double, double> init_excited_distr() {

double part_ground = 1.0*exp(-0.0/T_neutral); // partition function for ground state

double part_triplet = 3.0*exp(-E_EXC_TH_1*EV_TO_J/(K_BOLTZMANN*T_neutral)); // partition function for triplet excited state

double part_singlet = 1.0*exp(-E_EXC_TH_2*EV_TO_J/(K_BOLTZMANN*T_neutral)); // partition function for singlet excited state

double part_func_total = part_ground + part_triplet + part_singlet; // total partition function

double n_triplet = ((part_triplet/part_func_total)*GAS_DENSITY); // denisty population of tripet state

double n_singlet = ((part_singlet/part_func_total)*GAS_DENSITY); // density population of singlet state

return {n_triplet, n_singlet};

}

void print_excitation_densities(void) {

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 += e1_density[p]; // triplet state density

total_e2 += e2_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){

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 n_min = 1e-15 * GAS_DENSITY;

auto exc = init_excited_distr();

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

e1_density[p] = exc.first;

e2_density[p] = exc.second;

}

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

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

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

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

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,rnd;

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,e_sc,e_ej;

// - Artem

// Determine cell p where the electron is

double c0 = xe * INV_DX;

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

// Local densities -- Artem

double e1_dens = e1_density[p];

double e2_dens = e2_density[p];

double ground_dens = GAS_DENSITY - e1_dens - e2_dens;

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

// choose the type of collision based on the cross sections

// take into account energy loss in inelastic collisions

// generate scattering and azimuth angles

// in case of ionization handle the 'new' electron

t0 = sigma[E_ELA][eindex] * ground_dens;

t1 = t0 + sigma[E_EXC_1][eindex] * ground_dens;

t2 = t1 + sigma[E_EXC_2][eindex] * ground_dens;

t3 = t2 + sigma[E_ION][eindex] * ground_dens;

t4 = t3 + sigma[E_SUPER_1][eindex] * e1_dens; // account for superelastic triplet - Artem

t5 = t4 + sigma[E_SUPER_2][eindex] * e2_dens; // account for superelastic singlet- Artem

rnd = R01(MTgen);

if (rnd < (t0/t5)){ // elastic scattering

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

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

} else if (rnd < (t1/t5)){ // excitation 1 (triplet)

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

energy = fabs(energy - E_EXC_TH_1 * 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

//add new triplet excited He atom density to this grid point, sample velocities from Maxwellian distribution - Artem

//no need if we calculate rates and densities outside

// e1_density[p] += WEIGHT / (ELECTRODE_AREA * DX);

} else if (rnd < (t2/t5)){ // excitation 2 (singlet)

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

energy = fabs(energy - E_EXC_TH_2 * 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

//add new singet excited He atom, sample velocities from Maxwellian distribution - Artem

//no need if we calculate rates and densities outside

// e2_density[p] += WEIGHT / (ELECTRODE_AREA * DX);

} else if (rnd < (t3/t5)) { // ionization

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

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

// e_ej = 10.0 * tan(R01(MTgen) * atan(energy/EV_TO_J / 20.0)) * EV_TO_J; // energy of the ejected electron

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

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(sqrt(e_sc / energy)); // scattering angle for scattered electron

// chi2 = acos(sqrt(e_ej / energy)); // scattering angle for ejected electrons

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

}

else if (rnd < (t4/t5)) { // accounting for superelastic collisions - triplet - Artem

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

energy = fabs(energy + E_EXC_TH_1 * EV_TO_J); // add energy for deexcitation

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

//substract excited He atom density from the grid point - Artem

//no need if we calculate rates and densities outside

// e1_density[p] = std::max(e1_density[p] - WEIGHT / (ELECTRODE_AREA * DX), 0.0);

} else { // account for supeelastic collision - singlet - Artem

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

energy = fabs(energy + E_EXC_TH_2 * EV_TO_J); // add energy for deexcitation

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

//substract excited He atom density from the grid point - Artem

//no need if we calculate rates and densities outside

// e2_density[p] = std::max(e2_density[p] - WEIGHT / (ELECTRODE_AREA * DX), 0.0);

}

// 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() {

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

sum_electron_density[p_temp] += WEIGHT / Volume;

sum_rate1f[p_temp] += sigma[E_EXC_1][energy_index] * velocity * WEIGHT;

sum_rate2f[p_temp] += sigma[E_EXC_2][energy_index] * velocity * WEIGHT;

sum_rate1b[p_temp] += sigma[E_SUPER_1][energy_index] * velocity * WEIGHT;

sum_rate2b[p_temp] += sigma[E_SUPER_2][energy_index] * velocity * WEIGHT;

}

// averaging the rates each RF cycle

void average_rates() {

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++){

avg_electron_density[p] = sum_electron_density[p] * inv_NT;

}

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

avg_rate1f[p] = sum_rate1f[p] * inv_NT * inv_Volume;

avg_rate1b[p] = sum_rate1b[p] * inv_NT * inv_Volume;

avg_rate2f[p] = sum_rate2f[p] * inv_NT * inv_Volume;

avg_rate2b[p] = sum_rate2b[p] * inv_NT * inv_Volume;

}

void update_excited_dens() {

xvector e1_density_bw;

xvector e2_density_bw;

xvector gas_dens_local; // local density of ground state atoms

double S0, S1, S2; // source terms for diffusion equation

double diff1, diff2; // diffusion terms

const double v_mean = sqrt(8.0*K_BOLTZMANN*T_neutral/(M_PI*HE_MASS)); // mean spead of He excited states, assuming T~Tg

bool converged = false;

double accuracy = 1.0E-4; // demanded accuracy

double alpha = 0.48;

double dt = alpha*DX*DX/Diff_coeff;

int max_iter = 1E8; // maximum iteration

std::cout << Diff_coeff << " " << dt << " " << DX << "\n";

// creating a helping array for backward step and also works with initial condition

std::memcpy(e1_density_bw, e1_density, sizeof(e1_density));

std::memcpy(e2_density_bw, e2_density, sizeof(e2_density));

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

gas_dens_local[i] = GAS_DENSITY;

}

double rcmb_coeff = (-1.0)*gamma_k*v_mean/(2.0*(2.0-gamma_k)*Diff_coeff);

for (int idx = 0; idx < max_iter; idx++) {

// swapping the arrays

std::swap(e1_density_bw, e1_density);

std::swap(e2_density_bw, e2_density);

for (int i = 1; i < N_G-1; i++){

// calculate source terms of diffusion equation, S0 for neutral atoms of He

// S1 for triplet excited state, S2 for singlet excited state

// triplet and singlet are diffusing, neutral atoms are not

S0 = dt*(-(avg_rate1f[i]+avg_rate2f[i])*gas_dens_local[i] + avg_rate1b[i]*e1_density_bw[i] + avg_rate2b[i]*e2_density_bw[i]);

S1 = dt*(avg_rate1f[i]*gas_dens_local[i] - avg_rate1b[i]*e1_density_bw[i]);

S2 = dt*(avg_rate2f[i]*gas_dens_local[i] - avg_rate2b[i]*e2_density_bw[i]);

diff1 = e1_density_bw[i]*(1.0 - 2.0*alpha) + alpha*(e1_density_bw[i-1] + e1_density_bw[i+1]);

diff2 = e2_density_bw[i]*(1.0 - 2.0*alpha) + alpha*(e2_density_bw[i-1] + e2_density_bw[i+1]);

e1_density[i] = e1_density_bw[i] + S1;

e2_density[i] = e2_density_bw[i] + S2;

// gas_dens_local[i] += S0;

// Prevent negative densities

// gas_dens_local[i] = std::max(gas_dens_local[i], 0.0);

e1_density[i] = std::max(e1_density[i], 0.0);

e2_density[i] = std::max(e2_density[i], 0.0);

}

// boundary conditions: now it's zero so no excited states density at the electrodes, yield is zero

e1_density[0] = e1_density[1]/(rcmb_coeff*DX + 1.0);

e1_density[N_G-1] = e1_density[N_G-2]/(rcmb_coeff*DX + 1.0);

e2_density[0] = e2_density[1]/(rcmb_coeff*DX + 1.0);

e2_density[N_G-1] = e2_density[N_G-2]/(rcmb_coeff*DX + 1.0);

//checking convergence

double rel1 = 0.0, rel2 = 0.0;

for (int i = 1; i < N_G-1; ++i) {

double abs1 = fabs(e1_density[i] - e1_density_bw[i]);

double abs2 = fabs(e2_density[i] - e2_density_bw[i]);

double base1 = std::max(e1_density[i], 1e-12);

double base2 = std::max(e2_density[i], 1e-12);

rel1 = std::max(rel1, abs1 / e1_density_bw[i]);

rel2 = std::max(rel2, abs2 / e2_density_bw[i]);

}

if (rel1 < accuracy && rel2 < accuracy) {

converged = true;

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

std::cout << "Relative changes are: " << rel1 << " " << rel2 << "\n";

break;

}

if (!converged) {

std::cerr << "Steady-state not reached after " << max_iter << " iterations.\n";

}

void do_one_cycle (void){

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 = 0.45 * L; // min. position for EEPF collection

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

int k, t, p, energy_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);

double local_neut_dens = GAS_DENSITY - e1_density[p] - e2_density[p];

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

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 < x_e[k]) && (x_e[k] < MAX_X)){

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;

}

// 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++; // with energy sampled from Maxwellian distribution

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

x_e[N_e-1] = 0.0;

}

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

vx_e[N_e-1] = RMB_e(MTgen); 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);

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

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

double sigma_super_e = sigma[E_SUPER_1][energy_index] * e1_density[e_crdnt] + sigma[E_SUPER_2][energy_index] * e2_density[e_crdnt];

double ground_dens_local = GAS_DENSITY - e1_density[e_crdnt] - e2_density[e_crdnt];

double sigma_ground = (sigma[E_ELA][energy_index] + sigma[E_EXC_1][energy_index] +

sigma[E_EXC_2][energy_index] + sigma[E_ION][energy_index]) * ground_dens_local;

nu = (sigma_ground + sigma_super_e) * velocity;

nu_avg += nu;

p_coll = 1 - exp(- nu * 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);

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

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] += e1_density[p]; // Artem

e2_xt [p][t_index] += e2_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();

// updating densities

update_excited_dens();

// reset accumulators

memset(sum_rate1f, 0, sizeof(sum_rate1f));

memset(sum_rate1b, 0, sizeof(sum_rate1b));

memset(sum_rate2f, 0, sizeof(sum_rate2f));

memset(sum_rate2b, 0, sizeof(sum_rate2b));

memset(sum_electron_density, 0, sizeof(sum_electron_density));

}

//calculate nu electron average:

nu_avg /= (N_T*N_e);

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

print_excitation_densities();

}

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

// save particle coordinates //

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

void save_particle_data(){

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 - Artem

fwrite(e1_density, sizeof(double), N_G, f);

fwrite(e2_density, 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(){

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 -- Artem

fread(e1_density, sizeof(double), N_G, f);

fread(e2_density, 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(void) {

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

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

fprintf(f, "%8.5f %12e %12e\n", p*DX, e1_density[p], e2_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 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[80];

strcpy(fname,"pot_xt.dat"); save_xt_1(pot_xt, fname);

strcpy(fname,"efield_xt.dat"); save_xt_1(efield_xt, fname);

strcpy(fname,"ne_xt.dat"); save_xt_1(ne_xt, fname);

strcpy(fname,"ni_xt.dat"); save_xt_1(ni_xt, fname);

strcpy(fname,"je_xt.dat"); save_xt_1(je_xt, fname);

strcpy(fname,"ji_xt.dat"); save_xt_1(ji_xt, fname);

strcpy(fname,"powere_xt.dat"); save_xt_1(powere_xt, fname);

strcpy(fname,"poweri_xt.dat"); save_xt_1(poweri_xt, fname);

strcpy(fname,"meanee_xt.dat"); save_xt_1(meanee_xt, fname);

strcpy(fname,"meanei_xt.dat"); save_xt_1(meanei_xt, fname);

strcpy(fname,"ioniz_xt.dat"); save_xt_1(ioniz_rate_xt, fname);

strcpy(fname,"e1_xt.dat"); save_xt_1(e1_xt, fname);

strcpy(fname,"e2_xt.dat"); save_xt_1(e2_xt, fname);

}

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

// simulation report including stability and accuracy conditions //

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

void check_and_save_info(void){

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(); // Artem

save_eepf();

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[]){

// printf(">> eduPIC: starting...\n");

// printf(">> eduPIC: **************************************************************************\n");

// printf(">> eduPIC: This program comes with ABSOLUTELY NO WARRANTY\n");

// printf(">> eduPIC: This is free software, you are welcome to use, modify and redistribute it\n");

// printf(">> eduPIC: according to the GNU General Public License, https://www.gnu.org/licenses/\n");

// printf(">> eduPIC: **************************************************************************\n");

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

}

set_electron_cross_sections_ar();

set_ion_cross_sections_ar();

calc_total_cross_sections();

auto exc = init_excited_distr();

printf("Excited state population - superparticles!!!\n");

printf("%lf %lf\n", exc.first, exc.second);

printf("Sanity check \n");

fflush(stdout);

// for (int i = 0; i < 10; ++i) {

// double E = i * DE_CS;

// printf("σ_super1(%.3f eV) = %e σ_super2(%.3f eV) = %e\n",

// E, sigma[E_SUPER_1][i],

// E, sigma[E_SUPER_2][i]);

// }

std::ofstream file0("rates.dat"); file0 << std::scientific << std::setprecision(6);

std::ofstream file1("excited_densities.dat"); file1 << std::scientific << std::setprecision(6);

//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); // seed initial electrons & ions

printf(">> eduPIC: running initializing cycle\n");

Time = 0;

do_one_cycle();

print_excitation_densities();

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

file1 << x << " " << e1_density[i] << " " << e2_density[i] << "\n";

}

cycles_done = 1;

} else {

no_of_cycles = arg1; // run number of cycles specified in command line

load_particle_data(); // read previous configuration from file

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

}

cycles_done += no_of_cycles;

}

fclose(datafile);

save_particle_data();

if (measurement_mode) {

check_and_save_info();

}

printf(">> eduPIC: simulation of %d cycle(s) is completed.\n",no_of_cycles);

file0.close();

}