Skip to content

Instantly share code, notes, and snippets.

@eteq

eteq/quantities.ipynb

Created February 8, 2018 18:50
Show Gist options
  • Save eteq/abba0705394c38effb3ab30663ae9531 to your computer and use it in GitHub Desktop.
Save eteq/abba0705394c38effb3ab30663ae9531 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
quantities.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Using Astropy Quantities for astrophysical calculations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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.\n",
"\n",
"For an in-depth discussion of `Quantity` objects, see the [astropy documentation section](http://docs.astropy.org/en/stable/units/quantity.html)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Preliminaries"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We start by loading standard libraries and set up plotting for ipython notebooks."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"# You shouldn't use the `seed` function in real science code, but we use it here for example purposes.\n",
"# It makes the \"random\" number generator always give the same numbers wherever you run it.\n",
"np.random.seed(12345)\n",
"\n",
"# Set up matplotlib\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It is conventional to load the astropy `units` module as the variable `u`, demonstrated below. This will make working with `Quantity` objects much easier. \n",
"\n",
"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`."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import astropy.units as u\n",
"from astropy.constants import G, h, k_B"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 1. Galaxy mass"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this first example, we will use `Quantity` objects to estimate a hypothetical galaxy's mass, given its half-light radius and radial velocities of stars in the galaxy."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Lets assume that we measured the half light radius of the galaxy to be 29 pc projected on the sky at the distance of the galaxy. This radius is often called the \"effective radius\", so we will store it as a `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."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"Reff = 29 * u.pc"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A completely equivalent (but more verbose) way of doing the same thing is to use the `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](http://docs.astropy.org/en/stable/api/astropy.units.quantity.Quantity.html)."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"Reff = u.Quantity(29, unit=u.pc)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can access the value and unit of a `Quantity` using the `value` and `unit` attributes."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Half light radius\n",
"value: 29.0\n",
"unit: pc\n"
]
}
],
"source": [
"print(\"\"\"Half light radius\n",
"value: {0}\n",
"unit: {1}\"\"\".format(Reff.value, Reff.unit))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The `value` and `unit` attributes can also be accessed within the print function."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Half light radius\n",
"value: 29.0\n",
"unit: pc\n"
]
}
],
"source": [
"print(\"\"\"Half light radius\n",
"value: {0.value}\n",
"unit: {0.unit}\"\"\".format(Reff))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Furthermore, we can convert the radius in parsecs to any other unit of length using the ``to()`` method. Here, we convert it to meters."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"8.95e+17 m\n"
]
}
],
"source": [
"print(\"{0:.3g}\".format(Reff.to(u.m)))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we will first create a synthetic dataset of radial velocity measurements, assuming a normal distribution with a mean velocity of 206 km/s and a velocity dispersion of 4.3 km/s."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"vmean = 206\n",
"sigin = 4.3\n",
"v = np.random.normal(vmean, sigin, 500)*u.km/u.s"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"First 10 radial velocity measurements: \n",
"[ 205.11975706 208.05945635 203.76641353 203.61035969 214.45285646\n",
" 211.99164508 206.39950387 207.21150846 209.30679704 211.35966937] km / s\n",
"[ 205119.75706422 208059.45635365 203766.41352526 203610.35969131\n",
" 214452.85646176 211991.64508178 206399.50387 207211.50845717\n",
" 209306.79704073 211359.66936646] m / s\n"
]
}
],
"source": [
"print(\"\"\"First 10 radial velocity measurements: \n",
"{0}\n",
"{1}\"\"\".format(v[:10], v.to(u.m/u.s)[:10]))"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0,0.5,'N')"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAFBlJREFUeJzt3X2wJXV95/H3JzyYSSI1jgzs7OBk\nsJZFzW548IbFJesqKKvCCptIZNU4EbJTWxV3sUytjg9lVZLd2iFWYh431qxoxioFKaJCRFHCQ6hY\nCAxPAo4EJKwBRmbiOghkFkS/+0f3DYfJvXPnPvQ9c+/v/ao6dbp/p/v29zfnnvuZ7j7961QVkqR2\n/di4C5AkjZdBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWrcweMuYH8cfvjhtX79\n+nGXIUlLyq233vp3VbV6puWWRBCsX7+ebdu2jbsMSVpSkvyf/VnOQ0OS1DiDQJIaZxBIUuMMAklq\nnEEgSY0zCCSpcQaBJDXOIJCkxhkEktS4JXFlsTSkUzZfy8O798xp3bUrV/DVTacucEXS4jII1LyH\nd+/hwc1nzGnd9ZuuXOBqpMXnoSFJapxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhrndQRaMPO5\nMAu8OEsaF4NAC2Y+F2aBF2dJ4+KhIUlq3KBBkGRlksuSfDPJ9iSvSLIqydVJ7uufXzBkDZKkfRt6\nj+D3gauq6iXAccB2YBNwTVUdA1zTz0uSxmSwIEhyGPBK4CKAqnq6qnYDZwFb+8W2AmcPVYMkaWZD\n7hG8GNgFfCLJ7Uk+luQngSOragdA/3zEgDVIkmYwZBAcDJwI/ElVnQA8ySwOAyXZmGRbkm27du0a\nqkZJat6QQfAQ8FBV3dTPX0YXDI8mWQPQP++cauWq2lJVE1U1sXr16gHLlKS2DRYEVfUd4G+THNs3\nnQZ8A7gC2NC3bQAuH6oGSdLMhr6g7L8An0pyKPAA8A668Lk0yfnAt4FzBq5BDZjv7Sallg0aBFV1\nBzAxxUunDbldtWe+VzVLLfPKYklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGucdyvQc47ww\na+3KFXO+S5kXhUlzZxDoOcZ5YZb3K5bGw0NDktQ4g0CSGmcQSFLjDAJJapxBIEmNMwgkqXEGgSQ1\nziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxgw5DneRB4HHgh8AzVTWRZBXwGWA98CDw\nS1X1vSHrkCRNbzH2CF5dVcdX1UQ/vwm4pqqOAa7p5yVJYzKOQ0NnAVv76a3A2WOoQZLUGzoICvhK\nkluTbOzbjqyqHQD98xED1yBJ2oehb1V5SlU9kuQI4Ook39zfFfvg2Aiwbt26oeqTpOYNukdQVY/0\nzzuBzwEnAY8mWQPQP++cZt0tVTVRVROrV68eskxJatpgQZDkJ5M8f3IaOB24G7gC2NAvtgG4fKga\nJEkzG/LQ0JHA55JMbufTVXVVkluAS5OcD3wbOGfAGiRJMxgsCKrqAeC4Kdq/C5w21HYlSbPjlcWS\n1DiDQJIaZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJapxBIEmN\nMwgkqXEGgSQ1ziCQpMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNe7goTeQ5CBgG/Bw\nVZ2Z5GjgEmAVcBvwy1X19NB1SENYu3IF6zddOed1v7rp1AWuSJq9wYMAuADYDhzWz18IfKSqLkny\nUeB84E8WoQ5pwc3nD/lcA0RaaIMeGkpyFHAG8LF+PsCpwGX9IluBs4esQZK0b0OfI/g94D3Aj/r5\nFwK7q+qZfv4hYO1UKybZmGRbkm27du0auExJatdgQZDkTGBnVd062jzFojXV+lW1paomqmpi9erV\ng9QoSRr2HMEpwBuTvAH4cbpzBL8HrExycL9XcBTwyIA1SJJmMNgeQVW9r6qOqqr1wLnAtVX1VuA6\n4E39YhuAy4eqQZI0s3FcR/Be4N1J7qc7Z3DRGGqQJPUW4+ujVNX1wPX99APASYuxXUnSzLyyWJIa\nZxBIUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjFmWICS2uUzZfy8O798xp\n3bUrVyxwNZIOdAbBMvTw7j08uPmMcZchaYnYZxAk+dA+Xq6q+q0FrkeStMhm2iN4coq2nwB+lW4I\naYNAkpa4fQZBVf3O5HSS5wMXAOcBlwC/M916kqSlY8ZzBElWAe8G3gpsBU6squ8NXZgkaXHMdI7g\nw8AvAFuAf1lVTyxKVZKkRTPTHsGvA08BHwQ+kGSyPXQniw8bsDZpWVu7cgXrN105r/W/uunUBaxI\nrZrpHIEXnEkDme8f8fmEiDTKP/SS1DiDQJIaZxBIUuMMAklq3GBBkOTHk9yc5M4k9yT5jb796CQ3\nJbkvyWeSHDpUDZKkmQ25R/AUcGpVHQccD7wuycnAhcBHquoY4HvA+QPWIEmawWBBUJ3JC9AO6R8F\nnApc1rdvBc4eqgZJ0swGPUeQ5KAkdwA7gauBbwG7q+qZfpGHgLXTrLsxybYk23bt2jVkmZLUtEGD\noKp+WFXHA0cBJwEvnWqxadbdUlUTVTWxevXqIcuUpKYtyreGqmo3cD1wMrAyyeQVzUcBjyxGDZKk\nqQ35raHVSVb20yuA1wDbgeuAN/WLbQAuH6oGSdLMhrxV5Rpga5KD6ALn0qr6QpJvAJck+e/A7cBF\nA9YgSZrBYEFQVV8HTpii/QG68wWSpAOAVxZLUuMMAklqnEEgSY0zCCSpcUN+a0jSgOZzq0tvc6lR\nBoG0RM3nD7m3udQoDw1JUuMMAklqnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGmcQSFLjDAJJ\napxBIEmNMwgkqXEGgSQ1ziCQpMYZBJLUOINAkho3WBAkeVGS65JsT3JPkgv69lVJrk5yX//8gqFq\nkCTNbMg9gmeAX6+qlwInA7+W5GXAJuCaqjoGuKaflySNyWBBUFU7quq2fvpxYDuwFjgL2NovthU4\ne6gaJEkzW5RzBEnWAycANwFHVtUO6MICOGIxapAkTW3wm9cn+Sngz4B3VdX3k+zvehuBjQDr1q0b\nrsAD1Cmbr+Xh3XvmtO7alSsWuBpJy9mgQZDkELoQ+FRVfbZvfjTJmqrakWQNsHOqdatqC7AFYGJi\nooas80D08O49PLj5jHGXIakBQ35rKMBFwPaq+t2Rl64ANvTTG4DLh6pBkjSzIfcITgF+GbgryR19\n2/uBzcClSc4Hvg2cM2ANkqQZDBYEVfVXwHQnBE4baruSpNnxymJJapxBIEmNMwgkqXEGgSQ1ziCQ\npMYZBJLUOINAkhpnEEhS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTGGQSS1DiDQJIaZxBIUuMMAklq\nnEEgSY0zCCSpcQaBJDXOIJCkxhkEktQ4g0CSGjdYECT5eJKdSe4eaVuV5Ook9/XPLxhq+5Kk/XPw\ngD/7T4E/Aj450rYJuKaqNifZ1M+/d8AaJC2wUzZfy8O798xp3bUrV/DVTacucEWar8GCoKpuSLJ+\nr+azgFf101uB6zEIpCXl4d17eHDzGXNad/2mKxe4Gi2ExT5HcGRV7QDon49Y5O1LkvYy5KGheUmy\nEdgIsG7dujFXM3vz2X2GbhdakhbDYgfBo0nWVNWOJGuAndMtWFVbgC0AExMTtVgFLpT57D5L0mJa\n7ENDVwAb+ukNwOWLvH1J0l6G/ProxcCNwLFJHkpyPrAZeG2S+4DX9vOSpDEa8ltD/3Gal04bapuS\npNnzymJJatwB+60hScNZu3LFnL/T7zfalh+DQGqQV/dqlIeGJKlxy36PwHFRJGnfln0QOC6KJO2b\nh4YkqXEGgSQ1ziCQpMYZBJLUOINAkhq37L81NB9efSmpBQbBPngNgaQWeGhIkhpnEEhS4wwCSWqc\nQSBJjTMIJKlxBoEkNc4gkKTGeR2BpCZ4b5LpGQSSmuC9SabnoSFJatxYgiDJ65Lcm+T+JJvGUYMk\nqbPoh4aSHAT8MfBa4CHgliRXVNU3FrsWSYtrPgM5LsS2x2EpnJsYxzmCk4D7q+oBgCSXAGcBBoG0\nzC3nE67TWQrnJsZxaGgt8Lcj8w/1bZKkMRjHHkGmaKt/tFCyEdjYzz6R5N45b/DC58weDvzdXH/W\nAcx+LS3LsV/LsU/Q92uvvyOzMq51gZ/en4XGEQQPAS8amT8KeGTvhapqC7BloTeeZFtVTSz0zx03\n+7W0LMd+Lcc+wfLt16hxHBq6BTgmydFJDgXOBa4YQx2SJMawR1BVzyR5J/Bl4CDg41V1z2LXIUnq\njOXK4qr6IvDFcWybAQ43HSDs19KyHPu1HPsEy7df/yBV/+g8rSSpIQ4xIUmNW3ZBkOTjSXYmuXuk\n7bgkNya5K8mfJzmsb1+fZE+SO/rHR8dX+fSSvCjJdUm2J7knyQV9+6okVye5r39+Qd+eJH/QD+Hx\n9SQnjrcHU5tDv16V5LGR9+tD4+3B1PbRr3P6+R8lmdhrnff179e9Sf7deCrft9n2axl8vj6c5Jv9\nZ+hzSVaOrHPAv1+zUlXL6gG8EjgRuHuk7Rbg3/bT5wG/1U+vH13uQH0Aa4AT++nnA38NvAz4bWBT\n374JuLCffgPwJbprNk4Gbhp3HxaoX68CvjDuuufRr5cCxwLXAxMjy78MuBN4HnA08C3goHH3YwH6\ntdQ/X6cDB/ftF478Hi6J92s2j2W3R1BVNwD/d6/mY4Eb+umrgV9c1KLmqap2VNVt/fTjwHa6q7HP\nArb2i20Fzu6nzwI+WZ2vASuTrFnksmc0h34tCdP1q6q2V9VUF0aeBVxSVU9V1d8A99MNxXJAmUO/\nloR99OsrVfVMv9jX6K55giXyfs3GsguCadwNvLGfPofnXtB2dJLbk/xlkn+z+KXNTpL1wAnATcCR\nVbUDul9m4Ih+sSU3jMd+9gvgFUnuTPKlJD+z6IXO0l79ms5Sf7/2ZSl/vkadR7eXDUvw/ZpJK0Fw\nHvBrSW6l2/V7um/fAayrqhOAdwOfnjx/cCBK8lPAnwHvqqrv72vRKdoO2K+HzaJftwE/XVXHAX8I\nfH4x6psr36/l8flK8gHgGeBTk01TrH7Avl/7o4kgqKpvVtXpVfVy4GK6Y3r0u3bf7adv7dv/+fgq\nnV6SQ+h+ST9VVZ/tmx+dPOTTP+/s2/drGI8DwWz6VVXfr6on+ukvAockOXwMZc9omn5NZ6m/X1Na\nBp8vkmwAzgTeWv0JApbQ+7W/mgiCJEf0zz8GfBD4aD+/Ot39EUjyYuAY4IFx1TmdJAEuArZX1e+O\nvHQFsKGf3gBcPtL+9v7bQycDj00eajmQzLZfSf5Jvw5JTqL7/f3u4lW8f/bRr+lcAZyb5HlJjqb7\nPbx5yBrnYrb9WuqfrySvA94LvLGq/n5klSXxfs3KuM9WL/SD7n/8O4Af0CX3+cAFdN8E+GtgM89e\nSPeLwD103wC4Dfj3465/mj79PN2u59eBO/rHG4AXAtcA9/XPq/rlQ3fzn28BdzHyTY4D6TGHfr1z\n5P36GvCvx92HWfbrP/S/k08BjwJfHlnnA/37dS/w+nH3YSH6tQw+X/fTnQuYbPvoUnq/ZvPwymJJ\nalwTh4YkSdMzCCSpcQaBJDXOIJCkxhkEktQ4g0BLRpLr9x7pMcm7kvyvGdZ7Yo7b+80krxnZzk/M\ncv0kuTbJYf1InHfPvNY+f94rkvzvaV47NMkNScZysyktbQaBlpKL6e5xPercvn3BVdWHquov+tl3\nAbMKArrvot9Z+x6GYTZeB1w11QtV9TTdNRdvXqBtqSEGgZaSy4AzkzwP/mGAsH8K/FU//9+S3NKP\nH/8be6/c/w/9w0nuTndvijePvPaevu3OJJv7tj9N8qYk/7XfznXpxq0/P8lHRtb9T0mmutL2rTx7\ntfdoHS/uB2L7uSS/kuTz6e6T8TdJ3pnk3f3rX0uyamTV04C/SPIzSW5ON8b/15Mc07/++X6b0qwY\nBFoyqhu35ma6/xlDtzfwmaqqJKfTXep/EnA88PIkr9zrR/xC/9pxwGuADydZk+T1dENd/6vqBrT7\n7b22+wd0Y8m8uqpeDVwCvLEfnwbgHcAnpij5FODW0YYkx9KNafOOqrqlb/4XwFv62v8H8PfVDdR2\nI/D2fr3DgR9U1WPAfwZ+v6qOByboruqFbpTdn5vu30+ajkGgpWb08NDoYaHT+8ftdMMZvIQuGEb9\nPHBxVf2wqh4F/pLuD+drgE9UP55MVe19P4vnqKongWvp9k5eAhxSVXdNseiq6sa3n7Sabg/hbVV1\nx0j7dVX1eFXtAh4D/rxvv4vu5i6T/ftKP30j8P4k76UbjXVPX9cPgaeTPH9f9Ut7Mwi01HweOC3d\n7TdXVH9DEbrxlf5nVR3fP/5ZVV2017pTDR882T7bsVY+BvwK0+8NADzTD3Q46TG6sWtO2Wu5p0am\nfzQy/yNg8uTv6+nPD1TVp+nur7EH+HKSU0fWfx7w/2bTEckg0JJS3TDU1wMf57knib8MnNePKU+S\ntZOjzo64AXhzkoOSrKa7renNdP/TPm/yW0F7HZef9DjdvSwm67iJbijitzD9yep7gRePzD9Ndwjq\n7UneMnNvO/3omD9LN/DZ5EieD/SHrK7oXyPJC4FdVfWD/f3ZEjz7vw1pKbkY+Cwj3yCqqq8keSlw\nYz9S9RPA23j2Hg0AnwNeQTcaZgHvqarvAFclOR7YluRp4IvA+/fa5hbgS0l29OcJAC4Fjq+q701T\n55V091m+f6TOJ5OcCVyd5Mn97O/Lgdvr2REi3wy8LckPgO8Av9m3v7qvXZoVRx+V5ijJF4CPVNU1\n07y+hu7e0a+d53Y+CNxfVZfMsNxngffVEr5/sMbDIJBmKclKukNKd1bVOTMs+0vAVQt4LcF02zkU\nOLeqPjnkdrQ8GQSS1DhPFktS4wwCSWqcQSBJjTMIJKlxBoEkNc4gkKTG/X920S2bbiQg/AAAAABJ\nRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x106fbbe10>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure()\n",
"plt.hist(v, bins='auto', histtype=\"step\")\n",
"plt.xlabel(\"Velocity (km/s)\")\n",
"plt.ylabel(\"N\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we calculate the velocity dispersion of the galaxy. This demonstrates how you can perform basic operations like subtraction and division with `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."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Velocity dispersion: 4.36 km / s\n"
]
}
],
"source": [
"sigma = np.sqrt(np.sum((v - np.mean(v))**2) / np.size(v))\n",
"print(\"Velocity dispersion: {0:.2f}\".format(sigma))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note how we needed to use `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."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"sigma_scalar = np.sqrt(np.sum((v - np.mean(v))**2) / len(v))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In general, you should only use `numpy` functions with `Quantity` objects, *not* the `math` equivalents, unless you are sure you understand the consequences."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now for the actual mass calculation. If a galaxy is pressure-supported (for example, an elliptical or dwarf spheroidal galaxy), its mass within the stellar extent can be estimated using a straightforward formula: $M_{1/2}=4\\sigma^2 R_{eff}/G$. There are caveats to the use of this formula for science - see Wolf et al. 2010 for details. For demonstrating `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."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$3.3085932 \\times 10^{13} \\; \\mathrm{\\frac{km^{2}\\,kg\\,pc}{m^{3}}}$"
],
"text/plain": [
"<Quantity 33085931653008.55 kg km2 pc / m3>"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M = 4*sigma**2*Reff/G\n",
"M"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The result is in a composite unit, so it's not really obvious it's a mass. However, it can be decomposed to cancel all of the length units ($km^2 pc/m^3$) using the decompose() method."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$1.0209252 \\times 10^{36} \\; \\mathrm{kg}$"
],
"text/plain": [
"<Quantity 1.0209251756364422e+36 kg>"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M.decompose()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can also easily express the mass in whatever form you like - solar masses are common in astronomy, or maybe you want the default SI and CGS units."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Galaxy mass\n",
"in solar units: 5.13e+05 solMass\n",
"SI units: 1.02e+36 kg\n",
"CGS units: 1.02e+39 g\n"
]
}
],
"source": [
"print(\"\"\"Galaxy mass\n",
"in solar units: {0:.3g}\n",
"SI units: {1:.3g}\n",
"CGS units: {2:.3g}\"\"\".format(M.to(u.Msun), M.si, M.cgs))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Or, if you want the log of the mass, you can just use ``np.log10`` as long as the logarithm's argument is dimensionless."
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$5.7104737 \\; \\mathrm{}$"
],
"text/plain": [
"<Quantity 5.710473687565016>"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.log10(M / u.Msun)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"However, you can't take the log of something with units, as that is not mathematically sensible."
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"tags": [
"raises-exception"
]
},
"outputs": [
{
"ename": "UnitTypeError",
"evalue": "Can only apply 'log10' function to dimensionless quantities",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mUnitConversionError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/quantity_helper.py\u001b[0m in \u001b[0;36mget_converter\u001b[0;34m(from_unit, to_unit)\u001b[0m\n\u001b[1;32m 21\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 22\u001b[0;31m \u001b[0mscale\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfrom_unit\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_to\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mto_unit\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 23\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mUnitsError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/core.py\u001b[0m in \u001b[0;36m_to\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m 943\u001b[0m raise UnitConversionError(\n\u001b[0;32m--> 944\u001b[0;31m \"'{0!r}' is not a scaled version of '{1!r}'\".format(self, other))\n\u001b[0m\u001b[1;32m 945\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mUnitConversionError\u001b[0m: 'Unit(\"kg km2 pc / m3\")' is not a scaled version of 'Unit(dimensionless)'",
"\nDuring handling of the above exception, another exception occurred:\n",
"\u001b[0;31mUnitConversionError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/quantity_helper.py\u001b[0m in \u001b[0;36mhelper_dimensionless_to_dimensionless\u001b[0;34m(f, unit)\u001b[0m\n\u001b[1;32m 107\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 108\u001b[0;31m return ([get_converter(unit, dimensionless_unscaled)],\n\u001b[0m\u001b[1;32m 109\u001b[0m dimensionless_unscaled)\n",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/quantity_helper.py\u001b[0m in \u001b[0;36mget_converter\u001b[0;34m(from_unit, to_unit)\u001b[0m\n\u001b[1;32m 24\u001b[0m return from_unit._apply_equivalencies(\n\u001b[0;32m---> 25\u001b[0;31m from_unit, to_unit, get_current_unit_registry().equivalencies)\n\u001b[0m\u001b[1;32m 26\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mAttributeError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/core.py\u001b[0m in \u001b[0;36m_apply_equivalencies\u001b[0;34m(self, unit, other, equivalencies)\u001b[0m\n\u001b[1;32m 880\u001b[0m \"{0} and {1} are not convertible\".format(\n\u001b[0;32m--> 881\u001b[0;31m unit_str, other_str))\n\u001b[0m\u001b[1;32m 882\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mUnitConversionError\u001b[0m: 'kg km2 pc / m3' (mass) and '' (dimensionless) are not convertible",
"\nDuring handling of the above exception, another exception occurred:\n",
"\u001b[0;31mUnitTypeError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-17-598955917a11>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlog10\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mM\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/quantity.py\u001b[0m in \u001b[0;36m__array_ufunc__\u001b[0;34m(self, function, method, *inputs, **kwargs)\u001b[0m\n\u001b[1;32m 619\u001b[0m \u001b[0;31m# consistent units between two inputs (e.g., in np.add) --\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 620\u001b[0m \u001b[0;31m# and the unit of the result (or tuple of units for nout > 1).\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 621\u001b[0;31m \u001b[0mconverters\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0munit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mconverters_and_unit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunction\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0minputs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 622\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 623\u001b[0m \u001b[0mout\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'out'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/quantity_helper.py\u001b[0m in \u001b[0;36mconverters_and_unit\u001b[0;34m(function, method, *args)\u001b[0m\n\u001b[1;32m 459\u001b[0m \u001b[0;31m# will have (this is a tuple of units if there are multiple outputs).\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 460\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mfunction\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mUFUNC_HELPERS\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 461\u001b[0;31m \u001b[0mconverters\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mresult_unit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mUFUNC_HELPERS\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mfunction\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunction\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0munits\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 462\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 463\u001b[0m raise TypeError(\"Unknown ufunc {0}. Please raise issue on \"\n",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/quantity_helper.py\u001b[0m in \u001b[0;36mhelper_dimensionless_to_dimensionless\u001b[0;34m(f, unit)\u001b[0m\n\u001b[1;32m 111\u001b[0m raise UnitTypeError(\"Can only apply '{0}' function to \"\n\u001b[1;32m 112\u001b[0m \u001b[0;34m\"dimensionless quantities\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 113\u001b[0;31m .format(f.__name__))\n\u001b[0m\u001b[1;32m 114\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 115\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mUnitTypeError\u001b[0m: Can only apply 'log10' function to dimensionless quantities"
]
}
],
"source": [
"np.log10(M)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercises"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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`.\n",
"\n",
"$$v = \\sqrt{\\frac{G M_{\\odot}}{r}}$$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There's a much easier way to figure out the velocity of the Earth using just two units or quantities. Do that and then compare to the Kepler's law answer (the easiest way is probably to compute the percentage difference, if any)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(Completely optional, but a good way to convince yourself of the value of Quantity:) Do the above calculations by hand - you can use a calculator (or python just for its arithmatic) but look up all the appropriate conversion factors and use paper-and-pencil approaches for keeping track of them all. Which one took longer?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2. Molecular cloud mass"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this second example, we will demonstrate how using `Quantity` objects can facilitate a full derivation of the total mass of a molecular cloud using radio observations of isotopes of Carbon Monoxide (CO)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Setting up the data cube"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's assume that we have mapped the inner part of a molecular cloud in the J=1-0 rotational transition of ${\\rm C}^{18}{\\rm O}$ and are interested in measuring its total mass. The measurement produced a data cube with RA and Dec as spatial coordiates and velocity as the third axis. Each voxel in this data cube represents the brightness temperature of the emission at that position and velocity. Furthermore, we will assume that we have an independent measurement of distance to the cloud $d=250$ pc and that the excitation temperature is known and constant throughout the cloud: $T_{ex}=25$ K."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [],
"source": [
"d = 250 * u.pc\n",
"Tex = 25 * u.K"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will generate a synthetic dataset, assuming the cloud follows a Gaussian distribution in each of RA, Dec and velocity. We start by creating a 100x100x300 numpy array, such that the first coordinate is right ascension, the second is declination, and the third is velocity. We use the `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."
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"# Cloud's center\n",
"cen_ra = 52.25 * u.deg\n",
"cen_dec = 0.25 * u.deg\n",
"cen_v = 15 * u.km/u.s\n",
"\n",
"# Cloud's size\n",
"sig_ra = 3 * u.arcmin\n",
"sig_dec = 4 * u.arcmin\n",
"sig_v = 3 * u.km/u.s\n",
"\n",
"#1D coordinate quantities\n",
"ra = np.linspace(52, 52.5, 100) * u.deg\n",
"dec = np.linspace(0, 0.5, 100) * u.deg\n",
"v = np.linspace(0, 30, 300) *u.km/u.s\n",
"\n",
"#this creates data cubes of size for each coordinate based on the dimensions of the other coordinates\n",
"ra_cube, dec_cube, v_cube = np.meshgrid(ra, dec, v)\n",
"\n",
"data_gauss = np.exp(-0.5*((ra_cube-cen_ra)/sig_ra)**2 + \n",
" -0.5*((dec_cube-cen_dec)/sig_dec)**2 + \n",
" -0.5*((v_cube-cen_v)/sig_v)**2 )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The units of the exponential are dimensionless, so we multiply the data cube by K to get brightness temperature units. Radio astronomers use a rather odd set of units [K km/s] as of integrated intensity (that is, summing all the emission from a line over velocity). As an aside for experts, we're setting up our artificial cube on the main-beam temperature scale (T$_{\\rm MB}$) which is the closest we can normally get to the actual brightness temperature of our source."
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"data = data_gauss * u.K"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will also need to know the width of each velocity bin and the size of each pixel, so we calculate that now."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"dra = 18.0 arcsec\n",
"ddec = 18.0 arcsec\n",
"dv = 0.1 km / s\n"
]
}
],
"source": [
"# Average pixel size\n",
"# This is only right if dec ~ 0, because of the cos(dec) factor.\n",
"dra = (ra.max() - ra.min()) / len(ra)\n",
"ddec = (dec.max() - dec.min()) / len(dec)\n",
"\n",
"#Average velocity bin width\n",
"dv = (v.max() - v.min()) / len(v)\n",
"print(\"\"\"dra = {0}\n",
"ddec = {1}\n",
"dv = {2}\"\"\".format(dra.to(u.arcsec), ddec.to(u.arcsec), dv))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We are interested in the integrated intensity over all of the velocity channels, so we will create a 2D quantity array by summing our data cube along the velocity axis (multiplying by the velocity width of a pixel)."
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\mathrm{\\frac{K\\,km}{s}}$"
],
"text/plain": [
"Unit(\"K km / s\")"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"intcloud = np.sum(data*dv, axis=2)\n",
"intcloud.unit"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can plot the 2D quantity using matplotlib's imshow function, by passing the quantity's value. Similarly, we can set the correct extent using the values of $x_i$ and $x_f$. Finally, we can set the colorbar label to have proper units."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAUMAAAEKCAYAAACIZDejAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJztnXu0JUV97z9fzswwMDASGEQErjxE\nvV6vohLQYIyKZo2ECzEhN6AxPkgwUXybBCNRQx5LzboqGpY6IRgTHygYzEhQIN4g6lVkeIg8FYmR\nAeIwIA95zDDM7/7RXXPq1Onq7r13995nz/l91tqrd3VXVdfus0/tb9XvV7+SmeE4jrPY2WHSDXAc\nx1kIeGfoOI6Dd4aO4ziAd4aO4ziAd4aO4ziAd4aO4zhAz52hpNWSbpJ0s6RTKq6/WtKdkq4uX7/X\nZ3scx3FyLOmrYkkzwBnAS4D1wOWS1prZ9UnWz5vZyX21w3Ecpw19KsPDgJvN7BYz2wycDRzb4/0c\nx3GGpjdlCOwD3Bql1wOHV+T7TUnPB34AvNXMbk0zSDoJOKlMPtsnOh2nX7bCRjPbc9jyq1evto0b\nN7bKe8UVV1xoZquHvVdX9NkZquJcuvbvy8DnzGyTpD8APgW8aF4hszXAGoAZyZZ33VLHcebwIPzn\nKOU3btzIunXrWuWVtGqUe3VFnyJrPbBflN4XuD3OYGZ3mdmmMvl3wLN7bI/jOGPDgC0tXwuDPpXh\n5cDBkg4AbgOOB14eZ5C0t5ndUSaPAW7osT2O44wNAx6edCMGorfO0My2SDoZuBCYAc4ys+sknQas\nM7O1wJskHUPx83A38Oq+2uM4zjgJynB66FMZYmYXABck594dvX8n8M4+2+A4ziTwztBxHAfvDB3H\ncQDvDB3HcbbhnaHjOIuercCmxlwLCe8MHcfpgekbJvvKNsdxeqIbp2tJT44iW10t6T5Jb+m6ta4M\nHcfpge6UoZndBBwC26Jh3Qac10nlEd4ZOkMx01O9j/ZUrzNuehsmHwn8yMxGWjtdhXeGjuP0wFYG\nWI63SlIc1WFNGZyliuOBz43SshzeGTq19KUAB72fK8ZppLUy3GhmhzZlkrSMIoZBL6vWvDN0HKcH\nehkmvxS40sx+2nXF4J2hw/jV3zDUtdFV40Kkl87wBHoaIoN3ho7j9EK3naGknSn2U3pdZ5UmeGfo\nOE4PdNsZmtmDwB6dVViBd4aLkC6GxX15628dokz6eXzYvBDw4K6O4zhM43I87wwXCYOowUmu0Wxz\n7yb1GH9WV4mTwpi2p++doeM4PeDK0Fkg9KUEx+WGU6cp0vbWKUWfT5wk3hk6jrPoGWg53oLAO8Pt\njDbKrUkJtqmj70ANgzhZD6MUXSH2jQ+THcdx8M7QmRhNSq1ODdaVzV0bt8U5VnJN84Bx23Iq0RXi\nOPDO0HGcRY8rQ8dxHLwzdMbOMMPYXJmq87l6FkKkm5yxJR76NhlXfLjcF25NdhzHKXFl6PRMnSob\nRMml56rKtsnT1KY6csaPYQI2VCnFtvX7Er6u8WGy4zgO3hk6vTKIIkzztpkPnMmcb1tfrt6UWJXl\n5vvSubw27jJtSOurqsvnEbvAO0PHcZyS6fo58c5wCuhLEeaU4EzmfVXZKprmD9sEYQiKrU6l5dRd\nlfJsUpyuELvGrcmO4zj4MNnplEF8CHOKsC5vqvKqyrStdxBr8tLofaq6cmqsqk3DKLZh5iRdIQ7D\n9HWGkwxq7DjOdkvoDNu8mpG0m6RzJd0o6QZJz+26xb0qQ0mrgdMpflzPNLP3ZfIdB5wD/KKZreuz\nTdNAl4qwTu01KcS6egaxLteRqq2Q3pqkq/wA07YMYnnOKcS4vCvEUelUGZ4OfNXMjpO0DNi5y8qh\nx85Q0gxwBsVep+uByyWtNbPrk3y7Am8CLuurLY7jjJvuDCiSVgLPB14NYGabgc2dVB7R5zD5MOBm\nM7ulbPzZwLEV+f4C+ADTZnpyHKeGgYbJqySti14nJZUdCNwJfFLSVZLOlLSi6xb3OUzeB7g1Sq8H\nDo8zSHomsJ+ZnS/pHbmKyodzEoB6aOi0Ube0rs2wNr22rCZvmzy5NrWJQN00TN6cpKvaMIpEaLOE\nzxmGgQwoG83s0JrrS4BnAW80s8sknQ6cAvzZaG2cS5/KsKrfsm0XpR2ADwFvb6rIzNaY2aFmdqh3\nho4zLXRmQFkPrDezMJV2LkXn2Cl9KsP1wH5Rel/g9ii9K/A04BJJAI8D1ko6ZjEaUYZd3tZWEcb1\nL03OhXSqAqvOtTWotCVnMAnpKmNGlVqMy47Sjhg3pIxCd641ZvZfkm6V9GQzuwk4Eri+qdyg9NkZ\nXg4cLOkA4DbgeODl4aKZ3QusCmlJlwDvWIwdoeNsf3TuZ/hG4DOlJfkW4DVdVg49doZmtkXSycCF\nFD+mZ5nZdZJOA9aZ2dq+7j1NDLPUbhDXl1QFVp1bmhyrlGF6rU55pmUDdXOGubnCKnUWzj3ScL/c\nubY0LdnzecY6ul2OZ2ZXA3XziiPTq5+hmV0AXJCce3cm7wv6bIvjOONmun4ifDneAqSNIsyVySnC\neAlcqvJSZRjnzeWpU55pmwJ1yvCRJL00OR9bjHOqdBirctW/ay6s2CiBZxcf07cczztDx3F6wDtD\npyWjbKhUt7QuPdapveXJtTRdVT7kqVOGgwR3zSnDh5PzcZseSfIM8izrwn3VnWvCLcwp3hk6juOU\neGfoOM6ix4O7OkNSF5GmzsE5ZzjJDW9hdpi8Y3ItPV93LTWsdDVMDkaQZUk6/rfalNwn/ZcbxJ1m\nkGGtG1IGwYfJjuM4BTZdM6jeGY6ZNpP9uTxVsQmb3GOWJUeYrwjDcUWSrsubc8KO25lzUYn/RVIn\n69Q4Eo5x+8O5B2imrXqrcssZZhmeG1Iipkw6e2foOE73GFP3i+Cd4YQZJAhDmznDcEzdZ+J5wFQB\n7pyk40Bx6bWcQoxdXwbZHe+R5JgqwnB8MCqTex5tVPfW5FgVwiu9T6BuX5YpE0H9Y8xfL7nA8c7Q\ncZzucWXo5GhSLXXBV+sCqTZZkdM5v/h9UH0rkvSuUd4mFVkV3CFWiVXEgiE3V/hgko7rbAobVvU/\nWDVfGaerLNyj7P/sc4dMnVz2ztBxnO5xZei0ZRBfuDTAapvlcjlLMcxXgCuSdJUyzKnIqiV8g1iT\nc4pweZKuq7/Ns2zaJa+KzS3yBNz3sALvDB3HWfQYE/llkPRY4Ajg8cBDwLUU8VMbW+OdYc+0nSts\nYyGumitrWnGSHmFW3YXjbuWxShmm53LKMK4/F9YrXW0CeatxqP/+pM643rT+1FIc3zO1HueCysbn\nZiqu5ciJoEU7d2j0sJlnHkkvpNgkanfgKmADxdfy14GDJJ0L/B8zuy9Xh3eGjuP0w3iV4VHA75vZ\nT9ILkpYAR1Ps4f7FXAXeGTqO0z1jNqCY2R/VXNsCfKmpDu8Mx8wwhpM6g0HqZJ0bJu/MLOnQNxzD\ncPkxUd7cEDo3XAZQzvel/OeIZ2/S4XFYYheGx3UGmqTaynTTPsxhyF71THND6zrckBIxmTnDNwOf\npPgKnQk8EzjFzC5qKtvnvsmO4yxWgjJs82qBpB9L+r6kqyXV7aD52nJe8FeBPSl20Xtfm3u4MlyA\nNBlO2rjWBGfr1G0mfr9bctw9ScOsSkzzLktvUGVByWyCosiCslMpDXcqJeGum+ZWV2WMGcQZOlV3\ny5PzS5PrVWVSpeiqrwX9DJNfaGYbG/KoPB4FfNLMvqdyY/YmvDN0HKd7Jrc2+QpJFwEHAO+UtCst\nf7+8M+yBQcJ0Vc0D5nbHq5qKaztXGDtdp2oxqL+gDPeI8oZzK8PNc5OIVev9ctvjxdFYw/tyknBZ\neXzsPWVVZZmqvZzTubw06EP8Pl32F47LknRc/0xyLQ0QUUWTi01dnu2O9h90VTL0XWNma5I8Blwk\nyYBPVFwPnAgcAtxiZg9K2oOWG857Z+g4TvcM5nS90cyaNog/wsxuL52qL5Z0o5ldOu+2hXP1lVH6\nLuCuNo3wznBMNFmqqtRkm53umpytU6svzJ8H3CM57hnlXbFjcjGdWKySnoMow2A+Dubke+bWsbL8\nGs9sYh65cFyxr2+63C9NVz3T1AKdOk63ccZ2qzKdSmAzu708bpB0HnAYMK8zHAW3JjuO0z1BGbZ5\nNSBpRTn3h6QVFJbia7tusivDCVO3EVQuiGlVuKzU+loXliudK0znDFfEMnLP5Jgqw9TxEObvGhUI\nki1WealjYSYCxIo7Z4s8+uDc6tJjLDxzQWOXJ9fjpubmCFNFGP/tFs08YFu6XY63F3BeaRReAnzW\nzL4aZyjnHL8FfAW4xMwG3prPO0PHcfqhozkCM7sFeEZDtucAzwNWA38u6S7gQuArZvaDNvfxznDM\n5IIM1G2oVGd5bjtnuDIqE96nc4YrQyXxpOFeyblwDIWqojsMMmd4f3IMyjDnqwisvL04PlJKuKDk\nNiXH+H0aGiydM2zz/OtWpDStill0jH853hbgkvKFpL2BlwJ/KemJwHfM7PV1dXhn6DhOP0zwF8HM\n7gDOAs6StAPw3KYy3hk6jtM9E4pnWEXpbvOtpnzeGXZIG2frYeoLQ7gqN5B0F7zU2To9wvzleMEm\nMm8oHL/fJ0nnDClxI9LNUFJrBswOjxOXmnnGl/gfq6xn99vnVpcGe4D5nz9cC7cJw+gqd6Uw/586\nX4/KoolxOGUf0DtDx3G6x7cKdVLaOltXRbpu2hMZ5u9OlzOgxMItdbpWUHnhGCvDx5XHvZJ0yFOl\nDFMjSCD8c8SbIAdlmC7rS9fcxcqwlIIqj7vdXRyD6ouVYbhVGhIsXcYYP9OgFptcnOIm5UaEi9b5\nekIbQkk6AHgjsD9R/2ZmxzSV7bUzlLQaOJ3i+3Ommb0vuf4HwBsoHtvPgZPM7Po+2+Q4zpiYzC/A\nl4C/B748aAt66wwlzQBnUITaXg9cLmlt0tl91sw+XuY/BvgghZ/Qdk+dYmxyqalyrUn3R86F8oJZ\n15oVoXDqdV2lDMPx8Umeqrhfg8wZhrnCnJpMo7DCPAm4opR9K8s890ZZm55LlWtNTpHXzR2mQnbR\nM7mtQh82s48MU7CxMyzN0s9gdrep68zspy3qPgy4uXSYRNLZwLHAts4w2ZxlBcUjdBxn2plcZ3i6\npPcAFxG5nJrZlfkiBdnOUNJBwJ8ALwZ+CNxJ8cP6JEkPAp8APlWzBd8+wK1Rej1weMV93gC8jWLa\n5kWZtpwEnASzkRunjSZn6zpLdGpNbjNnmMZPqFqON8+snAZjgFkFmM4ZhvQ2j+1YW4UK47heMCsJ\n75k9tWLr3KypqTXIsHgiMHXULqXgruXcYd2+z23mDGv8vSubGL935+uSyRlQ/ifwSoq+JPRNRqZv\nialThn8JfAx4nZnNUWxlGJ2Xlzf9VKZ8Vb81T/mZ2RnAGZJeDpwKvKoizxpgDcBMEc/McZyFzmTm\nDF8GHGhmA6+MznaGZnZCzbUNwIcb6l4P7Bel9wVur8l/NkXnu+jJBXOtWo6XKpx0WV6Vn+G296lZ\nuWrOMLxPFWHwO1y2e5IBYFV53JG5hFFLFLl95X+VDb57btYgqUKRKmV4z9zjzmUVVZ81ncZMFWHd\ncrxcsFenhskNk79H8Y3eMGjBNnOGv1Fx+l7g+2WnmONy4ODS1H0bcDyFmozrPtjMflgmf41iOO44\nzvbAZDrDvYAbJV3O3DnDTlxrTqRY1/fvZfoFwHco5g5PM7N/qipkZlsknUwROWIGOMvMrpN0GrDO\nzNYCJ0t6McWP7c+oGCI7jjOFTG453nuGLdimM9wK/PdgQZa0F8Vw9nCKSLOVnSGAmV0AXJCce3f0\n/s1DtHm7oC5qTZNrxyCRrlPnawDtUr5Jw9ekR5gdOqej4WWPLd88sTzuxSxNBpQ95udddnNZ/4a5\nWdMhcVU7y88RPteKn89mzT2P9LlV7cvcxrUpLbPonKvrmIwy3NnMvhKfKP2Zv95UsE2k6/0TV5oN\nwJPM7G58+sRxnCqCNbnNq1v+TNI2y7GkP6Fw6WukjTL8hqTzgXPK9HHApWX47XvyxRYPowRoqCvb\ntBwsPpdTOqkiAvJhsOdFbmC+UWWnICf3L49PKI/7RoWCZSanDGO36JBnS1l/mWfP0gW1ShkGW0to\nd+JHtLxCGeaeT90zbXJ7is8Powi364ANkzOgHAOcL+mPKBZwPKU810ibzvANwG9QRJEVhSvNF0t3\nmxcO1VzHcbZ/JjBnYGYby9Vs/wZcARyXugbmaOwMzczK/QXuNbN/k7QzsAuzv9lOwjC7bM1k3sfp\nKsftZcmxae5wTiL1yK6KWj1vHjEowFQRBl8bmHWtySnDjcynVIaUsm636+fet8prPN1/Zfn8u+bm\nCtPnVvX8m/4Og4zw4u/EophXHLMylHQ/c/2YlwEHAsdJMjNbWV1yljauNb9Psfpjd+Agim/9x4Ej\nh2m04ziLgPGH/d+1OVc9bYfJhwGXlTf9YbkCxRmCQVRjmrfK6Tqd38o5Yc/E5tLUEzv10K5SYUvC\n5GG6Li8ownjOMFiLc8owPR9fK1XjktIZu2qNXa7d5TH+rEtL+ZZzsq6yENdZ+pvwgA0RHUvgMvjL\nOuA2Mzu629rbdYabzGxzuU0fkpbgARUcx6mjn7XJbwZuYO7+Zp3RpjP8uqQ/BXaS9BLg9RSxwpwx\nUWfdbFo6ti0AQd1eAWlUh9gpcdv7MHkX5gP3yKTjc7swl58znzR4Qyhb3m/F3fk2NUVfAJaV/5BN\nSxzrrMnOEHQ8TJa0L8Uqtb+iCOwySNllbdYqt/l7n0IRseb7wOsonKhPHaQxjuMsQh5t+WrHh4E/\nJjP4lvRnmfOPoQjn1Ugba/JW4O/Kl9MRdT5sg/i3peeyq1Xqlq3kTNEQRT0IynCXJJ0qxvh9OjdY\nNVcYIjBk6k/VX127K5aTbJszzRzJpKuu5eZwoXov5UXNYMvxVpUeK4E1ZaQqACQdDWwwsyskvSBT\nxy9L+isze1dU7nEUy4G/2KYRdfEMv0/N3KCZPb3NDRzHWaS0/2XYaGaH1lw/AjhG0lEUP4srJX3a\nzH4nynMMcK6kD5rZ2yQdDHwF+Bsz+0SbRtQpw2CteUN5DGuQX8HcLX0cx3Hm0qEBxczeCbwToFSG\n70g6QszsYUkvA84uo+o/F3iLmZ3X9j518Qz/s7z5EWZ2RHTpFEnfAk5re5PFTG5YNkwddddyriI7\npBnjk+lYuipqwZKQOYxTd0nSicdzZd6UOG9aPqk/3H9pNOaa5zdUHivmF5qeS5tnOwyLPuL1+J2u\ng1HluxRzi98ADgjnzeyDTXW0sSavkPQ8M/tmedNfYq5tz3EcZz49LLUxs0uASyouxZ6oH6k410jb\neIZnlVYZo1hl/9pBbuJ0wyBKZV7eKq/ipg2agdmvSHpcnklXXUvPx3lz9SXHmcgzoqn9FT4SIz07\nZ3DGvwLlz0eto401+QrgGZJWAjKze5vKOI7jjHmYfCpwhpn9LHP9RRSxDs/P1VFnTf4din2Nt8K8\nbT3D7nl7h+Gz4zjONsYf6fr7FKG7HgauZHY3z4OBQyii2Px1XQV1ynAP4CpJV1CEwgmVPxH4FYpF\npKeM+AEcx9keMWDg/elGuJ3ZvwD/UrrUHAHsDdwHfBo4ycweaqqjzpp8uqS/pdhv9Ajg6RSbyN8A\nvNLMfjL6R3AGYZBRx7y8Wysubk2OlZ7DWzLHhzPp+H369Xq4Im+uvuRYtUlxrv0VimSkZ+cMx2Ti\nGf6QITeWq50zNLNHgYvLl+M4TjsmF+l6aNpYk50RSL8Pwyz+r/tO5URSepxTSaqkHskcAbaUmZcE\n5RaCLYT0A0k6fU/F+fh6Wj6pP9w/blPazvQhRJ+16bm0ebbDsOiX501ud7yh8c7QcZx+mMAvgaTd\ny83qBsajFDmO0z1hmNxd1Jq2XCbpHElHKQRhbUmbsP9/DXzAzO4p078AvN3MPIzXCKTDtKU111Kq\nzqfDsnQEXDvM3Jyk41FsWIW+MsQbDMPYe5Jj1aKkXDzDeA+UtJ6k/nD/uE25dlcM89ORdG74OshU\nRFWZKRsR9k8/wV3b8CTgxRQLQz4q6fPAP5jZD5oKtlGGLw0dIUDp1HjUsC11HGeRMAFlaAUXm9kJ\nwO8BrwK+K+nrkp5bV7bNnOGMpB3NbBOApJ2AHUdutdOaukn/VOmkeba5esW/0uFkUFtBfT2QHOP3\n25RhUHV3lccQh7AqVmEa2Trc8K7oXFpfSN/T3KZUNaZKMTqVez5VSnGe8ckZnAkZUCTtAfwO8Erg\np8AbgbUUjtfnAAfkyrbpDD8NfE3SJyk+4msp9k52HMfJMxlT+rcpwg3+upmtj86vk/TxuoJt1iZ/\nQNI1FONwAX9hZheO0trFTPixbBMMIDdHtbUiTzjmptEejdTSTFBUqbIK6XhH7PB+z9JAt+Sn5Ylc\nNOu4wjb7Jt9WHsvd8Cjr33L33PvHbcq1uzw+WjFnmE4zZl2QmP8/PIjAcTVZMjnXmlPN7AvxCUm/\nZWbnmNn76wq2da25AdgSNpGXtKuZ+SbyjuNUM+bleBGnAF9Izr2TYohci28i3wPxD2Jb36WqecA0\nXTW/tTk5pgbhKjfnFam/dKq+4p+51Ni7Kow8gqW46isU5gpzyjAOfBTq+8+56fS+VWo1HBO/7Sr3\n73BMFWI6pxi/b/o7DMKiVIxj/NCSXkph2N1H0keiSyuZXdtZi28i7zhO94x/Od7tFBvMH0MRWCZw\nP/DWNhX4JvIdEP7mwwQFrfIzTK/VWZPDuZyrXaUyTK2xQWEFFRb77wdjcYgZvKKM5LbTj8sT4SsU\n3yFnYU73SIZtc4SU9T1U1n9n0pa4TalaTKzhVcow93zqnukw/p6DsF0v1Rt/cNfvAd+T9Bkza6UE\nU3wTecdx+mG8w+QvmNn/pgg7GIs1UbgfNu7m2aYzPIUi9H+8ifyZQ7TXiQg/mmFOscpC3LiqhGYl\nmE4PAuxeTukphOtN5+di4RYUWVCGQezts6E4Lgs/wnGhYGlO3VE3lcfYmlxakTffPSc5TxHG1aft\nLD+HlZ8r/qy551EXmyL33Ov8D7drlTcM4x8mv7k8Hl2bq4ZWm8hL+hLwJTO7sym/4zhOl8vxJC0H\nLqX4dV0CnGtm75lzO7M7yrcbgYfKfutJwFMo9k9uJGvsVMF7JW0EbgRuknSnpHcP8CFWS7pJ0s2S\n5kXFlvQ2SddLukbS1yQ9oW3djuMscLpbjrcJeJGZPYNiJclqSc/J5L0UWC5pH+BrwGuAf2hzkzpl\n+BaKCNe/aGb/ASDpQOBjkt5qZh+qq1jSDHAG8BIKf4nLJa01s+ujbFcBh5rZg5L+EPgA8NttGr49\n07R0LB6m5ZysM/7Ic96vCB4uqeEk3mAxxF8II97USrRHWWhlPI4dwIByX/lpwmq8YE8Jw+U7k2Pc\nzlBN+TnqPmtueJw6YUPzEkenBR06XZuZMeuvtbR85Yy4KvuTE4GPlotGrmpznzo3uN8FTggdYdmo\nWyjW/f1ui7oPA242s1vMbDNwNnBsnMHM/t3Mwvf1O8C+bRrtOM4U0GGgBkkzkq4GNgAXm9ll+ax6\nLvAK4F/Lc60Wl9RlWmpmG9OTZnanpCpPkJR9gFuj9Hrg8Jr8J5IZ20s6icLxm4EClC0gchGvq1RH\n+nBT95nYsT89l4u9UOWzvCJ1qQnqrEoZhkYFZRgavM1CE8mA5XfPLZM2NvZ9SVVpUIA/TdKxa01o\nZ+JiU+Uz3jamQ9UzzQmbmq1WPMJ1YDADyipJ66L0GjNbM6e6YguSQyTtBpwn6Wlmdm1FXW+mWHFy\nnpldV45m/71NI+o6w7rFNG0W2lT1W5XSttyW9FCKXffmFyoezBqAmblmc8dxFirth8kbzezQNhnN\n7B5JlwCrgXmdoZldSjFvGNK3AG9qU3ddZ/gMSfdVnBfVq/NT1gP7Rel9KbzE51YmvRh4F/ArIUzY\nYqAuYMMgrh1BxYQHl7qQhPOxu0n4o+5aFt42d7giOcKsuksnVMKNg+SKV9jtnJTNlYH5Dt9BCf5X\ncqyaMyzv+UBZb/hc8Wdtei5VKrBNmK+URbncro5urcl7Ao+UHeFOFEFjKoMulBbkdwD7E/VvZvai\npvvUbRU6zIKKmMuBgyUdQBGa5Hjg5XEGSc8EPgGsNrMNI97PcZyFQrd+hnsDnyqNsjsAXzCz8zN5\nz6GInXDmoC3obUMoM9si6WTgQgoBdFY5hj8NWGdma4G/oVjxf0653O8nZnZMX22aBEEx5CxVqfM1\nzKrFnEJpM2eYi8UAsyIuCMCdS6WloPdj3R8aljYq3CBUHM8zhvJt5gxzSwFThVihDMO2P4lRuXLO\nMH0udXOGTYp8kHnBRa0YO+oMzewa4Jkts28xs48Nc59ed8czswsoVqzE594dvX9xn/d3HGdCTC6e\n4ZclvR44j9nZENrsmOdbhXbIKAEb6urLBWOAvCJMfe7qpumCkNsjqK9Y0aWSNr1hnTJMy6aqMi4f\njsFSnPoXxsqwfJ+6G6ZxG+L3qTW5ZqeA2iAOXbBoLM2T+aCvKo9/FJ0z4MCmgt4ZOo7TPRPaHc/M\nsnucNOH7JjuO0wsT2ByPMhL/qZLWlOmDJbUK3uDKcMw0OV/HeXZI0lXOvk1Ra4IBIY4fs2NyDKPi\npWUlK6vCcaRj9TDu3BbnsOIGy5I6wpg0dqBKLTw5Q0rUpvvKNqS+1/clx7j6pqg1dcvx6p5/rsxi\nZ/xBa7bxSYrgrr9UptdTWJhz1udtuDJ0HKcXtrZ8dcxBZvYByt85M3uIlgvXXBlOmCrn69QQk90L\nmVlVl+5H90CSjm0i4dyyzHEmskCsCK4tuQgQwcE6VoaDGFDS9XKpMizl3wORmkx8rue51FRtpBdu\nk36MOteapojXi9ptpoEJKsPNpWO2AUg6iLljkSzeGTqO0wsT+rF4L/BVYD9Jn6GIvPWaNgW9MxwT\nbZ2v4zy5CNexiszNGQYl+ECShlklFerZITnOaVf5m7oyKMScS03sqD2IMkzrS+YOQ4SvqjgN6Zxh\n+Fyxa03bOcMq15qmXfLqWOyqcSuT2SnUzC6SdAXwHIrh8ZurAs5U4Z2h4zi9MIkfBElfM7MjmQ3f\nFZ+rxTvDHoiVQ84Bu82i/zSYoXnyAAAPGElEQVRaVpqG2V/fdO4wpB9M0vH7NCpXOkcZ1/9weXK3\nUqItS63AsTLMBXdILdJxg8v6NpdKNN3mJI4DkVt+90ByhLyzdTpXWLU7Xm6OsOp6k1pcbFbmcc8Z\nllsD7EwRDuwXmDWarAQe36YO7wwdx+mFMf8AvI4iOv/jKVxrQmd4H0XE/Ua8M1yApH6GqSKMhdVM\nci4on2AZDiqpKhBE0/2r7hnqX1EquJ3LYywMlTY8qcwi6dkUjDZ1P4RZRZiqx7rgrmkorza74w0T\nmMEpGPfSZDM7HThd0hvN7KPD1OGdoeM4nTOh1XiY2Ucl/RLz4xn+Y1NZ7wzHTJNVuSpvTinG1zYn\neR5O0nGZVLBlBBzQ7F5YFfVradnwmUQaVCnbXICJnEKsOpceY2ty27nCNhtCtWGxW5FjJqGkJf0T\ncBBwddQEA7wzdBxn/EzQ6fpQ4KnljnoD4cvxHMfphQktx7sWeNwwBV0Z9kxTjMM2v545QwrMd8QO\nx7phcs7JuspQ0HaYHLvupG0Zpv50lV7sLpMOi1OXmiqf7ibDSZ1rzSgGlcVqdJmgMlwFXC/pu8wN\n7toYQd87Q8dxOmeCneF7hy3oneGEGMaQUre8Kac8hzHU1Cm3NDRYUIRxtK6mTbWrInXnFGJ6P5hv\nKKlTkWm727jWhDa54WR4JmhN/vqwZb0zdBynF8b5AyHpfqr3ZRdgZrayqQ7vDMdEl3OHVUEd0l/h\ncJ+q2EVNbagKHhvmBtNQYenSvrr60/vE9ecUaKoUYf7cYG7JHczfHzk3V1j1/H2ucHi6HCZL2o/C\nNeZxFH+ONaWT9ez9zHatKjsI3hk6jtMLHf4gbAHebmZXStoVuELSxWZ2fXe38M5w4tTNHbYJJppa\nmEM9bcInNe0PHNeTOlenc4VVy/3aWJPT+dCcQowVbjqP2CZv6mSdU33xuabd8XyeME+Xy/HM7A7g\njvL9/ZJuAPYBvDN0HGfhM4AyXCVpXZReY2ZrqjJK2p9iQ/nLRmhaJd4Zjpk2eyvn8lR9udoG0BzW\nah2UYBroICjCNAxYm3tVLX1Ll8ltStJV84DptfR8/D43Z1hlOc79E/tcYXsGtCZvNLNDmzJJ2gX4\nIvAWM7uvKf+geGfoOE7ndO1nKGkpRUf4GTP75w6r3oZ3ho7j9EJXc4aSBPw9cIOZfbCjaufhneEC\nIf7iVLnQNDFI3vRLWmcwCEOdMBwOw+aZ5HxXw+SQzjlHx+/TPKmxpKp8ep86A0rOlcYNJ810rAyP\nAF4JfF/S1eW5PzWzC7q7hXeGjuP0RFedoZl9k5Z7H4+Cd4YToo0hJUeVMknV2CDhiFIFFDtQp4aS\nNDBElTLMtaGq3TllmDN0VF1L87TZ6W6UmIVVuOFkLpNajjcK3hk6jtM5EwzUMDTeGS5AUkfs9EtV\npcKGUTZNoapg/pxh6mSdc7CO89S1MXfvVOVVOYLn5gGr9jPJzRVW4XOF3TBtz8s7Q8dxOseVoTMw\ndXOHTQpxEKp+pdNzVapph+RcOmeY2wivLTmrbqoQq5Rh2qY6Zdh0n7o9kHMKZ9r+2cfJuHfH64Je\nw/5LWi3pJkk3Szql4vrzJV0paYuk4/psi+M44+XRlq+FQm/KUNIMxebNLwHWA5dLWptEmvgJ8Grg\nHX21Y1qIvxSpyhokmMOw94zrr1OGOSU46pxhmq5Thk15qlReG0WY4opweNyaPJfDgJvN7BYASWcD\nxxJFmjCzH5fXpk1RO45Tg88ZzmUf4NYovR44fJiKJJ0EnARj8LycAqpU5ChfvKCsQl1dKcM2NCnD\nqnnMQTZqaqsIp+0fdxqYtmfaZ2dY1W8NvJcpQBnOZw3AjDRUHY7jjI9pNKD02RmuB/aL0vsCt/d4\nP8dxFhCuDGe5HDhY0gHAbcDxwMt7vN92Q87dpsqQkuZtE4svHSqmw+OqYXjOpSbNNyg5N5a6YWzb\nfY0HrTfNmzJt/9yTZBqVYW+uNWa2BTgZuBC4AfiCmV0n6TRJxwBI+kVJ64HfAj4h6bq+2uM4zvgw\nirnoNq+FQq9O12WInQuSc++O3l9OMXx2KmhSiDDfIXsUg0rV/XKO2WmbYjeKQXbHy13L3beqfJ2x\nJVdmkKV2rgiHY9qUoa9AcRync9y1xumFYZbspQpxGLVXV98oSmqQsoPM/7kiXDh4Z+g4jlPiw2Sn\nN8YV1KHKmjxIWLFB7p3LWzdnOEietvW3KeO0x5fjOY7j4MNkZ0wME9ShqyV8bSzcozCIkmvjU9mm\nnqYyznBM2zPtNYSX4ziLk+B03ebVhKSzJG2QdG1PzQVcGU49bZVa1aqVQBfBHrpiGOU2jBJsU9YZ\njQ6f7T8Afwv8Y3dVzsc7Q8dxOqfLOUMzu1TS/h1Vl8U7Q8dxOmdAa/IqSeui9JoyUtVY8c5wO6HO\n7Qaql/ClZQPDusl0QZv7jeLUPch9nNEYwKi20cwO7a8l7fDO0HGcznHXGmfiNClEqDeuxHVUMS5j\nyyD1D+LWM23/oNPMtD1rd61xHKdzOnat+RzwbeDJktZLOrGPNrsy3E6pmwdMaVKKdfVOAleC00GH\n1uQTOqqqFu8MHcfpnK342mRngVK3hC+ljfLqa36li2V9rgYXBtP2d/DO0HGczpnGPVC8M1yEDDKf\nmGMhfdGnTYEsFqbt7+KdoeM4neN+ho7jOHhwV2dKaeNkPWmmTWU40/c3887QcZzOcQOKs92R+3Xv\nSzFOm5pw8kzb39I7Q8dxOseVobNomLZffWf8TNt3xDtDx3E6x63JjuM4uJ+h4zgO4J2h4zjONtyA\n4jjOoseVoeM4TokrQ8dxFj0GbJ50IwbEO0PHcTrHna4dx3FKpm3OsNfd8SStlnSTpJslnVJxfUdJ\nny+vXyZp/z7b4zjOeAgGlDavhUJvnaGkGeAM4KXAU4ETJD01yXYi8DMzeyLwIeD9fbXHcZzx0tVW\noeOiT2V4GHCzmd1iZpuBs4FjkzzHAp8q358LHClJPbbJcZwxEJbjtXktFPqcM9wHuDVKrwcOz+Ux\nsy2S7gX2ADbGmSSdBJxUJjc9CNf20uJ+WEXyeRYw09RWmK72TlNbAZ48SuGtcOEDxWduw4J4Ln12\nhlUKz4bIg5mtAdYASFpnZoeO3rzxME3tnaa2wnS1d5raCkV7RylvZqu7asu46HOYvB7YL0rvC9ye\nyyNpCfAY4O4e2+Q4jlNJn53h5cDBkg6QtAw4Hlib5FkLvKp8fxzwf81snjJ0HMfpm96GyeUc4MnA\nhRRR4s8ys+sknQasM7O1wN8D/yTpZgpFeHyLqtf01eaemKb2TlNbYbraO01thelr78jIhZjjOE7P\nTteO4zjTgneGjuM4LLC1yZJ+DNxPsUpni5kdKulvgP9FEQTjR8BrzOyeNmUXalvL8jPAOuA2Mzu6\nz7aO0l5Jy4FLgR0pvi/nmtl7Fmhb9wP+EXgcxeKGNWZ2ep9tHaW9ZdmzgKOBDWb2tAXe1tXA6RQ2\ngDPN7H19t3esmNmCeQE/BlYl534VWFK+fz/w/rZlF2pby+tvAz4LnL+Q20vhC7pL+X4pcBnwnAXa\n1r2BZ5XvdwV+ADx1oT7b8trzgWcB1y7w78EMRUd5ILAM+N44nu04Xwt+mGxmF5nZljL5HQp/xQVJ\n27ZK2hf4NeDMcbWtijbttYKfl8ml5WvsVreWbb3DzK4s398P3ECxymnstP0umNmlTNi3tmVb2yyv\nnWoWWmdowEWSriiX4KW8FvjKkGW7ZpS2fhj4Y8a7Tn3o9kqakXQ1sAG42Mwu67GdMNqzBaCMgPRM\nCiXbNyO3d4wM29aq5bUT+aHpiwU1ZwgcYWa3S3oscLGkG8tfTiS9C9gCfGbQsguprZLC/NAVkl7Q\nY/s6aS+AmT0KHCJpN+A8SU8zsz7Xh4/yPUDSLsAXgbeY2X09trOT9o6ZYdvaaunsNLOglKGZ3V4e\nNwDnUUhzJL2KYpL5FVZOYLQtuwDbegRwTDmRfTbwIkmf7rOtI7Y3ruMe4BKg13Wno7RV0lKKjvAz\nZvbPfbazi/aOmxHa2mZ57XQz6UnL8AJWALtG7/8fxT/dauB6YM9Byy7Etib1vIAxGFBGfLZ7AruV\n73cCvgEcvUDbKgpr8oen4Xsb1bE/YzCgjPhslwC3AAcwa0D5H+N6zmP5W066AdHDPrB8wN8DrgPe\nVZ6/mWKu4ury9fHy/OOBC+rKLsS2JvWMqzMc5dk+HbgKuIYidNq7F3Bbn0cxdLsmynfUQm1vmf4c\ncAdFaL/1wIkLuK1HUVjof9T3/9gkXr4cz3EchwU2Z+g4jjMpvDN0HMfBO0PHcRzAO0PHcRzAO0PH\ncRzAO0MnQdKjkq6WdK2kL5erTuLrb5X0sKTH1NSxt6TzM9cukTRURCFJR0v682HKOk4T3hk6KQ+Z\n2SFWhJO6G3hDcv0Eiv1tXlZTx9uAv+uhbf9KsXpn5x7qdhY53hk6dXybaDG+pIOAXYBTKTrFHL8J\nfLUss5OksyVdI+nzFKtYQn2/Kunbkq6UdE65phhJR0m6UdI3JX0kqEwrnGIvoVg25jid4p2hU0kZ\nfPZI5u5oeALFiolvAE8uF/un5Q4AfmZmm8pTfwg8aGZPB/4KeHaZbxVFp/piM3sWRaDbt5XBZD8B\nvNTMnkexHDBmHfDL3XxKx5nFO0MnZacyXNddwO7AxdG144GzzWwr8M/Ab1WU3xu4M0o/H/g0gJld\nQ7FUDuA5wFOBb5X3exXwBOApwC1m9h9lvs8l9W+gWCbmOJ2y0EJ4OZPnITM7pDSQnE8xZ/gRSU8H\nDqYI+wTFYv1bgDPS8sDy5FzVmk9RxEacM9yW9MyG9i0v7+E4neLK0KnEzO4F3gS8owyLdQLwXjPb\nv3w9HthH0hOSoj+giMISuBR4BYCkp1EEfoAiovIRkp5YXttZ0pOAG4EDy+CsAL+d1P8kioARjtMp\n3hk6WczsKooIJ8eXr/OSLOeV5+MyDwA/Cp0c8DFgF0nXUET3/m6Z707g1cDnymvfAZ5iZg8Brwe+\nKumbwE+Be6NbvJDCquw4neJRa5zOkfQy4NlmduqQ5Xcxs5+rGI+fAfzQzD4kaS/gs2Z2ZJftdRxw\nZej0gJmdR7EL27D8fmlUuQ54DIV1GeC/AW8frXWOU40rQ8dxHFwZOo7jAN4ZOo7jAN4ZOo7jAN4Z\nOo7jAN4ZOo7jAPD/AXcTmQBfRNcWAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x10f605630>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#Note that we display RA in the convential way by going from max to min\n",
"plt.imshow(intcloud.value, \n",
" origin='lower', \n",
" extent=[ra.value.max(), ra.value.min(), dec.value.min(), dec.value.max()], \n",
" cmap='hot', \n",
" interpolation='nearest', \n",
" aspect='equal')\n",
"plt.colorbar().set_label(\"Intensity ({})\".format(intcloud.unit))\n",
"plt.xlabel(\"RA (deg)\")\n",
"plt.ylabel(\"Dec (deg)\");"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Measuring The Column Density of CO"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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",
"\n",
"$N=C \\frac{\\int T_B(V) dV}{1-e^{-B}}$\n",
"\n",
"where the constants $C$ and $B$ are given by:\n",
"\n",
"$C=3.0\\times10^{14} \\left(\\frac{\\nu}{\\nu_{13}}\\right)^2 \\frac{A_{13}}{A} {\\rm K^{-1} cm^{-2} \\, km^{-1} \\, s}$\n",
"\n",
"$B=\\frac{h\\nu}{k_B T}$\n",
"\n",
"(Rohlfs & Wilson \"Tools of Radio Astronomy\"). "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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. \n",
"\n",
"First, we look up the wavelength for these emission lines and store them as quantities. "
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"lambda13 = 2.60076 * u.mm\n",
"lambda18 = 2.73079 * u.mm"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since the wavelength and frequency of light are related using the speed of light, we can convert between them. However, doing so just using the to() method fails, as units of length and frequency are not convertible:"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"tags": [
"raises-exception"
]
},
"outputs": [
{
"ename": "UnitConversionError",
"evalue": "'mm' (length) and 'Hz' (frequency) are not convertible",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mUnitConversionError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-25-b4a9b54d7f21>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mnu13\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlambda13\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mHz\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/quantity.py\u001b[0m in \u001b[0;36mto\u001b[0;34m(self, unit, equivalencies)\u001b[0m\n\u001b[1;32m 842\u001b[0m \u001b[0;31m# and don't want to slow down this method (esp. the scalar case).\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 843\u001b[0m \u001b[0munit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mUnit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0munit\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 844\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_new_view\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_to_value\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0munit\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mequivalencies\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0munit\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 845\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 846\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mto_value\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0munit\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mequivalencies\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/quantity.py\u001b[0m in \u001b[0;36m_to_value\u001b[0;34m(self, unit, equivalencies)\u001b[0m\n\u001b[1;32m 814\u001b[0m \u001b[0mequivalencies\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_equivalencies\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 815\u001b[0m return self.unit.to(unit, self.view(np.ndarray),\n\u001b[0;32m--> 816\u001b[0;31m equivalencies=equivalencies)\n\u001b[0m\u001b[1;32m 817\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 818\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mto\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0munit\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mequivalencies\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/core.py\u001b[0m in \u001b[0;36mto\u001b[0;34m(self, other, value, equivalencies)\u001b[0m\n\u001b[1;32m 975\u001b[0m \u001b[0mIf\u001b[0m \u001b[0munits\u001b[0m \u001b[0mare\u001b[0m \u001b[0minconsistent\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 976\u001b[0m \"\"\"\n\u001b[0;32m--> 977\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_get_converter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mequivalencies\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mequivalencies\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvalue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 978\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 979\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0min_units\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mother\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1.0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mequivalencies\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/core.py\u001b[0m in \u001b[0;36m_get_converter\u001b[0;34m(self, other, equivalencies)\u001b[0m\n\u001b[1;32m 909\u001b[0m \u001b[0;32mpass\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 910\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 911\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mexc\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 912\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 913\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_to\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/core.py\u001b[0m in \u001b[0;36m_get_converter\u001b[0;34m(self, other, equivalencies)\u001b[0m\n\u001b[1;32m 895\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 896\u001b[0m return self._apply_equivalencies(\n\u001b[0;32m--> 897\u001b[0;31m self, other, self._normalize_equivalencies(equivalencies))\n\u001b[0m\u001b[1;32m 898\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mUnitsError\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mexc\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 899\u001b[0m \u001b[0;31m# Last hope: maybe other knows how to do it?\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/core.py\u001b[0m in \u001b[0;36m_apply_equivalencies\u001b[0;34m(self, unit, other, equivalencies)\u001b[0m\n\u001b[1;32m 879\u001b[0m raise UnitConversionError(\n\u001b[1;32m 880\u001b[0m \"{0} and {1} are not convertible\".format(\n\u001b[0;32m--> 881\u001b[0;31m unit_str, other_str))\n\u001b[0m\u001b[1;32m 882\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 883\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_get_converter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mother\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mequivalencies\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mUnitConversionError\u001b[0m: 'mm' (length) and 'Hz' (frequency) are not convertible"
]
}
],
"source": [
"nu13 = lambda13.to(u.Hz)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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](http://docs.astropy.org/en/stable/units/equivalencies.html).\n",
"\n",
"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:"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": [
"nu13 = lambda13.to(u.Hz, equivalencies=u.spectral())\n",
"nu18 = lambda18.to(u.Hz, equivalencies=u.spectral())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next, we look up Einstein coefficients (in units of s$^{-1}$), and calculate the ratios in constant $C$. Note how the ratios of frequency and Einstein coefficient units are dimensionless, so the unit of $C$ is unchanged."
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$2.1792458 \\times 10^{14} \\; \\mathrm{\\frac{s}{K\\,km\\,cm^{2}}}$"
],
"text/plain": [
"<Quantity 217924582157305.03 s / (cm2 K km)>"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A13 = 7.4e-8 / u.s\n",
"A18 = 8.8e-8 / u.s\n",
"\n",
"C = 3e14 * (nu18/nu13)**3 * (A13/A18) / (u.K * u.cm**2 * u.km *(1/u.s))\n",
"C"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we move on to calculate the constant $B$. This is given by the ratio of $\\frac{h\\nu}{k_B T}$, where $h$ is Planck's constant, $k_B$ is the Boltzmann's constant, $\\nu$ is the emission frequency, and $T$ is the excitation temperature. The constants were imported from `astropy.constants`, and the other two values are already calculated, so here we just take the ratio."
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"B = h * nu18 / (k_B * Tex)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The units of $B$ are Hz s, which can be decomposed to a dimensionless unit if you actually care about it's value. Usually this is not necessary, though. Quantities are at their best if you just use them without worrying about intermediate units, and only convert at the very end when you want a final answer."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.21074888275227613 Hz s\n",
"0.21074888275227613\n"
]
}
],
"source": [
"print('{0}\\n{1}'.format(B, B.decompose()))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"At this point we have all the ingredients to calculate the number density of $\\rm CO$ molecules in this cloud. We already integrated (summed) over the velocity channels above to show the integrated intensity map, but we'll do it again here for clarity. This gives us the column density of CO for each spatial pixel in our map. We can then print out the peak column column density."
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Peak CO column density: \n"
]
},
{
"data": {
"text/latex": [
"$8.5782066 \\times 10^{15} \\; \\mathrm{\\frac{1}{cm^{2}}}$"
],
"text/plain": [
"<Quantity 8578206553813961.0 1 / cm2>"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"NCO = C * np.sum(data*dv, axis=2) / (1 - np.exp(-B))\n",
"print(\"Peak CO column density: \")\n",
"np.max(NCO)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### CO to Total Mass"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We are using CO as a tracer for the much more numerous H$_2$, the quantity we are actually trying to infer. Since most of the mass is in H$_2$, we calculate its column density by multiplying the CO column density with the (known/assumed) H$_2$/CO ratio."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Peak H2 column density: \n"
]
},
{
"data": {
"text/latex": [
"$5.0611419 \\times 10^{22} \\; \\mathrm{\\frac{1}{cm^{2}}}$"
],
"text/plain": [
"<Quantity 5.061141866750237e+22 1 / cm2>"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"H2_CO_ratio = 5.9e6\n",
"NH2 = NCO * H2_CO_ratio\n",
"print(\"Peak H2 column density: \")\n",
"np.max(NH2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"That's a peak column density of roughly 50 magnitudes of visual extinction (assuming the conversion between N$_{\\rm H_2}$ and A$_V$ from Bohlin et al. 1978), which seems reasonable for a molecular cloud."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We obtain the mass column density by multiplying the number column density by the mass of an individual H$_2$ molecule. "
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [],
"source": [
"mH2 = 2 * 1.008 * u.Dalton #aka atomic mass unit/amu\n",
"rho = NH2 * mH2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A final step in going from the column density to mass is summing up over the area area. If we do this in the straightforward way of length x width of a pixel, this area is then in units of ${\\rm deg}^2$."
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"2.5e-05 deg2\n"
]
}
],
"source": [
"dap = dra * ddec\n",
"print(dap)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now comes an important subtlety: in the small angle approximation, multiplying the pixel area with the square of distance yields the cross-sectional area of the cloud that the pixel covers, in *physical* units, rather than angular units. So it is tempting to just multiply the area and the square of the distance."
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1.5625 deg2 pc2\n"
]
}
],
"source": [
"da = dap * d**2 # don't actually do it this way - use the version below instead!\n",
"print(da)"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$0.00047596472 \\; \\mathrm{pc^{2}}$"
],
"text/plain": [
"<Quantity 0.0004759647184167322 pc2>"
]
},
"execution_count": 35,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dap.to(u.steradian).value * d**2"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"But this is **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:"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"tags": [
"raises-exception"
]
},
"outputs": [
{
"ename": "UnitConversionError",
"evalue": "'deg2 pc2' and 'cm2' (area) are not convertible",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mUnitConversionError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-36-d7c4d4dcf9cc>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mda\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mu\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcm\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/quantity.py\u001b[0m in \u001b[0;36mto\u001b[0;34m(self, unit, equivalencies)\u001b[0m\n\u001b[1;32m 842\u001b[0m \u001b[0;31m# and don't want to slow down this method (esp. the scalar case).\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 843\u001b[0m \u001b[0munit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mUnit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0munit\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 844\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_new_view\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_to_value\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0munit\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mequivalencies\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0munit\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 845\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 846\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mto_value\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0munit\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mequivalencies\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/quantity.py\u001b[0m in \u001b[0;36m_to_value\u001b[0;34m(self, unit, equivalencies)\u001b[0m\n\u001b[1;32m 814\u001b[0m \u001b[0mequivalencies\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_equivalencies\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 815\u001b[0m return self.unit.to(unit, self.view(np.ndarray),\n\u001b[0;32m--> 816\u001b[0;31m equivalencies=equivalencies)\n\u001b[0m\u001b[1;32m 817\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 818\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mto\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0munit\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mequivalencies\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/core.py\u001b[0m in \u001b[0;36mto\u001b[0;34m(self, other, value, equivalencies)\u001b[0m\n\u001b[1;32m 975\u001b[0m \u001b[0mIf\u001b[0m \u001b[0munits\u001b[0m \u001b[0mare\u001b[0m \u001b[0minconsistent\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 976\u001b[0m \"\"\"\n\u001b[0;32m--> 977\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_get_converter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mother\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mequivalencies\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mequivalencies\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mvalue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 978\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 979\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0min_units\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mother\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvalue\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1.0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mequivalencies\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/core.py\u001b[0m in \u001b[0;36m_get_converter\u001b[0;34m(self, other, equivalencies)\u001b[0m\n\u001b[1;32m 909\u001b[0m \u001b[0;32mpass\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 910\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 911\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mexc\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 912\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 913\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_to\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mother\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/core.py\u001b[0m in \u001b[0;36m_get_converter\u001b[0;34m(self, other, equivalencies)\u001b[0m\n\u001b[1;32m 895\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 896\u001b[0m return self._apply_equivalencies(\n\u001b[0;32m--> 897\u001b[0;31m self, other, self._normalize_equivalencies(equivalencies))\n\u001b[0m\u001b[1;32m 898\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mUnitsError\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mexc\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 899\u001b[0m \u001b[0;31m# Last hope: maybe other knows how to do it?\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/core.py\u001b[0m in \u001b[0;36m_apply_equivalencies\u001b[0;34m(self, unit, other, equivalencies)\u001b[0m\n\u001b[1;32m 879\u001b[0m raise UnitConversionError(\n\u001b[1;32m 880\u001b[0m \"{0} and {1} are not convertible\".format(\n\u001b[0;32m--> 881\u001b[0;31m unit_str, other_str))\n\u001b[0m\u001b[1;32m 882\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 883\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_get_converter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mother\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mequivalencies\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mUnitConversionError\u001b[0m: 'deg2 pc2' and 'cm2' (area) are not convertible"
]
}
],
"source": [
"da.to(u.cm**2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The solution is to use the `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."
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$0.00047596472 \\; \\mathrm{pc^{2}}$"
],
"text/plain": [
"<Quantity 0.00047596471841673214 pc2>"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"da = (dap * d**2).to(u.pc**2, equivalencies=u.dimensionless_angles())\n",
"da"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$4.5318534 \\times 10^{33} \\; \\mathrm{cm^{2}}$"
],
"text/plain": [
"<Quantity 4.531853390818705e+33 cm2>"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"da.to(u.cm**2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, multiplying the column density with the pixel area and summing over all the pixels gives us the cloud mass."
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$317.63786 \\; \\mathrm{M_{\\odot}}$"
],
"text/plain": [
"<Quantity 317.63786094118024 solMass>"
]
},
"execution_count": 39,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"M = np.sum(rho * da)\n",
"M.decompose().to(u.solMass)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercises"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The astro material was pretty heavy on that one, so lets focus on some associated statistics using `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?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Similarly, compute the standard deviation and variance (if you don't know the relevant functions, look it up in the numpy docs or just type np.<tab> and a code cell). Do they have the units you expect?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 3. Using Quantities with Functions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`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](http://en.wikipedia.org/wiki/Mars_Climate_Orbiter#Cause_of_failure). Below, we provide a simple example."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Suppose you are working on an instrument, and the bigwig funding it asks for a function to give an analytic estimate of the response function. You determine from some tests it's basically a Lorentzian, but with a different scale along the two axes. Your first thought might be to do this:"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [],
"source": [
"def response_func(xinarcsec, yinarcsec):\n",
" xscale = 0.9\n",
" yscale = 0.85\n",
" xfactor = 1 / (1 + xinarcsec/xscale)\n",
" yfactor = 1 / (1 + yinarcsec/yscale)\n",
" \n",
" return xfactor * yfactor"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"You meant the inputs to be in arcsec, but you send that to your hapless collaborator, and they don't look closely and think the inputs are instead supposed to be in arcmin. So they do:"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"0.19640564826700893"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"response_func(1.0, 1.2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And now they tell all their friends how terrible the instrument is, because it's supposed to have arcsecond resolution, but your function clearly shows it can only resolve an arcmin at best. But you can solve this by requiring they pass in `Quantity` objects. The new function could simply be:"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [],
"source": [
"def response_func(x, y):\n",
" xscale = 0.9 * u.arcsec\n",
" yscale = 0.85 * u.arcsec\n",
" xfactor = 1 / (1 + x/xscale)\n",
" yfactor = 1 / (1 + y/yscale)\n",
" \n",
" return xfactor * yfactor"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And your collaborator now has to pay attention. If they just blindly put in a number they get an error:"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {
"tags": [
"raises-exception"
]
},
"outputs": [
{
"ename": "UnitsError",
"evalue": "Can only apply 'add' function to dimensionless quantities when other argument is not a quantity (unless the latter is all zero/infinity/nan)",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mUnitsError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-43-5d7d1ca80126>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mresponse_func\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1.0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1.2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;32m<ipython-input-42-c1a22f103934>\u001b[0m in \u001b[0;36mresponse_func\u001b[0;34m(x, y)\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0mxscale\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0.9\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marcsec\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 3\u001b[0m \u001b[0myscale\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0.85\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mu\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marcsec\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0mxfactor\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m1\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mxscale\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 5\u001b[0m \u001b[0myfactor\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m1\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0myscale\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/quantity.py\u001b[0m in \u001b[0;36m__array_ufunc__\u001b[0;34m(self, function, method, *inputs, **kwargs)\u001b[0m\n\u001b[1;32m 619\u001b[0m \u001b[0;31m# consistent units between two inputs (e.g., in np.add) --\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 620\u001b[0m \u001b[0;31m# and the unit of the result (or tuple of units for nout > 1).\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 621\u001b[0;31m \u001b[0mconverters\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0munit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mconverters_and_unit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunction\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0minputs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 622\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 623\u001b[0m \u001b[0mout\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mkwargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'out'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/astroconda3/lib/python3.6/site-packages/astropy/units/quantity_helper.py\u001b[0m in \u001b[0;36mconverters_and_unit\u001b[0;34m(function, method, *args)\u001b[0m\n\u001b[1;32m 481\u001b[0m \u001b[0;34m\"argument is not a quantity (unless the \"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 482\u001b[0m \u001b[0;34m\"latter is all zero/infinity/nan)\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 483\u001b[0;31m .format(function.__name__))\n\u001b[0m\u001b[1;32m 484\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 485\u001b[0m \u001b[0;31m# _can_have_arbitrary_unit failed: arg could not be compared\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mUnitsError\u001b[0m: Can only apply 'add' function to dimensionless quantities when other argument is not a quantity (unless the latter is all zero/infinity/nan)"
]
}
],
"source": [
"response_func(1.0, 1.2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Which is their cue to provide the units explicitly:"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$0.0001724307 \\; \\mathrm{}$"
],
"text/plain": [
"<Quantity 0.00017243069807384764>"
]
},
"execution_count": 44,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"response_func(1.0*u.arcmin, 1.2*u.arcmin)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The funding agency is impressed at the resolution you achieved, and your instrument is saved. You now go on to win the Nobel Prize due to discoveries the instrument makes. And it was all because you used `Quantity` as the input of code you shared."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercise"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Write a function that computes the Keplerian velocity you worked out in section 1 (using `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. "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"astropy-tutorials": {
"author": "Ana Bonaca <ana.bonaca@yale.edu>, Erik Tollerud <erik.tollerud@yale.edu>, Jonathan Foster <jonathan.b.foster@yale.edu>",
"date": "April 2014",
"description": "Demonstrates use of the astropy.units, astropy.constants, and the Quantity object for calculating the mass of a galaxy from its velocity dispersion and determining masses of molecular clouds from CO intensity maps. Includes use of matplotlib for making a histogram and an image with a colorbar. Also includes good practices for using quantities in functions you might distribute to other people.",
"link_name": "Using Astropy Quantities for astrophysical calculations",
"name": "",
"published": true
},
"kernelspec": {
"display_name": "Python 3",
"language": "python3",
"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.6.3"
}
},
"nbformat": 4,
"nbformat_minor": 1
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment