Particles and Formulary

We’ll start this notebook with a quick refresher on astropy.units, before going through examples from plasmapy.particles ⚛️ and plasmapy.formulary 🧮. After an interlude about type hint annotations and decorators, we’ll discuss how to create a formulary function.

Let’s start with some preliminary imports & settings. To execute a cell in a Jupyter notebook, press Shift + Enter. If you have to restart the notebook, please re-execute the following cell.

If using Google Colab, click Run anyway when prompted. (If prompted again, select Restart runtime when the installation finishes.)

import sys

if 'google.colab' in str(get_ipython()):
    if 'plasmapy' not in sys.modules:
        !pip install plasmapy==2024.7.0

import astropy.units as u
import numpy as np
from plasmapy.particles import *
from plasmapy.particles.particle_class import valid_categories
from plasmapy.formulary import *
from plasmapy.utils.decorators import validate_quantities

Quick refresher on astropy.units

We’ll be using astropy.units throughout the rest of this notebook, so we’ll begin with a brief reminder o fhow to use it. We typically import this subpackage as u.

We can create a physical quantity by multiplying or dividing a number or array with a unit.

60 * u.km

$60 \; \mathrm{km}$

This operation creates a Quantity object: a number, sequence, or array that has been assigned a physical unit. We can create Quantity objects with compound units.

V = 88 * u.imperial.mile / u.hour
print(V)

88.0 mi / h

We can use Quantity.to() to convert a Quantity to different units.

V.to("m/s")

$39.33952 \; \mathrm{\frac{m}{s}}$

Particles

The plasmapy.particles subpackage contains functions to access basic particle data and classes to represent particles.

from plasmapy.particles import *

Particle properties

There are several functions that provide information about different particles that might be present in a plasma. The input of these functions is a representation of a particle, such as a string for the atomic symbol or the element name.

atomic_number("Fe")

26

We can provide a number that represents the atomic number.

element_name(26)

'iron'

We can provide standard symbols or the names of particles.

is_stable("e-")

True

charge_number("proton")

1

We can represent isotopes with the atomic symbol followed by a hyphen and the mass number. Let’s use half_life to return the half-life of a radioactive particle in seconds as a Quantity.

half_life("C-14")

$1.8082505 \times 10^{11} \; \mathrm{s}$

We typically represent an ion in a string by putting together the atomic symbol or isotope symbol, a space, the charge number, and the sign of the charge.

charge_number("Fe-56 13+")

13

Functions in plasmapy.particles are quite flexible in terms of string inputs representing particles. An input is particle-like if it can be used to represent a physical particle.

particle_mass("iron-56 +13")

$9.2870305 \times 10^{-26} \; \mathrm{kg}$

particle_mass("iron-56+++++++++++++")

$9.2870305 \times 10^{-26} \; \mathrm{kg}$

Most of these functions take additional arguments, with Z representing the charge number of an ion and mass_numb representing the mass number of an isotope. These arguments are often keyword-only to avoid ambiguity.

particle_mass("Fe", Z=13, mass_numb=56)

$9.2870305 \times 10^{-26} \; \mathrm{kg}$

Particle objects

Up until now, we have been using functions that accept representations of particles and then return particle properties. With the Particle class, we can create objects that represent physical particles.

proton = Particle("p+")
electron = Particle("electron")
iron56_nucleus = Particle("Fe", Z=26, mass_numb=56)

Particle properties can be accessed via attributes of the Particle class.

proton.mass

$1.6726219 \times 10^{-27} \; \mathrm{kg}$

electron.charge

$-1.6021766 \times 10^{-19} \; \mathrm{C}$

electron.charge_number

-1

iron56_nucleus.nuclear_binding_energy

$7.8868678 \times 10^{-11} \; \mathrm{J}$

Antiparticles

We can get antiparticles of fundamental particles by using the antiparticle attribute of a Particle.

electron.antiparticle

Particle("e+")

We can also use the tilde (~) operator on a Particle to get its antiparticle.

