irfpy.util.energyspectrum
¶
Generic purpose energy spectrum implementation
Code author: Yoshifumi Futaana
DifferentialFluxGrid
represents the differential flux on spherical grids.
CountRateGrid
represents the count rate on spherical grids.
The conversion in between, gfactor information should be available.
For more serious data analysis, the response of the sensor is also needed.
These “instrument model” is not straightforwardly obtained.
However, the simplest implementation is anyway available in irfpy.util.fluxtools
.
A counter part (velocity distribution function) is implemented
in irfpy.util.vdf
module. irfpy.util.fluxtools
will
be used for conversion.
Differential flux is represented as a unit of
/cm2 s eV sr
in plasma physics.
The differential flux, j is a function of
the energy and directions.
Directions are defined by the polar coordinate system, namely, the angle \(\theta\) as the angle from zaxis and the angle \(\phi\) as the angle from x axis as in the projection into xy plane.

Representation of the differential flux on grids in polar coordinates. 

Represent the count rate (counts per sec) on grids in polar frame. 

Create 
 class irfpy.util.energyspectrum.DifferentialFluxGrid(dflux, energy_grid, theta_grid, phi_grid, mass=1.67262192369e27, energy_grid_size='log', theta_grid_size='linear', phi_grid_size='linear')[source]¶
Bases:
RegularGridVolumeContainer
Representation of the differential flux on grids in polar coordinates.
The natural log of differential flux can be assigned at a grid of (E, \(\theta\), \(\phi\)). Also, moment calculation can be possible.
 Parameters:
dflux – Differential flux. Should be a shape of (ne, nt, np). Unit of /cm2 s eV sr.
energy_grid – Energy grid in eV. Ascending order.
theta_grid – Theta (polar angle) angle array in radians. Ascending order.
phi_grid – Phi (azimuth) angle array in radians. Ascending order.
mass – Mass in kg. It is used for calculating the moment values.
energy_grid_size – Specify the energy grid size. It is used to calculate the moment values. Use ‘linear’, ‘log’, or np.array with the same shape as energy_grid.
theta_grid_size – Same as energy_grid_size but for theta direction.
phi_grid_size – Same as energy_grid_size, but for phi direction.
Note
(Note on 190118): For integrators, we use simple approach (see
integral_flux()
). One can have better answer by integrating \(\sin\theta d\theta\) properly. Will see how we should develop. UNIT = 'cm^2 s^1 eV^1 sr^1'¶
Unit of the differential flux
 dflux¶
Differential flux
 energy_grid¶
Energy grid
 ne¶
Number of energy grid
 theta_grid¶
Theta grid (radians)
 nt¶
Number of theta grid
 phi_grid¶
Phi grid (radians)
 np¶
Number of phi grid
 energy_grid_size¶
Width of energy grid (i.e. dE)
 theta_grid_size¶
Width of theta grid (i.e. dTheta)
 phi_grid_size¶
Width of phi grid (i.e. dPhi)
 mass¶
Mass of the particles
 density()[source]¶
Calculate the density by integrating over the valid counts.
 Returns:
Density in cm^3
Calculating density from differential flux is formulated as
\[n = \int\int\int \sqrt{\frac{m}{2E}} J(E, \theta, \phi) \sin \theta dE d\theta d\phi\]This formulation can be obtained simply by definition
\[n = \int\int\int f(v_x, v_y, v_z) dv_x dv_y dv_z\]with the variables expressed in (E, theta, phi):
\[dv_x dv_y dv_z = \sqrt{\frac{2E}{m^3}}\sin\theta\]and the relation
\[j = \frac{2E}{m^2} f\]
 integral_flux()[source]¶
 Returns:
Integral flux. Unit is /cm2 s.
The integral is, in MKSA unit, expressed by
\[\begin{split}\vec{F} &= \int\int\int \vec{j}dE \sin\theta d\theta d\phi \\ &= \vec{j_0} \sin\theta_0 \Delta E\Delta\theta \Delta\phi\end{split}\]
 velocity_vector()[source]¶
Returns the velocity vector in km/s
 Returns:
