''' Calculate mass ring number from energy and mass in Vex IMA Bible.
The way of calculating the mass ring number from energy and mass
is described in the ImaBible.pdf file.
This module implements that.
.. autosummary::
rm
dm
massline
Users may use :func:`rm` function to calculate the mass ring number
from energy and mass together with PACC index.
The mass ring number looks to have 0-based index.
It may say that the range of mass ring number could be (-0.5, 31.5).
>>> print('{:.2f}'.format(rm(100, 32, 3))) # Mass ring number for 100 eV/q, 32 amu/q particle with PACC=3.
8.49
:func:`dm` will give the width, although it is almost always 1.2.
Users also want to know the iso-mass line. This will be done by
:func:`massline` function.
>>> energies, masses = massline(32, 6, enestep=np.logspace(1, 4, 96))
>>> print('{e:.2f} {m:.2f}'.format(e=energies[15], m=masses[15]))
29.76 7.03
See also the example in :mod:`snippet_mima.mima_masscontour`.
*Details*: The :attr:`table` contains the table needed for calculation.
It is 8-element list, but only used for 0, 3, and 6th element.
Each corresponds to the table for each PACC.
The table is (5, 12) shaped. Contents is in :func:`_load_table4` function.
>>> print(_table[0].shape)
(5, 12)
>>> print(_table[4])
None
'''
import numpy as np
from scipy.interpolate import InterpolatedUnivariateSpline
def _load_table4():
''' Load the table 4. This is constant over the mission. V6.2, 04 Dec 2006.
>>> tbl = _load_table4()
>>> print(tbl.shape)
(15, 12)
:returns: A table. (15, 12 shaped).
[0:5, :] is for low pacc, [5:10, :] is for mid pac,
and [10:15, :] is for high pac.
For column, refer to the table below.
Column
::
0: PaccIndex
1: Pacc
2: Elimit
3: #
4: GfitP0
5: GfitP1
6: GfitD0
7: GfitD1
8: Kpacc
9: Kmass
10: Dmass
11: Omass
'''
tbl = np.array([[ 0.00000000e+00, 3.90000000e+01, 3.01000000e+02,
0.00000000e+00, 8.32400000e-04, 5.30200000e-01,
7.90500000e-05, 2.61100000e-01, 1.00000000e+00,
1.20000000e+00, -5.00000000e+01, 1.00000000e+00],
[ 0.00000000e+00, 3.90000000e+01, 3.01000000e+02,
1.00000000e+00, 1.09500000e-02, 2.10600000e+00,
1.04000000e-03, 2.79700000e-01, 9.30100000e-16,
2.20000000e+00, -3.00000000e+01, 2.00000000e+00],
[ 0.00000000e+00, 3.90000000e+01, 3.01000000e+02,
2.00000000e+00, 1.44000000e-01, -8.26800000e-03,
1.36800000e-02, -1.07300000e-02, -8.17500000e-16,
1.70000000e+01, -3.00000000e-01, 1.40000000e+01],
[ 0.00000000e+00, 3.90000000e+01, 3.01000000e+02,
3.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
1.90000000e+01, 0.00000000e+00, 1.60000000e+01],
[ 0.00000000e+00, 3.90000000e+01, 3.01000000e+02,
4.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
3.00000000e+01, 9.00000000e+00, 3.20000000e+01],
[ 3.00000000e+00, 1.90200000e+03, 3.01000000e+02,
0.00000000e+00, 1.11400000e-02, -6.66000000e-03,
2.08400000e-03, 8.74700000e-01, 9.65600000e-01,
1.00000000e+00, -5.00000000e+01, 1.00000000e+00],
[ 3.00000000e+00, 1.90200000e+03, 3.01000000e+02,
1.00000000e+00, 6.34300000e-02, 2.14900000e+00,
1.18600000e-02, 2.87800000e-01, 1.20900000e-01,
2.00000000e+00, -3.00000000e+01, 2.00000000e+00],
[ 3.00000000e+00, 1.90200000e+03, 3.01000000e+02,
2.00000000e+00, 3.61200000e-01, -1.80000000e-02,
6.75600000e-02, -1.40900000e-02, -8.73700000e-02,
1.40000000e+01, -0.00000000e+00, 1.40000000e+01],
[ 3.00000000e+00, 1.90200000e+03, 3.01000000e+02,
3.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
1.60000000e+01, 0.00000000e+00, 1.60000000e+01],
[ 3.00000000e+00, 1.90200000e+03, 3.01000000e+02,
4.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
3.20000000e+01, 9.00000000e+00, 3.20000000e+01],
[ 6.00000000e+00, 3.61500000e+03, 1.28000000e+02,
0.00000000e+00, 2.80100000e+01, 2.85800000e+00,
-5.48100000e-01, 1.11500000e+00, 4.09400000e-02,
4.00000000e-01, -5.00000000e+01, 1.00000000e+00],
[ 6.00000000e+00, 3.61500000e+03, 1.28000000e+02,
1.00000000e+00, -1.63100000e+00, 1.02600000e+00,
2.43600000e-01, 1.15300000e-01, 4.86700000e-01,
1.00000000e+00, -3.00000000e+01, 2.00000000e+00],
[ 6.00000000e+00, 3.61500000e+03, 1.28000000e+02,
2.00000000e+00, 8.10300000e-03, -5.92400000e-03,
-4.88500000e-03, -2.61400000e-03, -1.12700000e-01,
3.00000000e+01, -3.00000000e-01, 1.40000000e+01],
[ 6.00000000e+00, 3.61500000e+03, 1.28000000e+02,
3.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
3.40000000e+01, 0.00000000e+00, 1.60000000e+01],
[ 6.00000000e+00, 3.61500000e+03, 1.28000000e+02,
4.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
8.20000000e+01, 9.00000000e+00, 3.20000000e+01]])
return tbl
def _load_table4_pacc(pi):
''' Return the table 4 contents according to PACC.
:param pi: Pacc index. 0, 3 or 6.
:returns: (5, 12) shaped array. See :func:`load_table4` for contents.
'''
tbl = _load_table4()
if pi == 0:
tbl = tbl[:5, :]
elif pi == 3:
tbl = tbl[5:10, :]
elif pi == 6:
tbl = tbl[10:, :]
else:
raise ValueError('PI={} not used.'.format(pi))
return tbl
_table = [_load_table4_pacc(0), None, None,
_load_table4_pacc(3), None, None,
_load_table4_pacc(6), None]
''' Table for Mass Ring calculation.
An list of numpy array with a shape of (5,12).
This is eight elements, each corresponds to PACC index.
>>> print(table[0].shape)
(5, 12)
>>> print(table[4])
None
'''
def _meff(m_per_q, pi):
''' Calculate Meff.
Meff = Lin(M/Q, Kmass_i(PI), Omass_i(PI)
Here Lin(X, Ytab, Xtab) is the linear interpolation of X on the tabulated Xtab and Ytab.
:param m_per_q: Mass per charge. Float or array.
:param pi: Pacc index. 0, 3 or 6.
:returns: Meff value
:rtype: np.array (that's why one need float() for the later examples)
>>> print(_meff(1, 0), _meff(1, 3), _meff(1, 6))
1.2 1.0 0.4
>>> print(_meff(15, 0), _meff(15, 3), _meff(15, 6))
18.0 15.0 32.0
>>> print(_meff(48, 0), _meff(48, 3), _meff(48, 6))
41.0 48.0 130.0
'''
tbl = _table[pi]
kmass = tbl[:, 9]
omass = tbl[:, 11]
f = InterpolatedUnivariateSpline(omass, kmass, k=1) # k=1 means linear
return_val = f(m_per_q)
return return_val
def _pacceff(m_per_q, pi):
r''' Calculate Pacceff.
.. math::
Pacc_{eff} = Pacc(PI) \cdot (Kpacc0(PI) + Kpacc1(PI) / Meff + Kpacc2(PI) / Meff ^ 2)
>>> peff = _pacceff(12.8, 0) # M/q = 12.8 for PACC=0
>>> print('{:.2f}'.format(float(peff)))
39.00
>>> peff = _pacceff(37.6, 3) # M/q = 37.6 for PACC=3
>>> print('{:.2f}'.format(float(peff)))
1842.57
>>> peff = _pacceff(3.75, 6) # M/q = 3.75 for PACC=6
>>> print('{:.2f}'.format(float(peff)))
469.56
'''
tbl = _table[pi]
pacc = tbl[0, 1] # PACC in the table (tbl[1, :] is all same)
kpacc0 = tbl[0, 8]
kpacc1 = tbl[1, 8]
kpacc2 = tbl[2, 8]
Meff = _meff(m_per_q, pi)
coeffic = kpacc0 + kpacc1 / Meff + kpacc2 / Meff / Meff
retval = pacc * coeffic
return retval
def _geff(e_per_q, m_per_q, pi):
r''' Geff calculation.
.. math::
G_{eff} = 10^3 / \sqrt{(E/Q + Pacc_{eff})\cdot M_{eff}}
:param e_per_q: Given by the user.
Technical note shows it as the function of energy index,
it is complicated to get the energy from the index,
so that this is specified by the user.
Pacceff and Meff is obtained from the other functions (
:func:`_pacceff` and :func:`_meff`).
>>> Geff = _geff(1200, 12.8, 0) # M/q = 12.8 for PACC=0
>>> print('{:.4f}'.format(float(Geff)))
7.2114
>>> Geff = _geff(2400, 37.6, 3) # M/q = 37.6 for PACC=3
>>> print('{:.4f}'.format(float(Geff)))
2.5038
>>> Geff = _geff(750, 3.75, 6) # M/q = 3.75 for PACC=6
>>> print('{:.4f}'.format(float(Geff)))
12.5222
'''
PaccEff = _pacceff(m_per_q, pi)
Meff = _meff(m_per_q, pi)
geff = 1e3 / np.sqrt((e_per_q + PaccEff) * Meff)
return geff
def _glim(m_per_q, pi):
r''' Glim
.. math::
G_{lim} = 10^3 / \sqrt{(Elimit(PI) + Pacc_{eff})\cdot(M_{eff})}
>>> Glim = _glim(12.8, 0) # M/q = 12.8 for PACC=0
>>> print('{:.3f}'.format(float(Glim)))
13.766
>>> Glim = _glim(37.6, 3) # M/q = 37.6 for PACC=3
>>> print('{:.3f}'.format(float(Glim)))
3.522
>>> Glim = _glim(3.75, 6) # M/q = 3.75 for PACC=6
>>> print('{:.3f}'.format(float(Glim)))
17.889
'''
## Elimit is obtained from table [0, 2]
tbl = _table[pi]
elimit = tbl[0, 2]
Meff = _meff(m_per_q, pi)
PaccEff = _pacceff(m_per_q, pi)
glim = 1e3 / np.sqrt((elimit + PaccEff) * Meff)
return glim
[docs]def rm(e_per_q, m_per_q, pi):
''' Return Rm.
>>> Rm = rm(1200, 12.8, 0) # M/q = 12.8 for PACC=0
:param e_per_q: Energy per charge (in eV/q)
:param m_per_q: Mass per charge (in amu/q)
:param pi: Pacc index, 0, 3 or 6.
:returns: Mass ring number.
>>> print('{:.2f}'.format(float(Rm)))
15.29
>>> Rm = rm(2400, 37.6, 3) # M/q = 37.6 for PACC=3
>>> print('{:.2f}'.format(float(Rm)))
5.26
>>> Rm = rm(750, 3.75, 6) # M/q = 3.75 for PACC=6
>>> print('{:.2f}'.format(float(Rm)))
14.78
>>> Rm = rm(100, 12.8, 0) # M/q = 12.8 for PACC=0
>>> print('{:.2f}'.format(float(Rm)))
67.50
>>> Rm = rm(250, 37.6, 3) # M/q = 37.6 for PACC=3
>>> print('{:.2f}'.format(float(Rm)))
7.45
>>> Rm = rm(350, 3.75, 6) # M/q = 3.75 for PACC=6
>>> print('{:.2f}'.format(float(Rm)))
17.15
'''
tbl = _table[pi]
elimit = tbl[0, 2]
if e_per_q < elimit:
return _rm_low(e_per_q, m_per_q, pi)
else:
return _rm_high(e_per_q, m_per_q, pi)
[docs]def massline(m_per_q, pi, enestep=None):
""" Return the iso-mass contour.
For the given mass and post acceleration,
the mass ring numbers as a function of energies are
returned.
:param m_per_q: Mass per charge (amu/q)
:param pi: Post acceleration index. 0, 3, or 6.
:keyword enestep: Energy step. Default is logspace(0, 4, 96).
:return: (enestep, massring) with each element np.array
with shape of given enestep.
"""
if enestep is None:
enestep = np.logspace(0, 4, 96)
ml = np.zeros_like(enestep)
for i in range(len(enestep)):
ml[i] = rm(enestep[i], m_per_q, pi)
return np.array(enestep), np.array(ml)
[docs]def dm(e_per_q, m_per_q, pi):
''' Return Dm.
:param e_per_q: Energy per charge (in eV/q)
:param m_per_q: Mass per charge (in amu/q)
:param pi: Pacc index, 0, 3 or 6.
:returns: Mass ring number.
>>> Dm = dm(1200, 12.8, 0) # M/q = 12.8 for PACC=0
>>> print('{:.2f}'.format(float(Dm)))
1.72
>>> Dm = dm(2400, 37.6, 3) # M/q = 37.6 for PACC=3
>>> print('{:.2f}'.format(float(Dm)))
1.51
>>> Dm = dm(750, 3.75, 6) # M/q = 3.75 for PACC=6
>>> print('{:.2f}'.format(float(Dm)))
2.15
>>> Dm = dm(100, 12.8, 0) # M/q = 12.8 for PACC=0
>>> print('{:.2f}'.format(float(Dm)))
6.36
>>> Dm = dm(250, 37.6, 3) # M/q = 37.6 for PACC=3
>>> print('{:.2f}'.format(float(Dm)))
0.90
>>> Dm = dm(350, 3.75, 6) # M/q = 3.75 for PACC=6
>>> print('{:.2f}'.format(float(Dm)))
2.27
'''
tbl = _table[pi]
elimit = tbl[0, 2]
if e_per_q < elimit:
return _dm_low(e_per_q, m_per_q, pi)
else:
return _dm_high(e_per_q, m_per_q, pi)
def _rm_high(e_per_q, m_per_q, pi):
r''' Return Rm when the energy is higher than limit.
.. math::
Rm = GfitP1_0(PI) + GfitP1_i(PI)\cdot G_{eff} + GfitP1_2(PI)\cdot G_{eff}^2
'''
tbl = _table[pi] # (5, 12) shape
gfp10 = tbl[0, 5]
gfp11 = tbl[1, 5]
gfp12 = tbl[2, 5]
Geff = _geff(e_per_q, m_per_q, pi)
Rm = gfp10 + gfp11 * Geff + gfp12 * Geff * Geff
return Rm
def _rm_low(e_per_q, m_per_q, pi):
r''' Return Rm for E<Elimit
.. math::
dM = GfitP0_0(PI) + GfitP0_1(PI)\cdot G_{lim} + GfitP0_2(PI)\cdot G_{lim}^2 - (GfitP1_0(PI) + GfitP1_1(PI)\cdot G_{lim} + GfitP1_2(PI)\cdot G_{lim}^2 ) \\
Rm = GfitP0_0(PI) + GfitP0_1(PI)\cdot G_{eff} + GfitP0_2(PI)\cdot G_{eff}^2 - dM
'''
tbl = _load_table4_pacc(pi) # (5, 12) array
gfp00 = tbl[0, 4]
gfp01 = tbl[1, 4]
gfp02 = tbl[2, 4]
gfp10 = tbl[0, 5]
gfp11 = tbl[1, 5]
gfp12 = tbl[2, 5]
Glim = _glim(m_per_q, pi)
Geff = _geff(e_per_q, m_per_q, pi)
dM = gfp00 + gfp01 * Glim + gfp02 * Glim * Glim - (gfp10 + gfp11 * Glim + gfp12 * Glim * Glim)
Rm = gfp00 + gfp01 * Geff + gfp02 * Geff * Geff - dM
return Rm
def _dm_high(e_per_q, m_per_q, pi):
r''' Dm for high energy part.
.. math::
Dm = GfitD1_0(PI) + GfitD1_1(PI)\cdot G_{eff} + GfitD1_2(PI)\cdot G_{eff}^2
'''
tbl = _table[pi]
gfd10 = tbl[0, 7]
gfd11 = tbl[1, 7]
gfd12 = tbl[2, 7]
Geff = _geff(e_per_q, m_per_q, pi)
dm = gfd10 + gfd11 * Geff + gfd12 * Geff * Geff
return dm
def _dm_low(e_per_q, m_per_q, pi):
r''' Dm for low energy part.
.. math::
Dm = GfitD0_0(PI) + GfitD0_1(PI)\cdot G_{eff} + GfitD0_2(PI)\cdot G_{eff}^2
'''
tbl = _table[pi]
gfd10 = tbl[0, 6]
gfd11 = tbl[1, 6]
gfd12 = tbl[2, 6]
Geff = _geff(e_per_q, m_per_q, pi)
dm = gfd10 + gfd11 * Geff + gfd12 * Geff * Geff
return dm