Skip to content

Commit

Permalink
Split ionization specializations into files according to rate used. C…
Browse files Browse the repository at this point in the history 
…orrect bug in BSI
  • Loading branch information
@GianlucaOberreit
GianlucaOberreit committed Sep 15, 2024
1 parent de467bd commit 4253299
Show file tree
Hide file tree
Showing 5 changed files with 137 additions and 138 deletions.
62 changes: 62 additions & 0 deletions src/Ionization/IonizationBSI.cpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,62 @@
#include "IonizationTunnel.h"
#include "IonizationTables.h"

template <>
IonizationTunnel<3>::IonizationTunnel(Params &params, Species *species) : Ionization(params, species)
{
DEBUG("Creating the Tunnel BSI Ionizaton class");

// Ionization potential & quantum numbers (all in atomic units 1 au = 27.2116 eV)
for (int Z = 0; Z < (int)atomic_number_; Z++)
{
DEBUG("Z : " << Z);

Potential[Z] = IonizationTables::ionization_energy(atomic_number_, Z) * eV_to_au;
Azimuthal_quantum_number[Z] = IonizationTables::azimuthal_atomic_number(atomic_number_, Z);

DEBUG("potential: " << Potential[Z] << " Az.q.num: " << Azimuthal_quantum_number[Z]);

double cst = ((double)Z + 1.0) * sqrt(2.0 / Potential[Z]);
alpha_tunnel[Z] = cst - 1.0;
beta_tunnel[Z] = pow(2, alpha_tunnel[Z]) * (8. * Azimuthal_quantum_number[Z] + 4.0) / (cst * tgamma(cst)) *
Potential[Z] * au_to_w0;
gamma_tunnel[Z] = 2.0 * pow(2.0 * Potential[Z], 1.5);
}
DEBUG("Finished Creating the BSI Tunnel Ionizaton class");
}

template <>
double IonizationTunnel<3>::ionizationRate(const int Z, const double E, const double oldZ)
{
auto normal = [this](const int Z, const double E) -> double {
double delta = gamma_tunnel[Z] / E;
return beta_tunnel[Z] * exp(-delta * one_third + alpha_tunnel[Z] * log(delta));
};

auto linear = [this](const int Z, const double E) -> double {
const double ratio_of_IPs = IH / IonizationTables::ionization_energy(atomic_number_, Z);
return au_to_w0 * (0.8 * E * pow(ratio_of_IPs, 0.5));
};
auto quadratic = [this](const int Z, const double E) -> double {
const double ratio_of_IPs = IH / IonizationTables::ionization_energy(atomic_number_, Z);
return au_to_w0 * (2.4 * (pow(E, 2)) * pow(ratio_of_IPs, 2));
};

double BSI_rate_quadratic = quadratic(oldZ+1, E);
double BSI_rate_linear = linear(oldZ+1, E);
double Tunnel_rate = normal(oldZ, E);

if (BSI_rate_quadratic >= BSI_rate_linear)
{
return quadratic(Z, E);
}
else if (std::min(Tunnel_rate, BSI_rate_quadratic) == BSI_rate_quadratic)
{
return linear(Z, E);
}
else
{
return normal(Z, E);
}

}
35 changes: 35 additions & 0 deletions src/Ionization/IonizationFullPPT.cpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,35 @@
#include "IonizationTunnel.h"
#include "IonizationTables.h"

template <>
IonizationTunnel<1>::IonizationTunnel(Params &params, Species *species) : Ionization(params, species)
{
DEBUG("Creating the Tunnel Ionizaton class");

// Ionization potential & quantum numbers (all in atomic units 1 au = 27.2116 eV)
Magnetic_quantum_number.resize(atomic_number_);

for (unsigned int Z = 0; Z < atomic_number_; Z++)
{
DEBUG("Z : " << Z);

Potential[Z] = IonizationTables::ionization_energy(atomic_number_, Z) * eV_to_au;
Azimuthal_quantum_number[Z] = IonizationTables::azimuthal_atomic_number(atomic_number_, Z);
Magnetic_quantum_number[Z] = IonizationTables::magnetic_atomic_number(atomic_number_, Z);

DEBUG("Potential: " << Potential[Z] << " Az.q.num: " << Azimuthal_quantum_number[Z]
<< " M.q.num: " << Magnetic_quantum_number[Z]);

double cst = ((double)Z + 1.0) * sqrt(2.0 / Potential[Z]);
double abs_m = abs(Magnetic_quantum_number[Z]);
alpha_tunnel[Z] = cst - 1.0 - abs_m;
beta_tunnel[Z] = pow(2, alpha_tunnel[Z]) * (8. * Azimuthal_quantum_number[Z] + 4.0) / (cst * tgamma(cst)) *
Potential[Z] * au_to_w0 * tgamma(Azimuthal_quantum_number[Z] + abs_m + 1) /
(tgamma(abs_m + 1) * tgamma(Azimuthal_quantum_number[Z] - abs_m + 1));
gamma_tunnel[Z] = 2.0 * pow(2.0 * Potential[Z], 1.5);
}

DEBUG("Finished Creating the Tunnel Ionizaton class");
}


