03-UnitsQuantities/Astropy_Units.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Astropy Units, Quantities, and Constants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Astropy includes a powerful framework for units that allows users to attach units to scalars and arrays.  These quantities can be manipulated or combined, while keeping track of the units.\n",
    "\n",
    "For more information about the features presented below, please see the\n",
    "[astropy.units](https://docs.astropy.org/en/stable/units/index.html) docs.\n",
    "\n",
    "Also note that this tutorial assumes you have little or no knowledge of the astropy units docs.  If you're moderately familiar with them and are interested in some more complex examples, you might instead prefer the Astropy tutorial on [\"Using Astropy Quantities for astrophysical calculations\"](http://learn.astropy.org/Quantities.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Table of contents:\n",
    "\n",
    "* <a href=\"#Representing-units\">Representing units</a>\n",
    "* <a href=\"#Composite-units\">Composite units</a>\n",
    "* <a href=\"#Quantity-objects\">Quantity objects</a>\n",
    "* <a href=\"#Quantity-attributes\">Quantity attributes</a>\n",
    "* <a href=\"#Quantity-arithmetic-operations\">Quantity arithmetic operations</a>\n",
    "* <a href=\"#Combining-Quantities\">Combining Quantities</a>\n",
    "* <a href=\"#Converting-units\">Converting units</a>\n",
    "* <a href=\"#Decomposing-units\">Decomposing units</a>\n",
    "* <a href=\"#Integration-with-Numpy-functions\">Integration with NumPy functions</a>\n",
    "* <a href=\"#Defining-new-units\">Defining new units</a>\n",
    "* <a href=\"#Using-physical-constants\">Using physical constants</a>\n",
    "* <a href=\"#Equivalencies\">Equivalencies</a>\n",
    "* <a href=\"#Performance-consideration\">Performance consideration</a>\n",
    "* <a href=\"#Putting-it-all-together:--a-concise-example\">Putting it all together: a concise example</a>\n",
    "* <a href=\"#Exercises\">Exercises</a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Representing units"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we need to import the astropy units subpackage (**`astropy.units`**).  Because we probably want to use units in many expressions, it is most concise to rename the subpackage as **`u`**.  This is the standard convention, but note that this will conflict with any variable called **`u`**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import astropy.units as u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will also import numpy here, which we use later.\n",
    "# \"np\" is the common naming standard for this import.\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Units can then be accessed as **`u.<unit>`**.  For example, the meter unit is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Units have docstrings, which give some explanatory text about them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.m.__doc__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.pc.__doc__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and a physical type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.m.physical_type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.s.physical_type"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Many units also have aliases:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.m.aliases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.meter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.arcsec.aliases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.arcsecond"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SI and cgs units are available by default, but Imperial units require the **`imperial`** prefix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "raises-exception"
    ]
   },
   "outputs": [],
   "source": [
    "# This is undefined.\n",
    "u.inch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use this instead.\n",
    "u.imperial.inch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please see the complete list of [available units](https://astropy.readthedocs.org/en/stable/units/index.html#module-astropy.units.si)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Composite units"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Composite units are created using Python numeric operators (e.g., \"`*`\" (multiplication), \"`/`\" (division), and \"`**`\" (power))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.km / u.s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.imperial.mile / u.h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(u.eV * u.Mpc) / u.Gyr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.cm**3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.m / u.kg / u.s**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ``Quantity`` objects\n",
    "The most useful feature of units is the ability to attach them to scalars or arrays, creating `Quantity` objects. A `Quantity` object contains both a value and a unit.  The most convenient way to create a `Quantity` object is by multiplying the value with its unit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "3.7 * u.au  # Quantity object"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An equivalent (but more verbose) way of doing the same thing is to use the `Quantity` object's initializer, as demonstrated below.  In general, the more concise 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](https://docs.astropy.org/en/stable/api/astropy.units.quantity.Quantity.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.Quantity(3.7, unit=u.au)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.Quantity(1 * u.cm, unit=u.m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Where quantities really shine is when you make an array `Quantity` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We want to create the composite unit first, hence the parenthesis.\n",
    "x = [1.2, 6.8, 3.7] * (u.pc / u.year)\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## `Quantity` attributes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The units and value of a `Quantity` can be accessed separately via the ``value`` and ``unit`` attributes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = 5. * u.Mpc\n",
    "q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q.value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q.unit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [1.2, 6.8, 3.7] * (u.pc / u.year)\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x.value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x.unit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## `Quantity` arithmetic operations\n",
    "\n",
    "\"`*`\" (multiplication), \"`/`\" (division), and \"`**`\" (power) operations can be performed on `Quantity` objects with `float`/`int` values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = 3.1 * u.km"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q * 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q / 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "1 / q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q ** 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combining Quantities\n",
    "\n",
    "Quantities can be combined using Python numeric operators:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q1 = 3. * (u.m / u.s)\n",
    "q1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q2 = 5. * (u.cm / u.s / u.g**2)\n",
    "q2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q1 * q2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q1 / q2  # Note the \"second\" unit cancelled out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sometimes, one more step needed for \"clean\" output unit.\n",
    "# Also see: Decomposing units (below)\n",
    "(q1 / q2).to(u.g ** 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q1 ** 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [1.2, 6.8, 3.7] * (u.pc / u.year)\n",
    "x * 3  # Element-wise multiplication."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When adding or subtracting quantities, the units must be **compatible** (but not necessarily identical)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add two quantities.\n",
    "(3 * u.m) + (5 * u.m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we add two distance quantities that do not have identical units:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(3 * u.km) + (5 * u.cm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "raises-exception"
    ]
   },
   "outputs": [],
   "source": [
    "# This will fail because the units are incompatible.\n",
    "(3 * u.km) + (5. * u.km / u.s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Converting units\n",
    "\n",
    "Units can be converted to other equivalent units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = 2.5 * u.year\n",
    "q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert year to seconds.\n",
    "q.to(u.s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert degree squared to steradian.\n",
    "(7. * u.deg**2).to(u.sr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert mph to kph.\n",
    "(55. * (u.imperial.mile / u.h)).to(u.km / u.h)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Important Note**:  Converting a unit (not a `Quantity`) gives only the scale factor:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.Msun.to(u.kg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To keep the units, use a `Quantity` (value and unit) object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(1. * u.Msun).to(u.kg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decomposing units"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The units of a `Quantity` object can be decomposed into a set of base units using the\n",
    "``decompose()`` method. By default, units will be decomposed to SI unit bases:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = 8. * (u.cm * u.pc / u.g / u.year**2)\n",
    "q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q.decompose()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To decompose into cgs unit bases:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q.decompose(u.cgs.bases)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.cgs.bases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.si.bases"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Units will not cancel out unless they are identical:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = 7 * u.m / (7 * u.km)\n",
    "q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But they will cancel by using the `decompose()` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = q.decompose()\n",
    "x  # This is a \"dimensionless\" Quantity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "repr(x.unit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integration with NumPy functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most [NumPy](https://www.numpy.org) functions understand `Quantity` objects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sin(30)  # np.sin assumes the input is in radians."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sin(30 * u.degree)  # Awesome!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = 100 * (u.kg * u.kg)\n",
    "np.sqrt(q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(10) * u.km\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.mean(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some NumPy \"ufuncs\" (universal functions) require dimensionless quantities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "raises-exception"
    ]
   },
   "outputs": [],
   "source": [
    "np.log10(4 * u.m)  # This is not possible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.log10(4 * u.m / (4 * u.km))  # Note that the units cancelled."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Care needs to be taken with dimensionless units.\n",
    "\n",
    "For example, passing ordinary values to an inverse trigonometric function gives a result without units:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.arcsin(1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`u.dimensionless_unscaled` creates a ``Quantity`` with a \"dimensionless unit\" and therefore gives a result *with* units:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.arcsin(1.0 * u.dimensionless_unscaled)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.arcsin(1.0 * u.dimensionless_unscaled).to(u.degree)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Important Note:**  In-place array operations do not work with units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "raises-exception"
    ]
   },
   "outputs": [],
   "source": [
    "a = np.arange(10.)\n",
    "a *= 1.0 * u.kg  # In-place operator will fail."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assign to a *new* array instead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = a * 1.0 * u.kg\n",
    "a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also, Quantities lose their units with some NumPy operations, for example:\n",
    "\n",
    "* np.append\n",
    "* np.dot\n",
    "* np.hstack\n",
    "* np.vstack\n",
    "* np.where\n",
    "* np.choose\n",
    "* np.vectorize\n",
    "\n",
    "See [Quantity Known Issues](https://docs.astropy.org/en/stable/known_issues.html#quantities-lose-their-units-with-some-operations) for more details."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining new units"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also define custom units for something that isn't built in to Astropy.\n",
    "\n",
    "Let's define a unit called **\"sol\"** that represents a Martian day."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sol = u.def_unit('sol', 1.0274912510 * u.day)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(1. * u.yr).to(sol)  # 1 Earth year in Martian sol units."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's define Mark Watney's favorite unit, the [**Pirate-Ninja**](https://en.wikipedia.org/wiki/List_of_humorous_units_of_measurement#Pirate_Ninja):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pirate_ninja = u.def_unit('☠️👤', 1.0 * (u.kW * u.hr / sol))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "5.2 * pirate_ninja"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Mars oxygenator power requirement for 6 people.\n",
    "(44.1 * pirate_ninja).to(u.W)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using physical constants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [astropy.constants](https://docs.astropy.org/en/stable/constants/index.html) module contains physical constants relevant for astronomy.  They are defined as ``Quantity`` objects using the ``astropy.units`` framework.\n",
    "\n",
    "Please see the complete list of [available physical constants](https://docs.astropy.org/en/stable/constants/index.html#module-astropy.constants).  Additions are welcome!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy import constants as const"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "const.G"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "const.c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "const.L_sun"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Constants are Quantities, thus they can be coverted to other units:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "const.L_sun.to(u.erg / u.s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also note that even constants are not constant. `astropy.constants` supports [different versions](https://docs.astropy.org/en/stable/constants/index.html#collections-of-constants-and-prior-versions)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import astropy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "astropy.physical_constants.get()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "astropy.astronomical_constants.get()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.constants import iau2012\n",
    "iau2012.L_sun"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "raises-exception"
    ]
   },
   "outputs": [],
   "source": [
    "# Changing the entire science state like this after\n",
    "# imports is currently not possible, see\n",
    "# https://github.com/astropy/astropy/issues/8781\n",
    "astropy.astronomical_constants.set('iau2012')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equivalencies\n",
    "\n",
    "Equivalencies can be used to convert quantities that are not strictly the same physical type, but in a specific context are interchangable.  A familiar physics example is the mass-energy equivalency: these are strictly different physical types, but it is often understood that you can convert between the two using $E=mc^2$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.constants import m_p  # Proton mass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "raises-exception"
    ]
   },
   "outputs": [],
   "source": [
    "# This raises an error because mass and energy are different units.\n",
    "(m_p).to(u.eV)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This succeeds using equivalencies.\n",
    "(m_p).to(u.MeV, u.mass_energy())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This concept extends further in `astropy.units` to include some common practical astronomy situations where the units have no direct physical connection, but it is often useful to have a \"quick shorthand.\"  For example, astronomical spectra are often given as a function of wavelength, frequency, or even energy of the photon.  Suppose you want to find the [Lyman-limit](https://en.wikipedia.org/wiki/Lyman_limit) wavelength:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "raises-exception"
    ]
   },
   "outputs": [],
   "source": [
    "# This raises an error.\n",
    "(13.6 * u.eV).to(u.Angstrom)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Normally, you can convert `u.eV` only to the following units:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.eV.find_equivalent_units()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But by using a spectral equivalency, you can also convert `u.eV` to the following units:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u.eV.find_equivalent_units(equivalencies=u.spectral())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now back to Lyman-limit.\n",
    "(13.6 * u.eV).to(u.Angstrom, equivalencies=u.spectral())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or if you remember the 21cm HI line, but cannot remember the frequency, you could do:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(21. * u.cm).to(u.GHz, equivalencies=u.spectral())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To go one step further, the units of a spectrum's *flux* are further complicated by being dependent on the units of the spectrum's \"X-axis\" (i.e., $f_{\\lambda}$ for flux per unit wavelength or $f_{\\nu}$ for flux per unit frequency).  `astropy.units` supports this use case, but it is necessary to supply the location in the spectrum where the conversion is done:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = 1e-18 * (u.erg / u.s / u.cm**2 / u.AA)\n",
    "q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q.to(u.uJy, equivalencies=u.spectral_density(1. * u.um))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There's a lot of flexibility with equivalencies, including a variety of other useful built-in equivalencies.  So if you want to know more, you might want to check out the [equivalencies narrative documentation](https://docs.astropy.org/en/stable/units/equivalencies.html) or the [astropy.units.equivalencies reference docs](https://docs.astropy.org/en/stable/units/index.html#module-astropy.units.equivalencies)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Domain-specific units"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some astropy modules define more specific units and equivalencies.\n",
    "\n",
    "For example, ``astropy.cosmology`` defines units for cosmological redshift and\n",
    "the Hubble parameter and equivalencies for both."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import astropy.cosmology.units as cu\n",
    "\n",
    "cu.__all__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = 1100 * cu.redshift\n",
    "repr(z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Redshift isn't a \"real\" unit (so the ``dimensionless_redshift`` equivalency\n",
    "is turned on by default) but it is very useful for converting between cosmological\n",
    "distance measures, for example the background temperature. We can do this\n",
    "in a specific cosmological context:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.cosmology import Planck18  # cosmological context\n",
    "\n",
    "z.to(u.K, equivalencies=cu.redshift_temperature(Planck18))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to know more about cosmology and units, you might want to check out the [astropy.cosmology.units  reference docs](https://docs.astropy.org/en/stable/cosmology/units.html#cosmological-units-and-equivalencies)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Performance consideration\n",
    "\n",
    "For more details about performance on big datasets, see [units performance tips](https://docs.astropy.org/en/stable/units/index.html#performance-tips).\n",
    "\n",
    "When working with a big data array, the `<<` operator can be used to speed up Quantity initialization. Let's consider this array:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "big_array = np.random.random(10_000_000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defining the Quantity the simplest way without parenthesis.\n",
    "%timeit big_array * u.m / u.s / u.kg / u.sr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defining the Quantity with the * operator and parenthesis.\n",
    "%timeit big_array * (u.m / u.s / u.kg / u.sr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Defining the Quantity with the << operator (parenthesis not needed).\n",
    "%timeit big_array << u.m / u.s / u.kg / u.sr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Replace the ``Quantity`` values below with your slow/fast timing to see how much things sped up on your machine. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "((89.4 * u.ms) / (67.1 * u.us)).to(u.dimensionless_unscaled)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Putting it all together:  a concise example\n",
    "\n",
    "Let's estimate the (circular) orbital speed of the Earth around the Sun using Kepler's Law:\n",
    "\n",
    "$$v = \\sqrt{\\frac{G M_{\\odot}}{r}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v = np.sqrt(const.G * 1 * u.M_sun / (1 * u.au))\n",
    "v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's a velocity unit... but it sure isn't obvious when you look at it!\n",
    "\n",
    "Let's use a variety of the available ``Quantity`` methods to get something more sensible:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v.decompose()  # Remember that the default uses SI bases"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v.decompose(u.cgs.bases)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "v.to(u.km / u.s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises\n",
    "\n",
    "## Exercise 1\n",
    "\n",
    "The *James Webb Space Telescope (JWST)* is in a halo orbit around the second Sun-Earth Lagrange (L2) point:\n",
    "\n",
    "&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;☀️ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; 🌎 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; L2 &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; *(not to scale)*\n",
    "\n",
    "\n",
    "L2 is located at a distance from the Earth (opposite the Sun) of approximately:\n",
    "\n",
    "$$ r \\approx R \\left(\\frac{M_{earth}}{3 M_{sun}}\\right) ^{(1/3)} $$\n",
    "\n",
    "where $R$ is the Sun-Earth distance.\n",
    "\n",
    "Calculate the Earth-L2 distance in kilometers and miles.\n",
    "\n",
    "*Hints*:\n",
    "\n",
    "* $M_{earth}$ and $M_{sun}$ are defined [constants](https://docs.astropy.org/en/stable/constants/#reference-api) \n",
    "\n",
    "* The mile unit is defined as ``u.imperial.mile`` (see [imperial units](https://docs.astropy.org/en/v0.2.1/units/index.html#module-astropy.units.imperial))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Answer here (km)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Answer here (mile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2\n",
    "\n",
    "The L2 point is about 1.5 million kilometers away from the Earth opposite the Sun.\n",
    "The total mass of the *James Webb Space Telescope (JWST)* is about 6500 kg.\n",
    "\n",
    "Using the value you obtained above for the Earth-L2 distance, calculate the gravitational force in Newtons between: \n",
    "\n",
    "* *JWST* (at L2) and the Earth\n",
    "* *JWST* (at L2) and the Sun\n",
    "\n",
    "*Hint*: the gravitational force between two masses separated by a distance *r* is:\n",
    "\n",
    "$$ F_g = \\frac{G m_1 m_2}{r^2} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Answer here (JWST and Earth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Answer here (JWST and Sun)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3\n",
    "\n",
    "Calculate the [Schwarzschild radius](https://en.wikipedia.org/wiki/Schwarzschild_radius) in units of solar radii of the Sgr A*, the Milky Way's supermassive black hole with $M = 4.31 \\times 10^6 M_\\odot$, given\n",
    "\n",
    "$$r_\\mathrm{s} = \\frac{2 G M}{c^2}$$\n",
    "\n",
    "Also calculate the angular size of the event horizon on the sky in microarcseconds, given the distance to the galactic center $d_{center} = 7.94$ kpc, given\n",
    "\n",
    "$$\\theta = \\mathrm{arctan}\\frac{2 r_\\mathrm{s}}{d_{center}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Answer radius here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Answer angular size here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4\n",
    "\n",
    "We can make a very rough estimate of the temperature of material in the vicinity of Sgr A* by assuming hydrostatic equilibrium, so that the thermal energy of the gas balances the gravitational force:\n",
    "\n",
    "$$ kT \\sim GM m_p /R $$\n",
    "\n",
    "where $m_p$ is the mass of a proton and $R$ is the distance away from the black hole. Using Astropy constants and the properties of Sgr A* described in Example 3, compute the temperature of the gas required to balance the black hole's gravitation pull at 1 million Schwarzschild radii derived above. Use this equation:\n",
    "\n",
    "$$ T = \\frac{G M m_p}{10^6 r_s k} $$\n",
    "\n",
    "Then use the [Astropy units temperature equivalencies](https://docs.astropy.org/en/stable/units/equivalencies.html#temperature-equivalency) to find the energy of that gas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Answer temperature here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Answer energy here"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "vscode": {
   "interpreter": {
    "hash": "69ffed3fe34942d91c870dd5e7b1cdc8962d4c36821e73c483e75fd1432e7bb3"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}