LearningPython.units-M31halo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using `astropy.units` and `astropy.constants` to check M31 mass estimates\n",
    "This notebook walks a check of the erratum for Lehner+ (2015) using the `astropy.units` functionality. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "##Set up matplotlib for display in the browser:\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "mpl.rc(\"savefig\", dpi=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "##Very most basic/standard imports:\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "##Automatic unit tracking\n",
    "import astropy.units as u\n",
    "##Costants\n",
    "import astropy.constants as c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The problem: the CGM mass in M31\n",
    "We made an error in Lehner+ (2015) [=LHW15] that required an erratum. The error was in our mass calculation, which went awry. Here we will recalculate the mass in M31's CGM to check our numbers in the erratum.\n",
    "\n",
    "The total mass of gas within an impact parameter is given:\n",
    "\n",
    "$M_g = \\pi \\rho^2 \\mu m_H N({\\rm Si}) ({\\rm Si/H})^{-1}_\\odot$.\n",
    "\n",
    "### Quantities \n",
    "#### Column density in M31\n",
    "We start with the surface density of Si in the inner 50 kpc of M31 from LHW15."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Column of Si ions (I+II+III+IV)\n",
    "sigmaSi = 7.4e13*u.cm**-2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Impact parameter [kpc] and $\\mu$:\n",
    "We'll define $\\rho$ and $\\mu$ next:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rho = 50.*u.kpc\n",
    "mu = 1.3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Solar abundance of Si:\n",
    "Now we're going to need the solar abundance of Si as well. \n",
    "\n",
    "ND users could pull this from the `pyND.solarabundance` module. Or from `pyigm`. For now we'll hard-code it at $10^{-4.49}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "solarSi = 10.**(-4.49)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Hydrogen mass...in units of solar masses!\n",
    "Next we grab the hydrogen nuclear mass using `astropy.constants`. We could do this inline below, but this demonstrates how to make use of the constants functionality explicitly.\n",
    "\n",
    "We'll print the mass to see what it looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Name   = Proton mass\n",
      "  Value  = 1.672621898e-27\n",
      "  Uncertainty  = 2.1e-35\n",
      "  Unit  = kg\n",
      "  Reference = CODATA 2014\n"
     ]
    }
   ],
   "source": [
    "m_H = c.m_p\n",
    "print(m_H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's a lot of info! We're also going to be interested in the mass in different units. Here's an example of transforming the units."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(<Quantity 1.672621777e-24 g>, <Quantity 8.408937594892163e-58 solMass>)\n"
     ]
    }
   ],
   "source": [
    "print(m_H.to('g'),m_H.to('Msun'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculating the total gas mass\n",
    "Last we calculate the total gas mass. Note the use of the `.to('')` functionality to make sure all the distance units are on the same scale and the final result is in solar masses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "187000856.71642455 solMass\n"
     ]
    }
   ],
   "source": [
    "Mg = np.pi*rho.to('cm')**2 * mu*m_H.to('Msun')*sigmaSi*(1./solarSi)\n",
    "print(Mg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# If the cm^2 units don't cancel, we can work around this:\n",
    "Mg=Mg.to('Msun')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Print the final gas mass:\n",
    "Now let's make a formatted printing of the result, including our chosen $\\rho$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gas mass within 50 kpc of M31: 1.87e+08 solMass\n"
     ]
    }
   ],
   "source": [
    "print('Gas mass within {0:0.0f} of M31: {1:0.2e}'.format(rho,Mg))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "```\n",
    "## Mass within $R_{\\rm vir}$\n",
    "The following recalculates the mass assuming a constant surface density of Si (see LHW15 for a more realistic version)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gas mass within 300 kpc of M31: 6.73e+09 solMass\n",
      "(Assumes constant surface density.)\n"
     ]
    }
   ],
   "source": [
    "rho = 300.*u.kpc\n",
    "Mg = np.pi*rho.to('cm')**2 *sigmaSi*mu*m_H.to('Msun')*(1./solarSi)\n",
    "\n",
    "print('Gas mass within {0:0.0f} of M31: {1:0.2e}'.format(rho,Mg))\n",
    "print('(Assumes constant surface density.)')"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}