~proton

Particle("p-")

Ionization and recombination

The recombine() and ionize() methods of a Particle representing an ion or neutral atom will return a different Particle with fewer or more electrons.

deuterium = Particle("D 0+")
deuterium.ionize()

Particle("D 1+")

When provided with a number, these methods tell how many bound electrons to add or remove.

alpha = Particle("alpha")
alpha.recombine(2)

Particle("He-4 0+")

Custom particles

Sometimes we want to use a particle with custom properties. For example, we might want to represent an average ion in a multi-species plasma or a dust particle. For that we can use CustomParticle.

cp = CustomParticle(9e-26 * u.kg, 2.18e-18 * u.C, symbol="Fe 13.6+")

Many of the attributes of CustomParticle are the same as in Particle.

cp.mass

$9 \times 10^{-26} \; \mathrm{kg}$

cp.charge

$2.18 \times 10^{-18} \; \mathrm{C}$

cp.symbol

'Fe 13.6+'

If we do not include one of the physical quantities, it gets set to numpy.nan (not a number) in the appropriate units.

CustomParticle(9.27e-26 * u.kg).charge

${\rm NaN} \; \mathrm{C}$

CustomParticle objects can provided to most of the commonly used functions in plasmapy.formulary, and we’re planning to improve interoperability in future releases of PlasmaPy.

Particle lists

The ParticleList class is a container for Particle and CustomParticle objects.

iron_ions = ParticleList(["Fe 12+", "Fe 13+", "Fe 14+"])

By using a ParticleList, we can access the properties of multiple particles at once.

iron_ions.mass

$[9.2721873 \times 10^{-26},~9.2720962 \times 10^{-26},~9.2720051 \times 10^{-26}] \; \mathrm{kg}$

iron_ions.charge

$[1.922612 \times 10^{-18},~2.0828296 \times 10^{-18},~2.2430473 \times 10^{-18}] \; \mathrm{C}$

iron_ions.symbols

['Fe 12+', 'Fe 13+', 'Fe 14+']

We can also create a ParticleList by adding Particle and/or CustomParticle objects together.

proton + electron

ParticleList(['p+', 'e-'])

We can also get an average particle.

iron_ions.average_particle()

CustomParticle(mass=9.272096197546505e-26 kg, charge=2.0828296241999997e-18 C)

ParticleList objects can also be provided to the most commonly used functions in plasmapy.formulary, with more complete interoperability expected in the future.

Particle Categorization

The categories attribute of a Particle provides a set of the categories that the particle belongs to as a set.

muon = Particle("muon")
muon.categories

{'charged', 'fermion', 'lepton', 'matter', 'unstable'}

The is_category() method of a Particle lets us determine if the particle belongs to one or more categories.

muon.is_category("lepton")

True

If we need to be more specific with is_category(), we can use its keyword arguments:

• require: for categories that a Particle must belong to

• exclude: for categories that the Particle cannot belong to

• any_of: for categories of which a Particle must belong to at least one of

electron.is_category(require="lepton", exclude="baryon", any_of={"boson", "fermion"})

True

All valid particle categories are included in plasmapy.particles.particle_class.valid_categories.

print(valid_categories)

{'baryon', 'boson', 'lepton', 'antilepton', 'antimatter', 'stable', 'post-transition metal', 'fermion', 'alkali metal', 'isotope', 'noble gas', 'neutrino', 'nonmetal', 'element', 'antineutrino', 'transition metal', 'uncharged', 'charged', 'electron', 'proton', 'matter', 'ion', 'alkaline earth metal', 'lanthanide', 'metal', 'antibaryon', 'custom', 'halogen', 'unstable', 'actinide', 'metalloid', 'neutron', 'positron'}

The is_category() method of ParticleList returns True if all particles match the criteria, and False otherwise.

particles = ParticleList(["e-", "p+", "n"])
particles.is_category(require="lepton")

False

If we want to do this operation for all of the particles, we can set the particlewise argument to True. We will then get a list of booleans.

