diff --git a/doc/Sphinx/Understand/ionization.rst b/doc/Sphinx/Understand/ionization.rst index 33734c1b7..1f97da880 100755 --- a/doc/Sphinx/Understand/ionization.rst +++ b/doc/Sphinx/Understand/ionization.rst @@ -72,6 +72,66 @@ over a period :math:`2\pi/\omega` leads to the well-known cycle-averaged ionizat \,I_p\,\left( \frac{2 (2 I_p)^{3/2}}{\vert E\vert} \right)^{2n^\star-\vert m \vert -3/2}\, \exp\!\left( -\frac{2 (2 I_p)^{3/2}}{3 \vert E\vert} \right)\,. + +In :program:`Smilei`, four models are available to compute the ionization rate of :eq:`ionizationRate1`. + +**Tunnelling Ionization (PPT-ADK) for :math:`m=0`** + +In the classical model, the ionization rate of :eq:`ionizationRate1` +is computed for :math:`\vert m \vert=0` only. +Indeed, as shown in [Ammosov1986]_, the ratio :math:`R` of the ionization rate +computed for :math:`\vert m\vert=0` by the rate computed for :math:`\vert m\vert=1` is: + +.. math:: + + R = \frac{\Gamma_{{\rm qs},\vert m \vert = 0}}{\Gamma_{{\rm qs},\vert m \vert = 1}} + = 2\frac{(2\,I_p)^{3/2}}{\vert E\vert} + \simeq 7.91\,10^{-3} \,\,\frac{(I_p[\rm eV])^{3/2}}{a_0\,\hbar\omega_0[\rm eV]}\,, + +where, in the practical units formulation, we have considered ionization +by a laser with normalized vector potential :math:`a_0=e\vert E\vert /(m_e c \omega_0)`, +and photon energy :math:`\hbar\omega_0` in eV. + + + +**Tunnelling Ionization (PPT-ADK) for :math:`|m|>0`** +In this model,d ependence on the magnetic quantum number :math:`m` is added. + +:math:`m` is attributed to eac electron in accordance with the following rulse: +1. Since :math:`\Gamma_z(m=0)>\Gamma_z(m=1)>\Gamma_z(m=2)>...` we assume that for electrons +with the same azimuthal quantum number :math:`l`, the states with the lowest value of +:math:`|m|` are ionized first. +2. Electrons with the same azimuthal qunatum number :math:`l` occupy the sub-shells in the +order of increasing :math:`|m|` and for the same :math:`|m|` in the order of increasing :math:`m`. + +With this algorithm, by knowing the atomic number A, we can assign a unique set of +quantum numbers :math:`nlm` to each electron on the atomic sub-shells and identify their extraction +order during successive ionization. + + +**Barrier Suppression Ionization (Tong&Ling)** +The formula proposed by Tong and Lin [Tong2005]_ extends the tunnelling ionization rate to the barrier-suppression +regime. This is achieved by introducing the empirical factor in :eq:`ionizationRate1`: + +.. math:: + + \Gamma_{Z^\star}^{TL} = \Gamma_{Z^\star} \times \exp (-2\alpha_{TL}n^{\star2}\frac{E}{(2I_p)^{3/2}}), + +where :math:`\alpha_TL` is an emprirical constant with value typically from 6 to 9. The actual value +should be guessed from empirical data. When such data is +not available, the formula can be used for qualitative analysis of the barrier-suppression +ionization (BSI), e.g. see [Ciappina2020]_. The module was tested to reproduce the results from this paper. + + +**Barrier Suppression IOnization (Ouatu)** +Ouatu implemented a different model of BSI initially proposed in [Kostyukov2018]_ which was used in the +study [Ouatu2022]_. + + + +To Be removed +"""""""""""""""""""" + In :program:`Smilei`, following [Nuter2011]_, the ionization rate of :eq:`ionizationRate1` is computed for :math:`\vert m \vert=0` only. Indeed, as shown in [Ammosov1986]_, the ratio :math:`R` of the ionization rate @@ -340,3 +400,11 @@ References .. [Schroeder2014] `C. B. Schroeder, J.-L. Vay, E. Esarey, S. S. Bulanov, C. Benedetti, L.-L. Yu, M. Chen, C. G. R. Geddes, and W. P. Leemans, Phys. Rev. ST Accel. Beams 17, 101301 `_ .. [Gibbon] P. Gibbon, Short Pulse Laser Interactions with Matter - An Introduction, Imperial College Press (2005) + +.. [Tong2005] `Tong X. M., Lin C. D., J. Phys. B: At. Mol. Opt. Phys. 38 2593 (2005) ` + +.. [Ciappina2020] `M. F. Ciappina, S. V. Popruzhenko., Laser Phys. Lett. 17 025301 (2020) ` + +.. [Kostyukov2018] `I. Yu. Kostyukov, A. A. Golovanov, Phys. Rev. A 98, 043407 (2018) ` + +.. [Ouatu2022] `I. Ouatu et al, Phys. Rev. E 106, 015205 (2022) ` diff --git a/ionization_module_tests_and_man/cmap_map.py b/ionization_module_tests_and_man/cmap_map.py new file mode 100644 index 000000000..5bb5f939c --- /dev/null +++ b/ionization_module_tests_and_man/cmap_map.py @@ -0,0 +1,33 @@ +import matplotlib +import numpy as np +import matplotlib.pyplot as plt + +def cmap_map(function, cmap): + """ Applies function (which should operate on vectors of shape 3: [r, g, b]), on colormap cmap. + This routine will break any discontinuous points in a colormap. + """ + cdict = cmap._segmentdata + step_dict = {} + # Firt get the list of points where the segments start or end + for key in ('red', 'green', 'blue'): + step_dict[key] = list(map(lambda x: x[0], cdict[key])) + step_list = sum(step_dict.values(), []) + step_list = np.array(list(set(step_list))) + # Then compute the LUT, and apply the function to the LUT + reduced_cmap = lambda step : np.array(cmap(step)[0:3]) + old_LUT = np.array(list(map(reduced_cmap, step_list))) + new_LUT = np.array(list(map(function, old_LUT))) + # Now try to make a minimal segment definition of the new LUT + cdict = {} + for i, key in enumerate(['red','green','blue']): + this_cdict = {} + for j, step in enumerate(step_list): + if step in step_dict[key]: + this_cdict[step] = new_LUT[j, i] + elif new_LUT[j,i] != old_LUT[j, i]: + this_cdict[step] = new_LUT[j, i] + colorvector = list(map(lambda x: x + (x[1], ), this_cdict.items())) + colorvector.sort() + cdict[key] = colorvector + + return matplotlib.colors.LinearSegmentedColormap('colormap',cdict,1024) \ No newline at end of file diff --git a/ionization_module_tests_and_man/input.py b/ionization_module_tests_and_man/input.py new file mode 100755 index 000000000..86cdc7563 --- /dev/null +++ b/ionization_module_tests_and_man/input.py @@ -0,0 +1,210 @@ +import numpy as np +# import sys +# sys.path.append('/home/ent/Smilei_ionization/simulations/') + + +l0 = 2.0*np.pi # laser wavelength +t0 = l0 # optical cicle +Lsim = [12.*l0] # length of the simulation +Tsim = 24.*t0 # duration of the simulation +resx = 128. # nb of cells in one laser wavelength +dx = dy = dz = l0/resx +cell_length = [dx] +dt = 0.95*dx/np.sqrt(1.) +rest = t0/dt +timesteps = int(Tsim/dt) + +Z = 18 +mass = 39.948 + + +# laser pulse input parameters +a0 = 100. +omega = 1. +N_cycles = 10. +xf = Lsim[0]/2. +tcenter = N_cycles*l0 + +Lmu = 0.8 # mu-m + +c = 299792458. # m/sec +omega_ref = 2*np.pi*c/(Lmu*1.e-6) # 1/sec +eps0 = 8.8541878128e-12 # F⋅m^−1 +charge_SI = 1.60217663e-19 # C +m_e_SI = 9.1093837e-31 # kg +N_r = eps0*m_e_SI*omega_ref**2/charge_SI**2 +print('Reference density = ',N_r*10**(-6), ' cm-3') + +I_S = 4.65e29 #W/cm^2 +l_C = 3.8615901e-13 #m +I_ref = 0.5*(a0*omega_ref*l_C/c)**2.*4.65e29 +print('Reference intensity = ',I_ref, ' W/cm^2') + + +# n0 = 5.e13/(N_r*1.e-6) # initial density 5.e13 cm^(-3) +# nppc = 16 # nb of particle per cells + + +n0 = 5.e13/(N_r*1.e-6) +thickness = 0.25*l0 +position = Lsim[0]/2.-thickness/2. +def nf(x): + return n0*np.heaviside(x-position, 1.)*np.heaviside(position+thickness-x, 1.) +nppc = 128 + +Main( + geometry = "1Dcartesian", + + interpolation_order = 2 , + + cell_length = cell_length, + grid_length = Lsim, + + number_of_patches = [1], + + timestep = dt, + simulation_time = Tsim, + + reference_angular_frequency_SI = omega_ref, + + EM_boundary_conditions = [ ["silver-muller", "silver-muller"] ], + +) + +def sin2_envelope(t, tcenter, N): + if (-np.pi*N <= t-tcenter) and (t-tcenter <= np.pi*N): + return np.sin((t-tcenter+np.pi*N)/2./N)**2 + else: + return 0. + + + +Laser( + box_side = "xmin", + omega = omega, + space_time_profile = [ lambda t: 0., + lambda t: a0*sin2_envelope(t, tcenter, N_cycles)\ + *np.cos(omega*((t-tcenter)-xf)) ], +) + + +Species( + name = 'atom_tunnel', + ionization_model = 'tunnel', + ionization_electrons = 'electron', + atomic_number = Z, + position_initialization = 'regular', + momentum_initialization = 'cold', + particles_per_cell = nppc, + mass = mass*1836.0, + charge = 0., + # number_density = n0, + number_density = nf, + boundary_conditions = [['remove','remove']] +) + + +Species( + name = 'atom_TL', + ionization_model = 'tunnel_TL', + ionization_tl_parameter = 6, + ionization_electrons = 'electron', + atomic_number = Z, + position_initialization = 'regular', + momentum_initialization = 'cold', + particles_per_cell = nppc, + mass = mass*1836.0, + charge = 0., + # number_density = n0, + number_density = nf, + boundary_conditions = [['remove','remove']] +) + +Species( + name = 'atom_BSI', + ionization_model = 'tunnel_BSI', + ionization_electrons = 'electron', + atomic_number = Z, + position_initialization = 'regular', + momentum_initialization = 'cold', + particles_per_cell = nppc, + mass = mass*1836.0, + charge = 0., + # number_density = n0, + number_density = nf, + boundary_conditions = [['remove','remove']] +) + +Species( + name = 'atom_full_PPT', + ionization_model = 'tunnel_full_PPT', + ionization_electrons = 'electron', + atomic_number = Z, + position_initialization = 'regular', + momentum_initialization = 'cold', + particles_per_cell = nppc, + mass = mass*1836.0, + charge = 0., + # number_density = n0, + number_density = nf, + boundary_conditions = [['remove','remove']] +) + +Species( + name = 'electron', + position_initialization = 'regular', + momentum_initialization = 'cold', + particles_per_cell = 0, + mass = 1.0, + charge = -1.0, + charge_density = 0.0, + boundary_conditions = [['remove','remove']], +# time_frozen = 2.*Tsim +) + + +# DiagScalar( +# every = rest +# ) + + +DiagParticleBinning( + deposited_quantity = "weight", + every = 1, + species = ["atom_tunnel"], + axes = [ + ["charge", -0.5, Z+0.5, Z+1] + ] +) + +DiagParticleBinning( + deposited_quantity = "weight", + every = 1, + species = ["atom_TL"], + axes = [ + ["charge", -0.5, Z+0.5, Z+1] + ] +) + +DiagParticleBinning( + deposited_quantity = "weight", + every = 1, + species = ["atom_BSI"], + axes = [ + ["charge", -0.5, Z+0.5, Z+1] + ] +) + +DiagParticleBinning( + deposited_quantity = "weight", + every = 1, + species = ["atom_full_PPT"], + axes = [ + ["charge", -0.5, Z+0.5, Z+1] + ] +) + +DiagFields( + every = 100, + fields = ['Ex','Ey','Ez','Bx','By','Bz'] +) \ No newline at end of file diff --git a/ionization_module_tests_and_man/manual/ionization_manual.pdf b/ionization_module_tests_and_man/manual/ionization_manual.pdf new file mode 100644 index 000000000..09a49ff57 Binary files /dev/null and b/ionization_module_tests_and_man/manual/ionization_manual.pdf differ diff --git a/ionization_module_tests_and_man/manual/ionization_manual.tex b/ionization_module_tests_and_man/manual/ionization_manual.tex new file mode 100644 index 000000000..6090c2ca0 --- /dev/null +++ b/ionization_module_tests_and_man/manual/ionization_manual.tex @@ -0,0 +1,141 @@ +\documentclass[prd, preprint, +aps, +amsmath, +amssymb, +onecolumn, +%twocolumn, +nofootinbib, +%showpacs, +%showkeys, +superscriptaddress, +%longbibliography +]{revtex4-2} +\usepackage{graphicx,dcolumn,bm} +\usepackage{amsfonts,amsmath,amssymb,amsthm,amscd} +\usepackage[english]{babel} +\hyphenation{ALPGEN} +\hyphenation{EVTGEN} +\hyphenation{PYTHIA} +%\usepackage{multirow} +\usepackage{physics} +\sloppy\raggedbottom +\usepackage{needspace} +\usepackage{xcolor} +\usepackage{slashed} +\usepackage[export]{adjustbox} +\usepackage{tabularx} +\usepackage{hyperref} + +\begin{document} + + \title{Ionization modules for SMILEI beyond the PPT-ADK formula} + + \author{A.~A. Mironov} + \affiliation{LULI, Sorbonne Université, CNRS, CEA, École Polytechnique, Institut Polytechnique de Paris, F-75255 Paris, France; \url{mironov.hep@gmail.com}} + + \maketitle + + \section{Ionization rates} + \subsection{Tunneling ionization (PPT-ADK)} + All the notations are the same as on the reference page in the SMILEI manual \url{https://smileipic.github.io/Smilei/Understand/ionization.html}: + \begin{gather} + \label{PPT_definition} + \Gamma_{Z^*}=A_{n^*,l^*} B_{l,|m|}\left( \frac{2(2I_p)^{3/2}}{|E|} \right)^{2n^*-|m|-1} \exp\left\lbrace - \frac{2(2I_p)^{3/2}}{3|E|}\right\rbrace,\\ + A_{n^*,l^*} = \frac{2^{2n^*}}{n^*\Gamma(n^*+l^*+1)\Gamma(n^*-l^*)},\\ + B_{l,|m|} = \frac{(2l+1)\Gamma(l+|m|+1)}{2^{|m|}\Gamma(|m|+1)\Gamma(l-|m|+1)}, + \end{gather} + where $\Gamma(x)$ is the Gamma function. + + \subsubsection{Representation in the code (as in the standard ionization module)} + + First, note that the effective orbital angular momentum $l^*=n^*-1$. This allows simplifying $A_{n^*,l^*}$. Second, we assume that the ionization is successive (outer shell electrons are extracted first) and that at the instant of ionization, each electron has the magnetic quantum number $m=0$. The ionization probability for an electron with $m=0$ is higher than for the one with $|m|>0$. The hypothesis is that electrons with $|m|>0$ have enough time to relax to the $m=0$ state once it is not occupied. With this, we set $m=0$ for all the electrons, and rewrite the PPT-ADK formula accordingly (note that in the code notations $Z=Z^*-1$): + \begin{gather} + \label{tunnel} + \Gamma_{Z}=\beta(Z)\exp\left\lbrace -\frac{\Delta}{3} + \alpha(Z)\log\Delta \right\rbrace,\\ + \alpha(Z) = c^*-1,\quad c^*=(Z+1)\sqrt{\frac{2}{I_p}} \equiv 2n^*,\\ + \beta(Z) = I_p\frac{2^{\alpha(Z)} \left[ 8l+4 \right]}{c^*\Gamma(c^*)} ,\\ + \label{gamma_Delta} + \Delta = \frac{\gamma(Z)}{E}, \quad \gamma(Z)=2(2I_p)^{3/2}. + \end{gather} + + \subsubsection{Bug fix in the standard model} + In \texttt{Smilei/src/Ionization/IonizationTunnel.cpp}, in the method \texttt{void IonizationTunnel::operator()}, array \texttt{Dnom[Z]} is not initialized. As a result, for high $Z$ atoms, the array can be uncontrollably allocated with random data, and the result of the operation \texttt{Dnom\_tunnel[k\_times+1] -= D\_sum} is not defined. As the value for \texttt{Dnom[k\_times+1]} is not defined anywhere before this line in the algorithm, the correction \texttt{Dnom\_tunnel[k\_times+1] = -D\_sum} resolves the bug. This correction is implemented in the standard ionization module and all new modules described below. However, it is recommended to review the code and initialize the array properly. + + \subsection{Tunnelling ionization (PPT-ADK) with account for $|m|>0$} + The model adds the dependence on the magnetic quantum number $m$ as in Eq.~\eqref{PPT_definition}. It is implemented in the code by changing the definitions for $\alpha(Z)$ and $\beta(Z)$ in Eq.~\eqref{tunnel} as follows: + \begin{gather} + \alpha(Z) = c^*-|m|-1,\\ + \beta(Z) = I_p\frac{2^{\alpha(Z)} \left[ 8l+4 \right]}{c^*\Gamma(c^*)} \times \frac{\Gamma(l+|m|+1)}{\Gamma(|m|+1)\Gamma(l-|m|+1)}. + \end{gather} + Note that $\beta(Z)$ depends both on $l^*$ and $l$. + + The model is implemented in file \texttt{Smilei/src/Ionization/IonizationTunnelFullPPT.cpp} and can be used in the namelist by putting \texttt{ionization\_model = 'tunnel\_full\_PPT'} in the corresponding atom \texttt{Species}. + + \subsubsection{Implementation details} + The magnetic quantum number $m$ is attributed to each electron in accordance with the following rules. + \begin{enumerate} + \item Since $\Gamma_Z(m=0)>\Gamma_Z(|m|=1)>\Gamma_Z(|m|=2)>\ldots$, we assume that electrons in the state with the lowest value of $|m|$ are ionized first. + + \clearpage + \item Electrons with the same azimuthal quantum number $l$ occupy the sub-shells in the order of increasing $|m|$ and for the same $|m|$ in the order of increasing $m$: + \begin{table}[!h] + \begin{tabular}{l | c | c | c | c | c | c | c | c | c | c | c} + Level order num & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & $\ldots$ \\ \hline + Magnetic q. num $m$ & 0 & 0 & -1 & -1 & 1 & 1 & -2 & -2 & 2 & 2 & $\ldots$ + \end{tabular} + \end{table} + + For example, for the atomic configuration with 3 $p$-electrons ($l=1$) in the outer shell, we choose $m=0,\, 0, -1$ for these electrons. + + \item The states with lower $|m|$ are ionized first. + + \end{enumerate} + With this algorithm, by knowing the atomic number $A$, we can assign a unique set of quantum numbers $nlm$ to each electron on the atomic sub-shells and identify their extraction order during successive ionization. + For example, the full configuration for ${}_{12}$Mg and the corresponding order for the ionization process is: + \begin{table}[h!] + \begin{tabular}{l | c | c | c | c | c | c | c | c | c | c | c | c} + Shell & $3s$ & $3s$ & $2p$ & $2p$ & $2p$ & $2p$ & $2p$ & $2p$ & $2s$ & $2s$ & $1s$ & $1s$ \\\hline + $n$ & 3 & 3 & 2 & 2 & 2 & 2 & 2 & 2 & 2 & 2 & 1 & 1 \\\hline + $l$ & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 0 \\\hline + $m$ & 0 & 0 & 0 & 0 & -1 & -1 & 1 & 1 & 0 & 0 & 0 & 0 \\\hline + Ionization order & 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 & 9 & 10 & 11 & 12 + \end{tabular} + \end{table} + + The values for $m$ for each atom with $A=1\ldots99$ are stored in \texttt{Smilei/src/Ionization/IonizationTables.h} in the table named \texttt{magneticQuantumNumber}, which is similar to the table for $l$ used in the tunnelling formula \eqref{tunnel} in the standard SMILEI package. The table for $m$ was generated based on the table for $l$ following the rules presented above. + + \textbf{Important note!} As a temporary solution for the $f$- ($l=3$) and higher shells, $m$ is always set to zero [namely, the rates are calculated with Eqs.~\eqref{tunnel}-\eqref{gamma_Delta}]. This is due to possible sub-level overlap for the $d$- and $f$-shells, that cannot be reproduced with the described algorithm for the table generation. This issue is to be resolved in the future. + + + \clearpage + \subsection{Tong \& Ling formula} + The formula proposed by Tong and Lin [\href{https://iopscience.iop.org/article/10.1088/0953-4075/38/15/001}{Tong X. M. and Lin C. D., J. Phys. B: At. Mol. Opt. Phys. 38 2593 (2005)}] extends the tunnelling ionization rate to the barrier-suppression regime. This is achieved by introducing the \textit{empirical} factor in Eq.~\eqref{PPT_definition}: + \begin{equation} + \Gamma_{Z^*}^{\text{TL}} = \Gamma_{Z^*}\times \exp\left( -2\alpha_{\text{TL}} n^{*2} \frac{E}{(2 I_p)^{3/2}}\right), + \end{equation} + where $\alpha_{\text{TL}}$ is an empirical constant. Typically, the value for it stays in the interval from 6 to 9. The actual value should be guessed from experimental data. When such data is not available, the formula can be used for qualitative analysis of the barrier-suppression ionization (BSI), e.g. see [\href{https://iopscience.iop.org/article/10.1088/1612-202X/ab6559}{M. F. Ciappina and S. V. Popruzhenko., Laser Phys. Lett. 17 025301 (2020)}]. The module was tested to reproduce the results from this paper. + + \subsubsection{Representation in the code} + In the code, the Tong-Ling formula is implemented as follows: + \begin{gather} + \label{tonglin} + \Gamma_{Z}^\text{TL}=\Gamma_{Z}\times \exp\left[ -E\lambda(Z) \right] ,\\ + \lambda(Z) = \frac{\alpha_{\text{TL}} c^{*2}}{\gamma(Z)}, + \end{gather} + where $\gamma(Z)$ is given by Eq.~\eqref{gamma_Delta}. + + The model is implemented in file \texttt{Smilei/src/Ionization/IonizationTunnelTL.cpp} and can be used in the namelist by putting \texttt{ionization\_model = 'tunnel\_TL'} in the corresponding atom \texttt{Species}. The $\alpha_{\text{TL}}$ can be set with the (newly added) parameter \texttt{ionization\_tl\_parameter = } in the corresponding \texttt{Species} + (set to 6 by default). + + + \clearpage + \subsection{BSI formula} + I. Ouatu implemented a different model of BSI initially proposed in [\href{https://journals.aps.org/pra/abstract/10.1103/PhysRevA.98.043407}{I. Yu. Kostyukov and A. A. Golovanov, Phys. Rev. A 98, 043407 (2018)}]. The original implementation by I. Ouatu is available on Github \url{https://github.com/iouatu/mySmilei}. It was used in the study [\href{https://journals.aps.org/pre/abstract/10.1103/PhysRevE.106.015205}{I. Ouatu et al, + Phys. Rev. E 106, 015205 (2022)}]. The module is included in this version of SMILEI without changes (except for some small corrections in the header file) and was tested to reproduce the results of this paper. + + The model uses the ionization rate Eq.~\eqref{tunnel} when applicable and switches to an empirical formula described in the original paper by Kostyukov \& Golovanov. There are no additional parameters introduced to tune the model. + + The model is implemented in file \texttt{Smilei/src/Ionization/IonizationTunnelBSI.cpp} and can be used in the namelist by putting \texttt{ionization\_model = 'tunnel\_BSI'} in the corresponding atom \texttt{Species}. In the implementation, the code relies on some additional functionality introduced in \texttt{Smilei/src/Ionization/continuity\_tool.cpp}. + + +\end{document} \ No newline at end of file diff --git a/ionization_module_tests_and_man/manual/ionization_manualNotes.bib b/ionization_module_tests_and_man/manual/ionization_manualNotes.bib new file mode 100644 index 000000000..8f3dc15db --- /dev/null +++ b/ionization_module_tests_and_man/manual/ionization_manualNotes.bib @@ -0,0 +1,2 @@ +@CONTROL{REVTEX42Control} +@CONTROL{apsrev42Control,author="08",editor="1",pages="0",title="0",year="1"} diff --git a/ionization_module_tests_and_man/run_test.py b/ionization_module_tests_and_man/run_test.py new file mode 100755 index 000000000..1d79a72f0 --- /dev/null +++ b/ionization_module_tests_and_man/run_test.py @@ -0,0 +1,164 @@ +# +# ANALYSIS OF TUNNEL IONISATION SIMULATION +# + +simulation_to_analyse = '.' + +import matplotlib as mpl +import matplotlib.pyplot as plt +import numpy as np +from cmap_map import cmap_map + +# IMPORT SMILEI's POST-PROCESSING TOOL +# ------------------------------------ + +import happi + +import subprocess + +### RUN SMILEI SIMULATION +subprocess.run(["../../smilei", "input.py"]) + + +# LOADING SIMULATION & VARIABLES FROM NAMELIST +# ------------------------------------------------------ + +S = happi.Open(simulation_to_analyse ,verbose=False) + +t0 = 2.*np.pi +Lv = 0. +Lp = S.namelist.Lsim[0] +dt = S.namelist.Main.timestep +Z_A = S.namelist.Z + +a0 = S.namelist.a0 +resx = S.namelist.resx +Tsim = S.namelist.Tsim +Lsim = S.namelist.Lsim + +print('- ref. ang. frequency w0 = '+str(S.namelist.Main.reference_angular_frequency_SI)) + +### PLOT FIELD + +def plot_field(): + Data = S.Field(0, "Ey").getData() + i = int(len(Data)/3*1.913) + x = np.linspace(0, Lsim, Data[i].size) + + fig = plt.figure() + ax = fig.add_subplot(111) + ax.plot(x/t0, Data[i], linewidth=2.5) + + plt.xlim(x[0]/t0,x[-1]/t0) + plt.xlabel(r'$x/\lambda$') + plt.ylabel('Ey [a0]') + + +# SIMULATION ANALYSIS & COMPARISON TO RATE EQUATIONS +# -------------------------------------------------- + + +# get corresponding time-steps +t = dt * np.array(S.ParticleBinning(0).get()['times']) + + +# Species diagnostic number -> ionization model name +diag_model_dict = {0: 'Tunneling', + 1: 'Tong-Lin', + 2: 'Barrier suppression (BSI)', + 3: 'Full PPT (m != 0)'} + + +### PLOT CHARGE STATES FOR ONE SPECIES + +def plot_charge_states(diag, compare_to=None): + ''' + diag = 0, 1, 2, or 3 --- diagnostic number + compare_to = None, 0, 1, 2, or 3 --- if not None, the function plots another diagnostic for comparison + ''' + + # assign a color to each charge state + cmap = mpl.colormaps.get_cmap('gist_ncar') + cmap_dark = cmap_map(lambda x: 0.8*x, cmap) + charge = np.linspace(0, Z_A+1) + + n = np.array( S.ParticleBinning(diag).getData() ) + n00 = n[0,0] + n = n/n00 + + fig = plt.figure() + ax = fig.add_subplot(111) + for Z in range(Z_A+1): + rgba = cmap(Z/float(Z_A+1)) + ax.plot(t/t0, n[:,Z], color=rgba, linewidth=2, label='Z = %d'%Z) + title = diag_model_dict[diag] + + if compare_to != None: + n2 = array( S.ParticleBinning(compare_to).getData() ) + n200 = n2[0,0] + n2 = n2/n200 + for Z in range(Z_A+1): + rgba = cmap_dark(Z/float(Z_A+1)) + ax.plot(t/t0, n2[:,Z], '--', color=rgba, linewidth=1, label='Z = %d'%Z) + title = 'Solid: ' + title + '; dashed: ' + diag_model_dict[compare_to] + ax.legend(loc='upper left', ncol=1 if compare_to==None else 2, handletextpad=0.1, labelspacing=0.1) + ax.set_xlabel(r'$t/T$') + ax.set_ylabel(r'$n_i/n_{total}$') + ax.set_title(title) + + +def compare_models(tmin=10.5, tmax=14): + ''' + tmin and tmax set the plot limits alon the x-axis + ''' + cmap = mpl.colormaps.get_cmap('gist_ncar') + charge = np.linspace(0, Z_A+1) + + ### plot the density of each charge state as a function of time + fig, *axes = plt.subplots(nrows=4, ncols=1) + for i in range(4): + ax = axes[0][i] + n = np.array( S.ParticleBinning(i).getData() ) + n00 = n[0,0] + n = n/n00 + for Z in range(Z_A+1): + rgba = cmap(Z/float(Z_A+1)) + ax.plot(t/t0, n[:,Z], color=rgba, linewidth=1.2, label='Z = %d'%Z) + ax.set_xlim([tmin, tmax]) + ax.grid() + ax.text(tmin+0.05, + 0.85, + diag_model_dict[i], + bbox=dict(boxstyle='round', facecolor='w', edgecolor='k', alpha=0.05)) + ax.set_ylabel(r'$n_i/n_{total}$') + if i<3: + ax.set_xticklabels([]) + if i==0: + ax.legend(loc='upper left', + ncol=Z/3, + handletextpad=0.1, + labelspacing=0.1, + bbox_to_anchor=(0., 1.6)) + axes[0][-1].set_xlabel('t/T') + + plt.subplots_adjust(hspace=0) + + fig.set_size_inches(7, 7) + + +if __name__ == '__main__': + plot_field() + + plot_charge_states(0) + plot_charge_states(1) + plot_charge_states(2) + plot_charge_states(3) + + compare_models(tmin=10.5, tmax=14) + + plt.show(block=True) + + + + + diff --git a/src/Ionization/Ionization.cpp b/src/Ionization/Ionization.cpp index debc48b71..8ba734f2c 100755 --- a/src/Ionization/Ionization.cpp +++ b/src/Ionization/Ionization.cpp @@ -1,71 +1,69 @@ +#include "Ionization.h" + #include -#include "Ionization.h" #include "Species.h" using namespace std; -Ionization::Ionization( Params ¶ms, Species *species ) +Ionization::Ionization(Params ¶ms, Species *species) { - reference_angular_frequency_SI = params.reference_angular_frequency_SI; - - dt = params.timestep; - invdt = 1./dt; - nDim_field = params.nDim_field; - nDim_particle = params.nDim_particle; - ionized_species_invmass = 1./species->mass_; - + + dt = params.timestep; + invdt = 1. / dt; + nDim_field = params.nDim_field; + nDim_particle = params.nDim_particle; + ionized_species_invmass = 1. / species->mass_; + // Normalization constant from Smilei normalization to/from atomic units - eV_to_au = 1.0 / 27.2116; - au_to_mec2 = 27.2116/510.998e3; - EC_to_au = 3.314742578e-15 * reference_angular_frequency_SI; // hbar omega / (me c^2 alpha^3) - au_to_w0 = 4.134137172e+16 / reference_angular_frequency_SI; // alpha^2 me c^2 / (hbar omega) - + eV_to_au = 1.0 / 27.2116; + au_to_mec2 = 27.2116 / 510.998e3; + EC_to_au = 3.314742578e-15 * reference_angular_frequency_SI; // hbar omega / (me c^2 alpha^3) + au_to_w0 = 4.134137172e+16 / reference_angular_frequency_SI; // alpha^2 me c^2 / (hbar omega) + #ifdef _OMPTASKS new_electrons_per_bin = new Particles[species->Nbins]; - ion_charge_per_bin_.resize( species->Nbins ); + ion_charge_per_bin_.resize(species->Nbins); #endif } +Ionization::~Ionization() {} -Ionization::~Ionization() +void Ionization::joinNewElectrons(unsigned int Nbins) { -} - - -void Ionization::joinNewElectrons( unsigned int Nbins ) -{ - // if tasks on bins are used for ionization, join the lists of new electrons + // if tasks on bins are used for ionization, join the lists of new electrons // created in each bin, to have the list of new electrons for this species and patch - + size_t start = new_electrons.size(); - + // Resize new_electrons size_t total_n_new = 0; - for( size_t ibin = 0 ; ibin < Nbins ; ibin++ ) { + for (size_t ibin = 0; ibin < Nbins; ibin++) { total_n_new += new_electrons_per_bin[ibin].size(); } - new_electrons.createParticles( total_n_new ); + new_electrons.createParticles(total_n_new); // Also resize ion_charge_ if necessary - if( save_ion_charge_ ) { - ion_charge_.resize( start + total_n_new ); + if (save_ion_charge_) { + ion_charge_.resize(start + total_n_new); } - + // Move each new_electrons_per_bin into new_electrons - for( size_t ibin = 0 ; ibin < Nbins ; ibin++ ) { + for (size_t ibin = 0; ibin < Nbins; ibin++) { size_t n_new = new_electrons_per_bin[ibin].size(); - for( size_t i=0; i #include -#include "Tools.h" -#include "Params.h" -#include "Patch.h" #include "Field.h" +#include "Params.h" #include "Particles.h" +#include "Patch.h" #include "Projector.h" +#include "Tools.h" +using namespace std; //! Class Ionization: generic class allowing to define Ionization physics class Ionization { - -public: + public: //! Constructor for Ionization - Ionization( Params ¶ms, Species *species ); + Ionization(Params ¶ms, Species *species); virtual ~Ionization(); - + //! Overloading of () operator - virtual void operator()( Particles *, unsigned int, unsigned int, std::vector *, Patch *, Projector *, int = 0 ) {}; + virtual void operator()(Particles *, unsigned int, unsigned int, std::vector *, Patch *, Projector *, int = 0) {}; //! method for envelope ionization - virtual void envelopeIonization( Particles *, unsigned int, unsigned int, std::vector *, std::vector *, std::vector *, std::vector *, Patch *, Projector *, int = 0, int = 0 ){}; - + virtual void envelopeIonization(Particles *, unsigned int, unsigned int, std::vector *, std::vector *, + std::vector *, std::vector *, Patch *, Projector *, int = 0, int = 0) {}; + // method for tunnel ionization using tasks - virtual void ionizationTunnelWithTasks( Particles *, unsigned int, unsigned int, std::vector *, Patch *, Projector *, int, int, double *, double *, double *, int = 0 ){}; + virtual void ionizationTunnelWithTasks(Particles *, unsigned int, unsigned int, std::vector *, Patch *, + Projector *, int, int, double *, double *, double *, int = 0) {}; // join the lists of electrons created through ionization when tasks are used void joinNewElectrons(unsigned int Nbins); Particles new_electrons; Particles *new_electrons_per_bin; - + //! Whether the initial charge (of the atom that was ionized) should be saved bool save_ion_charge_ = false; //! Temporarily contains the initial charge of the atom that was ionized std::vector ion_charge_; std::vector > ion_charge_per_bin_; - -protected: + protected: double eV_to_au; double au_to_mec2; double EC_to_au; double au_to_w0; - + double reference_angular_frequency_SI; double dt; double invdt; unsigned int nDim_field; unsigned int nDim_particle; double ionized_species_invmass; - -private: - + private: }; #endif diff --git a/src/Ionization/IonizationBSI.cpp b/src/Ionization/IonizationBSI.cpp new file mode 100644 index 000000000..abbac99f5 --- /dev/null +++ b/src/Ionization/IonizationBSI.cpp @@ -0,0 +1,71 @@ +#include "IonizationTables.h" +#include "IonizationTunnel.h" + +template <> +IonizationTunnel<3>::IonizationTunnel(Params ¶ms, Species *species) : Ionization(params, species) +{ + DEBUG("Creating the Tunnel BSI Ionizaton class"); + + atomic_number_ = species->atomic_number_; + Potential.resize(atomic_number_); + Azimuthal_quantum_number.resize(atomic_number_); + + one_third = 1.0 / 3.0; + + alpha_tunnel.resize(atomic_number_); + beta_tunnel.resize(atomic_number_); + gamma_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); + } + DEBUG("Finished Creating the BSI Tunnel Ionizaton class"); +} + +template <> +double IonizationTunnel<3>::ionizationRate1(const int Z, const double E) +{ + double ratio_of_IPs = IH / IonizationTables::ionization_energy(atomic_number_, Z); + + double BSI_rate_quadratic = 2.4 * (pow(E, 2))*pow(ratio_of_IPs, 2) * au_to_w0; + double BSI_rate_linear = 0.8 * E * pow(ratio_of_IPs, 0.5) * au_to_w0; + double delta = gamma_tunnel[Z] / E; + double Tunnel_rate = beta_tunnel[Z] * exp(-delta / 3.0 + alpha_tunnel[Z] * log(delta)); + + if (std::min(Tunnel_rate, BSI_rate_quadratic) == BSI_rate_quadratic) { + rate_formula = 1; + return BSI_rate_quadratic; + } else if (BSI_rate_quadratic >= BSI_rate_linear) { + rate_formula = 2; + return BSI_rate_quadratic; + } else { + rate_formula = 0; + return Tunnel_rate; + } +} + +template <> +double IonizationTunnel<3>::ionizationRate2(const int newZ, const double E) +{ + double ratio_of_IPs_newZ = IH / IonizationTables::ionization_energy(atomic_number_, newZ); + double delta = gamma_tunnel[newZ] / E; + if (rate_formula == 1) { + return au_to_w0 * (2.4 * (pow(E, 2))*pow(ratio_of_IPs_newZ, 2)); + } else if (rate_formula == 2) { + return au_to_w0 * (0.8 * E * pow(ratio_of_IPs_newZ, 0.5)); + } else { + return beta_tunnel[newZ] * exp(-delta * one_third + alpha_tunnel[newZ] * log(delta)); + } +} diff --git a/src/Ionization/IonizationFactory.h b/src/Ionization/IonizationFactory.h old mode 100755 new mode 100644 index 87c985eff..a7039bec0 --- a/src/Ionization/IonizationFactory.h +++ b/src/Ionization/IonizationFactory.h @@ -2,64 +2,78 @@ #define IonizationFactory_H #include "Ionization.h" -#include "IonizationTunnel.h" #include "IonizationFromRate.h" #include "IonizationTunnelEnvelopeAveraged.h" - #include "Params.h" - -#include "Tools.h" - #include "Species.h" +#include "IonizationTunnel.h" +#include "Tools.h" //! this class create and associate the right ionization model to species class IonizationFactory { -public: - static Ionization *create( Params ¶ms, Species *species ) + public: + static Ionization *create(Params ¶ms, Species *species) { Ionization *Ionize = NULL; - std::string model=species->ionization_model_; - - if( model == "tunnel" ) { - - if( species->max_charge_ > ( int )species->atomic_number_ ) { - ERROR( "Charge > atomic_number for species " << species->name_ ); - } + std::string model = species->ionization_model_; - if( params.Laser_Envelope_model ) { - ERROR( "The ionization model for species interacting with envelope is tunnel_envelope_averaged" ); - } - - Ionize = new IonizationTunnel( params, species ); + if( model == "tunnel" ) { + checkMaxCharge(species); + checkNotLaserEnvelopeModel(params); + Ionize = new IonizationTunnel<0>( params, species ); // The original model included in Smilei } else if( model == "tunnel_envelope_averaged" ) { - if( species->max_charge_ > ( int )species->atomic_number_ ) { - ERROR( "Charge > atomic_number for species " << species->name_ ); - } - if( species->particles->is_test ) { - ERROR( "Cannot ionize test species " << species->name_ ); - } - - Ionize = new IonizationTunnelEnvelopeAveraged( params, species ); - + checkMaxCharge(species); + checkTestParticle(species); if ( !params.Laser_Envelope_model ) { ERROR( "The ionization model tunnel_envelope_averaged needs a laser envelope"); } + Ionize = new IonizationTunnelEnvelopeAveraged( params, species ); + } else if( model == "from_rate" ) { - - if( species->max_charge_ > ( int )species->maximum_charge_state_ ) { + if ( species->max_charge_ > ( int ) species->maximum_charge_state_ ) { ERROR( "For species '" << species->name_ << ": charge > maximum_charge_state" ); } - + Ionize = new IonizationFromRate( params, species ); - + + } else if (model == "tunnel_full_PPT") { + checkMaxCharge(species); + checkNotLaserEnvelopeModel(params); + Ionize = new IonizationTunnel<1>(params, species); // FullPPT + } else if (model == "tunnel_TL") { + checkMaxCharge(species); + checkNotLaserEnvelopeModel(params); + Ionize = new IonizationTunnel<2>(params, species); // Tong&Ling + } else if (model == "tunnel_BSI") { + checkMaxCharge(species); + checkNotLaserEnvelopeModel(params); + Ionize = new IonizationTunnel<3>(params, species); // BSI } - + return Ionize; } + + private: + static void checkMaxCharge(const Species *species) { + if ( species->max_charge_ > ( int )species->atomic_number_ ) { + ERROR( "Charge > atomic_number for species " << species->name_ ); + } + } + + static void checkTestParticle(const Species *species) { + if( species->particles->is_test ) { + ERROR( "Cannot ionize test species " << species->name_ ); + } + } + static void checkNotLaserEnvelopeModel(const Params ¶ms) { + if ( params.Laser_Envelope_model ) { + ERROR( "The ionization model for species interacting with envelope is tunnel_envelope_averaged" ); + } + } }; #endif diff --git a/src/Ionization/IonizationFullPPT.cpp b/src/Ionization/IonizationFullPPT.cpp new file mode 100644 index 000000000..40d6d81cb --- /dev/null +++ b/src/Ionization/IonizationFullPPT.cpp @@ -0,0 +1,42 @@ +#include "IonizationTables.h" +#include "IonizationTunnel.h" + +template <> +IonizationTunnel<1>::IonizationTunnel(Params ¶ms, Species *species) : Ionization(params, species) +{ + DEBUG("Creating the Tunnel Ionizaton class"); + + atomic_number_ = species->atomic_number_; + Potential.resize(atomic_number_); + Azimuthal_quantum_number.resize(atomic_number_); + + one_third = 1.0 / 3.0; + + alpha_tunnel.resize(atomic_number_); + beta_tunnel.resize(atomic_number_); + gamma_tunnel.resize(atomic_number_); + + // 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"); +} diff --git a/src/Ionization/IonizationTL.cpp b/src/Ionization/IonizationTL.cpp new file mode 100644 index 000000000..bf43c20a6 --- /dev/null +++ b/src/Ionization/IonizationTL.cpp @@ -0,0 +1,54 @@ +#include "IonizationTables.h" +#include "IonizationTunnel.h" + +template <> +IonizationTunnel<2>::IonizationTunnel(Params ¶ms, Species *species) : Ionization(params, species) +{ + DEBUG("Creating the Tong-Lin Tunnel Ionizaton class"); + + atomic_number_ = species->atomic_number_; + Potential.resize(atomic_number_); + Azimuthal_quantum_number.resize(atomic_number_); + + one_third = 1.0 / 3.0; + + alpha_tunnel.resize(atomic_number_); + beta_tunnel.resize(atomic_number_); + gamma_tunnel.resize(atomic_number_); + + 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>::ionizationRate1(const int Z, const double E) +{ + const double delta = gamma_tunnel[Z] / E; + return beta_tunnel[Z] * exp(-delta * one_third + alpha_tunnel[Z] * log(delta) - E * lambda_tunnel[Z]); +} + +template <> +double IonizationTunnel<2>::ionizationRate2(const int Z, const double E) +{ + return ionizationRate1(Z, E); +} diff --git a/src/Ionization/IonizationTables.cpp b/src/Ionization/IonizationTables.cpp index 98df43bdb..860f50ae3 100755 --- a/src/Ionization/IonizationTables.cpp +++ b/src/Ionization/IonizationTables.cpp @@ -12,6 +12,12 @@ int IonizationTables::azimuthal_atomic_number( int atomic_number, int Zstar ) return ( int )azimuthalQuantumNumber[( atomic_number*( atomic_number-1 ) )/2 + Zstar]; } +// Gets the magnetic atomic number of a given ion +int IonizationTables::magnetic_atomic_number( int atomic_number, int Zstar ) +{ + return ( int )magneticQuantumNumber[( atomic_number*( atomic_number-1 ) )/2 + Zstar]; +} + // Gets the k-th binding energy in any neutral or ionized atom with atomic number Z and charge Zstar // We use the formula by Carlson et al., At. Data Nucl. Data Tables 2, 63 (1970) double IonizationTables::binding_energy( int atomic_number, int Zstar, int k ) diff --git a/src/Ionization/IonizationTables.h b/src/Ionization/IonizationTables.h index c93fe18fb..f623a0af6 100755 --- a/src/Ionization/IonizationTables.h +++ b/src/Ionization/IonizationTables.h @@ -1197,6 +1197,109 @@ const unsigned char azimuthalQuantumNumber[100/2*101] = { 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0 }; +/* Added by Arseny Mironov. Note that for l=3 m is always set to zero. This is a temporary solution. */ +const signed char magneticQuantumNumber[100/2*101] = { + /* Z=1 */ 0, + /* Z=2 */ 0, 0, + /* Z=3 */ 0, 0, 0, + /* Z=4 */ 0, 0, 0, 0, + /* Z=5 */ 0, 0, 0, 0, 0, + /* Z=6 */ 0, 0, 0, 0, 0, 0, + /* Z=7 */ 0, 0,-1, 0, 0, 0, 0, + /* Z=8 */ 0, 0,-1,-1, 0, 0, 0, 0, + /* Z=9 */ 0, 0,-1,-1, 1, 0, 0, 0, 0, + /* Z=10 */ 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=11 */ 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=12 */ 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=13 */ 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=14 */ 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=15 */ 0, 0,-1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=16 */ 0, 0,-1,-1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=17 */ 0, 0,-1,-1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=18 */ 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=19 */ 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=20 */ 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=21 */ 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=22 */ 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=23 */ 0, 0, 0,-1,-1, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=24 */ 0, 0, 0,-1,-1, 1, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=25 */ 0, 0, 0, 0,-1,-1, 1, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=26 */ 0, 0, 0, 0,-1,-1, 1, 1, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=27 */ 0, 0, 0,-1,-1, 1, 1,-2,-2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=28 */ 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=29 */ 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=30 */ 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=31 */ 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=32 */ 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=33 */ 0, 0,-1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=34 */ 0, 0,-1,-1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=35 */ 0, 0,-1,-1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=36 */ 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=37 */ 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=38 */ 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=39 */ 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=40 */ 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=41 */ 0, 0, 0,-1,-1, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=42 */ 0, 0, 0,-1,-1, 1, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=43 */ 0, 0, 0, 0,-1,-1, 1, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=44 */ 0, 0, 0,-1,-1, 1, 1,-2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=45 */ 0, 0, 0,-1,-1, 1, 1,-2,-2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=46 */ 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=47 */ 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=48 */ 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=49 */ 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=50 */ 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=51 */ 0, 0,-1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=52 */ 0, 0,-1,-1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=53 */ 0, 0,-1,-1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=54 */ 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=55 */ 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=56 */ 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=57 */ 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=58 */ 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=59 */ 0, 0, 0, 0, 0, 0, 0,-1,-1, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=60 */ 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=61 */ 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=62 */ 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=63 */ 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=64 */ 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=65 */ 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=66 */ 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=67 */ 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=68 */ 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=69 */ 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=70 */ 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=71 */ 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=72 */ 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=73 */ 0, 0, 0, 0,-1, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=74 */ 0, 0, 0, 0,-1,-1, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=75 */ 0, 0, 0, 0,-1,-1, 1, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=76 */ 0, 0, 0, 0, 0,-1,-1, 1, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=77 */ 0, 0, 0, 0,-1,-1, 1, 1,-2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=78 */ 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=79 */ 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=80 */ 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=81 */ 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=82 */ 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=83 */ 0, 0,-1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=84 */ 0, 0,-1,-1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=85 */ 0, 0,-1,-1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=86 */ 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=87 */ 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=88 */ 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=89 */ 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=90 */ 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=91 */ 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=92 */ 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0,-1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=93 */ 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0,-1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=94 */ 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0,-1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=95 */ 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0,-1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=96 */ 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=97 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=98 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0, + /* Z=99 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0, 0,-1,-1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1,-2,-2, 2, 2, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0,-1,-1, 1, 1, 0, 0, 0, 0 +}; + // Gets the ionization energy of a given ion double ionization_energy( int atomic_number, int Zstar ); @@ -1204,10 +1307,13 @@ double ionization_energy( int atomic_number, int Zstar ); // Gets the azimuthal atomic number of a given ion int azimuthal_atomic_number( int atomic_number, int Zstar ); +// Gets the magnetic atomic number of a given ion +int magnetic_atomic_number( int atomic_number, int Zstar ); + // Gets the k-th binding energy in any neutral or ionized atom with atomic number Z and charge Zstar // We use the formula by Carlson et al., At. Data Nucl. Data Tables 2, 63 (1970) double binding_energy( int atomic_number, int Zstar, int k ); }; -#endif \ No newline at end of file +#endif diff --git a/src/Ionization/IonizationTunnel.cpp b/src/Ionization/IonizationTunnel.cpp old mode 100755 new mode 100644 index abadcc5bd..ddbd51f64 --- a/src/Ionization/IonizationTunnel.cpp +++ b/src/Ionization/IonizationTunnel.cpp @@ -1,316 +1,38 @@ #include "IonizationTunnel.h" -#include "IonizationTables.h" - -#include - -#include "Particles.h" -#include "Species.h" - -using namespace std; - +#include "IonizationTables.h" -IonizationTunnel::IonizationTunnel( Params ¶ms, Species *species ) : Ionization( params, species ) +// Classic: 1 +template <> +IonizationTunnel<0>::IonizationTunnel(Params ¶ms, Species *species) : Ionization(params, species) { - DEBUG( "Creating the Tunnel Ionizaton class" ); - - atomic_number_ = species->atomic_number_; - - // Ionization potential & quantum numbers (all in atomic units 1 au = 27.2116 eV) - Potential.resize( atomic_number_ ); - Azimuthal_quantum_number.resize( atomic_number_ ); - for( int Zstar=0; Zstar<( int )atomic_number_; Zstar++ ) { - Potential [Zstar] = IonizationTables::ionization_energy( atomic_number_, Zstar ) * eV_to_au; - Azimuthal_quantum_number[Zstar] = IonizationTables::azimuthal_atomic_number( atomic_number_, Zstar ); - } - - for( unsigned int i=0; iatomic_number_; + Potential.resize(atomic_number_); + Azimuthal_quantum_number.resize(atomic_number_); + one_third = 1.0 / 3.0; + alpha_tunnel.resize(atomic_number_); + beta_tunnel.resize(atomic_number_); + gamma_tunnel.resize(atomic_number_); -void IonizationTunnel::operator()( Particles *particles, unsigned int ipart_min, unsigned int ipart_max, vector *Epart, Patch *patch, Projector *Proj, int ipart_ref ) -{ + // Ionization potential & quantum numbers (all in atomic units 1 au = 27.2116 eV) + for (unsigned int Z = 0; Z < atomic_number_; Z++) { + DEBUG("Z : " << Z); - unsigned int Z, Zp1, newZ, k_times; - double TotalIonizPot, E, invE, factorJion, delta, ran_p, Mult, D_sum, P_sum, Pint_tunnel; - vector IonizRate_tunnel( atomic_number_ ), Dnom_tunnel( atomic_number_ , 0.); - LocalFields Jion; - double factorJion_0 = au_to_mec2 * EC_to_au*EC_to_au * invdt; - - int nparts = Epart->size()/3; - double *Ex = &( ( *Epart )[0*nparts] ); - double *Ey = &( ( *Epart )[1*nparts] ); - double *Ez = &( ( *Epart )[2*nparts] ); - - for( unsigned int ipart=ipart_min ; ipartcharge( ipart ) ); - - // If ion already fully ionized then skip - if( Z==atomic_number_ ) { - continue; - } - - // Absolute value of the electric field normalized in atomic units - E = EC_to_au * sqrt( pow( *( Ex+ipart-ipart_ref ), 2 ) - +pow( *( Ey+ipart-ipart_ref ), 2 ) - +pow( *( Ez+ipart-ipart_ref ), 2 ) ); - if( E<1e-10 ) { - continue; - } - - // -------------------------------- - // Start of the Monte-Carlo routine - // -------------------------------- - - invE = 1./E; - factorJion = factorJion_0 * invE*invE; - delta = gamma_tunnel[Z]*invE; - ran_p = patch->rand_->uniform(); - IonizRate_tunnel[Z] = beta_tunnel[Z] * exp( -delta*one_third + alpha_tunnel[Z]*log( delta ) ); - - // Total ionization potential (used to compute the ionization current) - TotalIonizPot = 0.0; - - // k_times will give the nb of ionization events - k_times = 0; - Zp1=Z+1; - - if( Zp1 == atomic_number_ ) { - // if ionization of the last electron: single ionization - // ----------------------------------------------------- - if( ran_p < 1.0 -exp( -IonizRate_tunnel[Z]*dt ) ) { - TotalIonizPot += Potential[Z]; - k_times = 1; - } - - } else { - // else : multiple ionization can occur in one time-step - // partial & final ionization are decoupled (see Nuter Phys. Plasmas) - // ------------------------------------------------------------------------- - - // initialization - Mult = 1.0; - Dnom_tunnel[0]=1.0; - Pint_tunnel = exp( -IonizRate_tunnel[Z]*dt ); // cummulative prob. - - //multiple ionization loop while Pint_tunnel < ran_p and still partial ionization - while( ( Pint_tunnel < ran_p ) and ( k_times < atomic_number_-Zp1 ) ) { - newZ = Zp1+k_times; - delta = gamma_tunnel[newZ]*invE; - IonizRate_tunnel[newZ] = beta_tunnel[newZ] - * exp( -delta*one_third+alpha_tunnel[newZ]*log( delta ) ); - D_sum = 0.0; - P_sum = 0.0; - Mult *= IonizRate_tunnel[Z+k_times]; - for( unsigned int i=0; iran_p ) && ( k_times==atomic_number_-Zp1 ) ) { - TotalIonizPot += Potential[atomic_number_-1]; - k_times++; - } - }//END Multiple ionization routine - - // Compute ionization current - if (patch->EMfields->Jx_ != NULL){ // For the moment ionization current is not accounted for in AM geometry - factorJion *= TotalIonizPot; - Jion.x = factorJion * *( Ex+ipart ); - Jion.y = factorJion * *( Ey+ipart ); - Jion.z = factorJion * *( Ez+ipart ); - - Proj->ionizationCurrents( patch->EMfields->Jx_, patch->EMfields->Jy_, patch->EMfields->Jz_, *particles, ipart, Jion ); - } - - // Creation of the new electrons - // (variable weights are used) - // ----------------------------- + Potential[Z] = IonizationTables::ionization_energy(atomic_number_, Z) * eV_to_au; + Azimuthal_quantum_number[Z] = IonizationTables::azimuthal_atomic_number(atomic_number_, Z); - if( k_times !=0 ) { - new_electrons.createParticle(); - int idNew = new_electrons.size() - 1; - for( unsigned int i=0; iposition( i, ipart ); - } - for( unsigned int i=0; i<3; i++ ) { - new_electrons.momentum( i, idNew ) = particles->momentum( i, ipart )*ionized_species_invmass; - } - new_electrons.weight( idNew )=double( k_times )*particles->weight( ipart ); - new_electrons.charge( idNew )=-1; - - if( save_ion_charge_ ) { - ion_charge_.push_back( particles->charge( ipart ) ); - } - - // Increase the charge of the particle - particles->charge( ipart ) += k_times; - } - - } // Loop on particles -} + DEBUG("Potential: " << Potential[Z] << " Az.q.num: " << Azimuthal_quantum_number[Z]); -void IonizationTunnel::ionizationTunnelWithTasks( Particles *particles, unsigned int ipart_min, unsigned int ipart_max, - vector *Epart, Patch *patch, Projector *Proj, int ibin, int bin_shift, - double *b_Jx, double *b_Jy, double *b_Jz, int ipart_ref ) -{ + 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); + } - unsigned int Z, Zp1, newZ, k_times; - double TotalIonizPot, E, invE, factorJion, delta, ran_p, Mult, D_sum, P_sum, Pint_tunnel; - vector IonizRate_tunnel( atomic_number_ ), Dnom_tunnel( atomic_number_ ); - LocalFields Jion; - double factorJion_0 = au_to_mec2 * EC_to_au*EC_to_au * invdt; - - int nparts = Epart->size()/3; - double *Ex = &( ( *Epart )[0*nparts] ); - double *Ey = &( ( *Epart )[1*nparts] ); - double *Ez = &( ( *Epart )[2*nparts] ); - - for( unsigned int ipart=ipart_min ; ipartcharge( ipart ) ); - - // If ion already fully ionized then skip - if( Z==atomic_number_ ) { - continue; - } - - // Absolute value of the electric field normalized in atomic units - E = EC_to_au * sqrt( pow( *( Ex+ipart-ipart_ref ), 2 ) - +pow( *( Ey+ipart-ipart_ref ), 2 ) - +pow( *( Ez+ipart-ipart_ref ), 2 ) ); - if( E<1e-10 ) { - continue; - } - - // -------------------------------- - // Start of the Monte-Carlo routine - // -------------------------------- - - invE = 1./E; - factorJion = factorJion_0 * invE*invE; - delta = gamma_tunnel[Z]*invE; - ran_p = patch->rand_->uniform(); - IonizRate_tunnel[Z] = beta_tunnel[Z] * exp( -delta*one_third + alpha_tunnel[Z]*log( delta ) ); - - // Total ionization potential (used to compute the ionization current) - TotalIonizPot = 0.0; - - // k_times will give the nb of ionization events - k_times = 0; - Zp1=Z+1; - - if( Zp1 == atomic_number_ ) { - // if ionization of the last electron: single ionization - // ----------------------------------------------------- - if( ran_p < 1.0 -exp( -IonizRate_tunnel[Z]*dt ) ) { - TotalIonizPot += Potential[Z]; - k_times = 1; - } - - } else { - // else : multiple ionization can occur in one time-step - // partial & final ionization are decoupled (see Nuter Phys. Plasmas) - // ------------------------------------------------------------------------- - - // initialization - Mult = 1.0; - Dnom_tunnel[0]=1.0; - Pint_tunnel = exp( -IonizRate_tunnel[Z]*dt ); // cummulative prob. - - //multiple ionization loop while Pint_tunnel < ran_p and still partial ionization - while( ( Pint_tunnel < ran_p ) and ( k_times < atomic_number_-Zp1 ) ) { - newZ = Zp1+k_times; - delta = gamma_tunnel[newZ]*invE; - IonizRate_tunnel[newZ] = beta_tunnel[newZ] - * exp( -delta*one_third+alpha_tunnel[newZ]*log( delta ) ); - D_sum = 0.0; - P_sum = 0.0; - Mult *= IonizRate_tunnel[Z+k_times]; - for( unsigned int i=0; iran_p ) && ( k_times==atomic_number_-Zp1 ) ) { - TotalIonizPot += Potential[atomic_number_-1]; - k_times++; - } - }//END Multiple ionization routine - - // Compute ionization current - if (b_Jx != NULL){ // For the moment ionization current is not accounted for in AM geometry - factorJion *= TotalIonizPot; - Jion.x = factorJion * *( Ex+ipart ); - Jion.y = factorJion * *( Ey+ipart ); - Jion.z = factorJion * *( Ez+ipart ); - - Proj->ionizationCurrentsForTasks( b_Jx, b_Jy, b_Jz, *particles, ipart, Jion, bin_shift ); - } - - // Creation of the new electrons - // (variable weights are used) - // ----------------------------- - if( k_times !=0 ) { - new_electrons_per_bin[ibin].createParticle(); - int idNew = new_electrons_per_bin[ibin].size() - 1;//cout<<"ibin "<position( i, ipart ); - } - for( unsigned int i=0; i<3; i++ ) { - new_electrons_per_bin[ibin].momentum( i, idNew ) = particles->momentum( i, ipart )*ionized_species_invmass; - } - new_electrons_per_bin[ibin].weight( idNew )=double( k_times )*particles->weight( ipart ); - new_electrons_per_bin[ibin].charge( idNew )=-1; - - if( save_ion_charge_ ) { - ion_charge_per_bin_[ibin].push_back( particles->charge( ipart ) ); - } - - // // Increase the charge of the particle - particles->charge( ipart ) += k_times; - } - - } // Loop on particles + DEBUG("Finished Creating the Tunnel Ionizaton class"); } diff --git a/src/Ionization/IonizationTunnel.h b/src/Ionization/IonizationTunnel.h old mode 100755 new mode 100644 index 788b45ccd..0d483c616 --- a/src/Ionization/IonizationTunnel.h +++ b/src/Ionization/IonizationTunnel.h @@ -2,35 +2,354 @@ #define IONIZATIONTUNNEL_H #include - +#include #include #include "Ionization.h" +#include "Particles.h" +#include "Species.h" #include "Tools.h" class Particles; -//! calculate the particle tunnel ionization +template class IonizationTunnel : public Ionization { + public: + inline IonizationTunnel(Params ¶ms, Species *species); + + inline void operator()(Particles *, unsigned int, unsigned int, std::vector *, Patch *, Projector *, + int ipart_ref = 0) override; -public: - //! Constructor for IonizationTunnel: with no input argument - IonizationTunnel( Params ¶ms, Species *species ); - - //! apply the Tunnel Ionization model to the species (with ionization current) - void operator()( Particles *, unsigned int, unsigned int, std::vector *, Patch *, Projector *, int ipart_ref = 0 ) override; //! method for tunnel ionization with tasks - void ionizationTunnelWithTasks( Particles *, unsigned int, unsigned int, std::vector *, Patch *, Projector *, int, int, double *b_Jx, double *b_Jy, double *b_Jz, int ipart_ref = 0 ) override; - -private: - unsigned int atomic_number_; - std::vector Potential; - std::vector Azimuthal_quantum_number; - + inline void ionizationTunnelWithTasks(Particles *, unsigned int, unsigned int, std::vector *, Patch *, Projector *, + int, int, double *b_Jx, double *b_Jy, double *b_Jz, int ipart_ref = 0) override; + + private: + inline double ionizationRate1(const int Z, const double E); + inline double ionizationRate2(const int Z, const double E); + double one_third; + unsigned int atomic_number_; + std::vector Potential, Azimuthal_quantum_number; std::vector alpha_tunnel, beta_tunnel, gamma_tunnel; + + // To be conditionally prepared + // FullPPT + std::vector Magnetic_quantum_number; + + // Tong&Ling + double ionization_tl_parameter_; + std::vector lambda_tunnel; + + // BSI + const double au_to_eV = 27.2116; + const double IH = 13.598434005136; + int rate_formula; }; +template +inline void IonizationTunnel::operator()(Particles *particles, unsigned int ipart_min, unsigned int ipart_max, + vector *Epart, Patch *patch, Projector *Proj, int ipart_ref) +{ + unsigned int Z, Zp1, newZ, k_times; + double TotalIonizPot, E, invE, factorJion, ran_p, Mult, D_sum, P_sum, Pint_tunnel; + vector IonizRate_tunnel(atomic_number_), Dnom_tunnel(atomic_number_, 0.); + LocalFields Jion; + double factorJion_0 = au_to_mec2 * EC_to_au * EC_to_au * invdt; + + int nparts = Epart->size() / 3; + double *Ex = &((*Epart)[0 * nparts]); + double *Ey = &((*Epart)[1 * nparts]); + double *Ez = &((*Epart)[2 * nparts]); + + for (unsigned int ipart = ipart_min; ipart < ipart_max; ipart++) { + // Current charge state of the ion + Z = (unsigned int)(particles->charge(ipart)); + + // If ion already fully ionized then skip + if (Z == atomic_number_) { + continue; + } + + // Absolute value of the electric field normalized in atomic units + E = EC_to_au * sqrt(pow(*(Ex + ipart - ipart_ref), 2) + pow(*(Ey + ipart - ipart_ref), 2) + + pow(*(Ez + ipart - ipart_ref), 2)); + if (E < 1e-10) { + continue; + } + + // -------------------------------- + // Start of the Monte-Carlo routine + // -------------------------------- + + invE = 1. / E; + factorJion = factorJion_0 * invE * invE; + ran_p = patch->rand_->uniform(); + IonizRate_tunnel[Z] = ionizationRate1(Z, E); + + // Total ionization potential (used to compute the ionization current) + TotalIonizPot = 0.0; + + // k_times will give the nb of ionization events + k_times = 0; + Zp1 = Z + 1; + + if (Zp1 == atomic_number_) { + // if ionization of the last electron: single ionization + // ----------------------------------------------------- + if (ran_p < 1.0 - exp(-IonizRate_tunnel[Z] * dt)) { + TotalIonizPot += Potential[Z]; + k_times = 1; + } + + } else { + // else : multiple ionization can occur in one time-step + // partial & final ionization are decoupled (see Nuter Phys. Plasmas) + // ------------------------------------------------------------------------- + + // initialization + Mult = 1.0; + Dnom_tunnel[0] = 1.0; + Pint_tunnel = exp(-IonizRate_tunnel[Z] * dt); // cummulative prob. + + // multiple ionization loop while Pint_tunnel < ran_p and still partial ionization + while ((Pint_tunnel < ran_p) and (k_times < atomic_number_ - Zp1)) { + newZ = Zp1 + k_times; + IonizRate_tunnel[newZ] = ionizationRate2(newZ, E); + D_sum = 0.0; + P_sum = 0.0; + Mult *= IonizRate_tunnel[Z + k_times]; + for (unsigned int i = 0; i < k_times + 1; i++) { + Dnom_tunnel[i] = Dnom_tunnel[i] / (IonizRate_tunnel[newZ] - IonizRate_tunnel[Z + i]); + D_sum += Dnom_tunnel[i]; + P_sum += exp(-IonizRate_tunnel[Z + i] * dt) * Dnom_tunnel[i]; + } + Dnom_tunnel[k_times + 1] -= D_sum; // bug fix + P_sum = P_sum + Dnom_tunnel[k_times + 1] * exp(-IonizRate_tunnel[newZ] * dt); + Pint_tunnel = Pint_tunnel + P_sum * Mult; + + TotalIonizPot += Potential[Z + k_times]; + k_times++; + } // END while + + // final ionization (of last electron) + if (((1.0 - Pint_tunnel) > ran_p) && (k_times == atomic_number_ - Zp1)) { + TotalIonizPot += Potential[atomic_number_ - 1]; + k_times++; + } + } // END Multiple ionization routine + + // Compute ionization current + if (patch->EMfields->Jx_ != NULL) { // For the moment ionization current is not accounted for in AM geometry + factorJion *= TotalIonizPot; + Jion.x = factorJion * *(Ex + ipart); + Jion.y = factorJion * *(Ey + ipart); + Jion.z = factorJion * *(Ez + ipart); + + Proj->ionizationCurrents(patch->EMfields->Jx_, patch->EMfields->Jy_, patch->EMfields->Jz_, *particles, ipart, Jion); + } + + // Creation of the new electrons + // (variable weights are used) + // ----------------------------- + + if (k_times != 0) { + new_electrons.createParticle(); + int idNew = new_electrons.size() - 1; + for (unsigned int i = 0; i < new_electrons.dimension(); i++) { + new_electrons.position(i, idNew) = particles->position(i, ipart); + } + for (unsigned int i = 0; i < 3; i++) { + new_electrons.momentum(i, idNew) = particles->momentum(i, ipart) * ionized_species_invmass; + } + new_electrons.weight(idNew) = double(k_times) * particles->weight(ipart); + new_electrons.charge(idNew) = -1; + + if (save_ion_charge_) { + ion_charge_.push_back(particles->charge(ipart)); + } + + // Increase the charge of the particle + particles->charge(ipart) += k_times; + } + + } // Loop on particles +} + +template +inline double IonizationTunnel::ionizationRate1(const int Z, const double E) +{ + double delta = gamma_tunnel[Z] / E; + return beta_tunnel[Z] * exp(-delta * one_third + alpha_tunnel[Z] * log(delta)); +} + +template +inline double IonizationTunnel::ionizationRate2(const int Z, const double E) +{ + double delta = gamma_tunnel[Z] / E; + return beta_tunnel[Z] * exp(-delta * one_third + alpha_tunnel[Z] * log(delta)); +} + +// IonizationTunnel : 0 +template <> +IonizationTunnel<0>::IonizationTunnel(Params ¶ms, Species *species); + +// IonizationTunnelFullPPT: 1 +template <> +IonizationTunnel<1>::IonizationTunnel(Params ¶ms, Species *species); + +// Tong&Ling: 2 +template <> +IonizationTunnel<2>::IonizationTunnel(Params ¶ms, Species *species); + +template <> +double IonizationTunnel<2>::ionizationRate1(const int Z, const double E); + +template <> +double IonizationTunnel<2>::ionizationRate2(const int Z, const double E); + +// BSI: 3 +template <> +IonizationTunnel<3>::IonizationTunnel(Params ¶ms, Species *species); + +template <> +double IonizationTunnel<3>::ionizationRate1(const int Z, const double E); + +template <> +double IonizationTunnel<3>::ionizationRate2(const int Z, const double E); + +template +void IonizationTunnel::ionizationTunnelWithTasks(Particles *particles, unsigned int ipart_min, + unsigned int ipart_max, vector *Epart, Patch *patch, + Projector *Proj, int ibin, int bin_shift, double *b_Jx, + double *b_Jy, double *b_Jz, int ipart_ref) +{ + unsigned int Z, Zp1, newZ, k_times; + double TotalIonizPot, E, invE, factorJion, delta, ran_p, Mult, D_sum, P_sum, Pint_tunnel; + vector IonizRate_tunnel(atomic_number_), Dnom_tunnel(atomic_number_); + LocalFields Jion; + double factorJion_0 = au_to_mec2 * EC_to_au * EC_to_au * invdt; + + int nparts = Epart->size() / 3; + double *Ex = &((*Epart)[0 * nparts]); + double *Ey = &((*Epart)[1 * nparts]); + double *Ez = &((*Epart)[2 * nparts]); + + for (unsigned int ipart = ipart_min; ipart < ipart_max; ipart++) { + // Current charge state of the ion + Z = (unsigned int)(particles->charge(ipart)); + + // If ion already fully ionized then skip + if (Z == atomic_number_) { + continue; + } + + // Absolute value of the electric field normalized in atomic units + E = EC_to_au * sqrt(pow(*(Ex + ipart - ipart_ref), 2) + pow(*(Ey + ipart - ipart_ref), 2) + + pow(*(Ez + ipart - ipart_ref), 2)); + if (E < 1e-10) { + continue; + } + + // -------------------------------- + // Start of the Monte-Carlo routine + // -------------------------------- + + invE = 1. / E; + factorJion = factorJion_0 * invE * invE; + delta = gamma_tunnel[Z] * invE; + ran_p = patch->rand_->uniform(); + IonizRate_tunnel[Z] = beta_tunnel[Z] * exp(-delta * one_third + alpha_tunnel[Z] * log(delta)); + + // Total ionization potential (used to compute the ionization current) + TotalIonizPot = 0.0; + + // k_times will give the nb of ionization events + k_times = 0; + Zp1 = Z + 1; + + if (Zp1 == atomic_number_) { + // if ionization of the last electron: single ionization + // ----------------------------------------------------- + if (ran_p < 1.0 - exp(-IonizRate_tunnel[Z] * dt)) { + TotalIonizPot += Potential[Z]; + k_times = 1; + } + + } else { + // else : multiple ionization can occur in one time-step + // partial & final ionization are decoupled (see Nuter Phys. Plasmas) + // ------------------------------------------------------------------------- + + // initialization + Mult = 1.0; + Dnom_tunnel[0] = 1.0; + Pint_tunnel = exp(-IonizRate_tunnel[Z] * dt); // cummulative prob. + + // multiple ionization loop while Pint_tunnel < ran_p and still partial ionization + while ((Pint_tunnel < ran_p) and (k_times < atomic_number_ - Zp1)) { + newZ = Zp1 + k_times; + delta = gamma_tunnel[newZ] * invE; + IonizRate_tunnel[newZ] = beta_tunnel[newZ] * exp(-delta * one_third + alpha_tunnel[newZ] * log(delta)); + D_sum = 0.0; + P_sum = 0.0; + Mult *= IonizRate_tunnel[Z + k_times]; + for (unsigned int i = 0; i < k_times + 1; i++) { + Dnom_tunnel[i] = Dnom_tunnel[i] / (IonizRate_tunnel[newZ] - IonizRate_tunnel[Z + i]); + D_sum += Dnom_tunnel[i]; + P_sum += exp(-IonizRate_tunnel[Z + i] * dt) * Dnom_tunnel[i]; + } + Dnom_tunnel[k_times + 1] -= D_sum; + P_sum = P_sum + Dnom_tunnel[k_times + 1] * exp(-IonizRate_tunnel[newZ] * dt); + Pint_tunnel = Pint_tunnel + P_sum * Mult; + + TotalIonizPot += Potential[Z + k_times]; + k_times++; + } // END while + + // final ionization (of last electron) + if (((1.0 - Pint_tunnel) > ran_p) && (k_times == atomic_number_ - Zp1)) { + TotalIonizPot += Potential[atomic_number_ - 1]; + k_times++; + } + } // END Multiple ionization routine + + // Compute ionization current + if (b_Jx != NULL) { // For the moment ionization current is not accounted for in AM geometry + factorJion *= TotalIonizPot; + Jion.x = factorJion * *(Ex + ipart); + Jion.y = factorJion * *(Ey + ipart); + Jion.z = factorJion * *(Ez + ipart); + + Proj->ionizationCurrentsForTasks(b_Jx, b_Jy, b_Jz, *particles, ipart, Jion, bin_shift); + } + + // Creation of the new electrons + // (variable weights are used) + // ----------------------------- + if (k_times != 0) { + new_electrons_per_bin[ibin].createParticle(); + int idNew = new_electrons_per_bin[ibin].size() - + 1; // cout<<"ibin "<position(i, ipart); + } + for (unsigned int i = 0; i < 3; i++) { + new_electrons_per_bin[ibin].momentum(i, idNew) = particles->momentum(i, ipart) * ionized_species_invmass; + } + new_electrons_per_bin[ibin].weight(idNew) = double(k_times) * particles->weight(ipart); + new_electrons_per_bin[ibin].charge(idNew) = -1; + + if (save_ion_charge_) { + ion_charge_per_bin_[ibin].push_back(particles->charge(ipart)); + } + + // // Increase the charge of the particle + particles->charge(ipart) += k_times; + } + + } // Loop on particles +} #endif diff --git a/src/Python/pyinit.py b/src/Python/pyinit.py index f5aeeb7e1..79c4797df 100755 --- a/src/Python/pyinit.py +++ b/src/Python/pyinit.py @@ -430,6 +430,7 @@ class Species(SmileiComponent): ionization_rate = None atomic_number = None maximum_charge_state = 0 + ionization_tl_parameter = 6 is_test = False relativistic_field_initialization = False keep_interpolated_fields = [] diff --git a/src/Species/Species.h b/src/Species/Species.h index d4af3bf9d..8796e3f0d 100755 --- a/src/Species/Species.h +++ b/src/Species/Species.h @@ -79,6 +79,9 @@ class Species //! maximum charge state unsigned int maximum_charge_state_; + //! alpha parameter in the Tong-Lin ionization model + double ionization_tl_parameter_; + //! user defined ionization rate profile PyObject *ionization_rate_; diff --git a/src/Species/SpeciesFactory.h b/src/Species/SpeciesFactory.h index edef5e334..8d5476d8a 100755 --- a/src/Species/SpeciesFactory.h +++ b/src/Species/SpeciesFactory.h @@ -739,6 +739,9 @@ class SpeciesFactory this_species->maximum_charge_state_ = 0; PyTools::extract( "maximum_charge_state", this_species->maximum_charge_state_, "Species", ispec); + + this_species->ionization_tl_parameter_ = 6; + PyTools::extract( "ionization_tl_parameter", this_species->ionization_tl_parameter_, "Species", ispec); std::string model; PyTools::extract( "ionization_model", model, "Species", ispec ); @@ -760,7 +763,7 @@ class SpeciesFactory LINK_NAMELIST + std::string("#species") ); } - if( model == "tunnel" ){ + if( (model == "tunnel") || (model == "tunnel_BSI") || (model == "tunnel_TL") || (model == "tunnel_full_PPT") ){ if (params.Laser_Envelope_model){ ERROR_NAMELIST("An envelope is present, so tunnel_envelope or tunnel_envelope_averaged ionization model should be selected for species "<