"""
Common distribution functions for plasmas, such as the Maxwellian or
Kappa distributions. Functionality is intended to include generation,
fitting and calculation.
"""
__all__ = [
"Maxwellian_1D",
"Maxwellian_velocity_2D",
"Maxwellian_velocity_3D",
"Maxwellian_speed_1D",
"Maxwellian_speed_2D",
"Maxwellian_speed_3D",
"kappa_velocity_1D",
"kappa_velocity_3D",
]
import astropy.units as u
import numpy as np
from scipy.special import gamma
from plasmapy.formulary.speeds import kappa_thermal_speed, thermal_speed
from plasmapy.particles import ParticleLike, particle_input
from plasmapy.utils._units_definitions import (
SPEED_DISTRIBUTION_UNITS_1D,
SPEED_DISTRIBUTION_UNITS_2D,
SPEED_DISTRIBUTION_UNITS_3D,
SPEED_UNITS,
)
def _v_drift_conversion(v_drift):
# Helper method to assign equivalent value in SPEED_UNITS and/or remove units
if isinstance(v_drift, u.Quantity):
v_drift = v_drift.to_value(SPEED_UNITS)
return v_drift
[docs]
@particle_input
def Maxwellian_1D(
v,
T,
particle: ParticleLike = "e-",
v_drift=0,
vTh=np.nan,
units="units",
*,
mass_numb=None,
Z=None,
):
r"""
Probability distribution function of velocity for a Maxwellian
distribution in 1D.
Returns the probability density function at the velocity ``v`` in m/s
to find a particle ``particle`` in a plasma of temperature ``T``
following the Maxwellian distribution function.
Parameters
----------
v : `~astropy.units.Quantity`
The velocity in units convertible to m/s.
T : `~astropy.units.Quantity`
The temperature in kelvin.
particle : `str`, optional
Representation of the particle species(e.g., ``'p+'`` for protons,
``'D+'`` for deuterium, or ``'He-4 +1'`` for singly ionized
helium-4), which defaults to electrons.
v_drift : `~astropy.units.Quantity`, optional
The drift velocity in units convertible to m/s.
vTh : `~astropy.units.Quantity`, optional
Thermal velocity (most probable velocity) in m/s. This is used for
optimization purposes to avoid re-calculating ``vTh``, for example
when integrating over velocity-space.
units : `str`, optional
Selects whether to run function with units and unit checks (when
equal to "units") or to run as unitless (when equal to "unitless").
The unitless version is substantially faster for intensive
computations.
mass_numb : integer, |keyword-only|, optional
The mass number associated with ``particle``.
Z : real number, |keyword-only|, optional
The charge number associated with ``particle``.
Returns
-------
f : `~astropy.units.Quantity`
Probability density in units of velocity\ :sup:`-1`\ , normalized so that
:math:`\int_{-∞}^{+∞} f(v) dv = 1`.
Raises
------
`TypeError`
The parameter arguments are not Quantities and
cannot be converted into Quantities.
`~astropy.units.UnitConversionError`
If the parameters are not in appropriate units.
`ValueError`
If the temperature is negative, or the particle mass or charge state
cannot be found.
Notes
-----
In one dimension, the Maxwellian distribution function for a particle of
mass m, velocity v, a drift velocity V and with temperature T is:
.. math::
f = \sqrt{\frac{m}{2 \pi k_B T}} e^{-\frac{m}{2 k_B T} (v-V)^2}
\equiv \frac{1}{\sqrt{\pi v_{Th}^2}} e^{-(v - v_{drift})^2 / v_{Th}^2}
where :math:`v_{Th} = \sqrt{2 k_B T / m}` is the thermal speed
Examples
--------
>>> import astropy.units as u
>>> v = 1 * u.m / u.s
>>> Maxwellian_1D(v=v, T=30000 * u.K, particle="e-", v_drift=0 * u.m / u.s)
<Quantity 5.9163...e-07 s / m>
"""
if units == "units":
# unit checks and conversions
# checking velocity units
v = v.to_value(SPEED_UNITS)
# Catching case where drift velocities have default values,
v_drift = _v_drift_conversion(v_drift)
# convert temperature to kelvin
T = T.to_value(u.K, equivalencies=u.temperature_energy())
if not np.isnan(vTh):
# check units of thermal velocity
vTh = vTh.to_value(SPEED_UNITS)
if np.isnan(vTh):
# get thermal speed
vTh = thermal_speed(
T << u.K, particle=particle, method="most_probable"
).to_value(SPEED_UNITS)
# Get thermal velocity squared
vThSq = vTh**2
# Get square of relative particle velocity
vSq = (v - v_drift) ** 2
# calculating distribution function
coeff = (vThSq * np.pi) ** (-1 / 2)
expTerm = np.exp(-vSq / vThSq)
distFunc = coeff * expTerm
if units == "units":
return distFunc << SPEED_DISTRIBUTION_UNITS_1D
elif units == "unitless":
return distFunc
else:
raise ValueError(f"Units must be either 'units' or 'unitless', got {units}).")
[docs]
@particle_input
def Maxwellian_velocity_2D(
vx,
vy,
T,
particle: ParticleLike = "e-",
vx_drift=0,
vy_drift=0,
vTh=np.nan,
units="units",
*,
mass_numb=None,
Z=None,
):
r"""
Probability distribution function of velocity for a Maxwellian
distribution in 2D.
Return the probability density function for finding a particle with
velocity components ``vx`` and ``vy`` in m/s in an equilibrium plasma of
temperature ``T`` which follows the 2D Maxwellian distribution function.
This function assumes Cartesian coordinates.
Parameters
----------
vx : `~astropy.units.Quantity`
The velocity in x-direction units convertible to m/s.
vy : `~astropy.units.Quantity`
The velocity in y-direction units convertible to m/s.
T : `~astropy.units.Quantity`
The temperature, preferably in kelvin.
particle : `str`, optional
Representation of the particle species [e.g., ``'p+'`` for protons,
``'D+'`` for deuterium, or ``'He-4 +1'`` for :math:`He_4^{+1}`
(singly ionized helium-4)], which defaults to electrons.
vx_drift : `~astropy.units.Quantity`, optional
The drift velocity in x-direction in units convertible to m/s.
vy_drift : `~astropy.units.Quantity`, optional
The drift velocity in y-direction in units convertible to m/s.
vTh : `~astropy.units.Quantity`, optional
Thermal velocity (most probable) in m/s. This is used for
optimization purposes to avoid re-calculating ``vTh``, for example
when integrating over velocity-space.
units : `str`, optional
Selects whether to run function with units and unit checks (when
equal to "units") or to run as unitless (when equal to "unitless").
The unitless version is substantially faster for intensive
computations.
mass_numb : integer, |keyword-only|, optional
The mass number associated with ``particle``.
Z : real number, |keyword-only|, optional
The charge number associated with ``particle``.
Returns
-------
f : `~astropy.units.Quantity`
Probability density in Velocity\ :sup:`-1`\ , normalized so that
:math:`\iiint_{0}^∞ f(\vec{v}) d\vec{v} = 1`.
Raises
------
TypeError
A parameter argument is not a `~astropy.units.Quantity` and
cannot be converted into a `~astropy.units.Quantity`.
~astropy.units.UnitConversionError
If the parameters is not in appropriate units.
ValueError
If the temperature is negative, or the particle mass or charge state
cannot be found.
Notes
-----
In 2D, the Maxwellian velocity distribution function describing
the distribution of particles with speed :math:`v` in a plasma with
temperature :math:`T` is given by:
.. math::
f = (\pi v_{Th}^2)^{-1} \exp \left [-(\vec{v} -
\vec{V}_{drift})^2 / v_{Th}^2 \right ]
where :math:`v_{Th} = \sqrt{2 k_B T / m}` is the thermal speed.
See Also
--------
Maxwellian_1D
Examples
--------
>>> import astropy.units as u
>>> v = 1 * u.m / u.s
>>> Maxwellian_velocity_2D(
... vx=v,
... vy=v,
... T=30000 * u.K,
... particle="e-",
... vx_drift=0 * u.m / u.s,
... vy_drift=0 * u.m / u.s,
... )
<Quantity 3.5002...e-13 s2 / m2>
"""
if units == "units":
# unit checks and conversions
# checking velocity units
vx = vx.to_value(SPEED_UNITS)
vy = vy.to_value(SPEED_UNITS)
# catching case where drift velocities have default values, they
# need to be assigned units
vx_drift = _v_drift_conversion(vx_drift)
vy_drift = _v_drift_conversion(vy_drift)
# convert temperature to kelvin
T = T.to_value(u.K, equivalencies=u.temperature_energy())
if not np.isnan(vTh):
# check units of thermal velocity
vTh = vTh.to_value(SPEED_UNITS)
if np.isnan(vTh):
# get thermal speed
vTh = thermal_speed(
T << u.K, particle=particle, method="most_probable"
).to_value(SPEED_UNITS)
# accounting for thermal velocity in 2D
vThSq = vTh**2
# Get square of relative particle velocity
vSq = (vx - vx_drift) ** 2 + (vy - vy_drift) ** 2
# calculating distribution function
coeff = (vThSq * np.pi) ** (-1)
expTerm = np.exp(-vSq / vThSq)
distFunc = coeff * expTerm
if units == "units":
return distFunc << SPEED_DISTRIBUTION_UNITS_2D
elif units == "unitless":
return distFunc
else:
raise ValueError(f"Units must be either 'units' or 'unitless', got {units}).")
[docs]
@particle_input
def Maxwellian_velocity_3D(
vx,
vy,
vz,
T,
particle: ParticleLike = "e-",
vx_drift=0,
vy_drift=0,
vz_drift=0,
vTh=np.nan,
units="units",
*,
mass_numb=None,
Z=None,
):
r"""
Probability distribution function of velocity for a Maxwellian
distribution in 3D.
Return the probability density function for finding a particle with
velocity components ``vx``, ``vy``, and ``vz`` in m/s in an equilibrium
plasma of temperature ``T`` which follows the 3D Maxwellian distribution
function. This function assumes Cartesian coordinates.
Parameters
----------
vx : `~astropy.units.Quantity`
The velocity in x-direction in units convertible to m/s.
vy : `~astropy.units.Quantity`
The velocity in y-direction units convertible to m/s.
vz : `~astropy.units.Quantity`
The velocity in z-direction units convertible to m/s.
T : `~astropy.units.Quantity`
The temperature, preferably in kelvin.
particle : `str`, optional
Representation of the particle species (e.g., ``'p+'`` for protons,
``'D+'`` for deuterium, or ``'He-4 +1'`` for
singly ionized helium-4), which defaults to electrons.
vx_drift : `~astropy.units.Quantity`, optional
The drift velocity in x-direction units convertible to m/s.
vy_drift : `~astropy.units.Quantity`, optional
The drift velocity in y-direction units convertible to m/s.
vz_drift : `~astropy.units.Quantity`, optional
The drift velocity in z-direction units convertible to m/s.
vTh : `~astropy.units.Quantity`, optional
Thermal velocity (most probable) in m/s. This is used for
optimization purposes to avoid re-calculating ``vTh``, for example
when integrating over velocity-space.
units : `str`, optional
Selects whether to run function with units and unit checks (when
equal to "units") or to run as unitless (when equal to "unitless").
The unitless version is substantially faster for intensive
computations.
mass_numb : integer, |keyword-only|, optional
The mass number associated with ``particle``.
Z : real number, |keyword-only|, optional
The charge number associated with ``particle``.
Returns
-------
f : `~astropy.units.Quantity`
Probability density in Velocity^-1, normalized so that
:math:`\iiint_{0}^∞ f(\vec{v}) d\vec{v} = 1`.
Raises
------
`TypeError`
A parameter argument is not a `~astropy.units.Quantity` and
cannot be converted into a `~astropy.units.Quantity`.
`~astropy.units.UnitConversionError`
If the parameters is not in appropriate units.
`ValueError`
If the temperature is negative, or the particle mass or charge state
cannot be found.
Notes
-----
In 3D, the Maxwellian speed distribution function describing
the distribution of particles with speed :math:`v` in a plasma with
temperature :math:`T` is given by:
.. math::
f = (\pi v_{Th}^2)^{-3/2} \exp \left [-(\vec{v} -
\vec{V}_{drift})^2 / v_{Th}^2 \right ]
where :math:`v_{Th} = \sqrt{2 k_B T / m}` is the thermal speed.
See Also
--------
Maxwellian_1D
Examples
--------
>>> import astropy.units as u
>>> v = 1 * u.m / u.s
>>> Maxwellian_velocity_3D(
... vx=v,
... vy=v,
... vz=v,
... T=30000 * u.K,
... particle="e-",
... vx_drift=0 * u.m / u.s,
... vy_drift=0 * u.m / u.s,
... vz_drift=0 * u.m / u.s,
... )
<Quantity 2.0708...e-19 s3 / m3>
"""
if units == "units":
# unit checks and conversions
# checking velocity units
vx = vx.to_value(SPEED_UNITS)
vy = vy.to_value(SPEED_UNITS)
vz = vz.to_value(SPEED_UNITS)
# catching case where drift velocities have default values, they
# need to be assigned units
vx_drift = _v_drift_conversion(vx_drift)
vy_drift = _v_drift_conversion(vy_drift)
vz_drift = _v_drift_conversion(vz_drift)
# convert temperature to kelvin
T = T.to_value(u.K, equivalencies=u.temperature_energy())
if not np.isnan(vTh):
# check units of thermal velocity
vTh = vTh.to_value(SPEED_UNITS)
if np.isnan(vTh):
# get thermal velocity and thermal velocity squared
vTh = thermal_speed(
T << u.K, particle=particle, method="most_probable"
).to_value(SPEED_UNITS)
# accounting for thermal velocity in 3D
vThSq = vTh**2
# Get square of relative particle velocity
vSq = (vx - vx_drift) ** 2 + (vy - vy_drift) ** 2 + (vz - vz_drift) ** 2
# calculating distribution function
coeff = (vThSq * np.pi) ** (-3 / 2)
expTerm = np.exp(-vSq / vThSq)
distFunc = coeff * expTerm
if units == "units":
return distFunc << SPEED_DISTRIBUTION_UNITS_3D
elif units == "unitless":
return distFunc
else:
raise ValueError(f"Units must be either 'units' or 'unitless', got {units}).")
[docs]
@particle_input
def Maxwellian_speed_1D(
v,
T,
particle: ParticleLike = "e-",
v_drift=0,
vTh=np.nan,
units="units",
*,
mass_numb=None,
Z=None,
):
r"""
Probability distribution function of speed for a Maxwellian distribution
in 1D.
Return the probability density function for finding a particle with
speed ``v`` in m/s in an equilibrium plasma of temperature ``T`` which
follows the Maxwellian distribution function.
Parameters
----------
v : `~astropy.units.Quantity`
The speed in units convertible to m/s.
T : `~astropy.units.Quantity`
The temperature, preferably in kelvin.
particle : `str`, optional
Representation of the particle species [e.g., ``'p+'`` for protons, ``'D+'``
for deuterium, or ``'He-4 +1'`` for :math:`He_4^{+1}`
(singly ionized helium-4)], which defaults to electrons.
v_drift : `~astropy.units.Quantity`
The drift speed in units convertible to m/s.
vTh : `~astropy.units.Quantity`, optional
Thermal velocity (most probable) in m/s. This is used for
optimization purposes to avoid re-calculating ``vTh``, for example
when integrating over velocity-space.
units : `str`, optional
Selects whether to run function with units and unit checks (when
equal to "units") or to run as unitless (when equal to "unitless").
The unitless version is substantially faster for intensive
computations.
mass_numb : integer, |keyword-only|, optional
The mass number associated with ``particle``.
Z : real number, |keyword-only|, optional
The charge number associated with ``particle``.
Returns
-------
f : `~astropy.units.Quantity`
Probability density in speed\ :sup:`-1`\ , normalized so that
:math:`\int_{0}^∞ f(v) dv = 1`.
Raises
------
`TypeError`
The parameter arguments are not Quantities and
cannot be converted into Quantities.
`~astropy.units.UnitConversionError`
If the parameters is not in appropriate units.
`ValueError`
If the temperature is negative, or the particle mass or charge state
cannot be found.
Notes
-----
In one dimension, the Maxwellian speed distribution function describing
the distribution of particles with speed :math:`v` in a plasma with
temperature :math:`T` is given by:
.. math::
f(v) = 2 \frac{1}{(π v_{Th}^2)^{1/2}} \exp(-(v - V_{drift})^2 / v_{Th}^2 )
where :math:`v_{Th} = \sqrt{2 k_B T / m}` is the thermal speed.
Examples
--------
>>> import astropy.units as u
>>> v = 1 * u.m / u.s
>>> Maxwellian_speed_1D(v=v, T=30000 * u.K, particle="e-", v_drift=0 * u.m / u.s)
<Quantity 1.1832...e-06 s / m>
"""
if units == "units":
# unit checks and conversions
# checking velocity units
v = v.to_value(SPEED_UNITS)
# Catching case where drift velocities have default values, they
# need to be assigned units
v_drift = _v_drift_conversion(v_drift)
# convert temperature to kelvin
T = T.to_value(u.K, equivalencies=u.temperature_energy())
if not np.isnan(vTh):
# check units of thermal velocity
vTh = vTh.to_value(SPEED_UNITS)
if np.isnan(vTh):
# get thermal velocity and thermal velocity squared
vTh = thermal_speed(
T << u.K, particle=particle, method="most_probable"
).to_value(SPEED_UNITS)
# Get thermal velocity squared
vThSq = vTh**2
# Get square of relative particle velocity
vSq = (v - v_drift) ** 2
# calculating distribution function
coeff = 2 * (vThSq * np.pi) ** (-1 / 2)
expTerm = np.exp(-vSq / vThSq)
distFunc = coeff * expTerm
if units == "units":
return distFunc << SPEED_DISTRIBUTION_UNITS_1D
elif units == "unitless":
return distFunc
else:
raise ValueError(f"Units must be either 'units' or 'unitless', got {units}).")
[docs]
@particle_input
def Maxwellian_speed_2D(
v,
T,
particle: ParticleLike = "e-",
v_drift=0,
vTh=np.nan,
units="units",
*,
mass_numb=None,
Z=None,
):
r"""
Probability distribution function of speed for a Maxwellian
distribution in 2D.
Return the probability density function of finding a particle with
speed components ``vx`` and ``vy`` in m/s in an equilibrium plasma
of temperature ``T`` which follows the 2D Maxwellian distribution
function. This function assumes Cartesian coordinates.
Parameters
----------
v: `~astropy.units.Quantity`
The speed in units convertible to m/s.
T: `~astropy.units.Quantity`
The temperature, preferably in kelvin.
particle: |particle-like|, optional
Representation of the particle species(e.g., ``'p+'`` for protons,
``'D+'`` for deuterium, or ``'He-4 +1'`` for singly ionized
helium-4), which defaults to electrons.
v_drift: `~astropy.units.Quantity`
The drift speed in units convertible to m/s.
vTh: `~astropy.units.Quantity`, optional
Thermal velocity (most probable) in m/s. This is used for
optimization purposes to avoid re-calculating ``vTh``, for
example when integrating over velocity-space.
units: `str`, optional
Selects whether to run function with units and unit checks (when
equal to "units") or to run as unitless (when equal to "unitless").
The unitless version is substantially faster for intensive
computations.
mass_numb : integer, |keyword-only|, optional
The mass number associated with ``particle``.
Z : real number, |keyword-only|, optional
The charge number associated with ``particle``.
Returns
-------
f : `~astropy.units.Quantity`
Probability density in \ :sup:`-1`\ , normalized so that:
:math:`\iiint_{0}^∞ f(\vec{v}) d\vec{v} = 1`.
Raises
------
`TypeError`
A parameter argument is not a `~astropy.units.Quantity` and
cannot be converted into a `~astropy.units.Quantity`.
`~astropy.units.UnitConversionError`
If the parameters is not in appropriate units.
`ValueError`
If the temperature is negative, or the particle mass or charge
state cannot be found.
Notes
-----
In 2D, the Maxwellian speed distribution function describing the
distribution of particles with speed :math:`v` in a plasma with
temperature :math:`T` is given by:
.. math::
f = 2 π v (π v_{Th}^2)^{-1} \exp(-v^2 / v_{Th}^2)
where :math:`v_{Th} = \sqrt{2 k_B T / m}` is the thermal speed.
See Also
--------
Maxwellian_speed_1D
Examples
--------
>>> import astropy.units as u
>>> v = 1 * u.m / u.s
>>> Maxwellian_speed_2D(v=v, T=30000 * u.K, particle="e-", v_drift=0 * u.m / u.s)
<Quantity 2.199...e-12 s / m>
"""
if v_drift != 0:
raise NotImplementedError("Non-zero drift speed is work in progress.")
if units == "units":
# unit checks and conversions
# checking velocity units
v = v.to_value(SPEED_UNITS)
# Catching case where drift velocity has default value, and
# needs to be assigned units
v_drift = _v_drift_conversion(v_drift)
# convert temperature to kelvin
T = T.to_value(u.K, equivalencies=u.temperature_energy())
if not np.isnan(vTh):
# check units of thermal velocity
vTh = vTh.to_value(SPEED_UNITS)
if np.isnan(vTh):
# get thermal velocity and thermal velocity squared
vTh = thermal_speed(
T << u.K, particle=particle, method="most_probable"
).to_value(SPEED_UNITS)
# getting square of thermal speed
vThSq = vTh**2
# get square of relative particle speed
vSq = (v - v_drift) ** 2
# calculating distribution function
coeff1 = (np.pi * vThSq) ** (-1)
coeff2 = 2 * np.pi * (v - v_drift)
expTerm = np.exp(-vSq / vThSq)
distFunc = coeff1 * coeff2 * expTerm
if units == "units":
return distFunc << SPEED_DISTRIBUTION_UNITS_1D
elif units == "unitless":
return distFunc
else:
raise ValueError(f"Units must be either 'units' or 'unitless', got {units}).")
[docs]
@particle_input
def Maxwellian_speed_3D(
v,
T,
particle: ParticleLike = "e-",
v_drift=0,
vTh=np.nan,
units="units",
*,
mass_numb=None,
Z=None,
):
r"""
Probability distribution function of speed for a Maxwellian
distribution in 3D.
Return the probability density function for finding a particle with
speed components ``vx``, ``vy``, and ``vz`` in m/s in an equilibrium
plasma of temperature ``T`` which follows the 3D Maxwellian
distribution function. This function assumes Cartesian coordinates.
Parameters
----------
v : `~astropy.units.Quantity`
The speed in units convertible to m/s.
T : `~astropy.units.Quantity`
The temperature, preferably in kelvin.
particle : `str`, optional
Representation of the particle species(e.g., ``'p+'`` for protons, ``'D+'``
for deuterium, or ``'He-4 +1'`` for :math:`He_4^{+1}`
(singly ionized helium-4), which defaults to electrons.
v_drift : `~astropy.units.Quantity`
The drift speed in units convertible to m/s.
vTh : `~astropy.units.Quantity`, optional
Thermal velocity (most probable) in m/s. This is used for
optimization purposes to avoid re-calculating vTh, for example
when integrating over velocity-space.
units : `str`, optional
Selects whether to run function with units and unit checks (when
equal to "units") or to run as unitless (when equal to "unitless").
The unitless version is substantially faster for intensive
computations.
mass_numb : integer, |keyword-only|, optional
The mass number associated with ``particle``.
Z : real number, |keyword-only|, optional
The charge number associated with ``particle``.
Returns
-------
f : `~astropy.units.Quantity`
Probability density in speed\ :sup:`-1`\ , normalized so that:
:math:`\iiint_0^∞ f(\vec{v}) d\vec{v} = 1`.
Raises
------
`TypeError`
A parameter argument is not a `~astropy.units.Quantity` and
cannot be converted into a `~astropy.units.Quantity`.
`~astropy.units.UnitConversionError`
If the parameters is not in appropriate units.
`ValueError`
If the temperature is negative, or the particle mass or charge state
cannot be found.
Notes
-----
In 3D, the Maxwellian speed distribution function describing
the distribution of particles with speed :math:`v` in a plasma with
temperature :math:`T` is given by:
.. math::
f = 4 π v^{2} (π v_{Th}^2)^{-3/2} \exp(-v^{2} / v_{Th}^2)
where :math:`v_{Th} = \sqrt{2 k_B T / m}` is the thermal speed.
See Also
--------
Maxwellian_speed_1D
Examples
--------
>>> import astropy.units as u
>>> v = 1 * u.m / u.s
>>> Maxwellian_speed_3D(v=v, T=30000 * u.K, particle="e-", v_drift=0 * u.m / u.s)
<Quantity 2.60235...e-18 s / m>
"""
if v_drift != 0:
raise NotImplementedError("Non-zero drift speed is work in progress.")
if units == "units":
# unit checks and conversions
# checking velocity units
v = v.to_value(SPEED_UNITS)
# Catching case where drift velocity has default value, and
# needs to be assigned units
v_drift = _v_drift_conversion(v_drift)
# convert temperature to kelvin
T = T.to_value(u.K, equivalencies=u.temperature_energy())
if not np.isnan(vTh):
# check units of thermal velocity
vTh = vTh.to_value(SPEED_UNITS)
if np.isnan(vTh):
# get thermal velocity and thermal velocity squared
vTh = thermal_speed(
T << u.K, particle=particle, method="most_probable"
).to_value(SPEED_UNITS)
# getting square of thermal speed
vThSq = vTh**2
# get square of relative particle speed
vSq = (v - v_drift) ** 2
# calculating distribution function
coeff1 = (np.pi * vThSq) ** (-3 / 2)
coeff2 = 4 * np.pi * vSq
expTerm = np.exp(-vSq / vThSq)
distFunc = coeff1 * coeff2 * expTerm
if units == "units":
return distFunc << SPEED_DISTRIBUTION_UNITS_1D
elif units == "unitless":
return distFunc
else:
raise ValueError(f"Units must be either 'units' or 'unitless', got {units}).")
[docs]
@particle_input
def kappa_velocity_1D(
v,
T,
kappa,
particle: ParticleLike = "e-",
v_drift=0,
vTh=np.nan,
units="units",
*,
mass_numb=None,
Z=None,
):
r"""
Return the probability density at the velocity ``v`` in m/s
to find a particle ``particle`` in a plasma of temperature ``T``
following the Kappa distribution function in 1D. The slope of the
tail of the Kappa distribution function is set by 'kappa', which
must be greater than :math:`1/2`.
Parameters
----------
v : `~astropy.units.Quantity`
The velocity in units convertible to m/s.
T : `~astropy.units.Quantity`
The temperature in kelvin.
kappa : `~astropy.units.Quantity`
The kappa parameter is a dimensionless number which sets the slope
of the energy spectrum of suprathermal particles forming the tail
of the Kappa velocity distribution function. Kappa must be greater
than :math:`3/2`.
particle : `str`, optional
Representation of the particle species(e.g., ``'p`` for protons, ``'D+'``
for deuterium, or ``'He-4 +1'`` for :math:`He_4^{+1}`
(singly ionized helium-4)), which defaults to electrons.
v_drift : `~astropy.units.Quantity`, optional
The drift velocity in units convertible to m/s.
vTh : `~astropy.units.Quantity`, optional
Thermal velocity (most probable) in m/s. This is used for
optimization purposes to avoid re-calculating ``vTh``, for example
when integrating over velocity-space.
units : `str`, optional
Selects whether to run function with units and unit checks (when
equal to ``"units"``) or to run as unitless (when equal to
``"unitless"``). The unitless version is substantially faster for
intensive computations.
mass_numb : integer, |keyword-only|, optional
The mass number associated with ``particle``.
Z : real number, |keyword-only|, optional
The charge number associated with ``particle``.
Returns
-------
f : `~astropy.units.Quantity`
Probability density in velocity\ :sup:`-1`\ , normalized so that
:math:`\int_{-∞}^{+∞} f(v) dv = 1`.
Raises
------
`TypeError`
A parameter argument is not a `~astropy.units.Quantity` and
cannot be converted into a `~astropy.units.Quantity`.
`~astropy.units.UnitConversionError`
If the parameters is not in appropriate units.
`ValueError`
If the temperature is negative, or the particle mass or charge state
cannot be found.
Notes
-----
In one dimension, the Kappa velocity distribution function describing
the distribution of particles with speed :math:`v` in a plasma with
temperature :math:`T` and suprathermal parameter :math:`κ` is
given by:
.. math::
f = A_κ \left(1 + \frac{(\vec{v} -
\vec{V_{drift}})^2}{κ v_{Th},κ^2}\right)^{-κ}
where :math:`v_{Th},κ` is the kappa thermal speed
and :math:`A_κ = \frac{1}{\sqrt{π} κ^{3/2} v_{Th},κ^2
\frac{Γ(κ + 1)}{Γ(κ - 1/2)}}`
is the normalization constant.
As :math:`κ → ∞`, the kappa distribution function converges to the
Maxwellian distribution function.
Examples
--------
>>> import astropy.units as u
>>> v = 1 * u.m / u.s
>>> kappa_velocity_1D(
... v=v,
... T=30000 * u.K,
... kappa=4,
... particle="e-",
... v_drift=0 * u.m / u.s,
... )
<Quantity 6.75549...e-07 s / m>
See Also
--------
kappa_velocity_3D
~plasmapy.formulary.speeds.kappa_thermal_speed
"""
# must have kappa > 3/2 for distribution function to be valid
if kappa <= 3 / 2:
raise ValueError(f"Must have κ > 3/2, instead of {kappa}.")
if units == "units":
# unit checks and conversions
# checking velocity units
v = v.to_value(SPEED_UNITS)
# catching case where drift velocities have default values
v_drift = _v_drift_conversion(v_drift)
# convert temperature to kelvin
T = T.to_value(u.K, equivalencies=u.temperature_energy())
if not np.isnan(vTh):
# check units of thermal velocity
vTh = vTh.to_value(SPEED_UNITS)
if np.isnan(vTh):
# get thermal velocity and thermal velocity squared
vTh = kappa_thermal_speed(T << u.K, kappa, particle=particle).to_value(
SPEED_UNITS
)
# Get thermal velocity squared and accounting for 1D instead of 3D
vThSq = vTh**2
# Get square of relative particle velocity
vSq = (v - v_drift) ** 2
# calculating distribution function
expTerm = (1 + vSq / (kappa * vThSq)) ** (-kappa)
coeff1 = 1 / (np.sqrt(np.pi) * kappa ** (3 / 2) * vTh)
coeff2 = gamma(kappa + 1) / (gamma(kappa - 1 / 2))
distFunc = coeff1 * coeff2 * expTerm
if units == "units":
return distFunc << SPEED_DISTRIBUTION_UNITS_1D
elif units == "unitless":
return distFunc
else:
raise ValueError(f"Units must be either 'units' or 'unitless', got {units}).")
[docs]
@particle_input
def kappa_velocity_3D(
vx,
vy,
vz,
T,
kappa,
particle: ParticleLike = "e-",
vx_drift=0,
vy_drift=0,
vz_drift=0,
vTh=np.nan,
units="units",
*,
mass_numb=None,
Z=None,
):
r"""
Return the probability density function for finding a particle with
velocity components ``v_x``, ``v_y``, and ``v_z``in m/s in a suprathermal
plasma of temperature ``T`` and parameter ``kappa`` which follows the
3D Kappa distribution function. This function assumes Cartesian
coordinates.
Parameters
----------
vx : `~astropy.units.Quantity`
The velocity in x-direction units convertible to m/s.
vy : `~astropy.units.Quantity`
The velocity in y-direction units convertible to m/s.
vz : `~astropy.units.Quantity`
The velocity in z-direction units convertible to m/s.
T : `~astropy.units.Quantity`
The temperature, preferably in kelvin.
kappa : `~astropy.units.Quantity`
The kappa parameter is a dimensionless number which sets the slope
of the energy spectrum of suprathermal particles forming the tail
of the Kappa velocity distribution function. ``kappa`` must be greater
than :math:`3/2`.
particle : `str`, optional
Representation of the particle species(e.g., 'p+' for protons, 'D+'
for deuterium, or 'He-4 +1' for :math:`He_4^{+1}` : singly ionized
helium-4)), which defaults to electrons.
vx_drift : `~astropy.units.Quantity`, optional
The drift velocity in x-direction units convertible to m/s.
vy_drift : `~astropy.units.Quantity`, optional
The drift velocity in y-direction units convertible to m/s.
vz_drift : `~astropy.units.Quantity`, optional
The drift velocity in z-direction units convertible to m/s.
vTh : `~astropy.units.Quantity`, optional
Thermal velocity (most probable) in m/s. This is used for
optimization purposes to avoid re-calculating ``vTh``, for example
when integrating over velocity-space.
units : `str`, optional
Selects whether to run function with units and unit checks (when
equal to "units") or to run as unitless (when equal to "unitless").
The unitless version is substantially faster for intensive
computations.
mass_numb : integer, |keyword-only|, optional
The mass number associated with ``particle``.
Z : real number, |keyword-only|, optional
The charge number associated with ``particle``.
Returns
-------
f : `~astropy.units.Quantity`
Probability density in units of inverse velocity, normalized so that:
:math:`\iiint_{0}^∞ f(\vec{v}) d\vec{v} = 1`
Raises
------
`TypeError`
If any of the parameters is not a `~astropy.units.Quantity` and
cannot be converted into one.
`~astropy.units.UnitConversionError`
If the parameters is not in appropriate units.
`ValueError`
If the temperature is negative, or the particle mass or charge state
cannot be found.
Notes
-----
In three dimensions, the Kappa velocity distribution function describing
the distribution of particles with speed :math:`v` in a plasma with
temperature :math:`T` and suprathermal parameter :math:`κ` is given by:
.. math::
f = A_κ \left(1 + \frac{(\vec{v} -
\vec{V_{drift}})^2}{κ v_{Th},κ^2}\right)^{-(κ + 1)}
where :math:`v_{Th},κ` is the kappa thermal speed
and :math:`A_κ = \frac{1}{2 π (κ v_{Th},κ^2)^{3/2}}
\frac{Γ(κ + 1)}{Γ(κ - 1/2) Γ(3/2)}` is the
normalization constant.
As :math:`κ → ∞`, the kappa distribution function converges to the
Maxwellian distribution function.
See Also
--------
kappa_velocity_1D
~plasmapy.formulary.speeds.kappa_thermal_speed
Examples
--------
>>> import astropy.units as u
>>> v = 1 * u.m / u.s
>>> kappa_velocity_3D(
... vx=v,
... vy=v,
... vz=v,
... T=30000 * u.K,
... kappa=4,
... particle="e-",
... vx_drift=0 * u.m / u.s,
... vy_drift=0 * u.m / u.s,
... vz_drift=0 * u.m / u.s,
... )
<Quantity 3.7833...e-19 s3 / m3>
"""
# must have kappa > 3/2 for distribution function to be valid
if kappa <= 3 / 2:
raise ValueError(f"Must have kappa > 3/2, instead of {kappa}.")
if units == "units":
# unit checks and conversions
# checking velocity units
vx = vx.to_value(SPEED_UNITS)
vy = vy.to_value(SPEED_UNITS)
vz = vz.to_value(SPEED_UNITS)
# Catching case where drift velocities have default values
vx_drift = _v_drift_conversion(vx_drift)
vy_drift = _v_drift_conversion(vy_drift)
vz_drift = _v_drift_conversion(vz_drift)
# convert temperature to kelvin
T = T.to_value(u.K, equivalencies=u.temperature_energy())
if not np.isnan(vTh):
# check units of thermal velocity
vTh = vTh.to_value(SPEED_UNITS)
if np.isnan(vTh):
# get thermal velocity and thermal velocity squared
vTh = kappa_thermal_speed(T << u.K, kappa, particle=particle).to_value(
SPEED_UNITS
)
# getting square of thermal velocity
vThSq = vTh**2
# Get square of relative particle velocity
vSq = (vx - vx_drift) ** 2 + (vy - vy_drift) ** 2 + (vz - vz_drift) ** 2
# calculating distribution function
expTerm = (1 + vSq / (kappa * vThSq)) ** (-(kappa + 1))
coeff1 = 1 / (2 * np.pi * (kappa * vThSq) ** (3 / 2))
coeff2 = gamma(kappa + 1) / (gamma(kappa - 1 / 2) * gamma(3 / 2))
distFunc = coeff1 * coeff2 * expTerm
if units == "units":
return distFunc << SPEED_DISTRIBUTION_UNITS_3D
elif units == "unitless":
return distFunc
else:
raise ValueError(f"Units must be either 'units' or 'unitless', got {units}).")