40 changes: 40 additions & 0 deletions src/Ionization/IonizationTL.cpp
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,40 @@
#include "IonizationTunnel.h"
#include "IonizationTables.h"

template <>
IonizationTunnel<2>::IonizationTunnel(Params &params, Species *species) : Ionization(params, species)
{
DEBUG("Creating the Tong-Lin Tunnel Ionizaton class");

ionization_tl_parameter_ =
species->ionization_tl_parameter_; // species->ionization_tl_parameter_ is double
// Varies from 6 to 9. This is the alpha parameter in Tong-Lin exponential, see Eq. (6) in [M F Ciappina and S V Popruzhenko 2020 Laser Phys. Lett. 17 025301 2020].
lambda_tunnel.resize(atomic_number_);

// Ionization potential & quantum numbers (all in atomic units 1 au = 27.2116 eV)
for (int Z = 0; Z < (int)atomic_number_; Z++)
{
DEBUG("Z : " << Z);
Potential[Z] = IonizationTables::ionization_energy(atomic_number_, Z) * eV_to_au;
Azimuthal_quantum_number[Z] = IonizationTables::azimuthal_atomic_number(atomic_number_, Z);

DEBUG("potential: " << Potential[Z] << " Az.q.num: " << Azimuthal_quantum_number[Z]);

double cst = ((double)Z + 1.0) * sqrt(2.0 / Potential[Z]);
alpha_tunnel[Z] = cst - 1.0;
beta_tunnel[Z] = pow(2, alpha_tunnel[Z]) * (8. * Azimuthal_quantum_number[Z] + 4.0) / (cst * tgamma(cst)) *
Potential[Z] * au_to_w0;
gamma_tunnel[Z] = 2.0 * pow(2.0 * Potential[Z], 1.5);
lambda_tunnel[Z] = ionization_tl_parameter_ * pow(cst, 2) / gamma_tunnel[Z];
}

DEBUG("Finished Creating the Tong-Lin Tunnel Ionizaton class");
}

template <>
double IonizationTunnel<2>::ionizationRate(const int Z, const double E, const double oldZ)
{
const double delta = gamma_tunnel[Z] / E;
return beta_tunnel[Z] * exp(-delta * one_third + alpha_tunnel[Z] * log(delta) - E * lambda_tunnel[Z]);
}

133 changes: 0 additions & 133 deletions src/Ionization/IonizationTunnel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,137 +28,4 @@ IonizationTunnel<0>::IonizationTunnel(Params &params, Species *species) : Ioniza
DEBUG("Finished Creating the Tunnel Ionizaton class");
}

// IonizationFullPPT: 1
template <>
IonizationTunnel<1>::IonizationTunnel(Params &params, Species *species) : Ionization(params, species)
{
DEBUG("Creating the Tunnel Ionizaton class");

// Ionization potential & quantum numbers (all in atomic units 1 au = 27.2116 eV)
Magnetic_quantum_number.resize(atomic_number_);

for (unsigned int Z = 0; Z < atomic_number_; Z++)
{
DEBUG("Z : " << Z);

Potential[Z] = IonizationTables::ionization_energy(atomic_number_, Z) * eV_to_au;
Azimuthal_quantum_number[Z] = IonizationTables::azimuthal_atomic_number(atomic_number_, Z);
Magnetic_quantum_number[Z] = IonizationTables::magnetic_atomic_number(atomic_number_, Z);

DEBUG("Potential: " << Potential[Z] << " Az.q.num: " << Azimuthal_quantum_number[Z]
<< " M.q.num: " << Magnetic_quantum_number[Z]);

double cst = ((double)Z + 1.0) * sqrt(2.0 / Potential[Z]);
double abs_m = abs(Magnetic_quantum_number[Z]);
alpha_tunnel[Z] = cst - 1.0 - abs_m;
beta_tunnel[Z] = pow(2, alpha_tunnel[Z]) * (8. * Azimuthal_quantum_number[Z] + 4.0) / (cst * tgamma(cst)) *
Potential[Z] * au_to_w0 * tgamma(Azimuthal_quantum_number[Z] + abs_m + 1) /
(tgamma(abs_m + 1) * tgamma(Azimuthal_quantum_number[Z] - abs_m + 1));
gamma_tunnel[Z] = 2.0 * pow(2.0 * Potential[Z], 1.5);
}

DEBUG("Finished Creating the Tunnel Ionizaton class");
}

