IONBEA| Fit Function | \[\begin{align*} \sigma (pe) = \begin{cases} 0.0 & \text{if } \text{alpha} < 0.207 \\ 5.867 \times 10^{-17} \; z_p^2 n_{eff}^4 \left( \text{alpha} - \frac{0.164}{\text{alpha}^2} + \frac{0.1875}{\text{alpha}^3+\text{alpha}^2} \right) & \text{if } 0.207 \leq \text{alpha} < 1.207 \\ 1.467 \times 10^{-16} \; z_p^2 n_{eff}^4 \frac{1.0 - \frac{0.15}{\text{alpha}^2 - 1.0}}{\text{alpha}^2} &\text{otherwise} \end{cases} \end{align*}\] |
| Comments | Python code requires NumPy imported as `np`. Please refer to the source code for the information on `alpha` and `neff`. |
| Arguments |
|
||||||||||||||||||||||||
| Return values |
|
||||||||||||||||||||||||
| Code | c
c###################################################################
c
subroutine heionen(n, l, mult, sumen, eion, kermsg)
c
c this subroutine passes the ionization energy for excited states
c of helium taken from tables a.4 and a.5 given by janev et al.
c (see doc=h-he-plasma.)
c
c the input subroutine parameters are
c
c n = principal quantum number of the excited electronic state
c
c l = orbital angular momentum quantum number of the exited state
c
c mult = the spin multiplicity (2s+1) of the state
c
c sumen = indicates choice of ionization energy to be returned.
c if sumen=0 the ionization energy for the specific
c state (quantified by n, l amd mult) is returned
c if sumen=1 the ionization energy taken as an average
c over angular momentum and toatal spin is returned
c
c the output subroutine parameters are
c
c eion = ionization energy
c
c kermsg = blank if no errors
c
c written by j. j. smith , iaea atomic and molecular data unit
c
c------------------------------------------------------------------------
c
character*(*) kermsg
integer n, l, mult, sumen, lp, multp
dimension enl(4,4,2), en(7)
data (((enl(i,j,k),i=1,4),j=1,4),k=1,2) /32*0.0/
c
c data enl(1,1,1)
c * / 24.588 /,
c
c---- n=2 excited states
c
c 1 enl(2,1,1),enl(2,1,2),enl(2,2,1),enl(2,2,2)
c * / 3.973, 4.769, 3.371, 3.625 /,
c
c---- n=3 excited states
c
c 2 enl(3,1,1),enl(3,1,2),enl(3,2,1),enl(3,2,2),enl(3,3,1),
c * enl(3,3,2)
c * / 1.699, 1.871, 1.502, 1.582, 1.515, 1.515 /,
c
c---- n=4 excited states
c
c 3 enl(4,1,1),enl(4,1,2),enl(4,2,1),enl(4,2,2),enl(4,3,1),
c * enl(4,3,2),enl(4,4,1),enl(4,4,2)
c * / 0.9155, 0.9951, 0.8470, 0.8812, 0.8529, 0.8530, 0.8520,
c * 0.8520 /
c
c----- ionization energies summed over angular momentum and
c----- total spin
c
data (en(k),k=1,2) /2*0.0/
data en(3),en(4),en(5),en(6),en(7)
* / 1.609 , 0.8811, 0.52, 0.37, 0.29/
c
enl(1,1,1)=24.588
c---- n=2 excited states
enl(2,1,1)=3.973
enl(2,1,2)=4.769
enl(2,2,1)=3.371
enl(2,2,2)=3.625
c---- n=3 excited states
enl(3,1,1)=1.699
enl(3,1,2)=1.871
enl(3,2,1)=1.502
enl(3,2,2)=1.582
enl(3,3,1)=1.515
enl(3,3,2)=1.515
c---- n=4 excited states
enl(4,1,1)=0.9155
enl(4,1,2)=0.9951
enl(4,2,1)=0.8470
enl(4,2,2)=0.8812
enl(4,3,1)=0.8529
enl(4,3,2)=0.8530
enl(4,4,1)=0.8520
enl(4,4,2)=0.8520
kermsg =' '
if (sumen .eq. 1) then
if (n .ge. 8) then
eion =13.58/(n*n)
return
else
eion =en(n)
if (eion .eq. 0.0) kermsg =
* 'ionization energy for n value not in table in heionen'
endif
else
lp = l+1
if (mult .eq. 1) then
multp = 1
else if (mult .eq. 3) then
multp = 2
else
kermsg = 'invalid spin multipclity as input in heionen'
endif
eion = enl(n,lp,multp)
if (eion .eq. 0.0) kermsg =
* 'ionization energy not in table in heionen'
endif
return
end
c###################################################################
c
subroutine ionbea(pe, pcf, kncf, pxs, kermsg)
c
c this is a subroutine to calculate cross sections (cm[2])
c versus energy (ev) for heavy particle impact ionization
c in the binary encounter approximation (bea). this routine is
c only valid for a hydrogenic core.
c
c pe = collision energy in ev
c
c the coefficient data array passed should contain
c
c pcf(1) = ant , the atomic number of target
c
c pcf(2) = zct, the charge of the core of the target
c
c pcf(3) = zp, the charge of the projectile ion
c
c pcf(4) = eion, binding energy of target atom. if this value is
c zero and z0 = 2 (he) the binding energy is determined
c from the routine heionen.
c
c pcf(5) = itrans, integer which defines the type of transition.
c for transitions defined only in terms of the initial
c and final principal quantum numbers (n,m), itrans=1.
c for tranitions between (nl,ml') states itrans=2
c
c pcf(6) = sumen, this indicates choice of ionization energy to be
c returned. if sumen=0 the ionization energy for the
c specific state (quantified by n, l amd mult) is
c returned. if sumen=1 the ionization energy taken as
c an average over angular momentum and total spin is
c returned. this only applies to helium as a target.
c
c pcf(7) = n, the principal quantum number of the initial state
c
c pcf(8) = lin, the orbital angular monentum of the inital state
c
c pcf(9) = mult, spin multiplicity of the initial state
c
c - warning- .
c
c the coefficient array pcf is updated by this routine to
c include energy independent constants. these coefficients can be
c used in subsequent calls for the same entry. the coefficeients
c added are:
c
c pcf(10) = eth,threshold enegry for ionization
c
c pcf(11) = neff, effective value of the principal quantum number
c used. (see doc=h-he-plasma for definition)
c
c kermsg = blank if no errors
c
c pxs = cross section in cm[2]
c
c
c------------------------------------------------------------------------
c
double precision pe, pcf, pxs
c
integer ant, n, zct, zp, itrans, sumen, lin, mult
real neff
dimension pcf(11)
character*(*) kermsg
c
data ry /1.358e+01/
c
ant = pcf(1)
zct = pcf(2)
zp = pcf(3)
eion = pcf(4)
itrans = pcf(5)
sumen = pcf(6)
c
if (itrans .ne. 1 .and. itrans .ne. 2) then
kermsg = ' error invalid fifth coefficient passed to #ionbea'
return
endif
c
if (sumen .ne. 0 .and. sumen .ne. 1) then
kermsg = ' error invalid sixth coefficient passed to #ionbea'
return
endif
c
n = pcf(7)
an = n
kermsg = ' '
if (kncf .le. 9) then
c
c first call to ionbea determine energy independent
c parameters and place in pcf for further use
c
c determine the binding threshold energy (eion)
c
if (eion .gt. 0.0 ) then
enion = eion
neff=n
else if (itrans .eq. 1) then
if (ant .eq. 1) then
enion = ry / (an*an)
neff=n
else if (ant .eq. 2) then
call heionen(n, 0, 0, 1, enion, kermsg)
neff = sqrt(ry/enion)
else
kermsg = 'target not covered by ionbea '
return
endif
c
else if (itrans .eq. 2 .and. ant .eq. 2) then
if (kncf .ne. 9) then
kermsg = ' invalid number of coefficients passed to ionbea'
return
endif
c
lin = pcf(8)
mult = pcf(9)
call heionen(n, lin, mult, sumen, enion, kermsg)
neff = sqrt(ry/enion)
endif
c
c place energy independent parameters in coefficient array and
c update kncf
c
pcf(10) = enion
pcf(11) = neff
kncf = 11
else if (kncf .eq. 11) then
enion = pcf(10)
neff = pcf(11)
c
else
kermsg = ' incorrect number of coefficients passed to ionbea'
return
endif
c
if(pe .lt. enion) then
pxs=0.0
return
endif
c
c determine alpha
c
alpha= (6.3246e-03/zp) * neff * dsqrt (pe)
c
c determine from the value of alpha which formula to use to
c calculate the cross section
c
alpsq = alpha * alpha
if (alpha .lt. 0.207) then
pxs = 0.0
else if (alpha .lt. 1.207) then
pxs = 5.867e-17 * (zp**2) * (neff**4) *
1 ( alpha - (0.164/alpsq) + 0.1875/(alpsq*(alpha+ 1.0)))
else
pxs = 1.467e-16 * (zp**2) * (neff**4) *
1 ( 1.0 - (0.15/(alpsq-1.0)))/alpsq
endif
c
return
end |
| Arguments |
|
||||||||||||||||
| Return values |
|
||||||||||||||||
| Code | def heionen(n, l, mult, sumen):
"""
This function passes the ionization energy for excited states of helium
taken from tables a.4 and a.5 given by janev et al.
n: principal quantum number of the excited electronic state
l: orbital angular momentum quantum number of the exited state
mult: the spin multiplicity (2s+1) of the state
sumen: indicates choice of ionization energy to be returned.
"""
enl = np.zeros((4, 4, 2))
en = np.zeros(7)
# Define the ionization energy values
enl[0, 0, 0] = 24.588
enl[1, 0, 0] = 3.973
enl[1, 0, 1] = 4.769
enl[1, 1, 0] = 3.371
enl[1, 1, 1] = 3.625
enl[2, 0, 0] = 1.699
enl[2, 0, 1] = 1.871
enl[2, 1, 0] = 1.502
enl[2, 1, 1] = 1.582
enl[2, 2, 0] = 1.515
enl[2, 2, 1] = 1.515
enl[3, 0, 0] = 0.9155
enl[3, 0, 1] = 0.9951
enl[3, 1, 0] = 0.8470
enl[3, 1, 1] = 0.8812
enl[3, 2, 0] = 0.8529
enl[3, 2, 1] = 0.8530
enl[3, 3, 0] = 0.8520
enl[3, 3, 1] = 0.8520
en[2:5] = [1.609, 0.8811, 0.52, 0.37, 0.29]
eion = 0.0
if sumen == 1:
if n >= 8:
eion = 13.58 / (n * n)
else:
try:
eion = en[n - 1]
except IndexError:
raise ValueError('Ionization energy for n value not in table in heionen')
else:
lp = l + 1
if mult == 1:
multp = 1
elif mult == 3:
multp = 2
else:
raise ValueError('Invalid spin multiplicity as input in heionen')
try:
eion = enl[n - 1, lp - 1, multp - 1]
except IndexError:
raise ValueError('Ionization energy not in table in heionen')
return eion
def ionbea(pe, pcf, kncf):
"""
This function provides cross section in cm^2 as a function of energy (ev)
for heavy particle impact ionization in the binary encounter approximation (bea).
The function is only valid for a hydrogenic core.
pe= collision energy (eV)
pcf is the coefficient data array, where
pcf[0] = the atomic number of target
pcf[1] = the charge of the core of the target
pcf[2] = the charge of the projectile ion
pcf[3] = binding energy of target atom. if this value is zero and z0 = 2 (he) the binding energy is
determined from the function heionen.
pcf[4] = type of transition
pcf[5] = this indicates choice of ionization energy to be returned.
If 0, the ionization energy for the specific state
(quantified by n, l amd mult) is returned.
if 1, the ionization energy taken as an average over angular
momentum and total spin is returned.
this only applies to helium as a target.
pcf[6] = the principal quantum number of the initial state
pcf[7] = the orbital angular monentum of the inital state
pcf[8] = spin multiplicity of the initial state
pcf[9] = threshold enegry for ionization
pcf[10] = effective value of the principal quantum number
"""
ant = int(pcf[0])
zct = int(pcf[1])
zp = int(pcf[2])
eion = pcf[3]
itrans = int(pcf[4])
sumen = int(pcf[5])
ry = 13.58
if itrans != 1 and itrans != 2:
raise ValueError('Error: Invalid fifth coefficient passed to #ionbea')
if sumen != 0 and sumen != 1:
raise ValueError('Error: Invalid sixth coefficient passed to #ionbea')
n = int(pcf[6])
an = n
if kncf <= 9:
if eion > 0.0:
enion = eion
neff = n
elif itrans == 1:
if ant == 1:
enion = ry / (an * an)
neff = n
elif ant == 2:
# Call heionen(n, 0, 0, 1, enion) here and calculate neff
neff = np.sqrt(ry / enion)
else:
raise ValueError('Target not covered by ionbea')
elif itrans == 2 and ant == 2:
if kncf != 9:
raise ValueError('Invalid number of coefficients passed to ionbea')
lin = int(pcf[7])
mult = int(pcf[8])
# Call heionen(n, lin, mult, sumen, enion) here and calculate neff
neff = np.sqrt(ry / enion)
pcf[9] = enion
pcf[10] = neff
kncf = 11
elif kncf == 11:
enion = pcf[9]
neff = pcf[10]
else:
raise ValueError('Incorrect number of coefficients passed to ionbea')
if pe < enion:
pxs = 0.0
return pxs
alpha = (6.3246e-03 / zp) * neff * np.sqrt(pe)
alpsq = alpha * alpha
if alpha < 0.207:
pxs = 0.0
elif alpha < 1.207:
pxs = 5.867e-17 * (zp ** 2) * (neff ** 4) * (alpha - (0.164 / alpsq) + 0.1875 / (alpsq * (alpha + 1.0)))
else:
pxs = 1.467e-16 * (zp ** 2) * (neff ** 4) * (1.0 - (0.15 / (alpsq - 1.0))) / alpsq
return pxs |