This page was generated by nbsphinx from docs/notebooks/getting_started/units.ipynb.
Interactive online version: .

Using Astropy Units

In scientific computing, we often represent physical quantities as numbers.

distance_in_miles = 50
time_in_hours = 2
velocity_in_mph = distance_in_miles / time_in_hours
print(velocity_in_mph)

25.0

Representing a physical quantity as a number has risks. We might unknowingly perform operations with different units, like time_in_seconds + time_in_hours. We might even accidentally perform operations with physically incompatible units, like length + time, without catching our mistake. We can avoid these problems by using a units package.

This notebook introduces astropy.units with an emphasis on the functionality needed to work with plasmapy.particles and plasmapy.formulary. We typically import this subpackage as u.

import astropy.units as u

Contents

1. Unit basics

2. Unit operations

3. Unit conversations

4. Detaching units and values

5. Equivalencies

6. Physical constants

7. Units in PlasmaPy

8. Optimizing unit operations

9. Physical Types

Unit basics

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

distance = 60 * u.km
print(distance)

60.0 km

This operation creates a Quantity: a number, sequence, or array that has been assigned a physical unit.

type(distance)

astropy.units.quantity.Quantity

We can also create an object by using the Quantity class itself.

time = u.Quantity(120, u.min)

We can create Quantity objects with compound units.

88 * u.imperial.mile / u.hour

$88 \; \mathrm{\frac{mi}{h}}$

We can even create Quantity objects that are explicitly dimensionless.

3 * u.dimensionless_unscaled

$3 \; \mathrm{}$

We can also create a Quantity based off of a NumPy array or a list.

import numpy as np

np.array([2.5, 3.2, 1.1]) * u.kg

$[2.5,~3.2,~1.1] \; \mathrm{kg}$

[2, 3, 4] * u.m / u.s

$[2,~3,~4] \; \mathrm{\frac{m}{s}}$

Unit operations

Operations between Quantity objects handle unit conversions automatically. We can add Quantity objects together as long as their units have the same physical type.

1 * u.m + 25 * u.cm

$1.25 \; \mathrm{m}$

Units get handled automatically during operations like multiplication, division, and exponentiation.

velocity = distance / time
print(velocity)

0.5 km / min

area = distance**2
print(area)

3600.0 km2

Attempting an operation between physically incompatible units gives us an error, which we can use to find bugs in our code.

3 * u.m + 3 * u.s

UnitConversionError: 's' (time) and 'm' (length) are not convertible


During handling of the above exception, another exception occurred:

UnitConversionError: Can only apply 'add' function to quantities with compatible dimensions

Quantity objects behave very similarly to NumPy arrays because Quantity is a subclass of numpy.ndarray.

balmer_series = [656.279, 486.135, 434.0472, 410.1734] * u.nm
Hα = balmer_series[0]
print(Hα)

656.279 nm

np.max(balmer_series)

$656.279 \; \mathrm{nm}$

Most frequently encountered NumPy and SciPy functions can be used with Quantity objects. However, Quantity objects lose their units with some operations.

Unit conversions

The to method allows us to convert a Quantity to different units of the same physical type. This method accepts strings that represent a unit (including compound units) or a unit object.

velocity.to("m/s")

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

velocity.to(u.m / u.s)

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

The si and cgs attributes convert the Quantity to SI or CGS units, respectively.

velocity.si

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

velocity.cgs

$833.33333 \; \mathrm{\frac{cm}{s}}$

Detaching units and values

The value attribute of a Quantity provides the number (as a NumPy scalar) or NumPy array without the unit.

time.value

120.0

The unit attribute of a Quantity provides the unit without the value.

time.unit

$\mathrm{min}$

Equivalencies

Plasma scientists often use the electron-volt (eV) as a unit of temperature. This is a shortcut for describing the thermal energy per particle, or more accurately the temperature multiplied by the Boltzmann constant, \(k_B\). Because an electron-volt is a unit of energy rather than temperature, we cannot directly convert electron-volts to kelvin.

u.eV.to("K")

UnitConversionError: 'eV' (energy/torque/work) and 'K' (temperature) are not convertible

