In this tutorial we present some examples showing how astropy's `Quantity`

object can make astrophysics calculations easier. The examples include calculating the mass of a galaxy from its velocity dispersion and determining masses of molecular clouds from CO intensity maps. We end with an example of good practices for using quantities in functions you might distribute to other people.

For an in-depth discussion of `Quantity`

objects, see the astropy documentation section.

We start by loading standard libraries and set up plotting for ipython notebooks.

In [1]:

```
#!/usr/bin/env python
from __future__ import print_function, division
import math
import numpy as np
import matplotlib.pyplot as plt
# You shouldn't use the `seed` function in real science code, but we use it here for example purposes.
# It makes the "random" number generator always give the same numbers wherever you run it.
np.random.seed(12345)
# Set up matplotlib and use a nicer set of plot parameters
%config InlineBackend.rc = {}
import matplotlib
matplotlib.rc_file("../../templates/matplotlibrc")
import matplotlib.pyplot as plt
%matplotlib inline
```

It is conventional to load the astropy `units`

module as the variable `u`

, demonstrated below. This will make working with `Quantity`

objects much easier.

Astropy also has a `constants`

module, where typical physical constants are available. The constants are stored as objects of a subclass of `Quantity`

, so they behave just like a `Quantity`

. Here, we'll only need the gravitational constant `G`

, Planck's constant `h`

, and Boltzmann's constant, `k_B`

.

In [2]:

```
import astropy.units as u
from astropy.constants import G, h, k_B
```

`Quantity`

objects to estimate a hypothetical galaxy's mass, given its half-light radius and radial velocities of stars in the galaxy.

`Quantity`

object with the name `Reff`

. The easiest way to create a `Quantity`

object is just by multiplying the value with its unit. Units are accessed as u."unit", in this case u.pc.

In [3]:

```
Reff = 29 * u.pc
```

`Quantity`

object's initializer, demonstrated below. In general, the simpler form (above) is preferred, as it is closer to how such a quantity would actually be written in text. The initalizer form has more options, though, which you can learn about from the astropy reference documentation on Quantity.

In [4]:

```
Reff = u.Quantity(29, unit=u.pc)
```

We can access the value and unit of a `Quantity`

using the `value`

and `unit`

attributes.

In [5]:

```
print("""Half light radius
value: {0}
unit: {1}""".format(Reff.value, Reff.unit))
```

The `value`

and `unit`

attributes can also be accessed within the print function.

In [6]:

```
print("""Half light radius
value: {0.value}
unit: {0.unit}""".format(Reff))
```

`to()`

method. Here, we convert it to meters.

In [7]:

```
print("{0:.3g}".format(Reff.to(u.m)))
```

In [8]:

```
vmean = 206
sigin = 4.3
v = np.random.normal(vmean, sigin, 500)*u.km/u.s
```

In [9]:

```
print("""First 10 radial velocity measurements:
{0}
{1}""".format(v[:10], v.to(u.m/u.s)[:10]))
```

In [10]:

```
plt.figure()
plt.hist(v, bins=20, histtype="step")
plt.xlabel("Velocity (km/s)")
plt.ylabel("N")
```

Out[10]:

`Quantity`

objects, and also use them in standard numpy functions such as `mean()`

and `size()`

. They retain their units through these operations just as you would expect them to.

In [11]:

```
sigma = np.sqrt(np.sum((v - np.mean(v))**2) / np.size(v))
print("Velocity dispersion: {0:.2f}".format(sigma))
```

`numpy`

square root function, because the resulting velocity dispersion quantity is a `numpy`

array. If we used the python standard `math`

library's `sqrt`

function instead, we get an error.

In [12]:

```
sigma_scalar = math.sqrt(np.sum((v - np.mean(v))**2) / np.size(v))
```

`numpy`

functions with `Quantity`

objects, *not* the `math`

equivalents, unless you are sure you understand the consequences.

`Quantity`

, just accept that this is often good enough. For the calculation we can just multiply the quantities together, and `astropy`

will keep track of the units.

In [13]:

```
M = 4*sigma**2*Reff/G
M
```

Out[13]:

In [14]:

```
M.decompose()
```

Out[14]:

In [15]:

```
print("""Galaxy mass
in solar units: {0:.3g}
SI units: {1:.3g}
CGS units: {2:.3g}""".format(M.to(u.Msun), M.si, M.cgs))
```

`np.log10`

as long as the logarithm's argument is dimensionless.

In [16]:

```
np.log10(M / u.Msun)
```

Out[16]:

