''' Empirical energy and angular spectra
The empirical energy and angular spectra module.
- Energy spectrum: :func:`j` (with unit support) and :func:`j2` (no unit support)
- Angular response: :class:`AngularAV13`
To avoid the uncertainty of the unit conversion,
the :mod:`irfpy.util.unitq` module is used.
(https://irfpy.irf.se/projects/util/api/api_irfpy.util.unitq.html)
'''
import numpy as np
import irfpy.util.unitq as u
import irfpy.util.constant as k
from irfpy.cena import empirical_lib
[docs]def temperature_ena(velocity):
''' Calculate the ENA temperature.
:param velocity: Velocity of the solar wind (unitful by irfpy.util.unitq)
:returns: Temperature of the reflected ENAs.
:rtype: ``float`` with unit of the order of ``u.K``
The ENA temperature model is based on
Futaana et al., 2012; 2013.
>>> kT_ena = temperature_ena(300 * u.km / u.s)
>>> print(kT_ena.nr(u.K)) # in Kelvin
912900.0
The result is converted to eV (temperature)
>>> print('{:.2f}'.format(kT_ena.nr(u.eV_asT))) # in eV (temperature)
78.67
'''
vel = velocity.nr(u.m / u.s)
T = (vel * 3.12 - 2.31e4) * u.K
return T
[docs]def guess_vsw(t_ena):
''' Calculate the solar wind velocity from the obtained temperatuer.
Inverse function of :func:`temperature_ena`.
See Futaana et al., 2013 for application.
:param t_ena: ENA temperature. Should have the unit of temperature.
:returns: The guessed velocity of the solar wind.
:rtype: ``float`` with unit equivalent to speed ``u.m / u.s``.
>>> print(guess_vsw(912900.0 * u.K).nr(u.km / u.s))
300.0
'''
T = t_ena.nr(u.K)
vel = (T + 2.31e4) / 3.12 * u.m / u.s
return vel
[docs]def flux_half_maxwell(nena, Tena, E0=0. * u.J, angular_coefficient=2 * np.pi):
r''' Calculate the flux of half-maxwell distribution function
:param nena: ENA full density with unit of the order of /u.cm ** 3
:param Tena: ENA temperature with unit of the order of u.K
:keyword E0: Limit in the lower bound of the energy with unit.
:keyword angular_coefficient: Coefficient of the part of :math:`\int\int f_ad^2\Omega`
In case of the isotropic distribution, it will be :math:`2\pi` which is a default.
>>> n = 1 / u.cm**3
>>> T = 100 * u.eVT
>>> print('%.2e' % flux_half_maxwell(n, T).nr(1/u.cm**2/u.s))
7.81e+06
>>> n = 4.28 / u.cm**3
>>> T = 135.7 * u.eVT
>>> print('%.2e' % flux_half_maxwell(n, T).nr(1/u.cm**2/u.s))
3.89e+07
'''
e_lowlimit = E0
nena_mksa = nena.nr(u.m**-3)
Tena_mksa = Tena.nr(u.K)
e_lowlimit_mksa = e_lowlimit.nr(u.J)
# Alias
n = nena_mksa
m = k.m_hydrogen
pi = np.pi
kT = k.kB * Tena_mksa
E0 = e_lowlimit_mksa
# The partial flux should be in MKSA
partial_flux = n * np.sqrt(2.0 / (m * pi * kT)) * (kT + E0) * np.exp(-E0 / kT)
partial_flux = partial_flux / u.m ** 2 / u.s
partial_flux *= (angular_coefficient / (2 * pi))
return partial_flux
[docs]def refrate():
''' Return the average reflection rate (Futaana et al., 2012).
'''
return 0.19
[docs]def j(Fsw, vel, ref=None):
r''' Get the empirical function of the differential flux.
See Futaana et al., 2012.
:param Fsw: Solar wind flux with unit equivalent to ``u.cm ** -2 * u.s ** -1``.
If :math:`SZA` (solar zenith angle) plays into the game, the geometric reduction should be compensated prior by multiplying the :math:`\cos SZA`.
:param vel: Solar wind velocity with unit of speed (e.g, ``u.m / u.s``).
:keyword ref: Reflection rate. Default is taken from :func:`refrate`.
:returns: An empirical function is returned.
The returned function has the following signature:
- Input: Energy with unit of energy (e.g., ``eV``)
- Ouptput: Full-differential flux (i.e., the unit of ``1 / u.cm ** 2 / u.s / u.ster / u.eV``).
The isotropic distrubition is assumed.
>>> Fsw = 3e8 / u.cm ** 2 / u.s # The solar wind flux
>>> Vsw = 400 * u.km / u.s # The solar wind speed
>>> j_empir = j(Fsw, Vsw) # The empirical function
>>> E = 30 * u.eV_asE # 30 eV of energy as input
>>> j30 = j_empir(E) # Differential flux at 30 eV
>>> print('%.1f' % j30.nr(1 / u.eV_asE / u.cm ** 2 / u.s / u.ster))
18383.8
'''
if ref is None:
ref = refrate()
def je_func(energy):
''' A function taking unitful energy to return the differential flux.
'''
kb = k.kB * u.J / u.K
m = k.m_hydrogen * u.kg
je = 0.5 * energy / (kb * temperature_ena(vel)) ** 2
je = je * ref * Fsw
je = je / np.pi
inexp = -energy / (kb * temperature_ena(vel))
je = je * np.exp(inexp.ns) # ns gives simplified number.
return je / u.ster
return je_func
[docs]def j2(Fsw__cm2s, Vsw_km_s, ref=None):
''' Equivalent to :func:`j` but with ``float`` version.
:param Fsw__cm2s: The flux of the solar wind in the unit of ``/cm2 /s``
:param Vsw_km_s: The speed of the solar wind in the unit of ``km/s``
>>> Fsw = 3e8
>>> Vsw = 400.
>>> espec = j2(Fsw, Vsw)
>>> j30 = espec(30) # 30 eV
>>> print('%.1f' % j30)
18383.8
'''
if ref is None:
ref = refrate()
def je_func(energy_eV):
kb = k.kB * u.J / u.K
m = k.m_hydrogen * u.kg
energy = energy_eV * u.eV_asE
Fsw = Fsw__cm2s / u.cm**2 / u.s
vel = Vsw_km_s * u.km / u.s
je = 0.5 * energy / (kb * temperature_ena(vel)) ** 2
je = je * ref * Fsw
je = je / np.pi
inexp = -energy / (kb * temperature_ena(vel))
je = je * np.exp(inexp.ns)
return je.nr(1 / u.cm**2 / u.s / u.eV_asE)
return je_func
[docs]def data_maxwellfit(counts, ch, initprm=None):
''' This is to calculate the best parameters for the obtained data.
:param counts: Count array with shape of (16,)
:type counts: ``np.ma.masked_array``
:param ch: Channel of CENA
:type ch: ``int``
:param initprm: Initial parameter. Tuple of three float.
:returns: Tuple of the fitted results. First element is a tuple of the density and temperature.
They are unitfull with :mod:`irfpy.util.units`.
The second element is a background.
The third element is the log-likeliness calculated during fitting.
This can be used for calculating AIC.
The foruth element is the fitted function. The function will take float energies (np.array)
in eV to return the differential flux in /cm2 sr s eV.
The Maxwell + BG fitting from the given count rate data.
>>> cnts = [float(v) for v in 'nan nan nan nan 0.041152 0.078189 0.164609 0.378601 0.460905 0.547325 0.242798 0.057613 nan nan nan nan'.split()]
>>> fit = data_maxwellfit(cnts, 3)
>>> n = fit[0][0]
>>> print('%.2f' % n.nr(u.cm**-3))
2.44
>>> T = fit[0][1]
>>> print('%.2f' % T.nr(u.eVT))
87.87
'''
if initprm is None:
initprm = (20, 10, 0.1)
maxbgfunc = empirical_lib.flux_function_maxwell_bg(ch)
retval = empirical_lib.run_fitting(counts, ch, maxbgfunc, initprm)
prm, s, f, likely = retval
bg = prm[2] ** 2
n, T = empirical_lib.maxwell_parameter_to_physical_values(prm[0], prm[1])
return (n * (u.cm ** -3), T * (u.eVT)), bg, likely, f
[docs]def maxwell_func(n, T):
''' Return unitful Maxwell function of differential flux from the given n and T.
>>> n = 2.2 * (u.cm ** -3)
>>> T = 85.3 * (u.eVT)
>>> j_mx = maxwell_func(n, T) # Return a function.
>>> j100 = j_mx(100 * u.eV_asE) # Differential flux of 100 eV
>>> print('%.1f' % j100.nr(1 / u.cm ** 2 / u.ster / u.eV_asE / u.s))
10743.2
The :func:`maxwell_func` can be used, instead of the function returned by
:func:`data_maxwellfit`, to plot the fitted differential flux.
'''
n_ul = n.nr(u.m ** -3)
T_ul = T.nr(u.K)
kT = k.kB * T_ul
m = k.m_hydrogen
pi = np.pi
def je(e):
''' Energy spectra for density of %f /cm3 and temperature of %f eV.
''' % (n.nr(u.cm ** -3), T.nr(u.eVT))
e_ul = e.nr(u.J)
inexp = -e_ul / kT
coeff1 = (m / 2. / pi / kT) ** (1.5)
j = 2. * n_ul / (m ** 2) * coeff1 * e_ul * np.exp(inexp)
return j * (1 / u.m ** 2 / u.ster / u.J / u.s)
return je
[docs]class AngularAS11:
r""" Implementation of angular response, *f_s*, in Schaufelberger 2011.
.. note::
It is better using :func:`AngularAV13` as it is more understandable definition
.. math::
J_{ENA}(SZA, \phi, \theta) = J_{SW}\cdot R_{SS} \cdot f_S(SZA, \phi, \theta)
This class is for *f_S*.
Note that J_ENA is in the unit of [#/cm2 sr s] and J_SW is [#/cm2 s]
"""
def __init__(self):
pass
[docs] def fs(self, sza_deg, phi_deg, theta_deg):
r""" Return the angular response, fs. The unit is [1/sr].
:param sza_deg: Solar zenith angle, in degrees
:param phi_deg: Scattering azimuth angle, in degrees. 0 is for the sun direction
:param theta_deg: Scattering polar angle, in degrees. 0 is for the zenith dirction
:return: Sattering function fs in 1/steradian
.. math::
f_S(SZA, \phi, \theta) = f_0(SZA) \cdot
f_1(SZA, \phi) \cdot
f_2(SZA, \phi) \cdot
f_3(SZA, \theta)
The form is in Schaufelberger et al (2011).
"""
sza0 = np.deg2rad(sza_deg)
phi0 = np.deg2rad(phi_deg)
theta0 = np.deg2rad(theta_deg)
f0 = self._f0(sza0)
f1 = self._f1(sza0, phi0)
f2 = self._f2(sza0, phi0)
f3 = self._f3(sza0, theta0)
return f0 * f1 * f2 * f3
def _f0(self, sza):
z0 = self._z0(sza)
return 0.009 * z0
def _f1(self, sza, phi):
z1 = self._z1(sza)
return z1 * np.cos(2 * phi) + (1 - z1)
def _f2(self, sza, phi):
z2 = self._z2(sza)
return z2 * np.cos(phi) + (1 - z2)
def _f3(self, sza, theta):
z3 = self._z3(sza) # sza here in radian, z3 in degrees
return (1 - z3 / 90) * np.sin(theta + np.deg2rad(z3)) + z3 / 90
def _z0(self, sza):
return 1.24 * np.cos(1.44 * sza - np.deg2rad(48.7))
def _z1(self, sza):
return 0.30 * sza + np.deg2rad(1.72)
def _z2(self, sza):
return 0.24 * np.cos(np.deg2rad(74.48) - 1.52 * sza)
def _z3(self, sza):
return 90 - 1.03 * np.rad2deg(sza)
[docs]class AngularAV13(AngularAS11):
""" Implementation of angular response, *f_s*, in Vorburger 2013.
.. note::
As of my understanding, the AS11 model and AV13 model
definitions are different.
One recommended way to use is to normalize the
angular response to re-distribute to the integrated
(energy-)differential flux.
"""
def __init__(self):
AngularAS11.__init__(self)
def _f0(self, sza):
# sza in radian
t1 = np.cos(sza)
t2 = (0.074 * (sza - 3.23)
* (2 * sza * np.cos(1.03 * sza) + sza * np.sin(1.03 * sza)
* np.pi - 4 * (sza - 1.53)) * (np.cos(1.52 * sza - 1.30) - 4.17))
return t1 / t2