r""" Generic purpose energy spectrum implementation
.. codeauthor:: Yoshifumi Futaana
:class:`DifferentialFluxGrid` represents the differential flux on spherical grids.
:class:`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 :mod:`irfpy.util.fluxtools`.
A counter part (velocity distribution function) is implemented
in :mod:`irfpy.util.vdf` module. :mod:`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 :math:`\theta` as the angle from
z-axis and the angle :math:`\phi` as the angle from
x axis as in the projection into x-y plane.
.. autosummary::
DifferentialFluxGrid
CountRateGrid
maxwellDifferentialFluxGrid
"""
import numpy as np
from irfpy.util.polarframe import RegularGridVolumeContainer
from irfpy.util import exception as ex
from scipy import constants as _c
[docs]class DifferentialFluxGrid(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.
"""
UNIT = "cm^-2 s^-1 eV^-1 sr^-1"
""" Unit of the differential flux """
def __init__(self, dflux, energy_grid, theta_grid, phi_grid, mass=_c.proton_mass,
energy_grid_size='log', theta_grid_size='linear', phi_grid_size='linear'):
r"""
:param dflux: Differential flux. Should be a shape of (ne, nt, np). Unit of /cm2 s eV sr.
:param energy_grid: Energy grid in eV. Ascending order.
:param theta_grid: Theta (polar angle) angle array in radians. Ascending order.
:param phi_grid: Phi (azimuth) angle array in radians. Ascending order.
:keyword mass: Mass in kg. It is used for calculating the moment values.
:keyword 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.
:keyword theta_grid_size: Same as energy_grid_size but for theta direction.
:keyword phi_grid_size: Same as energy_grid_size, but for phi direction.
.. note::
(Note on 190118):
For integrators, we use simple approach (see :meth:`integral_flux`).
One can have better answer by integrating :math:`\sin\theta d\theta` properly.
Will see how we should develop.
"""
RegularGridVolumeContainer.__init__(self, np.ma.masked_array(dflux), energy_grid, theta_grid, phi_grid,
rgridsize=energy_grid_size, tgridsize=theta_grid_size, pgridsize=phi_grid_size)
self.dflux = self.data # For alias.
""" Differential flux"""
self.energy_grid = self.rgrid # Alias
""" Energy grid"""
self.ne = self.energy_grid.shape[0]
""" Number of energy grid"""
self.theta_grid = self.tgrid # Alias
""" Theta grid (radians)"""
self.nt = self.theta_grid.shape[0]
""" Number of theta grid"""
self.phi_grid = self.pgrid # Alias
""" Phi grid (radians)"""
self.np = self.phi_grid.shape[0]
""" Number of phi grid"""
self.energy_grid_size = self.rgridsize # Alias
""" Width of energy grid (i.e. dE) """
self.theta_grid_size = self.tgridsize # Alias
""" Width of theta grid (i.e. dTheta)"""
self.phi_grid_size = self.pgridsize # Alias
""" Width of phi grid (i.e. dPhi) """
if self.dflux.shape != (self.ne, self.nt, self.np):
raise ex.IrfpyException(
"Given dflux information is not consistent with grid information: {} != {}, {}, {}".format(
self.dflux.shape, self.ne, self.nt, self.np))
self.mass = mass
""" Mass of the particles"""
self._to3d()
def _to3d(self):
""" Calculate the 1-D vector (table) to 3-D array in MKSA.
"""
qe = 1.60217657e-19 # Elementary charge
j = self.dflux # in cm^-2 sr^-1 s^-1 eV^-1
self._j2 = j * 10000 / qe # in m^-2 sr^-1 s^-1 J^-1
e = self.energy_grid # in eV!, (nE) shape
eJ = e * qe # Now in Joule = kg m^2 s^-2, (nE) shape
# eJ2 = np.vstack([eJ] * self.nt).T # (nE, nT) shape
# self._eJ2 = np.dstack([eJ2] * self.np) # (nE, nT, nP) shape
# t = np.vstack([self.theta_grid] * self.ne) # (nE, nT) shape
# t2 = np.dstack([t] * self.np) # (nE, nT, nP) shape
# self._t2 = t2
# self._sint2 = np.sin(self._t2)
self._eJ2, self._t2, self._p2 = np.meshgrid(eJ, self.theta_grid, self.phi_grid, indexing='ij')
de = self.energy_grid_size # (nE,) shape
deJ = de * qe # Now in Joule.
dt = self.theta_grid_size # (nT,) shape
dp = self.phi_grid_size # (nP,) shape
de3d, dt3d, dp3d = np.meshgrid(deJ, dt, dp, indexing='ij') # (nT, nE, nP) shape for each.
self._de2, self._dt2, self._dp2 = de3d, dt3d, dp3d
[docs] def density(self):
r""" Calculate the density by integrating over the valid counts.
:returns: Density in cm^-3
Calculating density from differential flux is formulated as
.. math::
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
.. math::
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):
.. math::
dv_x dv_y dv_z = \sqrt{\frac{2E}{m^3}}\sin\theta
and the relation
.. math::
j = \frac{2E}{m^2} f
"""
mass = self.mass
eJ2 = self._eJ2 # in MKSA / (nE, nT, nP) array
sqrt_part = np.sqrt(mass / 2 / eJ2)
j2 = self._j2 # in MKSA / (nE, nT, nP) array
sintheta = np.sin(self._t2)
de3d, dt3d, dp3d = self._de2, self._dt2, self._dp2
dn = sqrt_part * j2 * sintheta * de3d * dt3d * dp3d
# The energy vector may contain nan, do not sum these:
np.ma.masked_invalid(dn)
n = np.nansum(dn.data) * 1e-6 ## m^-3 to cm^-3
return n
[docs] def integral_flux(self):
r"""
:return: Integral flux. Unit is /cm2 s.
The integral is, in MKSA unit, expressed by
.. math::
\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
"""
j = self._j2 # in MKSA
e = self._eJ2 # in MKSA
t = self._t2
p = self._p2
cost = np.cos(t)
sint = np.sin(t)
cosp = np.cos(p)
sinp = np.sin(p)
de = self._de2 # in MKSA
dt = self._dt2
dp = self._dp2
jx = j * sint * cosp
jy = j * sint * sinp
jz = j * cost
dfx = jx * de * dt * dp * sint
# The energy vector may contain nan, do not sum these:
np.ma.masked_invalid(dfx)
fx = np.nansum(dfx.data) # in MKSA
dfy = jy * de * dt * dp * sint
np.ma.masked_invalid(dfy)
fy = np.nansum(dfy.data)
dfz = jz * de * dt * dp * sint
np.ma.masked_invalid(dfz)
fz = np.nansum(dfz.data)
f = np.array([fx, fy, fz])
f = f * 1e-4 # To make /m2 s to /cm2 s
return f
[docs] def velocity_vector(self):
""" Returns the velocity vector in km/s
:returns: (vx, vy, vz) in km/s.
The velocity vector is
.. math::
F / n
where F is the :meth:`integral_flux` and n is the denisyt :meth:`density`.
"""
flux = self.integral_flux() # /cm2 s
density = self.density() # /cm3
vel = flux / density # cm / s
return vel / 1e5 # in km/s
[docs] def pressure_tensor(self):
r""" Return the pressure tensor in Pa.
The pressure tensor is defined by
.. math::
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 :math:`v_i = v a_i`, where *a_i* is the unit vector of the velocity vector,
we can rewrite the pressure as
.. math::
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 :math:`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 :math:`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
"""
m = self.mass # in kg
e = self._eJ2
t = self._t2
p = self._p2
de = self._de2
dt = self._dt2
dp = self._dp2
ax = np.cos(p) * np.sin(t)
ay = np.sin(p) * np.sin(t)
az = np.cos(t)
sqrt_part = np.sqrt(2 * e / m) # MKSA
sint = np.sin(t)
j = self._j2 # MKSA
coeff = sqrt_part * sint * j * de * dt * dp
pxx = (ax * ax * coeff).sum()
pyy = (ay * ay * coeff).sum()
pzz = (az * az * coeff).sum()
pxy = (ax * ay * coeff).sum()
pyz = (ay * az * coeff).sum()
pzx = (az * ax * coeff).sum()
ptensor = np.array([[pxx, pxy, pzx],
[pxy, pyy, pyz],
[pzx, pyz, pzz]]) * m
n = self.density() * 1e6 # In m^-3
v = self.velocity_vector() * 1e3 # in m/s
Vx, Vy, Vz = v
ptensor2 = np.array([[Vx * Vx, Vx * Vy, Vx * Vz],
[Vy * Vx, Vy * Vy, Vy * Vz],
[Vz * Vx, Vz * Vy, Vz * Vz]]) * m * n
return ptensor - ptensor2
[docs] def plot3d(self, ax=None, dynamic_range_scale=6, **kwds):
""" Provide 3D quick look for the energy spectrum.
:param ax: Axis object in matplotlib. This should be produced by
projection="3d" option. **plt.subplot(111, projection="3d")**
:param dynamic_range_scale: Specify the dynamic range. (Max is determined by
the maxmum of the data, and the min is deterimned by this value.
:param kwds: Other keyword to ax.scatter object
:return: 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*.
"""
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
if ax is None:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
elist = np.log10(self.energy_grid.copy())
emin = int(elist.min()) - 1
emax = int(elist.max()) + 1
elist = elist - emin
tlist = self.theta_grid
plist = self.phi_grid
e, t, p = np.meshgrid(elist, tlist, plist, indexing='ij')
ex = e * np.sin(t) * np.cos(p)
ey = e * np.sin(t) * np.sin(p)
ez = e * np.cos(t)
jlist = np.log10(self.dflux)
jmax = int(jlist.max()) + 1
jmin = jmax - dynamic_range_scale # Dynamic range to be 6 as default.
clist = jlist.clip(jmin, jmax)
clist2 = clist.copy()
clist = (clist - jmin) / (jmax - jmin) * 255.
slist = clist / 255. * 50.
# sct = ax.scatter(ex, ey, ez, c=clist, s=slist, **kwds)
sct = ax.scatter(ex, ey, ez, vmin=jmin, vmax=jmax, c=clist2, s=slist, **kwds)
ax.set_xlabel("log(Ex) - {}(offset) [eV]".format(emin))
ax.set_ylabel("log(Ey) - {}(offset) [eV]".format(emin))
ax.set_zlabel("log(Ez) - {}(offset) [eV]".format(emin))
return ax, sct
[docs]class LogDifferentialFluxGrid(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.
"""
UNIT = "cm^-2 s^-1 eV^-1 sr^-1"
""" Unit of the differential flux """
def __init__(self, dflux, energy_grid, theta_grid, phi_grid, mass=_c.proton_mass,
energy_grid_size='log', theta_grid_size='linear', phi_grid_size='linear'):
""" Initilizer of differential flux.
:param dflux: LN of Differential flux. Should be a shape of (ne, nt, np). Unit of /cm2 s eV sr.
:param energy_grid: Energy grid in eV. Ascending order.
:param theta_grid: Theta (polar angle) angle array in radians. Ascending order.
:param phi_grid: Phi (azimuth) angle array in radians. Ascending order.
:keyword mass: Mass in kg. It is used for calculating the moment values.
:keyword 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.
:keyword theta_grid_size: Same as energy_grid_size but for theta direction.
:keyword phi_grid_size: Same as energy_grid_size, but for phi direction.
"""
RegularGridVolumeContainer.__init__(self, np.ma.masked_array(dflux), energy_grid, theta_grid, phi_grid,
rgridsize=energy_grid_size, tgridsize=theta_grid_size, pgridsize=phi_grid_size)
self.lndflux = self.data # For alias.
""" Differential flux"""
self.energy_grid = self.rgrid # Alias
""" Energy grid"""
self.ne = self.energy_grid.shape[0]
""" Number of energy grid"""
self.theta_grid = self.tgrid # Alias
""" Theta grid (radians)"""
self.nt = self.theta_grid.shape[0]
""" Number of theta grid"""
self.phi_grid = self.pgrid # Alias
""" Phi grid (radians)"""
self.np = self.phi_grid.shape[0]
""" Number of phi grid"""
self.energy_grid_size = self.rgridsize # Alias
""" Width of energy grid (i.e. dE) """
self.theta_grid_size = self.tgridsize # Alias
""" Width of theta grid (i.e. dTheta)"""
self.phi_grid_size = self.pgridsize # Alias
""" Width of phi grid (i.e. dPhi) """
if self.lndflux.shape != (self.ne, self.nt, self.np):
raise ex.IrfpyException(
"Given dflux information is not consistent with grid information: {} != {}, {}, {}".format(
self.lndflux.shape, self.ne, self.nt, self.np))
self.mass = mass
""" Mass of the particles"""
self._to3d()
self.linearlized = self._linearlize()
def _linearlize(self):
d = DifferentialFluxGrid(np.exp(self.lndflux), self.energy_grid, self.theta_grid, self.phi_grid, mass=self.mass,
energy_grid_size=self.energy_grid_size, theta_grid_size=self.theta_grid_size, phi_grid_size=self.phi_grid_size)
return d
def _to3d(self):
""" Calculate the 1-D vector (table) to 3-D array in MKSA.
"""
qe = 1.60217657e-19 # Elementary charge
j = self.lndflux # in cm^-2 sr^-1 s^-1 eV^-1
self._j2 = j + np.log(10000 / qe) # in m^-2 sr^-1 s^-1 J^-1
e = self.energy_grid # in eV!, (nE) shape
eJ = e * qe # Now in Joule = kg m^2 s^-2, (nE) shape
self._eJ2, self._t2, self._p2 = np.meshgrid(eJ, self.theta_grid, self.phi_grid, indexing='ij')
de = self.energy_grid_size # (nE,) shape
deJ = de * qe # Now in Joule.
dt = self.theta_grid_size # (nT,) shape
dp = self.phi_grid_size # (nP,) shape
de3d, dt3d, dp3d = np.meshgrid(deJ, dt, dp, indexing='ij') # (nT, nE, nP) shape for each.
self._de2, self._dt2, self._dp2 = de3d, dt3d, dp3d
[docs] def density(self):
""" Calculate the density in cm^-3.
Refer to :meth:`DifferentialFluxGrid.density`.
"""
return self.linearlized.density()
[docs] def integral_flux(self):
r""" Calculate the integral flux in /cm2 s.
Refer to :meth:`DifferentialFluxGrid.integral_flux`
"""
return self.linearlized.integral_flux()
[docs] def velocity_vector(self):
""" Calculate the velocity vector in km/s
Refer to :meth:`DifferentialFluxGrid.velocity_vector`.
"""
return self.linearlized.velocity_vector()
[docs] def pressure_tensor(self):
r""" 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 :math:`nkT = nmv_{th}^2 = 1e6 * 1.67e-27 * 30.96e3^2 = 1.60e-12` [Pa].
Thus, the temperature is :math:`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
"""
return self.linearlized.pressure_tensor()
[docs]class CountRateGrid(RegularGridVolumeContainer):
""" Represent the count rate (counts per sec) on grids in polar frame.
"""
UNIT= "/s"
def __init__(self, cps, energy_grid, theta_grid, phi_grid,
energy_grid_size='log', theta_grid_size='linear', phi_grid_size='linear'):
""" Initilizer of differential flux.
:param cps: Count rate, unit of s^-1
:param energy_grid: Energy grid in eV. Ascending order.
:param theta_grid: Theta (polar angle) angle array in radians. Ascending order.
:param phi_grid: Phi (azimuth) angle array in radians. Ascending order.
:keyword mass: Mass in kg.
:keyword 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.
:keyword theta_grid_size: Same as energy_grid_size but for theta direction.
:keyword phi_grid_size: Same as energy_grid_size, but for phi direction.
"""
RegularGridVolumeContainer.__init__(self, np.ma.masked_array(cps), energy_grid, theta_grid, phi_grid,
rgridsize=energy_grid_size, tgridsize=theta_grid_size, pgridsize=phi_grid_size)
self.cps = self.data # Alias
""" Data. Array of (ne, nt, np) shape.
"""
self.energy_grid = self.rgrid # Alias
""" Energy grid. (ne,) shape array.
"""
self.ne = self.energy_grid.shape[0]
""" Number of energy grid"""
self.theta_grid = self.tgrid # Alias
""" Theta grid. (nt,) shape array"""
self.nt = self.theta_grid.shape[0]
""" Number of theta grid."""
self.phi_grid = self.pgrid # Alias
""" Phi grid. (np,) shape array"""
self.np = self.phi_grid.shape[0]
""" Number of phi grid. """
self.energy_grid_size = self.rgridsize # Alias
""" Energy grid size, dE. (ne,) shape array """
self.theta_grid_size = self.tgridsize # Alias
r""" Theta grid size, d:math:`\theta`. (nt,) shape array """
self.phi_grid_size = self.pgridsize # Alias
r""" Phi grid size, d:math:`\phi`. (np,) shape array """
if self.cps.shape != (self.ne, self.nt, self.np):
raise ex.IrfpyException(
"Given dflux information is not consistent with grid information: {} != {}, {}, {}".format(
self.cps.shape, self.ne, self.nt, self.np))
[docs]def maxwellDifferentialFluxGrid(n, vb, vth, elist, tlist, plist, m=_c.proton_mass):
r""" Create :class:`DifferentialFluxGrid` object with Maxwell distribution.
:param n: Density in m^-3
:param vb: Bulk velocity in m/s. (3,) shaped array.
:param vth: Thermal velocity in m/s. Scalar.
:param elist: List of energy grid. It should be the center values. Unit should be ``eV``
:param tlist: List of center of the theta grid values (radians)
:param plist: List of center of the phi grid (radians)
:keyword m: Mass for the particle in kg. Default is proton mass.
:return: The :class:`DifferentialFluxGrid` object
This function creates the differential flux object, :class:`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
.. math::
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 :math:`v_th=\sqrt{kT / m}`.
.. versionchanged:: 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.
"""
e = _c.elementary_charge
elist_c = np.array(elist)
kT = m * (vth ** 2)
elist_J = elist_c * e
vlist = np.sqrt(2 * elist_J / m)
_n = np.newaxis # For alias
vx = vlist[:, _n, _n] * np.sin(tlist[_n, :, _n]) * np.cos(plist[_n, _n, :])
vy = vlist[:, _n, _n] * np.sin(tlist[_n, :, _n]) * np.sin(plist[_n, _n, :])
vz = vlist[:, _n, _n] * np.cos(tlist[_n, :, _n])
coeff0 = 2 * n / (m ** 2)
coeff1 = m / (2 * np.pi * kT)
coeff1 = np.power(coeff1, 1.5)
coeff2 = elist_J[:, _n, _n]
incoeff0 = m / (2 * kT)
incoeff1 = (vx - vb[0]) ** 2 + (vy - vb[1]) ** 2 + (vz - vb[2]) ** 2
jval = coeff0 * coeff1 * coeff2 * np.exp(-incoeff0 * incoeff1) * 1e-4 * e
dgrid = DifferentialFluxGrid(jval, elist_c, tlist, plist, mass=m)
return dgrid
[docs]def maxwellLogDifferentialFluxGrid(n, vb, vth, elist, tlist, plist, m=_c.proton_mass):
r""" Create :class:`LogDifferentialFluxGrid` object with Maxwell distribution.
:param n: Density in m^-3
:param vb: Bulk velocity in m/s. (3,) shaped array.
:param vth: Thermal velocity in m/s. Scalar.
:param elist: List of energy grid. It should be the center values. Unit should be ``eV``
:param tlist: List of center of the theta grid values (radians)
:param plist: List of center of the phi grid (radians)
:keyword m: Mass for the particle in kg. Default is proton mass.
:return: The :class:`DifferentialFluxGrid` object
This function creates the differential flux object, :class:`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
.. math::
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 :math:`v_th=\sqrt{kT / m}`.
"""
e = _c.elementary_charge
elist_c = np.array(elist)
kT = m * (vth ** 2)
log_kT = np.log(kT)
elist_J = elist_c * e
vlist = np.sqrt(2 * elist_J / m)
_n = np.newaxis # For alias
vx = vlist[:, _n, _n] * np.sin(tlist[_n, :, _n]) * np.cos(plist[_n, _n, :])
vy = vlist[:, _n, _n] * np.sin(tlist[_n, :, _n]) * np.sin(plist[_n, _n, :])
vz = vlist[:, _n, _n] * np.cos(tlist[_n, :, _n])
log_m = np.log(m)
log_2 = np.log(2)
log_pi = np.log(np.pi)
coeff0 = np.log(2) + np.log(n) - 2 * log_m
coeff1 = (log_m - log_2 - log_pi - log_kT) * 1.5
coeff2 = np.log(elist_J)[:, _n, _n]
incoeff0 = m / (2 * kT)
incoeff1 = (vx - vb[0]) ** 2 + (vy - vb[1]) ** 2 + (vz - vb[2]) ** 2
jval = coeff0 + coeff1 + coeff2 -incoeff0 * incoeff1 + np.log(1e-4) + np.log(e)
dgrid = LogDifferentialFluxGrid(jval, elist_c, tlist, plist, mass=m)
return dgrid
[docs]def guess_grid_size_phi(pgrid, descend=False):
""" 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)) # doctest: +NORMALIZE_WHITESPACE
[10. 10. 10. 10. 10.]
>>> pgrid = [20, 10, 0, 350, 340] # Descending order case,
>>> print(guess_grid_size_phi(pgrid, descend=True)) # doctest: +NORMALIZE_WHITESPACE
[10. 10. 10. 10. 10.]
If you miss the descending order keyword, the value should be messy...
>>> print(guess_grid_size_phi(pgrid)) # doctest: +NORMALIZE_WHITESPACE
[350. 350. 350. 350. 350.]
"""
ppgrid = np.array(pgrid, dtype=None)
nppgrid = len(ppgrid)
if descend:
ppgrid = ppgrid[::-1] # Now, ppgrid is assumed to ascending
for idx in range(nppgrid - 1):
while ppgrid[idx + 1] - ppgrid[idx] < 0: # If the next argument is smaller
ppgrid[idx + 1:] += 360
upper = np.zeros_like(ppgrid, dtype=float)
lower = np.zeros_like(ppgrid, dtype=float)
upper[:-1] = (ppgrid[:-1] + ppgrid[1:]) / 2.0
lower[1:] = (ppgrid[:-1] + ppgrid[1:]) / 2.0
meansize = np.median((upper[1:-1] - lower[1:-1])) / 2.0
if meansize > 180:
import warnings
warnings.warn("guess_grid_size_phi: Guessed size too big ({:.3f}). May be descend option needed? Continue processing though.".format(meansize))
upper[-1] = ppgrid[-1] + meansize
lower[0] = ppgrid[0] - meansize
size = upper - lower
return size
[docs]def guess_grid_size_theta(tgrid):
""" 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
"""
ttgrid = np.array(tgrid)
upper = np.zeros_like(ttgrid, dtype=float)
lower = np.zeros_like(ttgrid, dtype=float)
### Convert to descend
if ttgrid[0] < ttgrid[1]:
descend = False
ttgrid = ttgrid[::-1]
else:
descend = True
upper[1:] = (ttgrid[:-1] + ttgrid[1:]) / 2.0
lower[:-1] = (ttgrid[:-1] + ttgrid[1:]) / 2.0
meansize = np.median((upper[1:-1] - lower[1:-1])) / 2.0
# if meansize < 0:
# upper, lower = lower, upper
# meansize = -meansize
upper[0] = ttgrid[0] + meansize
lower[-1] = ttgrid[-1] - meansize
size = upper - lower
if not descend:
size = size[::-1]
return size
[docs]def guess_gridsize_energy(egrid):
""" 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
"""
upper = np.zeros_like(egrid, dtype=float)
lower = np.zeros_like(egrid, dtype=float)
logegrid = np.log10(egrid)
### Convert to descend
if logegrid[0] < logegrid[1]:
descend = False
logegrid = logegrid[::-1]
else:
descend = True
upper[1:] = (logegrid[:-1] + logegrid[1:]) / 2.0
lower[:-1] = (logegrid[:-1] + logegrid[1:]) / 2.0
meansize = np.median((upper[1:-1] - lower[1:-1])) / 2.0
# if meansize < 0:
# upper, lower = lower, upper
# meansize = -meansize
upper[0] = logegrid[0] + meansize
lower[-1] = logegrid[-1] - meansize
upper = np.power(10, upper)
lower = np.power(10, lower)
size = upper - lower
if not descend:
size = size[::-1]
return size