However, you can't take the log of something with units, as that is not mathematically sensible.

In [17]:

```
np.log10(M)
```

Use `Quantity`

and Kepler's law in the form given below to determine the (circular) orbital speed of the Earth around the sun in km/s. You should not have to look up an constants or conversion factors to do this calculation - it's all in `astropy.units`

and `astropy.constants`

.

In [18]:

```
```

In [18]:

```
```

In [18]:

```
```

`Quantity`

objects can facilitate a full derivation of the total mass of a molecular cloud using radio observations of isotopes of Carbon Monoxide (CO).

In [18]:

```
d = 250 * u.pc
Tex = 25 * u.K
```

`numpy.meshgrid`

function to create data cubes for each of the three coordinates, and then use them in the formula for a Gaussian to generate an array with the synthetic data cube. In this cube, the cloud is positioned at the center of the cube, with $\sigma$ and the center in each dimension shown below. Note in particular that the $\sigma$ for RA and Dec have different units from the center, but `astropy`

automatically does the relevant conversions before computing the exponential.

In [19]:

```
# Cloud's center
cen_ra = 52.25 * u.deg
cen_dec = 0.25 * u.deg
cen_v = 15 * u.km/u.s
# Cloud's size
sig_ra = 3 * u.arcmin
sig_dec = 4 * u.arcmin
sig_v = 3 * u.km/u.s
#1D coordinate quantities
ra = np.linspace(52, 52.5, 100) * u.deg
dec = np.linspace(0, 0.5, 100) * u.deg
v = np.linspace(0, 30, 300) *u.km/u.s
#this creates data cubes of size for each coordinate based on the dimensions of the other coordinates
ra_cube, dec_cube, v_cube = np.meshgrid(ra, dec, v)
data_gauss = np.exp(-0.5*((ra_cube-cen_ra)/sig_ra)**2 +
-0.5*((dec_cube-cen_dec)/sig_dec)**2 +
-0.5*((v_cube-cen_v)/sig_v)**2 )
```

In [20]:

```
data = data_gauss * u.K
```

In [21]:

```
# Average pixel size
# This is only right if dec ~ 0, because of the cos(dec) factor.
dra = (ra.max() - ra.min()) / len(ra)
ddec = (dec.max() - dec.min()) / len(dec)
#Average velocity bin width
dv = (v.max() - v.min()) / len(v)
print("""dra = {0}
ddec = {1}
dv = {2}""".format(dra.to(u.arcsec), ddec.to(u.arcsec), dv))
```

In [22]:

```
intcloud = np.sum(data*dv, axis=2)
intcloud.unit
```

Out[22]:

In [23]:

```
#Note that we display RA in the convential way by going from max to min
plt.imshow(intcloud.value,
origin='lower',
extent=[ra.value.max(), ra.value.min(), dec.value.min(), dec.value.max()],
cmap='hot',
interpolation='nearest',
aspect='equal')
plt.colorbar().set_label("Intensity ({})".format(intcloud.unit))
plt.xlabel("RA (deg)")
plt.ylabel("Dec (deg)");
```

In order to calculate the mass of the molecular cloud, we need to measure its column density. A number of assumptions are required for the following calculation; the most important are that the emission is optically thin (typically true for ${\rm C}^{18}{\rm O}$) and that conditions of local thermodynamic equilibrium hold along the line of sight. In the case where the temperature is large compared to the separation in energy levels for a molecule and the source fills the main beam of the telescope, the total column density for ${\rm C}^{13}{\rm O}$ is

$N=C \frac{\int T_B(V) dV}{1-e^{-B}}$

where the constants $C$ and $B$ are given by:

$C=3.0\times10^{14} \left(\frac{\nu}{\nu_{13}}\right)^2 \frac{A_{13}}{A} {\rm K^{-1} cm^{-2} \, km^{-1} \, s}$

$B=\frac{h\nu}{k_B T}$

(Rohlfs & Wilson "Tools of Radio Astronomy").

Here we have given an expression for $C$ scaled to the values for ${\rm C}^{13}{\rm O}$ ($\nu_{13}$ and $A_{13}$). In order to use this relation for ${\rm C}^{18}{\rm O}$, we need to rescale the frequencies ${\nu}$ and Einstein coefficients $A$. $C$ is in funny mixed units, but that's okay. We'll define it as a `Quantities`

object and not have to worry about it.

First, we look up the wavelength for these emission lines and store them as quantities.

In [24]:

```
lambda13 = 2.60076 * u.mm
lambda18 = 2.73079 * u.mm
```

In [25]:

```
nu13 = lambda13.to(u.Hz)
```

Fortunately, `astropy`

comes to the rescue by providing a feature called "unit equivalencies". Equivalencies provide a way to convert between two physically different units that are not normally equivalent, but in a certain context have a one-to-one mapping. For more on equivalencies, see the equivalencies section of astropy's documentation.

In this case, calling the `astropy.units.spectral()`

function provides the equivalencies necessary to handle conversions between wavelength and frequency. To use it, provide the equivalencies to the `equivalencies`

keyword of the `to()`

call:

In [26]:

```
nu13 = lambda13.to(u.Hz, equivalencies=u.spectral())
nu18 = lambda18.to(u.Hz, equivalencies=u.spectral())
```

In [27]:

```
A13 = 7.4e-8 / u.s
A18 = 8.8e-8 / u.s
C = 3e14 * (nu18/nu13)**3 * (A13/A18) / (u.K * u.cm**2 * u.km *(1/u.s))
C
```

Out[27]:

`astropy.constants`

, and the other two values are already calculated, so here we just take the ratio.

In [28]:

```
B = h * nu18 / (k_B * Tex)
```

In [29]:

```
print('{0}\n{1}'.format(B, B.decompose()))
```

In [30]:

```
NCO = C * np.sum(data*dv, axis=2) / (1 - np.exp(-B))
print("Peak CO column density: ")
np.max(NCO)
```

Out[30]:

In [31]:

```
H2_CO_ratio = 5.9e6
NH2 = NCO * H2_CO_ratio
print("Peak H2 column density: ")
np.max(NH2)
```

Out[31]:

In [32]:

```
mH2 = 2 * 1.008 * u.Dalton #aka atomic mass unit/amu
rho = NH2 * mH2
```

In [33]:

```
dap = dra * ddec
print(dap)
```

*physical* units, rather than angular units. So it is tempting to just multiply the area and the square of the distance.

In [34]:

```
da = dap * d**2 # don't actually do it this way - use the version below instead!
print(da)
```

In [35]:

```
dap.to(u.steradian).value * d**2
```

Out[35]:

**wrong**, because `astropy.units`

treats angles (and solid angles) as actual physical units, while the small-angle approximation assumes angles are dimensionless. So if you, e.g., try to convert to a different area unit, it will fail:

In [36]:

```
da.to(u.cm**2)
```

`dimensionless_angles`

equivalency, which allows angles to be treated as dimensionless. This makes it so that they will automatically convert to radians and become dimensionless when a conversion is needed.

In [37]:

```
da = (dap * d**2).to(u.pc**2, equivalencies=u.dimensionless_angles())
da
```

Out[37]:

In [38]:

```
da.to(u.cm**2)
```

Out[38]:

In [39]:

```
M = np.sum(rho * da)
M.decompose().to(u.solMass)
```

Out[39]:

`Quantity`

's array capabililities. Compute the median and mean of the `data`

with the `np.mean`

and `np.median`

functions. Why are their values so different?

In [40]:

```
```

In [40]:

```
```

`Quantity`

is also a useful tool if you plan to share some of your code, either with collaborators or the wider community. By writing functions that take `Quantity`

objects instead of raw numbers or arrays, you can write code that is agnostic to the input unit. In this way, you may even be able to prevent the destruction of Mars orbiters. Below, we provide a simple example.

In [40]:

```
def response_func(xinarcsec, yinarcsec):
xscale = 0.9
yscale = 0.85
xfactor = 1 / (1 + xinarcsec/xscale)
yfactor = 1 / (1 + yinarcsec/yscale)
return xfactor * yfactor
```

In [41]:

```
response_func(1.0, 1.2)
```

Out[41]:

`Quantity`

objects. The new function could simply be:

In [42]:

```
def response_func(x, y):
xscale = 0.9 * u.arcsec
yscale = 0.85 * u.arcsec
xfactor = 1 / (1 + x/xscale)
yfactor = 1 / (1 + y/yscale)
return xfactor * yfactor
```

In [43]:

```
response_func(1.0, 1.2)
```

Which is their cue to provide the units explicitly:

In [44]:

```
response_func(1.0*u.arcmin, 1.2*u.arcmin)
```

Out[44]:

`Quantity`

as the input of code you shared.

`Quantity`

input and outputs, of course), but allowing for an arbitrary mass and orbital radius. Try it with some reasonable numbers for satellites orbiting the Earth, a moon of Jupiter, or an extrasolar planet. Feel free to use wikipedia or similar for the masses and distances.

In [45]:

```
```