Source code for irfpy.util.fluxtools

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