""" Tools for flux (differential, velocity distribution)
.. codeauthor:: Yoshifumi Futaana
This module provide tools for flux manipulations
(differential flux, count rate, velocity distribution function)
Containers of differential flux and count rate in a spherical coordinates are
implemented in :mod:`irfpy.util.energyspectrum`.
The velocity distribution functions are implemented
in :mod:`irfpy.util.vdf` module.
.. autosummary::
irfpy.util.energyspectrum.DifferentialFluxGrid
irfpy.util.energyspectrum.CountRateGrid
irfpy.util.vdf.VelocityDistributionGrid
vdf2dflux
dflux2vdf
cnt2dflx_simple_convfact
.. codeauthor:: Yoshifumi Futaana
"""
import logging
_logger = logging.getLogger(__name__)
import numpy as np
from irfpy.util import physics as ph
from irfpy.util import vdf
from irfpy.util import energyspectrum
from irfpy.util import nptools
[docs]def vdf2dflux(vdfun):
r"""Conversion from velocity distribution function to differential flux.
:param vdfun: Velocity distribution function on polar grid.
:type vdfun: :class:`irfpy.util.vdf.VelocityDistributionGrid`
:return: Differential flux object.
:rtype: :class:`irfpy.util.energyspectrum.DifferentialFluxGrid`
In MKSA, the conversion from VDF to Dflx is
.. math::
J = \frac{2 f E}{m^2}
Limitation: So far, the energy/theta/phi_grid_size is limited to default
(see :class:`irfpy.util.energyspectrum.DifferentialFluxGrid`.
Only the central value is used for conversion.
>>> vgrid = np.linspace(0.01, 1000e3, 96) # 96 velocity space
>>> tgrid = np.linspace(30, 150, 32) # 32 polar steps
>>> pgrid = np.linspace(-180, 180, 65)[:-1] # 64 azimuth steps
>>> maxwell = vdf.maxwellVelocityDistributionGrid(1e6, [-400e3, 0, 0], 80e3, vgrid, np.deg2rad(tgrid), np.deg2rad(pgrid))
>>> maxwell_dflx = vdf2dflux(maxwell)
>>> print('%.4f' % maxwell_dflx.density())
0.9999
>>> v = maxwell_dflx.velocity_vector()
>>> print('%.2f %.2f %.2f' % (v[0], v[1], v[2]))
-400.02 -0.00 0.00
"""
v = vdfun.rgrid
t = vdfun.tgrid
p = vdfun.pgrid
m = vdfun.mass
fval = vdfun.data
qe = 1.602e-19
ene = (v ** 2) * m / 2 # in MKSA, Joule
ene3, t3, p3 = np.meshgrid(ene, t, p, indexing='ij')
jval = fval * 2 * ene3 / (m ** 2) # in MKSA. /m^2 sr J s
jval = jval * qe * 1e-4 # in /cm2 sr eV s
ene_eV = ene / qe
return energyspectrum.DifferentialFluxGrid(jval, ene_eV, t, p, mass=m)
[docs]def dflux2vdf(dflx):
r""" Conversion from differential flux to velocity distribution function.
:param dflx: Differential flux object in :class:`irfpy.util.energyspectrum.DifferentialFluxGrid` object.
:return: Velocity distribution function in :class:`irfpy.util.vdf.VelocityDistributionGrid` class.
>>> vgrid = np.linspace(0.01, 1000e3, 96) # 96 velocity space
>>> tgrid = np.linspace(30, 150, 32) # 32 polar steps
>>> pgrid = np.linspace(-180, 180, 65)[:-1] # 64 azimuth steps
>>> maxwell = vdf.maxwellVelocityDistributionGrid(1e6, [-400e3, 0, 0], 80e3, vgrid, np.deg2rad(tgrid), np.deg2rad(pgrid))
>>> maxwell_dflx = vdf2dflux(maxwell)
>>> maxwell_vdf = dflux2vdf(maxwell_dflx) # Back again.
>>> print('%.5f' % ((maxwell.data - maxwell_vdf.data) ** 2).max())
0.00000
"""
ene = dflx.rgrid # in eV
t = dflx.tgrid
p = dflx.pgrid
m = dflx.mass
qe = 1.602e-19
ene_Joule = ene * qe
j = dflx.data # in cm-2 sr-1 s-1 eV-1
j_mksa = j * 10000 / qe # in m-2 sr-1 s-1 J-1
ene_Joule, t3, p3 = np.meshgrid(ene_Joule, t, p, indexing='ij')
fval = j_mksa / 2 / ene_Joule * (m ** 2) # in MKSA
vgrid = np.sqrt(ene_Joule[:, 0, 0] * 2 / m)
return vdf.VelocityDistributionGrid(fval, vgrid, t, p, mass=m)
[docs]def cnt2dflx_simple_convfact(count_rate, gfactor):
r""" Conversion from count rate to differential flux using g-factor.
Conversion from count rate to differential flux using g-factor (table)
is implemented.
.. math::
J = \frac{C}{G \times E}
where *C* is the count rate [/s], G is the geometric factor
including efficiency [cm2 sr eV/eV], and *E* is the
energy [eV]. The obtained differential flux is
J, with the unit of [/cm2 sr eV s].
The usage is identical to :func:`dflx2cnt_simple_convfact`,
which provides inverse function.
:param count_rate: Calculated count rate.
:type count_rate: :class:`irfpy.util.energyspectrum.CountRateGrid`
:param gfactor: Geometric factor, in the same shape as ``count_rate.data`` *(iene, ipol, iaz)*, or broadcast-able.
If gfactor depends on the energy, and gfactor has the shape of (iene), which is not broadcastable.
In such case, you may need to make the gfactor array by using :func:`irfpy.util.nptools.broadcast_r`.
See the sample below.
:returns: Differential flux object.
:rtype: :class:`irfpy.util.energyspectrum.DifferentialFluxGrid`
"""
estep = count_rate.energy_grid # Energy grid, (nE) shape
cps = count_rate.cps # Differential flux (nE, nPol, nAz) shape
estep3d = nptools.broadcast_r(estep, cps.shape) # Energy grid, (nE, nPol, nAz) shape
dflx = cps / estep3d / gfactor
dflx = energyspectrum.DifferentialFluxGrid(dflx,
count_rate.energy_grid, count_rate.theta_grid, count_rate.phi_grid,
energy_grid_size=count_rate.energy_grid_size,
theta_grid_size=count_rate.theta_grid_size,
phi_grid_size=count_rate.phi_grid_size)
return dflx
[docs]def dflx2cnt_simple_convfact(diff_flux, gfactor):
r""" Differential flux is converted to count rate.
With a simple conversion of "by definition".
.. math::
C = J \times G \times E
where *J* is the differential flux [/cm2 s sr eV], *G* is the geometric factor including efficiency [cm2 sr eV/eV],
*E* is the energy [eV] (not energy range).
The returned count rate will be [/s].
:param diff_flux: Differential flux object.
:type diff_flux: :class:`irfpy.util.energyspectrum.DifferentialFluxGrid`
:param gfactor: Geometric factor, in the same shape as ``dflx.data`` *(iene, ipol, iaz)*, or broadcast-able.
If gfactor depends on the energy, and gfactor has the shape of (iene), which is not broadcastable.
In such case, you may need to make the gfactor array by using :func:`irfpy.util.nptools.broadcast_r`.
See the sample below.
:return: Calculated count rate.
:rtype: :class:`irfpy.util.energyspectrum.CountRateGrid`
For comprehensive example, I first prepare the energy array for a "virtual" sensor:
>>> elist = [10, 20, 40, 80, 160, 320, 640, 1280, 2560, 5120] # 10 energy steps
Then, angle information of the sensor are defined.
>>> tlist = np.radians([15, 45, 75, 105, 135, 165]) # Polar angle, 6 channel
>>> plist = np.radians([0, 45, 90, 135, 180, 225, 270, 315]) # Azimuth, 8 channels
As an exercise, I prepare the Maxwell distribution
>>> n = 1e6 # 1 /cc
>>> vb = [200e3, 200e3, 0] # 200 km/s vx and vy
>>> vth = 70e3 # 70 km/s thermal velocity
>>> diff_flux = energyspectrum.maxwellDifferentialFluxGrid(n, vb, vth, elist, tlist, plist)
>>> print('{:.2f}'.format( diff_flux.data.max()))
68765.73
Now the ``diff_flux`` is an object of :class:`irfpy.util.energyspectrum.DifferentialFluxGrid`
Giving the g-factor of 1e-5 cm2 sr eV/eV, you may get the count rate object
(:class:`irfpy.util.energyspectrum.CountRateGrid` object).
>>> cps = dflx2cnt_simple_convfact(diff_flux, 1e-5)
>>> print('{:.1f}'.format(cps.cps.max()))
440.1
If the g-factor has energy dependence, and you have g-factor array, you may need to make a
g-factor array to be *(nEne, nPol, nAz)* shape.
>>> gfactor = np.array([1, 2, 3, 4, 5, 5, 4, 3, 2, 1]) * 1e-5
>>> gfactor_3d = nptools.broadcast_r(gfactor, [10, 6, 8])
>>> print(gfactor_3d.shape)
(10, 6, 8)
>>> cps = dflx2cnt_simple_convfact(diff_flux, gfactor_3d)
>>> print('{:.1f}'.format(cps.cps.max()))
1760.4
The inversion function is also prepared as :func:`cnt2dflx_simple_convfact`.
>>> dflx2 = cnt2dflx_simple_convfact(cps, gfactor_3d)
The results should be (almost) identical.
>>> diff = dflx2.dflux - diff_flux.dflux
>>> print('{:.10f}'.format((diff ** 2).max()))
0.0000000000
"""
estep = diff_flux.energy_grid # Energy grid, (nE) shape
dflx = diff_flux.dflux # Differential flux (nE, nPol, nAz) shape
estep3d = nptools.broadcast_r(estep, dflx.shape) # Energy grid, (nE, nPol, nAz) shape
nene, npol, naz = dflx.shape
cnts = dflx * estep3d * gfactor
cps = energyspectrum.CountRateGrid(cnts, diff_flux.energy_grid, diff_flux.theta_grid, diff_flux.phi_grid,
energy_grid_size=diff_flux.energy_grid_size,
theta_grid_size=diff_flux.theta_grid_size,
phi_grid_size=diff_flux.phi_grid_size)
return cps
[docs]def fit_differential_flux_by_maxwell(diff_flux, log=False, initial_params=None, method='L-BFGS-B'):
""" Fit the energy spectrum by Maxwellian using least-square method.
:param diff_flux: Differential flux.
:type diff_flux: :class:`irfpy.util.energyspectrum.DifferentialFluxGrid`
:keyword log: If specified, the fitting is done in the log-space (not yet implemented)
:keyword initial_params: You may specify this as a tuple.
(density_cc, vx_km/s, vy_km/s, vz_km/s, vth_km/s)
If ``None`` is given, the initial parameters are calculated from
moment calculations with proton assumed.
(:meth:`irfpy.util.energyspectrum.DifferentialFluxGrid.density()`,
:meth:`irfpy.util.energyspectrum.DifferentialFluxGrid.velocity_vector()`,
:meth:`irfpy.util.energyspectrum.DifferentialFluxGrid.pressure_tensor()`)
:keyword method: Method for ``scipy.optimize.minimize``. Default is 'L-BFGS-B'.
:author: Yoshifumi Futaana
>>> estep = [10, 20, 40, 80, 160, 320, 640, 1280, 2560, 5120, 10240, 20480]
>>> tstep = np.deg2rad([22.5, 45, 67.5, 90, 112.5, 135, 157.5])
>>> pstep = np.deg2rad([0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330])
>>> maxwell_flux = energyspectrum.maxwellDifferentialFluxGrid(
... 10e6, [400e3, 100e3, -50e3], 200e3,
... estep, tstep, pstep)
>>> fit = fit_differential_flux_by_maxwell(maxwell_flux)
>>> print('{:.1f} {:.1f} {:.1f} {:.1f} {:.1f}'.format(fit[0], fit[1], fit[2], fit[3], fit[4]))
10.0 400.0 100.0 -50.0 200.0
"""
if log:
raise NotImplementedError('log is not yet supported. Sorry.')
### First, initial_params are obtained from "moment calculation".
if initial_params is None:
n = diff_flux.density() # Density in cm^-3
v = diff_flux.velocity_vector() # Velocity in km/s
p = diff_flux.pressure_tensor()
p = (p[0, 0] + p[1, 1] + p[2, 2]) / 3. # Average of the disgnostic component
temperature = ph.pressure_2_temperature(p * 1e9, n) # pressure in nPa, density in cm^-3, returned in K.
vth = ph.temperature_2_thermal_velocity(1.67e-27, temperature) # Mass in kg, tempature in K, vth in m/s
vth = vth / 1000 # Now in km/s
initial_params = (n, v[0], v[1], v[2], vth)
### Observed flux is in diff_flux
energy_grid = diff_flux.energy_grid
theta_grid = diff_flux.theta_grid
phi_grid = diff_flux.phi_grid
energy_grid_size = diff_flux.energy_grid_size
theta_grid_size = diff_flux.theta_grid_size
phi_grid_size = diff_flux.phi_grid_size
dfobs = diff_flux.dflux
def log_likelihood(prm):
""" Return the negative of log_likelihood.
"""
### Predicuted flux is produced by energyspectrum.DifferentialFluxGrid
diff_flux_p = energyspectrum.maxwellDifferentialFluxGrid(
abs(prm[0]) * 1e6, np.array([prm[1], prm[2], prm[3]]) * 1e3 , abs(prm[4]) * 1e3,
energy_grid, theta_grid, phi_grid,)
dfpred = np.clip(diff_flux_p.dflux, 0, diff_flux_p.dflux.max())
logL = ((dfpred - dfobs) ** 2).sum()
# _logger.debug('Fitting ongoing: {} {}'.format(str(prm), str(logL)))
return logL
from scipy.optimize import minimize
res = minimize(log_likelihood, np.array(initial_params) * 1.5,
bounds=[(0, 200), (-5000, 5000), (-5000, 5000), (-5000, 5000), (-5000, 5000)], method=method)
return res.x