Interactive online version: .

# Using Astropy Units

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

[2]:

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.

[3]:

import astropy.units as u


## Unit basics

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

[4]:

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.

[5]:

type(distance)

[5]:

astropy.units.quantity.Quantity


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

[6]:

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


We can create Quantity objects with compound units.

[7]:

88 * u.imperial.mile / u.hour

[7]:

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

We can even create Quantity objects that are explicitly dimensionless.

[8]:

3 * u.dimensionless_unscaled

[8]:

$3 \; \mathrm{}$

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

[9]:

import numpy as np

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

[9]:

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

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

[10]:

$[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.

[11]:

1 * u.m + 25 * u.cm

[11]:

$1.25 \; \mathrm{m}$

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

[12]:

velocity = distance / time
print(velocity)

0.5 km / min

[13]:

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.

[14]:

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.

[15]:

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

656.279 nm

[16]:

np.max(balmer_series)

[16]:

$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.

[17]:

velocity.to("m/s")

[17]:

$8.3333333 \; \mathrm{\frac{m}{s}}$
[18]:

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

[18]:

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

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

[19]:

velocity.si

[19]:

$8.3333333 \; \mathrm{\frac{m}{s}}$
[20]:

velocity.cgs

[20]:

$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.

[21]:

time.value

[21]:

120.0


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

[22]:

time.unit

[22]:

$\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.

[23]:

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.

[24]:

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

[24]:

$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.

[25]:

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

[25]:

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

## Physical constants

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

[26]:

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().

[27]:

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.

[28]:

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)


[29]:

2 * e.si

[29]:

$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.

[30]:

from plasmapy.particles import Particle, ParticleList

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

[31]:

alpha.charge

[31]:

$3.2043533 \times 10^{-19} \; \mathrm{C}$
[32]:

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

[32]:

$[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.

[33]:

from plasmapy.formulary import Alfven_speed, Debye_length, gyrofrequency

[34]:

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

[35]:

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

[35]:

$6895.6934 \; \mathrm{\frac{km}{s}}$
[36]:

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

[36]:

$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π$$.

[37]:

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

[37]:

$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).

[38]:

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

[38]:

$2.1822556 \; \mathrm{m}$
[39]:

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

[39]:

$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.

[40]:

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().

[41]:

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

[41]:

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

[42]:

u.get_physical_type("number density")

[42]:

PhysicalType('number density')


These physical type objects can be used for dimensional analysis.

[43]:

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

energy flux/irradiance