irfpy.util.maxwell

Module collecting maxwell distribution function related functions.

In plasma physics, Maxwell distribution function appears everywhere. In this module, Maxwell distribution functions are implemented.

Code author: Yoshifumi Futaana

mkfunc(n0, v0, vth[, multiarray])

Returns a function of the Maxwellian.

mkfunc_u(n0_u, v0_u, vth_u)

Unitful version of mkfunc().

mk1d(A, B, C)

Returns a callable of 1-D expression of Maxwellian distribution.

mkslice(n0, v0, vth)

Returns a callable of 1-D slice of the distribution function.

fit1d(x_array, y_array)

Fitting by 1-dimensional Maxwell distiribution function.

getRandom(vel0, vth[, n])

Get a random distributed numpy.array of the 3-D Maxwellian.

getRandom1D(v0, vth[, n])

Get a random distributed numpy.array of values for 1-D Maxwellian.

flux_integrated(n0, v0, vth[, vzmin])

Integrate the maxwellian, with a specific minimum range.

flux_precip(n0, vz0, vth)

Return the precipitating flux along +z.

irfpy.util.energyspectrum.maxwellDifferentialFluxGrid(n, ...)

Create DifferentialFluxGrid object with Maxwell distribution.

irfpy.util.vdf.maxwellVelocityDistributionGrid(n, ...)

Create VelocityDistributionGrid object with Maxwell distribution.

The theory starts from:

\[f(\vec{v})=n_0\left(\frac{m}{2\pi k}\right)^\frac{3}{2} \frac{1}{\sqrt{\mathrm{det} \textbf{T}}} \exp-\frac{m}{2k}(\vec{v}-\vec{v_b})^T \textbf{T}^{-1} (\vec{v}-\vec{v_b})\]

When the temperature is isotropic, the above temperature tensor can be treated as a scalar, becoming

\[f(\vec{v})=n_0\left(\frac{m}{2\pi kT}\right)^\frac{3}{2} \exp-\frac{m}{2kT}(\vec{v}-\vec{v_b})^2\]

Using the definition of the thermal speed,

\[v_\mathrm{th}^2 = \frac{kT}{m}\]
\[f(v_x,v_y,v_z) = n_0 ( 2\pi v_\mathrm{th}^2 ) ^{-\frac{3}{2}} \exp -\frac{ (v_x-v_{x0})^2 + (v_y-v_{y0})^2 + (v_z-v_{z0})^2 } {2 v_\mathrm{th}^2}\]

Cookbook

Calculate Maxwellian VDF

Making 3-D Maxwellian VDF with plasma parameters 5/cc, vx = -400km/s and vth=40km/s

>>> fsw = mkfunc(5e6, [-400e3, 0, 0], 40e3)

You can calculate the value of VDF as

>>> print('%.2e' % fsw([-400e3, 0, 0]))    # The peak
4.96e-09

Be careful with the unit system. Always use MKSA for mkfunc(). Alternatively you can use mkfunc_u() to handle the unit.

>>> import numpy as np
>>> fsw_u = mkfunc_u(5 / u.cc, np.array([-400, 0, 0]) * u.km/u.s, 40 * u.km / u.s)
>>> fsw_peak = fsw_u(np.array([-400, 0, 0]) * u.km / u.s)
>>> print("{:.2e}".format(fsw_peak))
4.96e-09 s**3/m**6

Maxwellian flux object

Using MaxwellDistributionG, the Maxwellian VDF and DFlux is mapped onto energy/velocity spaces.

Fitting by Maxwellian

Fitting of Maxwellian from a dataset \(x=\{x_i\}, y=\{y_i\}\) can be calculated via fit1d() function.

x = some_array
y = some_array
(A,B,C), success = fit1d(x,y)           # A * ( 2pi B^2)  exp - (v-C)^2/B^2
ffit = mk1d(A,B,C)
yfit = ffit(x)
plot(x,y,'-', x,yfit, 'o')

Creating a random value set

Sometimes, you want to get set of particles whose distribution function follows the given Maxwellian.

vels = getRandom([-400e3, 0, 0], 40e3, n=100)

Calculating integral flux

The net flux of Maxwellian gas but with a given lower velocity threshold can be calculated analitycally. This is useful for example to calculate the Langmuir probe or Retarding Potential Analyzer performance, and to estimate the outflowing flux (not including inward flux).

\[F_{v_z>v_{zmin}} = \int_{-\infty}^{\infty}dv_x\int_{-\infty}^{\infty}dv_y\int_{v_{zmin}}^{\infty}d_z[v_zf(v_x, v_y, v_z)]\]

where the precipitating direction is assumed (without losing generality) as \(+v_z\) direction.

The integral will result in

\[F = \frac{n_0v_{th}}{\sqrt{2\pi}}\exp-\frac{(v_{zmin}-v_{z0})^2}{2v_{th}^2} + \frac{n_0v_{z0}}{2}\mathrm{erfc}\frac{v_{zmin}-v_{z0}}{\sqrt{2}{v_{th}}}\]

This is implemented in flux_integrated().

Some useful formula

Some formula useful for Maxwellian are summarized here.

\[I^{(n)}(v) \equiv \int x^n\exp\Bigl(-\frac{v^2}{2\sigma^2}\Bigr) dv\]
\[\begin{split}I^{(0)}(v) = \int \exp\Bigl(-\frac{v^2}{2\sigma^2}\Bigr) dv &=& \sigma\sqrt\frac{\pi}{2} \mathrm{erf}\Bigl(\frac{v}{\sqrt{2}\sigma}\Bigr) + C \\ I^{(1)}(v) = \int v \exp\Bigl(-\frac{v^2}{2\sigma^2}\Bigr) dv &=& -\sigma^2\exp\Bigl(-\frac{v^2}{2\sigma^2}\Bigr) + C \\ I^{(2)}(v) = \int v^2 \exp\Bigl(-\frac{v^2}{2\sigma^2}\Bigr) dv &=& \sqrt\frac{\pi}{2} \sigma^3 \mathrm{erf}\Bigl(\frac{v}{\sqrt{2}\sigma}\Bigr) - \sigma^2 v \exp\Bigl(-\frac{v^2}{2\sigma^2}\Bigr) + C \\ I^{(3)}(v) = \int v^3 \exp\Bigl(-\frac{v^2}{2\sigma^2}\Bigr) dv &=& -\sigma^2(2\sigma^2+v^2)\exp\Bigl(-\frac{v^2}{2\sigma^2}\Bigr) + C \\ I^{(4)}(v) = \int v^4 \exp\Bigl(-\frac{v^2}{2\sigma^2}\Bigr) dv &=& 3\sqrt\frac{\pi}{2} \sigma^5 \mathrm{erf}\Bigl(\frac{v}{\sqrt{2}\sigma}\Bigr) - \sigma^2 v (3\sigma^2+v^2) \exp\Bigl(-\frac{v^2}{2\sigma^2}\Bigr) + C ~~~~~~ \textrm{To be checked}\\\end{split}\]

Further expression:

\[I^{(n+2)}(v) = \sigma^2(n+1)\cdot I^{(n)}(v) - v^{n+1}\sigma^2\exp\Bigl(-\frac{v^2}{2\sigma^2}\Bigr)\]

or

\[\int v^{n+2} \exp\Bigl(-\frac{v^2}{2\sigma^2}\Bigr) dv = \sigma^2(n+1)\int v^n \exp\Bigl(-\frac{v^2}{2\sigma^2}\Bigr) dv - v^{n+1}\sigma^2\exp\Bigl(-\frac{v^2}{2\sigma^2}\Bigr)\]

Here \(\mathrm{erf}\) is the error function, which can be defined as

\[\mathrm{erf}(x)=\frac{2}{\sqrt\pi}\int_0^x \exp(-t^2) dt\]

By definition, \(\mathrm{erf}(\pm\infty)=\pm 1\), and \(\mathrm{erf}(0)=0\). In python, one can get it by either math.erf() or scipy.special.erf() function.

Sometimes, one defines \(\mathrm{erfc}(x) \equiv 1 - \mathrm{erf}(x)\), also one can use math.erfc() or scipy.special.erfc().

Definite integration:

Integral

\(v \in [-\infty, \infty]\)

\(v\in [0, \infty]\)

\(I^{(0)}\)

\(\sqrt{2\pi}\sigma\)

\(\sqrt\frac{\pi}{2}\sigma\)

\(I^{(1)}\)

\(0\)

\(\sigma^2\)

\(I^{(2)}\)

\(\sqrt{2\pi}\sigma^3\)

\(\sqrt\frac{\pi}{2}\sigma^3\)

\(I^{(3)}\)

\(0\)

\(2\sigma^4\)

\(I^{(4)}\) (TBC)

\(3\sqrt{2\pi}\sigma^5\)

\(3\sqrt\frac{\pi}{2}\sigma^5\)
irfpy.util.maxwell.mkfunc(n0, v0, vth, multiarray=False)[source]

Returns a function of the Maxwellian.

Parameters:

  • n0 (Float) – Density (scalar)

  • v0 (array-like) – 3-D array (list, tuple, numpy.array, or whatever) of initial velocity, (vx0, vy0, vz0)

  • vth (Float) – Thermal velocity (scalar)

Returns:

A Maxwell distribution function, which get 3-element array (list, tuple, numpy.array) to return the distribution function.

\[f([v_x,v_y,v_z]) = \frac{n_0}{(2\pi v_{th}^2)^{3/2}} \exp -\frac{(v_x-v_{x0})^2 + (v_y-v_{y0})^2 + (v_z-v_{z0})^2 }{2 v_{th}^2}\]

Return type:

A function.

For n=5 /cc, 400 km/s bulk velocity with 40 km/s thermal velocity maxwellian is formed from the following. Use MKSA system.

>>> fsw = mkfunc(5e6, [-400e3, 0, 0], 40e3)
>>> print('%.2e' % fsw([-400e3, 0, 0]))    # The peak
4.96e-09

You may give (3, N) shaped array or (3, N, M)… shaped array.

>>> fval = fsw([[-400e3, -400e3], [0, -40e3], [0, 0]])
>>> from irfpy.util.with_context import printoptions
>>> with printoptions(precision=3):
...     print(fval)       
[4.960e-09   3.009e-09]

Note again that the expression of the Maxwellian does not include mass. Thus, this function can be used regardless the particles.

class irfpy.util.maxwell.MaxwellDistributionTT(n0, v0, temperature_tensor, tol=1e-08, m=1.67262192369e-27)[source]

Bases: object

Create the function for the Maxwell distribution with temperature tensor as argument.

Parameters:

  • n0 – Density in MKSA (/m^3)

  • v0 – Bulk velocity vector in MKSA (m/s). (3,) shape numpy array.

  • temperature_tensor – Temperature tensor in Kelvin. (3, 3) shape array, but must be symmetric. If it is not symmetric, NotSymmetricException is raised.

  • tol – Tolerance if the temperature tensor is symmetric or not. If the tensor is judged as non-symmetric, the NotSymmetricException is raised

Returns:

A function, fv(vx, vy, vz). All the units are assumed MKSA.

\[f(\vec{v}) = n_0 \left(\frac{m}{2\pi k}\right)^{3/2}\frac{1}{\sqrt{|T|}}\exp-\frac{m}{2k}(\vec{v}-\vec{v_0})^T\mathbf{T}^{-1}(\vec{v}-\vec{v_0}))\]
>>> f = MaxwellDistributionTT(1e6, [300e3, 0e3, 0e3], [[15000, 0, 0], [0, 5000, 0,], [0, 0, 30000]])
>>> print(f)
f(\vec{v}) = n_0 \left(\frac{m}{2\pi k}\right)^{3/2}\frac{1}{\sqrt{|T|}}\exp-\frac{m}{2k}(\vec{v}-\vec{v_0})^T\mathbf{T}^{-1}(\vec{v}-\vec{v_0}))
n_0 = 1000000.0
v_0 = 300000.0 0.0 0.0
T = 15000 0 0
    0 5000 0
    0 0 30000
m = 1.67262192369e-27

>>> print('{:.2e}'.format(f([300e3, 0, 0])))
5.64e-08
>>> print('{:.2e}'.format(f([350e3, 0, 0])))
2.33e-12
>>> print('{:.2e}'.format(f([300e3, 50e3, 0])))
3.96e-21
>>> print('{:.2e}'.format(f([300e3, 0, 50e3])))
3.63e-10
class irfpy.util.maxwell.LnMaxwellDistributionTT(n0, v0, temperature_tensor, tol=1e-08, m=1.67262192369e-27)[source]

Bases: object

Create the function for the log-natural Maxwell distribution with temperature tensor as argument.

Parameters:

  • n0 – Density in MKSA (/m^3)

  • v0 – Bulk velocity vector in MKSA (m/s). (3,) shape numpy array.

  • temperature_tensor – Temperature tensor in Kelvin. (3, 3) shape array, but must be symmetric. If it is not symmetric, NotSymmetricException is raised.

  • tol – Tolerance if the temperature tensor is symmetric or not. If the tensor is judged as non-symmetric, the NotSymmetricException is raised

Returns:

A function, lnfv(vx, vy, vz); it returns the log of VDF. All the units are assumed MKSA.

\[\begin{split}\ln f(\vec{v}) &= \ln\left(n_0 \left(\frac{m}{2\pi k}\right)^{3/2}\frac{1}{\sqrt{|T|}}\exp-\frac{m}{2\pi k}(\vec{v}-\vec{v_0})^T\mathbf{T}^{-1}(\vec{v}-\vec{v_0}))\right) \\ &= \ln n_0 + \frac{3}{2}(\ln m-\ln 2\pi k) - \frac{1}{2}\ln|T| - \frac{m}{2k}(\vec{v}-\vec{v_b})^T\mathbf{T}(\vec{v}-\vec{v_b})\end{split}\]

Note

This is much faster than MaxwellDistributionTT (~1000 times). However, np.exp costs significant amount time.

>>> f = LnMaxwellDistributionTT(1e6, [300e3, 0e3, 0e3], [[15000, 0, 0], [0, 5000, 0,], [0, 0, 30000]])
>>> print(f)
\ln f(\vec{v}) = \ln\left(n_0 \left(\frac{m}{2\pi k}\right)^{3/2}\frac{1}{\sqrt{|T|}}\exp-\frac{m}{2k}(\vec{v}-\vec{v_0})^T\mathbf{T}^{-1}(\vec{v}-\vec{v_0}))\right)
n_0 = 1000000.0
v_0 = 300000.0 0.0 0.0
T = 15000 0 0
    0 5000 0
    0 0 30000
m = 1.67262192369e-27

>>> import numpy as np
>>> print('{:.2e}'.format(np.exp(f([300e3, 0, 0]))))
5.64e-08
>>> print('{:.2e}'.format(np.exp(f([350e3, 0, 0]))))
2.33e-12
>>> print('{:.2e}'.format(np.exp(f([300e3, 50e3, 0]))))
3.96e-21
>>> print('{:.2e}'.format(np.exp(f([300e3, 0, 50e3]))))
3.63e-10
exception irfpy.util.maxwell.NotSymmetricException(value)[source]

Bases: IrfpyException

Exception raised when the matrix/tensor is not symmetric.

irfpy.util.maxwell.is_symmetric_tensor(tensor2d, tol=1e-08, raises=False)[source]

Check if the given tensor is symmetric or not.

Parameters:

  • tensor2d – 2-D tensor to be checked.

  • tol – Tolerance

  • raises – Instead of returning False, an exception (NotSymmetricException) is raised.

Returns:

If all the elements in tensor2d is smaller than tol, the tensor2d is considered as symmetric and return True. Otherwise, either False is returned, or NotSymmetricException is raised, depending on the raises keyword.

>>> t = [[1, 3, 2],
...      [3, 1, 4],
...      [2, 4, 3.0]]
>>> is_symmetric_tensor(t)
True

>>> t = [[1, 3, 2],
...      [3, 4, 5],
...      [1.999, 5, 3]]
>>> is_symmetric_tensor(t)
False

>>> is_symmetric_tensor(t, raises=True)    # doctests: +IGNORE_EXCEPTION_DETAIL
Traceback (most recent call last):
irfpy.util.maxwell.NotSymmetricException: 'Given tensor/matrix is not symmetric: [[1.    3.    2.   ]\n [3.    4.    5.   ]\n [1.999 5.    3.   ]]'

>>> is_symmetric_tensor(t, tol=0.01)
True
irfpy.util.maxwell.mkfunc_u(n0_u, v0_u, vth_u)[source]

Unitful version of mkfunc().

>>> f_u = mkfunc_u(30 / u.cc, [250, -15, 0] * u.km / u.s, 35 * u.km / u.s)
>>> print('%.3e' % f_u([230, 0, 0] * u.km / u.s).nr(u.nit('s**3/cm**6')))
3.442e-20

If you use non-unitful version, the results should be the same.

>>> f = mkfunc(30e6, [250e3, -15e3, 0], 35e3)
>>> print('%.3e' % f([230e3, 0, 0]))
3.442e-08
class irfpy.util.maxwell.MaxwellDistributionG(m=UnitConstant('proton_mass', 1.672621637e-27 * kg, 'm_p'), n=array(10.) * 1 / cc, vsw=array([400., 50., -30.]) * km / s, vth=array(60.) * km / s)[source]

Bases: object

Maxwell Distribution in Gridded space

This class derives the differential flux and the velocity distribution function, and save these values on a gridded bin.

Note

This class / functionality has been replaced by irfpy.util.energyspectrum.maxwellDifferentialFluxGrid() and :func:` irfpy.util.vdf.maxwellVelocityDistributionGrid` functions.

First, create object with a set of parameter.

>>> maxwellobj = MaxwellDistributionG(m=k.proton_mass, n=10 / u.cc,
...             vsw = [-400, 0, 50] * u.km / u.s,
...             vth = 150 * u.km / u.s)

As seen in the sample, parameters are unitted. From profiler point of view, the overhead in the computational resource (time) is very small by using unitq.

Second, set the grid if you want. Only E-t-p contour grid is supported. By default, grid is defined as 96 step for energy direction (100 - 20000 eV, logspace) 181 steps for theta, and 361 steps for phi.

>>> maxwellobj.set_etp_grid([100, 300, 1000, 3000, 10000],
...                 [-90, -60, -30, -15, 0, 15, 30, 60, 90],
...                 [-180, -120, -60, 0, 60, 120, 180])

Third, calculate.

>>> maxwellobj.do_calc()

Now you can refer to j and fval.

>>> print('%.3e (1/cm**2 /sr /eV / s)' % maxwellobj.j.max())
3.164e+05 (1/cm**2 /sr /eV / s)

>>> u.pprint(maxwellobj.fval.max(), fmt='%.3e')
1.724e-10 s**3/m**6

Other attributes that you may refer to is

  • m_u: Mass specified

  • n_u: Density specified

  • vsw_u: Bulk velocity specified

  • vth_u: Thermal velocity specified.

  • dive: Number of energy bin

  • divt: Number of theta (polar) bin

  • divp: Number of phi (azimuth) bin

  • j: Differential flux

  • fval: VDF

Initialize by the solar wind parameter

logger = <Logger irfpy.util.maxwell.MaxwellDistributionG (DEBUG)>
set_etp_grid(egrid1d=array([100., 105.7356327, 111.80024023, 118.21269138, 124.99293716, 132.16207294, 139.74240402, 147.75751505, 156.2323434, 165.19325679, 174.66813525, 184.68645794, 195.27939482, 206.47990365, 218.32283253, 230.84502831, 244.08545125, 258.08529622, 272.88812088, 288.53998118, 305.08957471, 322.58839213, 341.09087745, 360.65459736, 381.3404204, 403.21270627, 426.33950611, 450.79277426, 476.64859204, 503.98740457, 532.89427097, 563.45912906, 595.77707514, 629.9486599, 666.08020125, 704.28411511, 744.67926515, 787.39133262, 832.5532074, 880.30540144, 930.79648594, 984.1835536, 1040.63270737, 1100.31957726, 1163.42986678, 1230.15993071, 1300.717386, 1375.32175778, 1454.20516231, 1537.61302918, 1625.80486494, 1719.05506048, 1817.65374473, 1921.90768735, 2032.14125321, 2148.69741152, 2271.93880297, 2402.24886796, 2540.03303966, 2685.72000538, 2839.76304035, 3002.641418, 3174.86190116, 3356.96031867, 3549.50323256, 3753.0897008, 3968.35314109, 4195.96330166, 4436.62834504, 4691.09705135, 4960.161148, 5244.65777298, 5545.47207942, 5863.53998959, 6199.85110685, 6555.45179453, 6931.44843155, 7329.01085465, 7749.37599811, 8193.85174221, 8663.82098245, 9160.74593214, 9686.17267175, 10241.73595928, 10829.16431641, 11450.28540651, 12107.03172099, 12801.44659186, 13535.69054917, 14312.04804301, 15132.93455118, 16000.90409437, 16918.65718254, 17889.04921698, 18915.0993743, 20000.]), tgrid1d=array([-90., -89., -88., -87., -86., -85., -84., -83., -82., -81., -80., -79., -78., -77., -76., -75., -74., -73., -72., -71., -70., -69., -68., -67., -66., -65., -64., -63., -62., -61., -60., -59., -58., -57., -56., -55., -54., -53., -52., -51., -50., -49., -48., -47., -46., -45., -44., -43., -42., -41., -40., -39., -38., -37., -36., -35., -34., -33., -32., -31., -30., -29., -28., -27., -26., -25., -24., -23., -22., -21., -20., -19., -18., -17., -16., -15., -14., -13., -12., -11., -10., -9., -8., -7., -6., -5., -4., -3., -2., -1., 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27., 28., 29., 30., 31., 32., 33., 34., 35., 36., 37., 38., 39., 40., 41., 42., 43., 44., 45., 46., 47., 48., 49., 50., 51., 52., 53., 54., 55., 56., 57., 58., 59., 60., 61., 62., 63., 64., 65., 66., 67., 68., 69., 70., 71., 72., 73., 74., 75., 76., 77., 78., 79., 80., 81., 82., 83., 84., 85., 86., 87., 88., 89., 90.]), pgrid1d=array([-180., -179., -178., -177., -176., -175., -174., -173., -172., -171., -170., -169., -168., -167., -166., -165., -164., -163., -162., -161., -160., -159., -158., -157., -156., -155., -154., -153., -152., -151., -150., -149., -148., -147., -146., -145., -144., -143., -142., -141., -140., -139., -138., -137., -136., -135., -134., -133., -132., -131., -130., -129., -128., -127., -126., -125., -124., -123., -122., -121., -120., -119., -118., -117., -116., -115., -114., -113., -112., -111., -110., -109., -108., -107., -106., -105., -104., -103., -102., -101., -100., -99., -98., -97., -96., -95., -94., -93., -92., -91., -90., -89., -88., -87., -86., -85., -84., -83., -82., -81., -80., -79., -78., -77., -76., -75., -74., -73., -72., -71., -70., -69., -68., -67., -66., -65., -64., -63., -62., -61., -60., -59., -58., -57., -56., -55., -54., -53., -52., -51., -50., -49., -48., -47., -46., -45., -44., -43., -42., -41., -40., -39., -38., -37., -36., -35., -34., -33., -32., -31., -30., -29., -28., -27., -26., -25., -24., -23., -22., -21., -20., -19., -18., -17., -16., -15., -14., -13., -12., -11., -10., -9., -8., -7., -6., -5., -4., -3., -2., -1., 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27., 28., 29., 30., 31., 32., 33., 34., 35., 36., 37., 38., 39., 40., 41., 42., 43., 44., 45., 46., 47., 48., 49., 50., 51., 52., 53., 54., 55., 56., 57., 58., 59., 60., 61., 62., 63., 64., 65., 66., 67., 68., 69., 70., 71., 72., 73., 74., 75., 76., 77., 78., 79., 80., 81., 82., 83., 84., 85., 86., 87., 88., 89., 90., 91., 92., 93., 94., 95., 96., 97., 98., 99., 100., 101., 102., 103., 104., 105., 106., 107., 108., 109., 110., 111., 112., 113., 114., 115., 116., 117., 118., 119., 120., 121., 122., 123., 124., 125., 126., 127., 128., 129., 130., 131., 132., 133., 134., 135., 136., 137., 138., 139., 140., 141., 142., 143., 144., 145., 146., 147., 148., 149., 150., 151., 152., 153., 154., 155., 156., 157., 158., 159., 160., 161., 162., 163., 164., 165., 166., 167., 168., 169., 170., 171., 172., 173., 174., 175., 176., 177., 178., 179., 180.]))[source]

Define the grid. Not the energy in eV, angles in degree.

get_vtp_grid3d()[source]

Calculate the 3-D grid and return

Returns:

(vgrid, tgrid, pgrid). Velocity, theta and phi (degress) on the grid. Each is 3-D array with shape of (dive, divt, divp).

get_egrid3d()[source]

Calculate 3-D grid of energy

Returns:

Energy on the grid in eV. Shape of (dive, divt, divp)

get_vxyz_array()[source]

Calculate corresponding velocity on the grid

Returns:

(vxarr, vyarr, vzarr). Each data is 3-D array with shape of (dive, divt, divp).

do_calc()[source]

Run the calculation

irfpy.util.maxwell.mk1d(A, B, C)[source]

Returns a callable of 1-D expression of Maxwellian distribution.

Returns a callable of 1-D expression of the Maxwellian distribution (standard distribution). The formulation is

\[f(v) = A \exp -\frac{(v-B)^2}{2 C^2}\]

See mkfunc() also.

irfpy.util.maxwell.mkslice(n0, v0, vth)[source]

Returns a callable of 1-D slice of the distribution function.

Maxwell distribution \(f(v_x ,v_y ,v_z)\) is 3-D, but one can slice it along the direction of the bulk velocity \(\mathbf{v_0} = (v_{x0} ,v_{y0} ,v_{z0})\). This method can be widely used for plasma data analysis because one can pick energy spectra of the direction along the bulk velocity.

The formulation is

\[f(v) = n_0 ( 2\pi v_{th}^2 ) ^{-3/2} \exp -\frac{(v-v_0)^2}{2 v_{th}^2}\]

Here, \(v\) is the scalar value, which indeed expresses the vector \(\mathbf{v} = v \mathbf{v_0}/|\mathbf{v_0}|\). Never consider such like \(v=|\mathbf{v}|\).

Parameters:

  • n0 – Density (scalar)

  • v0 – Velocity (scalar)

  • vth – Thermal velocity (scalar)

Retturns:

A Maxwell distribution function, which get scalar velocity to return the distribution function.

irfpy.util.maxwell.fit1d(x_array, y_array)[source]

Fitting by 1-dimensional Maxwell distiribution function.

Parameters:

  • x_array (A sequence of flaot) – A sequence of x values.

  • y_array (A sequence of flaot) – A sequence of y values.

Returns:

((A, B, C), success_code). (A, B, C) is defined below. success_code is the same as defied in optimize.leastsq().

Fitting by least square algorithm. Initial values for iteration is estimated automatically. In order to reduce confusion in the sqrt term in the coefficient, the coefficient it iself is one parameter, A. So users must normalize the coefficiency to convert it to density or probablility functions.

The returned values is the best fit parameters of (A,B,C) which is formulated as

\[y_{array} = A \exp(-\frac{(x-B)^2}{2C^2})\]
irfpy.util.maxwell.getRandom(vel0, vth, n=10000)[source]

Get a random distributed numpy.array of the 3-D Maxwellian.

vel0 should be 3-d array of the bulk velocity vector. vth should be the scalar. The returned value is a numpy.array with a shape of [3, n]

irfpy.util.maxwell.getRandom1D(v0, vth, n=10000)[source]

Get a random distributed numpy.array of values for 1-D Maxwellian.

Get a random distributed numpy.array of the value following 1-D Maxwellian using scipy. Indeed this function is just a wrapper to scipy.norm functionality. The thermal velocity is same as the scale parameter of the normal distribution, and the bulk velocity is same as the location parameter.

irfpy.util.maxwell.flux_integrated(n0, v0, vth, vzmin=-inf)[source]

Integrate the maxwellian, with a specific minimum range.

Parameters:

  • n0 – Density, scalar.

  • v0 – Velocity, vector. (vx0, vy0, vz0)

  • vth – Thermal velocity. Scalar.

Keyword:

vzmin: If given, the lower end of integral (in z axis) is given.

From the Maxwellian velocity distribution,

\[f(v_x,v_y,v_z) = n_0 ( 2\pi v_\mathrm{th}^2 ) ^{-\frac{3}{2}} \exp -\frac{ (v_x-v_{x0})^2 + (v_y-v_{y0})^2 + (v_z-v_{z0})^2 } {2 v_\mathrm{th}^2}\]

the following integral is returned.

\[F = \int_{-\infty}^{\infty} dv_x \int_{-\infty}^{\infty} dv_y \int_{v_{zmin}}^{\infty} dv_z \bigl[v_z f(v_x, v_y, v_z)\bigr]\]

A mathematics says that the value is

\[F = \frac{n_0v_{th}}{\sqrt{2\pi}}\exp-\frac{(v_{zmin}-v_{z0})^2}{2v_{th}^2} + \frac{n_0v_{z0}}{2}\mathrm{erfc}\frac{v_{zmin}-v_{z0}}{\sqrt{2}{v_{th}}}\]

This integral is useful for calculating the flux but only for those with the given speed.

Example 1: For the default parameter, the returned value gives just n0 v0z, since the integral is taken over the whole velocity space.

>>> n0 = 30
>>> v0 = (0, 10, 20)
>>> vth = 15
>>> print(flux_integrated(n0, v0, vth))
600.0

Example 2: If the lower edge of the vz integral range is higher than the core component, the flux approaches to 0.

>>> n0 = 30
>>> v0 = (0, 15, 40)
>>> vth = 20
>>> vzmin = 140    # This is vz0 + 5 vth.
>>> f = flux_integrated(n0, v0, vth, vzmin=vzmin)
>>> print('{:.2e}'.format(f))
1.24e-03

Example 3: If no bulk velocity, and integral range from 0, then

>>> n0 = 1
>>> v0 = (0, 0, 0)
>>> vth = 5
>>> vzmin = 0
>>> f = flux_integrated(n0, v0, vth, vzmin=vzmin)
>>> print('{:.4f}'.format(f))
1.9947
irfpy.util.maxwell.flux_precip(n0, vz0, vth)[source]

Return the precipitating flux along +z.

This function is to calculate the precipitation flux with a non-zero temperature.

This is equivalent to flux_integrated(n0, (0, 0, vz0), vth, vzmin=0).

Todo

More documentation

Parameters:

  • n0 – Density

  • vz0 – Velocity along the precipitating direction. Usually, the direction of nadir or direction of the magnetic field according to the problem of interest.

  • vth – Thermal velocity.

