{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Astropy Units and Quantities\n",
    "\n",
    "## Exercise Solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 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",
    "         ☀️                                  🌎        L2              *(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](http://docs.astropy.org/en/stable/constants/#reference-api) \n",
    "\n",
    "* The mile unit is defined as ``u.imperial.mile`` (see [imperial units](http://docs.astropy.org/en/v0.2.1/units/index.html#module-astropy.units.imperial))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import astropy.units as u\n",
    "from astropy import constants as const"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "d_earth = u.au * (1. * const.M_earth / (3. * const.M_sun))**(1./3)\n",
    "d_earth.to(u.km)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_earth.to(u.imperial.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": [
    "m_JWST = 6500. * u.kg\n",
    "F = (const.G * const.M_earth * m_JWST) / d_earth**2\n",
    "F.to(u.N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_sun = d_earth + 1. * u.au\n",
    "F = (const.G * const.M_sun * m_JWST) / d_sun**2\n",
    "F.to(u.N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3\n",
    "\n",
    "Calculate the 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": [
    "# Schwarzschild radius:\n",
    "bh_mass = 4.31e6 * u.Msun\n",
    "r_s = 2 * const.G * bh_mass / const.c**2\n",
    "print(\"Schwarzschild radius = \", r_s.to(u.Rsun))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Size on the sky given small angle approximation\n",
    "import numpy as np\n",
    "sgr_a_distance = 7940 * u.pc\n",
    "angular_diameter = np.arctan(2 * r_s / sgr_a_distance)\n",
    "angular_diameter.to(u.uarcsec)"
   ]
  },
  {
   "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": [
    "# Temperature calculation\n",
    "radius = 1e6 * r_s\n",
    "temperature = (const.G * bh_mass * const.m_p) / (radius * const.k_B)\n",
    "temperature.to('K')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Energy calculation\n",
    "temperature.to('keV', equivalencies=u.temperature_energy())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}