# Source code for irfpy.vima.massring

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

''' Load the table 4.  This is constant over the mission. V6.2, 04 Dec 2006.

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

''' 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.
'''
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 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