(vx, vy, vz) in km/s.
The velocity vector is
\[F / n\]where F is the
integral_flux()
and n is the denisytdensity()
.
 pressure_tensor()[source]¶
Return the pressure tensor in Pa.
The pressure tensor is defined by
\[P_{i,j} = m\int\int\int v_i v_j j(E, \theta, \phi) \sin\theta \sqrt{m/2E} dE d\theta d\phi  nmV_i V_j\]where V_i is the ith component bulk velocity and n for density. If we write \(v_i = v a_i\), where a_i is the unit vector of the velocity vector, we can rewrite the pressure as
\[P_{i, j} = m\int\int\int a_i a_j \sqrt{2E/m} \sin\theta j(E, \theta, \phi) dE d\theta d\phi  nmV_iV_j\]Example follows.
>>> n = 1e6 # The density in m^3 >>> vb = [100e3, 300e3, 500e3] # The bulkvelocity in m/s >>> vth = 30.96e3 # 30.96e3 m/s is 10 eV temperautre for proton.
In this case, the thermal pressure is \(nkT = nmv_{th}^2 = 10^6 \times 1.67\times 10^{27} \times (30.96\times 10^3)^2 = 1.60\times 10^{12}\) [Pa]. Thus, the temperature is \(T = 1.60\times 10^{12} / nk = 1.60\times 10^{12} / 10^{6} / 1.38\times 10^{23} = 1.16\times 10^5 \sim 10 eV\). Yes. It makes sense.
>>> elist = np.logspace(2, 5, 50) # fro 10 to 1000 eV with 50 grids >>> tlist_b = np.linspace(0, np.pi, 181) # From 0 to pi with 181 separation >>> tlist = (tlist_b[1:] + tlist_b[:1]) / 2. # The central grid is calculated from the boundary. 180 grid >>> plist_b = np.linspace(np.pi, np.pi, 361) >>> plist = (plist_b[1:] + plist_b[:1]) / 2. # The central grid for phi. 360 grid.
>>> maxwell = maxwellDifferentialFluxGrid(n, vb, vth, elist, tlist, plist) >>> print('{:.2f}'.format(maxwell.density())) # This should be 1 /cm^3 1.00 >>> print('{v[0]:.2f} {v[1]:.2f} {v[2]:.2f}'.format(v=maxwell.velocity_vector())) 100.00 300.00 500.00
>>> ptensor = maxwell.pressure_tensor() >>> print('{:.2e}'.format(ptensor[0][0])) # Diagonal component should be 1.60e12 Pa 1.60e12
 plot3d(ax=None, dynamic_range_scale=6, **kwds)[source]¶
Provide 3D quick look for the energy spectrum.
 Parameters:
ax – Axis object in matplotlib. This should be produced by projection=”3d” option. plt.subplot(111, projection=”3d”)
dynamic_range_scale – Specify the dynamic range. (Max is determined by the maxmum of the data, and the min is deterimned by this value.
kwds – Other keyword to ax.scatter object
 Returns:
The axis object and the scatter plot object.
To make the color bar, one has to use colorbar command with scatter plot object as argument. See example in plot3d_energyspectrum.py.
 class irfpy.util.energyspectrum.LogDifferentialFluxGrid(dflux, energy_grid, theta_grid, phi_grid, mass=1.67262192369e27, energy_grid_size='log', theta_grid_size='linear', phi_grid_size='linear')[source]¶
Bases:
RegularGridVolumeContainer
Representation of the differential flux on grids in polar coordinates, while it is lnscale.
The natural log of differential flux can be assigned at a grid of (E, \(\theta\), \(\phi\)). Also, moment calculation can be possible.
Initilizer of differential flux.
 Parameters:
dflux – LN of Differential flux. Should be a shape of (ne, nt, np). Unit of /cm2 s eV sr.
energy_grid – Energy grid in eV. Ascending order.
theta_grid – Theta (polar angle) angle array in radians. Ascending order.
phi_grid – Phi (azimuth) angle array in radians. Ascending order.
mass – Mass in kg. It is used for calculating the moment values.
energy_grid_size – Specify the energy grid size. It is used to calculate the moment values. Use ‘linear’, ‘log’, or np.array with the same shape as energy_grid.
theta_grid_size – Same as energy_grid_size but for theta direction.
phi_grid_size – Same as energy_grid_size, but for phi direction.
 UNIT = 'cm^2 s^1 eV^1 sr^1'¶
Unit of the differential flux
 lndflux¶
Differential flux
 energy_grid¶
Energy grid
 ne¶
Number of energy grid
 theta_grid¶
Theta grid (radians)
 nt¶
Number of theta grid
 phi_grid¶
Phi grid (radians)
 np¶
Number of phi grid
 energy_grid_size¶
Width of energy grid (i.e. dE)
 theta_grid_size¶
Width of theta grid (i.e. dTheta)
 phi_grid_size¶
Width of phi grid (i.e. dPhi)
 mass¶
Mass of the particles
 density()[source]¶
Calculate the density in cm^3.
Refer to
DifferentialFluxGrid.density()
.
 integral_flux()[source]¶
Calculate the integral flux in /cm2 s.
Refer to
DifferentialFluxGrid.integral_flux()
 velocity_vector()[source]¶
Calculate the velocity vector in km/s
Refer to
DifferentialFluxGrid.velocity_vector()
.
 pressure_tensor()[source]¶
Return the pressure tensor in Pa.
Example follows. >>> n = 1e6 # The density in m^3 >>> vb = [100e3, 300e3, 500e3] # The bulkvelocity in m/s >>> vth = 30.96e3 # 30.96e3 m/s is 10 eV temperautre for proton.
In this case, the thermal pressure is \(nkT = nmv_{th}^2 = 1e6 * 1.67e27 * 30.96e3^2 = 1.60e12\) [Pa]. Thus, the temperature is \(T = 1.60e12 / nk = 1.60e12 / 1e6 / 1.38e23 = 1.16e5 ~ 10 eV\). Yes. It makes sense.
>>> elist = np.logspace(2, 5, 50) # fro 10 to 1000 eV with 50 grids >>> tlist_b = np.linspace(0, np.pi, 181) # From 0 to pi with 181 separation >>> tlist = (tlist_b[1:] + tlist_b[:1]) / 2. # The central grid is calculated from the boundary. 180 grid >>> plist_b = np.linspace(np.pi, np.pi, 361) >>> plist = (plist_b[1:] + plist_b[:1]) / 2. # The central grid for phi. 360 grid.
>>> maxwell = maxwellLogDifferentialFluxGrid(n, vb, vth, elist, tlist, plist) >>> print('{:.2f}'.format(maxwell.density())) # This should be 1 /cm^3 1.00 >>> print('{v[0]:.2f} {v[1]:.2f} {v[2]:.2f}'.format(v=maxwell.velocity_vector())) 100.00 300.00 500.00
>>> ptensor = maxwell.pressure_tensor() >>> print('{:.2e}'.format(ptensor[0][0])) # Diagonal component should be 1.60e12 Pa 1.60e12
 class irfpy.util.energyspectrum.CountRateGrid(cps, energy_grid, theta_grid, phi_grid, energy_grid_size='log', theta_grid_size='linear', phi_grid_size='linear')[source]¶
Bases:
RegularGridVolumeContainer
Represent the count rate (counts per sec) on grids in polar frame.
Initilizer of differential flux.
 Parameters:
cps – Count rate, unit of s^1
energy_grid – Energy grid in eV. Ascending order.
theta_grid – Theta (polar angle) angle array in radians. Ascending order.
phi_grid – Phi (azimuth) angle array in radians. Ascending order.
mass – Mass in kg.
energy_grid_size – Specify the energy grid size. It is used to calculate the moment values. Use ‘linear’, ‘log’, or np.array with the same shape as energy_grid.
theta_grid_size – Same as energy_grid_size but for theta direction.
phi_grid_size – Same as energy_grid_size, but for phi direction.
 UNIT = '/s'¶
 cps¶
Data. Array of (ne, nt, np) shape.
 energy_grid¶
Energy grid. (ne,) shape array.
 ne¶
Number of energy grid
 theta_grid¶
Theta grid. (nt,) shape array
 nt¶
Number of theta grid.
 phi_grid¶
Phi grid. (np,) shape array
 np¶
Number of phi grid.
 energy_grid_size¶
Energy grid size, dE. (ne,) shape array
 theta_grid_size¶
Theta grid size, d:math:theta. (nt,) shape array
 phi_grid_size¶
Phi grid size, d:math:phi. (np,) shape array
 irfpy.util.energyspectrum.maxwellDifferentialFluxGrid(n, vb, vth, elist, tlist, plist, m=1.67262192369e27)[source]¶
Create
DifferentialFluxGrid
object with Maxwell distribution. Parameters:
n – Density in m^3
vb – Bulk velocity in m/s. (3,) shaped array.
vth – Thermal velocity in m/s. Scalar.
elist – List of energy grid. It should be the center values. Unit should be
eV
tlist – List of center of the theta grid values (radians)
plist – List of center of the phi grid (radians)
m – Mass for the particle in kg. Default is proton mass.
 Returns:
The
DifferentialFluxGrid
object
This function creates the differential flux object,
DifferentialFluxGrid
, representing the Maxwell distribution. Given by the centers of the grid points, the differential fluxJ(E, theta, phi)
is expressed in the unit ofcm^2 s^1 sr^1 eV^1
.For example, prepare the following characteristics.
>>> n = 1e7 # 10 /cm^3 >>> vb = [150e3, 1, 1] # 150 km/s, mostly vx. >>> vth = 30.96e3 # T=10 eV for proton >>> elist = np.logspace(1, 3, 50) # fro 10 to 1000 eV with 50 grids >>> tlist_b = np.linspace(0, np.pi, 181) # From 0 to pi with 181 separation >>> tlist = (tlist_b[1:] + tlist_b[:1]) / 2. # The central grid is calculated from the boundary. 180 grid >>> plist_b = np.linspace(np.pi, np.pi, 361) >>> plist = (plist_b[1:] + plist_b[:1]) / 2. # The central grid for phi. 360 grid.
Then, create the differential flux grid.
>>> df_mx = maxwellDifferentialFluxGrid(n, vb, vth, elist, tlist, plist)
You can get the differential flux matrix.
>>> df = df_mx.dflux >>> print('{:.4e}'.format(df.max())) # Maximum differential flux in /cm2 s sr eV 4.9750e+06
You can get the moment by integration.
>>> print('{:.4e}'.format(df_mx.density())) # in /cm3 1.0003e+01 >>> print('{v[0]:.4e} {v[1]:.4e} {v[2]:.4e}'.format(v=df_mx.velocity_vector())) # in km/s 1.5001e+02 1.0000e03 1.0000e03
Theory: The differential flux j for a Maxwellian is
\[j(E, \theta, \phi) = \frac{2n_0}{m^2}\left(\frac{m}{2\pi kT}\right)^\frac{3}{2}E \exp\left(\frac{m}{2kT}\left((v_xV_x)^2+(v_yV_y)^2+(v_zV_z)^2\right)\right)\]in MKSA. Note \(v_th=\sqrt{kT / m}\).
Changed in version v4.4.13a2: The algorithm to calculate the differential flux has been changed for better precision and better performance. The results may change from the previous version.
 irfpy.util.energyspectrum.maxwellLogDifferentialFluxGrid(n, vb, vth, elist, tlist, plist, m=1.67262192369e27)[source]¶
Create
LogDifferentialFluxGrid
object with Maxwell distribution. Parameters:
n – Density in m^3
vb – Bulk velocity in m/s. (3,) shaped array.
vth – Thermal velocity in m/s. Scalar.
elist – List of energy grid. It should be the center values. Unit should be
eV
tlist – List of center of the theta grid values (radians)
plist – List of center of the phi grid (radians)
m – Mass for the particle in kg. Default is proton mass.
 Returns:
The
DifferentialFluxGrid
object
This function creates the differential flux object,
DifferentialFluxGrid
, representing the Maxwell distribution. Given by the centers of the grid points, the natural log of the differential fluxlog J(E, theta, phi)
is expressed in the unit ofcm^2 s^1 sr^1 eV^1
.For example, prepare the following characteristics.
>>> n = 1e7 # 10 /cm^3 >>> vb = [150e3, 1, 1] # 150 km/s, mostly vx. >>> vth = 30.96e3 # T=10 eV for proton >>> elist = np.logspace(1, 3, 50) # fro 10 to 1000 eV with 50 grids >>> tlist_b = np.linspace(0, np.pi, 181) # From 0 to pi with 181 separation >>> tlist = (tlist_b[1:] + tlist_b[:1]) / 2. # The central grid is calculated from the boundary. 180 grid >>> plist_b = np.linspace(np.pi, np.pi, 361) >>> plist = (plist_b[1:] + plist_b[:1]) / 2. # The central grid for phi. 360 grid.
Then, create the differential flux grid.
>>> df_mx = maxwellLogDifferentialFluxGrid(n, vb, vth, elist, tlist, plist)
You can get the differential flux matrix, but in log value.
>>> df = df_mx.lndflux >>> print('{:.2f}'.format(df.max())) # Maximum differential flux in /cm2 s sr eV, but log value 15.42 >>> print('{:.4e}'.format(np.exp(df.max()))) # maximum differntial flux in /cm2 s sr eV. 4.9750e+06
>>> print('{:.1f}'.format(df.min())) # Minimum value exists, though it is extremely small. 162.7
You can get the moment by integration.
>>> print('{:.4e}'.format(df_mx.density())) # in /cm3 1.0003e+01 >>> print('{v[0]:.4e} {v[1]:.4e} {v[2]:.4e}'.format(v=df_mx.velocity_vector())) # in km/s 1.5001e+02 1.0000e03 1.0000e03
Theory: The differential flux j for a Maxwellian is
\[j(E, \theta, \phi) = \frac{2n_0}{m^2}\left(\frac{m}{2\pi kT}\right)^\frac{3}{2}E \exp\left(\frac{m}{2kT}\left((v_xV_x)^2+(v_yV_y)^2+(v_zV_z)^2\right)\right)\]in MKSA. Note \(v_th=\sqrt{kT / m}\).
 irfpy.util.energyspectrum.guess_grid_size_phi(pgrid, descend=False)[source]¶
Guess the grid size for the phi table.
The table is considered as periodic with respective to 360. This means that the given table is in degrees. In case the table is in descending order, one have to specify it via the
descend
keyword.>>> pgrid = [340, 350, 0, 10, 20] >>> print(guess_grid_size_phi(pgrid)) [10. 10. 10. 10. 10.]
>>> pgrid = [20, 10, 0, 350, 340] # Descending order case, >>> print(guess_grid_size_phi(pgrid, descend=True)) [10. 10. 10. 10. 10.]
If you miss the descending order keyword, the value should be messy…
>>> print(guess_grid_size_phi(pgrid)) [350. 350. 350. 350. 350.]
 irfpy.util.energyspectrum.guess_grid_size_theta(tgrid)[source]¶
Guess the grid size for the theta table.
The table must be sorted (either descending or ascending) by user’s responsibility.
>>> tgrid = [20, 0, 20, 40, 60] >>> tgridsize = guess_grid_size_theta(tgrid) >>> for s in tgridsize: ... print('{:.2f}'.format(s)) 20.00 20.00 20.00 20.00 20.00
 irfpy.util.energyspectrum.guess_gridsize_energy(egrid)[source]¶
Guess the grid size for the energy table.
Energy table must be sorted (either descending or ascending order) by user’s responsibility.
For descending case:
>>> egrid = [1, 10, 100, 1000, 10000] >>> egridsize = guess_gridsize_energy(egrid) >>> for egridsize_i in egridsize: ... print('{:.2e}'.format(egridsize_i)) 2.85e+00 2.85e+01 2.85e+02 2.85e+03 2.85e+04
For ascending case,
>>> egrid = [10000, 1000, 100, 10, 1] >>> egridsize = guess_gridsize_energy(egrid) >>> for egridsize_i in egridsize: ... print('{:.2e}'.format(egridsize_i)) 2.85e+04 2.85e+03 2.85e+02 2.85e+01 2.85e+00