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, g-factor information should be available. For more serious data analysis, the response of the sensor is also needed. These “instrument model” is not straightforward-ly 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 z-axis and the angle \(\phi\) as the angle from x axis as in the projection into x-y plane.

DifferentialFluxGrid(dflux, energy_grid, ...)

Representation of the differential flux on grids in polar coordinates.

CountRateGrid(cps, energy_grid, theta_grid, ...)

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

maxwellDifferentialFluxGrid(n, vb, vth, ...)

Create DifferentialFluxGrid object with Maxwell distribution.
class irfpy.util.energyspectrum.DifferentialFluxGrid(dflux, energy_grid, theta_grid, phi_grid, mass=1.67262192369e-27, 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 denisyt density().

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 i-th 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.60e-12 Pa
1.60e-12
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.67262192369e-27, 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 ln-scale.

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.67e-27 * 30.96e3^2 = 1.60e-12\) [Pa]. Thus, the temperature is \(T = 1.60e-12 / nk = 1.60e-12 / 1e6 / 1.38e-23 = 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.60e-12 Pa
1.60e-12
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.67262192369e-27)[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 flux J(E, theta, phi) is expressed in the unit of cm^-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.0000e-03 1.0000e-03

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_x-V_x)^2+(v_y-V_y)^2+(v_z-V_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.67262192369e-27)[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 flux log J(E, theta, phi) is expressed in the unit of cm^-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.0000e-03 1.0000e-03

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_x-V_x)^2+(v_y-V_y)^2+(v_z-V_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