// Tong&Ling: 2
template <>
IonizationTunnel<2>::IonizationTunnel(Params &params, Species *species) : Ionization(params, species)
{
DEBUG("Creating the Tong-Lin Tunnel Ionizaton class");

ionization_tl_parameter_ =
species->ionization_tl_parameter_; // species->ionization_tl_parameter_ is double
// Varies from 6 to 9. This is the alpha parameter in Tong-Lin exponential, see Eq. (6) in [M F Ciappina and S V Popruzhenko 2020 Laser Phys. Lett. 17 025301 2020].
lambda_tunnel.resize(atomic_number_);

// Ionization potential & quantum numbers (all in atomic units 1 au = 27.2116 eV)
for (int Z = 0; Z < (int)atomic_number_; Z++)
{
DEBUG("Z : " << Z);
Potential[Z] = IonizationTables::ionization_energy(atomic_number_, Z) * eV_to_au;
Azimuthal_quantum_number[Z] = IonizationTables::azimuthal_atomic_number(atomic_number_, Z);

DEBUG("potential: " << Potential[Z] << " Az.q.num: " << Azimuthal_quantum_number[Z]);

double cst = ((double)Z + 1.0) * sqrt(2.0 / Potential[Z]);
alpha_tunnel[Z] = cst - 1.0;
beta_tunnel[Z] = pow(2, alpha_tunnel[Z]) * (8. * Azimuthal_quantum_number[Z] + 4.0) / (cst * tgamma(cst)) *
Potential[Z] * au_to_w0;
gamma_tunnel[Z] = 2.0 * pow(2.0 * Potential[Z], 1.5);
lambda_tunnel[Z] = ionization_tl_parameter_ * pow(cst, 2) / gamma_tunnel[Z];
}

DEBUG("Finished Creating the Tong-Lin Tunnel Ionizaton class");
}

template <>
double IonizationTunnel<2>::ionizationRate(const int Z, const double E, const double oldZ)
{
const double delta = gamma_tunnel[Z] / E;
return beta_tunnel[Z] * exp(-delta * one_third + alpha_tunnel[Z] * log(delta) - E * lambda_tunnel[Z]);
}

// BSI: 3
template <>
IonizationTunnel<3>::IonizationTunnel(Params &params, Species *species) : Ionization(params, species)
{
DEBUG("Creating the Tunnel BSI Ionizaton class");

// Ionization potential & quantum numbers (all in atomic units 1 au = 27.2116 eV)
for (int Z = 0; Z < (int)atomic_number_; Z++)
{
DEBUG("Z : " << Z);

Potential[Z] = IonizationTables::ionization_energy(atomic_number_, Z) * eV_to_au;
Azimuthal_quantum_number[Z] = IonizationTables::azimuthal_atomic_number(atomic_number_, Z);

DEBUG("potential: " << Potential[Z] << " Az.q.num: " << Azimuthal_quantum_number[Z]);

double cst = ((double)Z + 1.0) * sqrt(2.0 / Potential[Z]);
alpha_tunnel[Z] = cst - 1.0;
beta_tunnel[Z] = pow(2, alpha_tunnel[Z]) * (8. * Azimuthal_quantum_number[Z] + 4.0) / (cst * tgamma(cst)) *
Potential[Z] * au_to_w0;
gamma_tunnel[Z] = 2.0 * pow(2.0 * Potential[Z], 1.5);
}
DEBUG("Finished Creating the BSI Tunnel Ionizaton class");
}

template <>
double IonizationTunnel<3>::ionizationRate(const int Z, const double E, const double oldZ)
{
auto normal = [this](const int Z, const double E) -> double {
double delta = gamma_tunnel[Z] / E;
return beta_tunnel[Z] * exp(-delta * one_third + alpha_tunnel[Z] * log(delta));
};

auto linear = [this](const int Z, const double E) -> double {
const double ratio_of_IPs = IH / IonizationTables::ionization_energy(atomic_number_, Z);
return au_to_w0 * (0.8 * E * pow(ratio_of_IPs, 0.5));
};
auto quadratic = [this](const int Z, const double E) -> double {
const double ratio_of_IPs = IH / IonizationTables::ionization_energy(atomic_number_, Z);
return au_to_w0 * (2.4 * (pow(E, 2)) * pow(ratio_of_IPs, 2));
};

double BSI_rate_quadratic = quadratic(oldZ, E);
double BSI_rate_linear = linear(oldZ, E);
double Tunnel_rate = normal(oldZ, E);

if (BSI_rate_quadratic >= BSI_rate_linear)
{
return quadratic(Z, E);
}
else if (std::min(Tunnel_rate, BSI_rate_quadratic) == BSI_rate_quadratic)
{
return linear(Z, E);
}
else
{
return normal(Z, E);
}

}




5 changes: 0 additions & 5 deletions src/Ionization/IonizationTunnel.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,6 @@
#include "Ionization.h"
#include "Tools.h"

// because defining in header
#include "IonizationTables.h"
#include "Particles.h"
#include "Species.h"

Expand Down Expand Up @@ -204,9 +202,6 @@ IonizationTunnel<3>::IonizationTunnel(Params &params, Species *species);
template <>
double IonizationTunnel<3>::ionizationRate(const int Z, const double E, const double oldZ);




template <int Model>
void IonizationTunnel<Model>::ionizationTunnelWithTasks( Particles *particles, unsigned int ipart_min, unsigned int ipart_max,
vector<double> *Epart, Patch *patch, Projector *Proj, int ibin, int bin_shift,
Expand Down

0 comments on commit 4253299

Please sign in to comment.