particles.is_category(require="lepton", particlewise=True)

[True, False, False]

Nuclear reactions

We can use plasmapy.particles to calculate the energy of a nuclear reaction using the > operator.

deuteron = Particle("D+")
triton = Particle("T+")
alpha = Particle("α")
neutron = Particle("n")

energy = deuteron + triton > alpha + neutron

energy.to("MeV")

$17.589253 \; \mathrm{MeV}$

If the nuclear reaction is invalid, then an exception is raised that states the reason why.

deuteron + triton > alpha + 56 * neutron

---------------------------------------------------------------------------
ParticleError                             Traceback (most recent call last)
Cell In[44], line 1
----> 1 deuteron + triton > alpha + 56 * neutron

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/plasmapy/particles/particle_collections.py:229, in ParticleList.__gt__(self, other)
    226 from plasmapy.particles.nuclear import nuclear_reaction_energy
    228 other_as_particle_list = self._cast_other_as_particle_list(other)
--> 229 return nuclear_reaction_energy(
    230     reactants=self.symbols, products=other_as_particle_list.symbols
    231 )

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/plasmapy/particles/nuclear.py:309, in nuclear_reaction_energy(*args, **kwargs)
    306         raise ParticleError(errmsg) from e
    308 if total_baryon_number(reactants) != total_baryon_number(products):
--> 309     raise ParticleError(
    310         f"The baryon number is not conserved for {reactants = } and {products = }."
    311     )
    313 if total_charge(reactants) != total_charge(products):
    314     raise ParticleError(
    315         f"Total charge is not conserved for {reactants = } and {products = }."
    316     )

