''' Calculation of flux, see Section 10 in :ref:`vexelsbible`.
The flux conversion is made using following functions.
:func:`raw2cps` and :func:`raw2cps_np`
Function to convert the raw count to count-per-second.
:func:`raw2def` and :func:`raw2def_np`
Function to convert the raw count to differential energy flux.
:func:`raw2dnf` and :func:`raw2dnf_np`
Function to convert the raw count to differential number flux.
:func:`raw2psd` and :func:`raw2psd_np`
Function to convert the raw count to phase space density.
A sample use of this module can be found in the :mod:`vels_flux_onecount`.
'''
import logging
import numpy as np
# Followings are in irfpy.vels.bible package.
from irfpy.vels.bible import times
from irfpy.vels.bible.gfactor import gfactor, gfactor1d
[docs]def anode_area_ratio():
'''The anode area ratio used to calculate the flux.
According to the :ref:`vexelsbible`, it is defined as the proportion
of the anode area that actually measures counts (:math:`A_A`).
The function always returns ``0.87``.
>>> print('%.3f' % anode_area_ratio())
0.870
'''
return 0.87
[docs]def raw2cps(raw):
''' Return "counts-per-sec".
:param raw: Raw count
:type raw: It is normally an ``integer``, but ``float`` is also allowed.
The float is used for different binning parameter and time compression.
>>> cnt = 1
>>> cps = raw2cps(1)
>>> print('%.3f' % cps)
35.556
'''
logger = logging.getLogger(__name__ + '.raw2cps')
acc = times.acctime()
logger.debug('ACC=%f' % acc)
return raw / acc
[docs]def raw2cps_np(raw_np):
''' Return the "count-per-sec" for numpy array.
A numpy array version of :meth:`raw2cps` function.
:param raw_np: ``numpy.array`` of the raw count data.
:type raw_np: ``numpy.array``. Typically 128x16 shaped array, but not restricted.
:returns: Count per sec as numpy array.
:rtype: ``numpy.array`` with the same shape of the ``raw_np``.
>>> cnts = np.array([0,1,2,3,4,5])
>>> cpss = raw2cps_np(cnts)
>>> print(cpss[0])
0.0
>>> print(cpss.shape)
(6,)
>>> print('%.1f' % cpss[5])
177.8
'''
acc = times.acctime()
return raw_np / acc
[docs]def raw2def(raw, energy, anode):
''' Return the differential energy flux (DEF). See equation (7) in section 10 in :ref:`vexelsbible`.
:param raw: Raw count
:param energy: Energy in eV. The value should be taken using :mod:`irfpy.vels.bible.energy` module.
:type energy: ``float``
:param anode: Anode number
:type anode: ``int``
:returns: The differential energy flux in the unit of ``eV/eV cm2 sr s``.
>>> cnt = 1
>>> je_1 = raw2def(1, 30, 5) # one count level for 30 eV, anode 5
>>> print('%.2e' % je_1)
5.56e+06
'''
aa = anode_area_ratio()
gf = gfactor(energy, anode)
acc = times.acctime()
logger = logging.getLogger(__name__ + '.raw2def')
# logger.setLevel(logging.INFO)
logger.debug('ACC=%f s' % acc)
logger.debug('GF=%g cm2 sr eV/eV' % gf)
logger.debug('AA=%f' % aa)
return raw / (gf * acc * aa)
def __vectorize_input(raw_np, energy_np, anode_np):
''' An internal helper function that checks the shape of the arrays for raw2*_np methods.
:param raw_np: Raw counts. (*nE*, *nD*)
:param energy_np: Energy in eV. (*nE*,) or (*nE*, *nD*)
:param anode_np: Anode number. (*nD, )
:returns: (*nE*, *nD*) shaped three matrix. raw, energy and anode matrixes.
Raised ValueError in case of argument inconsistency.
'''
rawshape = raw_np.shape
if len(rawshape) != 2:
emsg = 'Dimension mismatch: len(raw.shape) is (%d) != (2)' % (len(rawshape))
raise ValueError(emsg)
ne = rawshape[0]
nd = rawshape[1]
anodeshape = anode_np.shape
if len(anodeshape) != 1:
emsg = 'Dimension mismatch: len(anode.shape) is (%d) != (1)' % (len(anodeshape))
raise ValueError(emsg)
if nd != anodeshape[0]:
emsg = 'Shape mismatch. args[0].shape[1] (=%d) and args[2].shape[0] (=%d) must be the same' % (
nd, anodeshape[0])
raise ValueError(emsg)
eneshape = energy_np.shape
if len(eneshape) == 1:
# If energy array is 1D, use meshgrid to get (nE, nD) array for energy and anode.
energy2_np, anode2_np = np.meshgrid(energy_np, anode_np)
energy2_np = energy2_np.T
anode2_np = anode2_np.T
elif len(eneshape) == 2:
# If energy array is 2D, use energy as it is, and anode should be repeated by ne times.
energy2_np = energy_np
anode2_np = anode_np.repeat(ne).reshape(nd, ne).T
else:
emsg = 'Dimension mismatch: len(ene.shape) is (%d) != (1 or 2)' % (len(rawshape))
raise ValueError(emsg)
eneshape = energy2_np.shape
if rawshape != eneshape:
emsg = 'Shape mismatch. args[0].shape (=%s) and args[1].shape (=%s) must be the same' % (
rawshape, eneshape)
raise ValueError(emsg)
return (raw_np, energy2_np, anode2_np)
[docs]def raw2def_np(raw_np, energy_np, anode_np):
''' Return the DEF. A vectorized version of :meth:`raw2def`.
:param raw_np: Raw count matrix. Shape should be (*nE*, *nD*) where nE is the
number of energy step, and nD is the number of the anode (direction).
:type raw_np: ``numpy.array`` of ``float`` with shape of (*nE*, *nD*)
:param energy_np: Energy step in eV. Shape should be (nE,) or (nE, nD)
:type energy_np: ``numpy.array`` of ``float`` with shape of (*nE*,) or (nE, nD)
:param anode_np: Anode number. Shape should be (*nD*,).
:type anode_np: ``numpy.array`` of ``int`` with shape of (*nD*,)
:returns: The differential energy flux in the unit of ``eV/eV cm2 sr s``.
Raises ``ValueError`` if the shapes are inconsistent.
>>> etbl = np.array([1., 10, 100, 1000, 10000]) # As a sample, 5 energy steps assumed.
>>> anod = np.array([0, 1, 2, 3]) # As a sample, 4 anode steps assumed.
>>> cnts = np.ones([5, 4]) # Count data should have shape of (*nE*, *nD*).
>>> deflux = raw2def_np(cnts, etbl, anod)
>>> print(deflux.shape)
(5, 4)
>>> print('%.3e' % deflux[3, 2]) # One count level of 1000 eV Anode=2
4.102e+06
>>> print('%.3e' % raw2def(1, 1000., 2)) # Scalar version for confirmation
4.102e+06
'''
rr, ee, dd = __vectorize_input(raw_np, energy_np, anode_np)
ne = rr.shape[0]
nd = rr.shape[1]
aa = anode_area_ratio()
acc = times.acctime()
gfmatrix = np.zeros([ne, nd])
for _d in range(nd):
gfmatrix[:, _d] = gfactor1d(ee[:, _d], anode_np[_d])
return rr / (gfmatrix * acc * aa)
[docs]def raw2dnf(raw, energy, anode):
''' Return the differential number flux (DNF).
Apart from the equation (8) in section 10 in :ref:`vexelsbible`,
I used here the unit in ``/eV cm2 sr s``.
:param raw: Raw count
:param energy: Energy in eV
:type energy: ``float``
:param anode: Anode number
:type anode: ``int``
:returns: The differential number flux in the unit of ``/eV cm2 sr s``.
>>> cnt = 1
>>> jn_1 = raw2dnf(1, 30, 5) # one count level for 30 eV, anode 5
>>> print('%.2e' % jn_1)
1.85e+05
'''
defflx = raw2def(raw, energy, anode)
logger = logging.getLogger(__name__ + '.raw2dnf')
logger.debug('DEF=%e' % defflx)
return defflx / energy
[docs]def raw2dnf_np(raw_np, energy_np, anode_np):
''' Return the DNF. A vectorized version of :meth:`raw2dnf`.
:param raw_np: Raw count matrix. Shape should be (*nE*, *nD*) where nE is the
number of energy step, and nD is the number of the anode (direction).
:type raw_np: ``numpy.array`` of ``float`` with shape of (*nE*, *nD*)
:param energy_np: Energy step in eV. Shape should be (nE,) or (nE, nD)
:type energy_np: ``numpy.array`` of ``float`` with shape of (*nE*,) or (nE, nD)
:param anode_np: Anode number. Shape should be (*nD*,).
:type anode_np: ``numpy.array`` of ``int`` with shape of (*nD*,)
:returns: The differential number flux in the unit of ``/eV cm2 sr s``.
Raises ``ValueError`` if the shapes are inconsistent.
>>> etbl = np.array([1., 10, 100, 1000, 10000]) # As a sample, 5 energy steps assumed.
>>> anod = np.arange(16) # As a sample, 16 anode steps assumed.
>>> cnts = np.ones([5, 16]) # Count data should have shape of (*nE*, *nD*).
>>> dnflux = raw2dnf_np(cnts, etbl, anod)
>>> print(dnflux.shape)
(5, 16)
>>> print('%.3e' % dnflux[2, 8]) # One count level of 100 eV Anode=8
5.086e+04
>>> print('%.3e' % raw2dnf(1, 100., 8)) # Scalar version for confirmation
5.086e+04
'''
rr, ee, dd = __vectorize_input(raw_np, energy_np, anode_np)
dnf_np = raw2def_np(raw_np, energy_np, anode_np) / ee
return dnf_np
[docs]def raw2psd(raw, energy, anode):
r''' Return the phase space density (PSD).
.. note::
I think the :ref:`vexelsbible` contains error in the equation 10.
According to the definition, earray is in ``eV`` and E is ``J/eV``.
Thus, :math:`earray \times E` should be in ``J`` = ``kg m2 / s2``.
To use the unit of ``cm``, the term should be changed to be ((1e2)^2 x earray x E).
So the right expressions is :math:`PSD [cm^{-6} s^3] = \frac{raw \times m_e^2}
{GF \times accutime \times 2 \times (1e4 \times earray \times E)^2 \times A_A}`
:param raw: Raw count
:param energy: Energy in eV
:type energy: ``float``
:param anode: Anode number
:type anode: ``int``
:returns: PSD in the unit of ``cm-6 s3``.
>>> cnt = 1
>>> f_1 = raw2psd(1, 30, 5) # one count level for 30 eV, anode 5
>>> print('%.2e' % f_1)
1.00e-27
'''
qe = 1.602e-19
me = 9.11e-31
me_qe = me / qe # Order of -12
me_qe2 = me_qe * me_qe # Order of -24
gf = gfactor(energy, anode) # Order of -6
acc = times.acctime()
aa = anode_area_ratio()
psd = (me_qe2 * raw) / (2 * 1e8 * gf * acc * energy * energy * aa)
return psd
[docs]def raw2psd_np(raw_np, energy_np, anode_np):
''' Return the PSD. A vectorized version of :meth:`raw2psd`.
:param raw_np: Raw count matrix. Shape should be (*nE*, *nD*) where nE is the
number of energy step, and nD is the number of the anode (direction).
:type raw_np: ``numpy.array`` of ``float`` with shape of (*nE*, *nD*)
:param energy_np: Energy step in eV. Shape should be (nE,)
:type energy_np: ``numpy.array`` of ``float`` with shape of (*nE*,)
:param anode_np: Anode number. Shape should be (*nD*,).
:type anode_np: ``numpy.array`` of ``int`` with shape of (*nD*,)
:returns: The differential number flux in the unit of ``s3/cm6``.
Raises ``ValueError`` if the shapes are inconsistent.
>>> etbl = np.array([30., 300, 3000]) # As a sample, 3 energy steps assumed.
>>> anod = np.arange(8) # As a sample, 8 anode steps assumed.
>>> cnts = np.ones([3, 8]) # Count data should have shape of (*nE*, *nD*).
>>> psd = raw2psd_np(cnts, etbl, anod)
>>> print(psd.shape)
(3, 8)
>>> print('%.3e' % psd[2, 5]) # One count level of 3000 eV Anode=5
9.917e-32
>>> print('%.3e' % raw2psd(1, 3000., 5)) # Scalar version for confirmation
9.917e-32
.. todo::
Too slow due to loops. Need implementation of vectorized version like raw2dnf.
'''
rr, ee, dd = __vectorize_input(raw_np, energy_np, anode_np)
ne = rr.shape[0]
nd = rr.shape[1]
psd_np = np.zeros_like(raw_np)
for id in range(nd):
for ie in range(ne):
psd_np[ie, id] = raw2psd(rr[ie, id], ee[ie, id], dd[ie, id])
return psd_np