To handle non-standard unit conversions, astropy.units allows the use of equivalencies. The conversion from eV to K can be done by using the temperature_energy() equivalency.

(1 * u.eV).to("K", equivalencies=u.temperature_energy())

$11604.518 \; \mathrm{K}$

Radians are treated dimensionlessly when the dimensionless_angles() equivalency is in effect. Note that this equivalency does not account for the multiplicative factor of \(2π\) that is used when converting between frequency and angular frequency.

(3.2 * u.rad / u.s).to("1 / s", equivalencies=u.dimensionless_angles())

$3.2 \; \mathrm{\frac{1}{s}}$

Physical constants

We can use astropy.constants to access the most commonly needed physical constants.

from astropy.constants import c, e, k_B

print(c)

  Name   = Speed of light in vacuum
  Value  = 299792458.0
  Uncertainty  = 0.0
  Unit  = m / s
  Reference = CODATA 2018

A Constant behaves very similarly to a Quantity. For example, we can use the Boltzmann constant to mimic the behavior of u.temperature_energy().

thermal_energy_per_particle = 0.6 * u.keV
temperature = thermal_energy_per_particle / k_B
print(temperature.to("MK"))

6.962710872930049 MK

Electromagnetic constants often need the unit system to be specified. Code within PlasmaPy uses SI units.

2 * e

TypeError: Constant 'e' does not have physically compatible units across all systems of units and cannot be combined with other values without specifying a system (eg. e.emu)

2 * e.si

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

Units in PlasmaPy

Now we can show some uses of astropy.units in PlasmaPy, starting with plasmapy.particles. Many of the attributes of Particle and ParticleList provide Quantity objects.

from plasmapy.particles import Particle, ParticleList

alpha = Particle("He-4 2+")

alpha.charge

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

ions = ParticleList(["O 1+", "O 2+", "O 3+"])
ions.mass

$[2.6566054 \times 10^{-26},~2.6565143 \times 10^{-26},~2.6564232 \times 10^{-26}] \; \mathrm{kg}$

Similarly, Quantity objects are the expected inputs and outputs of most functions in plasmapy.formulary. We can use them to calculate some plasma parameters for a typical region of the solar corona.

from plasmapy.formulary import Alfven_speed, Debye_length, gyrofrequency

B = 0.01 * u.T
n = 1e15 * u.m**-3
proton = Particle("p+")

Alfven_speed(B=B, density=n, ion=proton).to("km /s")

$6895.6934 \; \mathrm{\frac{km}{s}}$

gyrofrequency(B=B, particle="e-")

$1.75882 \times 10^{9} \; \mathrm{\frac{rad}{s}}$

The to_hz keyword provides the frequency in hertz rather than radians per second, and accounts for the factor of \(2π\).

gyrofrequency(B=B, particle="e-", to_hz=True)

$2.799249 \times 10^{8} \; \mathrm{Hz}$

Formulary functions perform calculations based on SI units, but accept input arguments in other units. Temperature can be given in units of temperature (e.g., kelvin) or energy (e.g., electron-volts).

Debye_length(T_e=1e6 * u.K, n_e=1e9 * u.m**-3)

$2.1822556 \; \mathrm{m}$

Debye_length(T_e=86.17 * u.eV, n_e=1e3 * u.cm**-3)

$2.1822134 \; \mathrm{m}$

Optimizing unit operations

Astropy’s documentation includes performance tips for using astropy.units in computationally intensive situations. For example, putting compound units in parentheses reduces the need to make multiple copies of the data.

volume = 0.62 * (u.barn * u.Mpc)

Physical types

A physical type corresponds to physical quantities with dimensionally compatible units. Astropy has functionality that represents different physical types. These physical type objects can be accessed using either the physical_type attribute of a unit or get_physical_type().

(u.m**2 / u.s).physical_type

PhysicalType({'diffusivity', 'kinematic viscosity'})

u.get_physical_type("number density")

PhysicalType('number density')

These physical type objects can be used for dimensional analysis.

energy_density = (u.J * u.m**-3).physical_type
velocity = u.get_physical_type("velocity")
print(energy_density * velocity)

energy flux/irradiance