ParticleError: The baryon number is not conserved for reactants = [Particle("D 1+"), Particle("T 1+")] and products = [Particle("He-4 2+"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n"), Particle("n")].

PlasmaPy formulary

The plasmapy.formulary subpackage contains a broad variety of formulas needed by plasma scientists across disciplines, in particular to calculate plasma parameters.

from plasmapy.formulary import *

Plasma beta in the solar corona

Plasma beta (\(β\)) is one of the most fundamental plasma parameters. \(β\) is the ratio of the plasma (gas) pressure to the magnetic pressure. How a plasma behaves depends strongly on \(β\). When \(β ≫ 1\), the magnetic field is not strong enough to exert much of a force on the plasma, so its motions start to resemble a gas. When \(β ≪ 1\), magnetic tension and pressure are the dominant macroscopic forces.

Let’s use plasmapy.formulary to calculate plasma β in different regions of the solar atmosphere and see what we can learn.

Solar corona

Let’s start by defining some plasma parameters for an active region in the solar corona.

B_corona = 50 * u.G
n_corona = 1e9 * u.cm ** -3
T_corona = 1 * u.MK

When we use these parameters in beta, we find that \(β\) is quite small so that the corona is magnetically dominated.

beta(T=T_corona, n=n_corona, B=B_corona)

$0.0013879798 \; \mathrm{}$

Solar photosphere

Let’s specify some characteristic plasma parameters for the solar photosphere, away from any sunspots.

T_photosphere = 5800 * u.K
B_photosphere = 400 * u.G
n_photosphere = 1e17 * u.cm ** -3

When we calculate for the photosphere, we find that it is an order of magnitude larger than 1, so plasma pressure forces are more important than magnetic tension and pressure.

beta(T_photosphere, n_photosphere, B_photosphere)

$12.578567 \; \mathrm{}$

Plasma parameters in Earth’s magnetosphere

The Magnetospheric Multiscale Mission (MMS) is a constellation of four identical spacecraft. The goal of MMS is to investigate the small-scale physics of magnetic reconnection in Earth’s magnetosphere. In order to do this, the spacecraft need to orbit in a tight configuration. But how tight does the tetrahedron have to be? Let’s use plasmapy.formulary to find out.

Physics background

Magnetic reconnection is the fundamental plasma process that converts stored magnetic energy into kinetic energy, thermal energy, and particle acceleration. Reconnection powers solar flares and is a key component of geomagnetic storms in Earth’s magnetosphere. Reconnection can also degrade confinement in fusion devices such as tokamaks.

The inertial length is the characteristic length scale for a particle to get accelerated or decelerated by electromagnetic forces in a plasma.

When the reconnection layer thickness is shorter than the ion inertial length, \(d_i ≡ c/ω_{pi}\), collisionless effects and the Hall effect enable reconnection to be fast (Zweibel & Yamada 2009). The inner electron diffusion region has a thickness of about the electron inertial length, \(d_e ≡ c/ω_{pi}\). (Here, \(ω_{pi}\) and \(ω_{pe}\) are the ion and electron plasma frequencies.)

Our goal: calculate \(d_i\) and \(d_e\) to get an idea of how far the MMS spacecraft should be separated from each other to investigate reconnection.

Length scales

Let’s choose some characteristic plasma parameters for the magnetosphere.

n = 1 * u.cm ** -3
B = 5 * u.nT
T = 30000 * u.K

Let’s calculate the ion inertial length, \(d_i\). On length scales shorter than \(d_i\), the Hall effect becomes important as the ions and electrons decouple from each other.

inertial_length(n=n, particle="p+").to("km")

$227.71077 \; \mathrm{km}$

The ion diffusion regions should therefore be a few hundred kilometers thick. Let’s calculate the electron inertial length next.

inertial_length(n=n, particle="e-").to("km")

$5.3140933 \; \mathrm{km}$

The electron diffusion region should therefore have a characteristic length scale of a few kilometers, which is significantly smaller than the ion diffusion region.

We can also calculate the gyroradii for different particles

gyroradius(B=B, particle=["p+", "e-"], T=T).to("km")

$[46.466051,~1.0843852] \; \mathrm{km}$

The four MMS spacecraft have separations of ten to hundreds of kilometers, and thus are well-positioned to investigate Hall physics during reconnection in the magnetosphere.

Frequencies

We can calculate some of the fundamental frequencies associated with magnetospheric plasma.

plasma_frequency(n=n, particle=["p+", "e-"])

$[1316.5493,~56414.602] \; \mathrm{\frac{rad}{s}}$

gyrofrequency(B=B, particle=["p+", "e-"])

$[0.47894166,~879.41001] \; \mathrm{\frac{rad}{s}}$

lower_hybrid_frequency(B=B, n_i=n, ion="p+")

$20.520326 \; \mathrm{\frac{rad}{s}}$

Most of the functions in plasmapy.formulary have descriptions of the plasma parameters that include the formula and physical interpretation. If we ever forget what the lower_hybrid_frequency represents, we can check out its documentation page!

Python interlude

Now that we’ve gone through how to use plasmapy.particles and plasmapy.formulary, we’ll move on to the process for adding a function to plasmapy.formulary. Before that, we need to introduce two capabilities for the Python language:

• Type hint annotations

• Decorators

Type hint annotations

In a statically typed language, variable types are explicitly declared. If s is declared to be a string, then s will always be a string. Type errors are found when code is compiled.

Python is a dynamically typed language. Variable types are determined at runtime rather than at compile time, and it’s possible for variables to even change types!

There are tradeoffs! ⚖️ Dynamically are more flexible, but at the cost of reduced type safety.

Type hint annotations provide a middle ground between statically vs. dynamically typed languages.

🏷 Let’s define a variable called name and provide it with a type hint annotation that says the variable should be a a str.

name: str = "Darth Fortran"

The type hint of str in the above cell didn’t actually do anything when we ran it. But when we read the code, it does tell us what type to expect name to be.

How about a list that starts empty but will eventually contain names?

names: list[str] = []

This type hint still doesn’t do anything at runtime, but it helps us read the code, and lets us use code quality tools that can help us find errors.

We can specify that a variable should be one of multiple types with |:

identifier: str | int

Type hints are particularly helpful when defining functions:

def stringify(number: int) -> str:
    return str(num)

Many Python packages, including PlasmaPy, make use of static type checking tools like mypy. These tools can help us find type errors such as in the following function before we even run the code:

def return_argument(x: int) -> str:
    return x

What if we want to specify that a function accepts a length and returns a volume?

def volume(d: u.Quantity[u.m]) -> u.Quantity[u.m**3]:
    return d**3

Decorators

In Python, functions are objects. This means that we can write a function that returns another function.

def return_function():

    print("Defining inner_function...")

    def inner_function():
        print("Calling inner_function!")

    return inner_function

function = return_function()

Defining inner_function...

function()

Calling inner_function!

Or we can pass a function as an argument to another function! We can use typing.Callable as the corresponding type hint annotation.

from typing import Callable

def apply_function(function: Callable, array):
    return function(array)

array = [1, 2, 3, 4, 5]

apply_function(max, array)

5

A function that modifies another function is a decorator.

Decorators in Python are a way to modify or enhance functions without changing their actual code. They wrap another function, potentially adding extra functionality before and after the original function runs.

def decorator(function: Callable):

    def decorated_function() -> None:

        print("Before calling the function.")
        result = function()
        print("After calling the function.")

        return result

    return decorated_function

def example_function():
    print("Inside original example_function!")

Let’s try it out!

modified_function = decorator(example_function)

example_function()

Inside original example_function!

In Python, we have syntactic sugar for decorators, using @:

@decorator
def function2():
    print("Inside function2!")

function2()

Before calling the function.
Inside function2!
After calling the function.

Fibonacci sequence

Let’s look at a recursive function that computes the \(n\)th Fibonacci number. If we define \(F_0 ≡ 0\) and \(F_1 ≡ 1\), then $ F_n = F_{n-1} + F_{n-2}$.

def fibonacci(n: int) -> int:
    print("Calculating Fibonacci number for n =", n)
    return n if n < 2 else fibonacci(n - 1) + fibonacci(n - 2)

fibonacci(2)

Calculating Fibonacci number for n = 2
Calculating Fibonacci number for n = 1
Calculating Fibonacci number for n = 0

1

fibonacci(4)

Calculating Fibonacci number for n = 4
Calculating Fibonacci number for n = 3
Calculating Fibonacci number for n = 2
Calculating Fibonacci number for n = 1
Calculating Fibonacci number for n = 0
Calculating Fibonacci number for n = 1
Calculating Fibonacci number for n = 2
Calculating Fibonacci number for n = 1
Calculating Fibonacci number for n = 0

3

[@functools.cache]: https://docs.python.org/3/library/functools.html#functools.cache

We’ll use [@functools.cache] to store the output of a function for a particular argument.

from functools import cache


@cache
def fibonacci_cached(n: int) -> int:
    print("Calculating Fibonacci number for n =", n)
    return n if n < 2 else fibonacci_cached(n - 1) + fibonacci_cached(n - 2)

fibonacci_cached(2)

Calculating Fibonacci number for n = 2
Calculating Fibonacci number for n = 1
Calculating Fibonacci number for n = 0

1

fibonacci_cached(5)

Calculating Fibonacci number for n = 5
Calculating Fibonacci number for n = 4
Calculating Fibonacci number for n = 3

5

Let’s try calling it again!

fibonacci_cached(5)

5

Writing a formulary function

[@particle_input]: https://docs.plasmapy.org/en/stable/api/plasmapy.particles.decorators.particle_input.html [@validate_quantities]: https://docs.plasmapy.org/en/stable/api/plasmapy.utils.decorators.validators.validate_quantities.html#validate-quantities [@check_relativistic]: https://docs.plasmapy.org/en/stable/api/plasmapy.utils.decorators.checks.check_relativistic.html#check-relativistic [astropy.units]: https://docs.astropy.org/en/stable/units/index.html [plasmapy.particles]: https://docs.plasmapy.org/en/stable/particles/index.html [plasmapy.formulary]: https://docs.plasmapy.org/en/stable/formulary/index.html

PlasmaPy has three decorators that help us with writing formulary functions.

• [@check_relativistic]

• [@validate_quantities]

• [@particle_input]

Let’s get some imports out of the way.

import astropy.units as u
from astropy import constants
from plasmapy.utils.decorators import check_relativistic
from plasmapy.particles import particle_input

Checking relativistic velocities

[@check_relativistic]: https://docs.plasmapy.org/en/stable/api/plasmapy.utils.decorators.checks.check_relativistic.html#check-relativistic [plasmapy.formulary]: https://docs.plasmapy.org/en/stable/formulary/index.html

Many functions that return velocities in [plasmapy.formulary] assume that relativistic velocities are unimportant. To check for this, we can use the [@check_relativistic] decorator.

@check_relativistic
def speed(distance, time):
    return distance / time

Let’s try this with a velocity approaching the speed of light.

speed(299792457 * u.m, 1 * u.s)

/home/docs/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/plasmapy/utils/decorators/checks.py:1405: RelativityWarning: speed is yielding a velocity that is 100.0% of the speed of light. Relativistic effects may be important.
  warnings.warn(

$2.9979246 \times 10^{8} \; \mathrm{\frac{m}{s}}$

Now let’s try this with a velocity exceeding the speed of light.

speed(299792459 * u.m, 1 * u.s)

---------------------------------------------------------------------------
RelativityError                           Traceback (most recent call last)
Cell In[86], line 1
----> 1 speed(299792459 * u.m, 1 * u.s)

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/plasmapy/utils/decorators/checks.py:1326, in check_relativistic.<locals>.decorator.<locals>.wrapper(*args, **kwargs)
   1322 @preserve_signature
   1323 @functools.wraps(f)
   1324 def wrapper(*args, **kwargs):
   1325     return_ = f(*args, **kwargs)
-> 1326     _check_relativistic(return_, f.__name__, betafrac=betafrac)
   1327     return return_

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/plasmapy/utils/decorators/checks.py:1400, in _check_relativistic(V, funcname, betafrac)
   1398     raise RelativityError(f"{funcname} is yielding an infinite velocity.")
   1399 elif beta >= 1:
-> 1400     raise RelativityError(
   1401         f"{funcname} is yielding a velocity that is {round(beta, 3)} "
   1402         f"times the speed of light."
   1403     )
   1404 elif beta >= betafrac:
   1405     warnings.warn(
   1406         f"{funcname} is yielding a velocity that is "
   1407         f"{round(beta * 100, 3)}% of the speed of "
   1408         f"light. Relativistic effects may be important.",
   1409         RelativityWarning,
   1410     )

RelativityError: speed is yielding a velocity that is 1.0 times the speed of light.

Validating quantities

Most plasmapy.formulary functions accept and return Quantity objects from astropy.units. For example, we might write a function to calculate the magnetic pressure:

\[p_B ≡ \frac{B^2}{2μ_o}\]

def magnetic_pressure(B):
    return B**2 / (2 * constants.mu0)

magnetic_pressure(1 * u.T)

$397887.36 \; \mathrm{\frac{A^{2}\,T^{2}}{N}}$

But the above function doesn’t prevent us from making mistakes simple mistakes…

magnetic_pressure(1 * u.B)

$397887.36 \; \mathrm{\frac{A^{2}\,byte^{2}}{N}}$

[@validate_quantities]: https://docs.plasmapy.org/en/latest/api/plasmapy.utils.decorators.validators.validate_quantities.html#validate-quantities

The [@validate_quantities] decorator makes sure that

@validate_quantities
def magnetic_pressure(B: u.Quantity[u.T]) -> u.Quantity[u.Pa]:
    return B**2 / (2 * constants.mu0)

magnetic_pressure(1 * u.T)

$397887.36 \; \mathrm{Pa}$

magnetic_pressure(1 * u.B)

---------------------------------------------------------------------------
UnitTypeError                             Traceback (most recent call last)
Cell In[92], line 1
----> 1 magnetic_pressure(1 * u.B)

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/plasmapy/utils/decorators/validators.py:199, in ValidateQuantities.__call__.<locals>.wrapper(*args, **kwargs)
    196         continue
    198     # validate argument & update for conversion
--> 199     arg = self._validate_quantity(
    200         bound_args.arguments[arg_name], arg_name, validations[arg_name]
    201     )
    202     bound_args.arguments[arg_name] = arg
    204 # call function

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/plasmapy/utils/decorators/validators.py:376, in ValidateQuantities._validate_quantity(self, arg, arg_name, arg_validations)
    374     arg = arg.to(unit, equivalencies=equiv)
    375 elif err is not None:
--> 376     raise err
    378 self._check_value(arg, arg_name, arg_validations)
    380 return arg

UnitTypeError: The argument 'B' to function magnetic_pressure() should be an astropy Quantity with the following unit: T

magnetic_pressure(1)

WARNING: UnitsWarning: The argument 'B' to function magnetic_pressure() has no specified units. Assuming units of T. To silence this warning, explicitly pass in an astropy Quantity (e.g. 5. * astropy.units.cm) (see http://docs.astropy.org/en/stable/units/) [plasmapy.utils.decorators.validators]

$397887.36 \; \mathrm{Pa}$

[@validate_quantities]: https://docs.plasmapy.org/en/latest/api/plasmapy.utils.decorators.validators.validate_quantities.html#validate-quantities [equivalency]: https://docs.astropy.org/en/stable/units/equivalencies.html [astropy.units]: https://docs.astropy.org/en/stable/units/

We can also provide arguments directly to [@validate_quantities], for example if we want to do a validation on the return value whilst using an [equivalency] from [astropy.units].

@validate_quantities(
    validations_on_return={"units": u.K, "equivalencies": u.temperature_energy()}
)
def get_temperature(T):
    return T

get_temperature(1 * u.eV)

$11604.518 \; \mathrm{K}$

The docstring for [@validate_quantities](https://docs.plasmapy.org/en/latest/api/plasmapy.utils.decorators.validators.validate_quantities.html#validate-quantities) contains examples for several more validations that it can do, such as:

• Allowing only non-negative values

• Restricting/allowing complex values

• Restricting/allowing numpy.nan values

• Allowing None instead of a Quantity

Particle inputs

Plasma parameters often depend on particle properties, especially the particle’s charge and mass.

Early on, we wrote formulary functions like this:

def get_particle_mass(particle_string):
    particle_object = Particle(particle_string)
    return particle_object.mass

get_particle_mass("p+")

$1.6726219 \times 10^{-27} \; \mathrm{kg}$

The above example isn’t particularly bad. But if we wanted to make a function compatible with ParticleList objects too, then it started to become a mess.

from collections.abc import Iterable


def get_particle_mass(particle):
    if isinstance(particle, str):
        particle_object = Particle(particle)
    elif isinstance(particle, Iterable):
        particle_object = ParticleList(particle)
    return particle_object.mass

get_particle_mass(["p+", "e-"])

$[1.6726219 \times 10^{-27},~9.1093837 \times 10^{-31}] \; \mathrm{kg}$

The @particle_input decorator transforms arguments annotated with ParticleLike into a Particle, CustomParticle, or ParticleList. This transformation occurs before the argument gets passed into the decorated function, so we end up with cleaner code!

@particle_input
def make_particle(particle: ParticleLike):
    return particle

Let’s try passing it a string representing a particle.

make_particle("p+")

Particle("p+")

[@particle_input]: https://docs.plasmapy.org/en/stable/api/plasmapy.particles.decorators.particle_input.html

The argument got converted into a Particle object by [@particle_input].

If we add Z and mass_numb as parameters, they’ll get incorporated into the particle.

@particle_input
def make_particle(particle: ParticleLike, Z=None, mass_numb=None):
    return particle

make_particle("He", Z=2, mass_numb=4)

Particle("He-4 2+")

[@particle_input]: https://docs.plasmapy.org/en/stable/api/plasmapy.particles.decorators.particle_input.html

If we name the parameter element, ion, or isotope, then [@particle_input] will make sure that the particle matches that identity.

@particle_input
def mass_number(isotope: ParticleLike):
    return isotope.mass_number

mass_number("Fe-56")

56

mass_number("e-")

---------------------------------------------------------------------------
InvalidIsotopeError                       Traceback (most recent call last)
Cell In[106], line 1
----> 1 mass_number("e-")

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/plasmapy/particles/decorators.py:990, in particle_input.<locals>.wrapper(callable__, instance, args, kwargs)
    983 @wrapt.decorator
    984 def wrapper(
    985     callable__: Callable[..., Any],
   (...)
    988     kwargs: MutableMapping[str, Any],
    989 ) -> Callable[..., Any]:
--> 990     bound_arguments = particle_validator.process_arguments(args, kwargs, instance)
    991     return callable__(  # type: ignore[no-any-return]
    992         *bound_arguments.args,
    993         **bound_arguments.kwargs,
    994     )

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/plasmapy/particles/decorators.py:673, in _ParticleInput.process_arguments(self, args, kwargs, instance)
    668 mass_numb = bound_arguments.arguments.pop("mass_numb", None)
    670 self.perform_pre_validations(Z, mass_numb)
    672 processed_kwargs = {
--> 673     parameter: self.process_argument(parameter, argument, Z, mass_numb)
    674     for parameter, argument in bound_arguments.arguments.items()
    675 }
    677 for parameter in processed_kwargs:
    678     bound_arguments.arguments[parameter] = processed_kwargs[parameter]

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/plasmapy/particles/decorators.py:597, in _ParticleInput.process_argument(self, parameter, argument, Z, mass_numb)
    595 self.verify_charge_categorization(particle)
    596 self.verify_particle_categorization(particle)
--> 597 self.verify_particle_name_criteria(parameter, particle)
    598 self.verify_allowed_types(particle)
    600 return particle

File ~/checkouts/readthedocs.org/user_builds/plasmapy-summer-school/envs/latest/lib/python3.12/site-packages/plasmapy/particles/decorators.py:492, in _ParticleInput.verify_particle_name_criteria(self, parameter, particle)
    489     meets_name_criteria = all(meets_name_criteria)  # type: ignore[arg-type]
    491 if not meets_name_criteria:
--> 492     raise exception(
    493         f"The argument {parameter} = {particle!r} to "
    494         f"{self.callable_.__name__} does not correspond to a "
    495         f"valid {parameter}."
    496     )

InvalidIsotopeError: The argument isotope = Particle("e-") to mass_number does not correspond to a valid isotope.

Putting it all together

[@validate_quantities]: https://docs.plasmapy.org/en/latest/api/plasmapy.utils.decorators.validators.validate_quantities.html#validate-quantities [@particle_input]: https://docs.plasmapy.org/en/stable/api/plasmapy.particles.decorators.particle_input.html [@check_relativistic]: https://docs.plasmapy.org/en/stable/api/plasmapy.utils.decorators.checks.check_relativistic.html#check-relativistic

Let’s use both [@particle_input], [@validate_quantities], and [@check_relativistic] to write an Alfvén speed function. The formula is

\[V_A ≡ \frac{B}{\sqrt{μ_0 ρ}}\]

where we approximate the mass density as \(ρ ≈ m_i n_i\).

@particle_input
@check_relativistic
@validate_quantities
def alfven_speed(
    B: u.Quantity[u.T],
    n: u.Quantity[u.m**-3],
    ion: ParticleLike,
) -> u.Quantity[u.m / u.s]:
    """Return the Alfvén speed."""
    return B / np.sqrt(constants.mu0 * n * ion.mass)

alfven_speed(B = 0.02 * u.T, n = 1e15 * u.m**-3, ion="p+")

$13795142 \; \mathrm{\frac{m}{s}}$