The calculation is based on the MKSA system. If you need non-MKSA system calculation, unit conversions are necessary.

>>> print('{:.4f}'.format(flux_precip(1, 0, 5)))
1.9947
irfpy.util.maxwell.differential_flux_u(density, bulk_velocity, thermal_velocity, m=1.67262192369e-27)[source]
class irfpy.util.maxwell.LnDifferentialFlux(density_mksa, bulk_velocity_mksa, thermal_velocity_mksa, m=1.67262192369e-27)[source]

Bases: object

Create a function to return the log natural of differential flux.

Parameters:

  • density_mksa – Density in m^-3

  • bulk_velocity_mksa – Bulk velocity vector. m/s

  • thermal_velocity_mksa – Scalar thermal velocity. m/s

  • m – Mass. kg.

Return a function (a callable object) that calculate the log differential flux for a maxwellian.

If you want to calculate the differential flux corresponding to the density 10/cm^3, Vx=150 km/s, Vy=Vz=1 m/s, and Vth=30 km/s (~10 eV for proton)

>>> dflux_func = LnDifferentialFlux(10e6, [150e3, 1, 1], 30e3)

The returned dflux_func is a function(-like) that can give you the differential flux. The unit for the dflux_func is

  • eV for energy (in/out)

  • degrees for angles (in)

    • polar angle, \(\theta\), from +z

    • azimutha angle, \(\phi\), from +x

  • cm for dimension (out)

  • steradian for solid angle (out)

  • /cm2 s sr eV for the differential flux (out)

Thus, if you want to know the differential flux for 120 eV from x-direction (\(\theta\) =90, \(\phi\) =0).

>>> lnj = dflux_func(120, 90, 0)    # 120 eV, theta=90, phi=0
>>> import numpy as np
>>> print('{:.2e}'.format(np.exp(lnj)))          # The returned flux is in /cm2 s sr eV
5.17e+06
e = 1.602176634e-19
log_e = -43.27775366589175
log_2 = 0.6931471805599453
log_2pi = 1.8378770664093453
log_1em4 = -9.210340371976182
class irfpy.util.maxwell.DifferentialFlux(density_mksa, bulk_velocity_mksa, thermal_velocity_mksa, m=1.67262192369e-27)[source]

Bases: object

Return a function that calculate the differential flux.

Parameters:

  • density_mksa – Density in m^-3

  • bulk_velocity_mksa – Bulk velocity vector. m/s

  • thermal_velocity_mksa – Scalar thermal velocity. m/s

  • m – Mass. kg

Returns:

The function to calculate the differentil flux, J.

J is the function to return the differential flux. Note that the unit of the energy is eV and that of the area is cm^2.

  • eV for energy (in/out)

  • degrees for angles (in)

    • polar angle, \(\theta\), from +z

    • azimutha angle, \(\phi\), from +x

  • cm for dimension (out)

  • steradian for solid angle (out)

  • /cm2 s sr eV for the differential flux (out)

>>> dflux_func = DifferentialFlux(10e6, [150e3, 1, 1], 30e3)

Thus, if you want to know the differential flux for 120 eV from x-direction (\(\theta\) =90, \(\phi\) =0).

>>> j = dflux_func(120, 90, 0)    # 120 eV, theta=90, phi=0
>>> print('{:.2e}'.format(j))          # The returned flux is in /cm2 s sr eV
5.17e+06
class irfpy.util.maxwell.DifferentialFlux_u(density, bulk_velocity, thermal_velocity, m=UnitConstant('proton_mass', 1.672621637e-27 * kg, 'm_p'))[source]

Bases: object

Return a function that calculate the differential flux, with unit support.

Parameters:

  • density –

  • bulk_velocity – Bulk velocity vector.

  • thermal_velocity – Scalar thermal velocity.

  • m – Mass

Returns:

The function to calculate the differentil flux, J. Unit support.

J is the function to return the differential flux.

>>> dflux_func = DifferentialFlux_u(10 / u.cc, [150 * u.km/u.s, 0.001 * u.km/u.s, 0.001 * u.km/u.s], 30 * u.km/u.s)

Thus, if you want to know the differential flux for 120 eV from x-direction (\(\theta\) =90, \(\phi\) =0).

>>> j = dflux_func(120 * u.eV, 90 * u.deg, 0 * u.deg)    # 120 eV, theta=90, phi=0
>>> print('{:.2e}'.format(j))          # The returned flux is in /cm2 s sr eV
5.17e+06 1/(cm**2*s*sr*eV)