intro_astroplan.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "![astropy logo](http://astropy.readthedocs.org/en/stable/_images/astropy_banner.svg)\n",
    "\n",
    "<br /><br />\n",
    "\n",
    "# [`astroplan`](https://astroplan.readthedocs.org/en/latest/) <a id='astroplan'></a>\n",
    "\n",
    "\n",
    "`astroplan` is an `astropy`-affiliated package that helps you calculate when objects are observable from observatories on the Earth (Morris et al. 2017, submitted).\n",
    "\n",
    "\n",
    "## Table of contents\n",
    "* [astroplan](#astroplan)\n",
    "  * [Intro to astroplan](#introastroplan)\n",
    "  * [Rise/set/transit times](#riseset)\n",
    "  * [Exoplanets and eclipsing binaries](#eclipsing)\n",
    "  * [Constraints](#constraints)\n",
    "  * [Plots](#plots)\n",
    "\n",
    "\n",
    "\n",
    "<br /><br />\n",
    "\n",
    "### Intro to astroplan <a id='introastroplan'></a>\n",
    "\n",
    "You can describe an observatory using the [`Observer`](https://astroplan.readthedocs.io/en/latest/api/astroplan.Observer.html#astroplan.Observer) object, which knows about some observatories:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "from astroplan import Observer\n",
    "\n",
    "# Use the `at_site` method to access famous observatories\n",
    "salt = Observer.at_site('SALT', )\n",
    "keck = Observer.at_site('Keck')\n",
    "\n",
    "from astropy.coordinates import EarthLocation\n",
    "print(\"Available observatories: \\n\\n{0}\"\n",
    "      .format(', '.join(EarthLocation.get_site_names())))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "You can also specify observatories not included in the database using an [`EarthLocation`](http://docs.astropy.org/en/stable/api/astropy.coordinates.EarthLocation.html) object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "from astropy.coordinates import EarthLocation\n",
    "import astropy.units as u\n",
    "\n",
    "latitude = -33.9249*u.deg\n",
    "longitude = 18.4241*u.deg\n",
    "elevation = 0*u.m\n",
    "\n",
    "# Specify the location of Cape Town:\n",
    "location = EarthLocation.from_geodetic(longitude, latitude, elevation)\n",
    "\n",
    "# Make an observer at Cape Town:\n",
    "cape_town = Observer(location=location, pressure=1*u.bar, temperature=20*u.deg_C)\n",
    "\n",
    "cape_town"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "The pressure and temperature will be used in calculations for atmospheric refraction when computing the altitude and azimuth of targets.\n",
    "\n",
    "We define a target in the sky with a `FixedTarget` object, which simply contains a `SkyCoord` and a name:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "from astroplan import FixedTarget\n",
    "from astropy.coordinates import SkyCoord\n",
    "\n",
    "coord = SkyCoord(ra='14h29m42.94853s', dec='-62d40m46.1631s', \n",
    "                 distance=4.224*u.lightyear)\n",
    "proxima = FixedTarget(coord=coord, name='Proxima Cen')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "We can use the `from_name` method to grab targets with their identifier, resolved via [Sesame](http://cds.u-strasbg.fr/cgi-bin/Sesame):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Targets are stored as `astroplan.FixedTarget` objects\n",
    "target_names = ['Proxima', 'LMC', 'SMC', 'alpha Mensa', 'Polaris']\n",
    "targets = [FixedTarget.from_name(target) for target in target_names]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Let's see which of these targets is \"up\" ($>0^\\circ$ altitude) now:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "from astropy.time import Time\n",
    "\n",
    "# Which targets are visible right now?\n",
    "cape_town.target_is_up(Time.now(), targets)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "The circumpolar souther stars are up, and the circumpolar northern stars are not up.\n",
    "\n",
    "# Rise, Set, Transit times <a id='riseset'></a>\n",
    "\n",
    "What time is the next sunset (i.e.: Sun at $0^\\circ$ altitude)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "sun_set = cape_town.sun_set_time(Time.now(), which='next')\n",
    "\n",
    "print(sun_set.iso)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "At what time is the nearest astronomical twilight?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "astro_twilight = cape_town.twilight_evening_astronomical(Time.now(), which='previous')\n",
    "\n",
    "print(astro_twilight.iso)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "When will Proxima Cen rise next above 10 deg altitude?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "proxima = FixedTarget.from_name(\"Proxima Cen\")\n",
    "\n",
    "proxima_rise = cape_town.target_rise_time(Time.now(), proxima, which='next', \n",
    "                                          horizon=10*u.deg)\n",
    "\n",
    "print(proxima_rise.iso)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "When was the last time Proxima transited the meridian?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "cape_town.target_meridian_transit_time(Time.now(), proxima, which='previous')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "What is the altitude of Proxima right now?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "cape_town.altaz(Time.now(), proxima)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "What is the altitude of Proxima in the next 24 hours?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# Make an array of times over the next 24 hours\n",
    "times = Time.now() + np.linspace(0, 24, 100)*u.hour\n",
    "\n",
    "# Compute the altitude of Proxima at ``times``\n",
    "altitudes = cape_town.altaz(times, proxima).alt\n",
    "\n",
    "ax = plt.gca()\n",
    "ax.plot_date(times.plot_date, altitudes.deg)\n",
    "ax.set(xlabel='Time UTC [MM-DD HH]', ylabel='Altitude [deg]')\n",
    "plt.setp(ax.get_xticklabels(), rotation=45, ha='right')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "<br />\n",
    "*** \n",
    "<br />\n",
    "\n",
    "### [Constraints](http://astroplan.readthedocs.io/en/latest/tutorials/constraints.html) <a id='constraints'></a>\n",
    "\n",
    "Often you will need to plan observations given various constraints, for example, you might need to observe: \n",
    "* Between civil twilights (sun at $<-6^\\circ$ altitude)\n",
    "* The target altitude must be $20^\\circ < $alt$ < 85^\\circ$ (or similar pointing limits for your telescope)\n",
    "* The target is separated from the moon by at least $50^\\circ$\n",
    "\n",
    "`astroplan` provides a framework for defining these observing constraints, and computing whether or not targets are observable given those constraints. Let's implement the above list of constraints: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "from astroplan import (AtNightConstraint, AltitudeConstraint, \n",
    "                       MoonSeparationConstraint, MoonIlluminationConstraint)\n",
    "\n",
    "constraints = [AtNightConstraint.twilight_civil(),\n",
    "               AltitudeConstraint(min=20*u.deg, max=85*u.deg), \n",
    "               MoonSeparationConstraint(min=2*u.deg)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "def plot_constraints(observability_grid, constraints, time_grid, name):\n",
    "    \"\"\"\n",
    "    Plot some constraints in a grid!\n",
    "    \"\"\"\n",
    "    extent = [-0.5, -0.5+len(time_grid), -0.5,  -0.5+len(constraints)]\n",
    "\n",
    "    fig, ax = plt.subplots()\n",
    "    ax.imshow(observability_grid, extent=extent)\n",
    "\n",
    "    ax.set_yticks(range(0, len(constraints)))\n",
    "    ax.set_yticklabels([c.__class__.__name__ for c in constraints])\n",
    "\n",
    "    ax.set_xticks(range(len(time_grid)))\n",
    "    ax.set_xticklabels([t.datetime.strftime(\"%H:%M\") for t in time_grid])\n",
    "\n",
    "    ax.set_xticks(np.arange(extent[0], extent[1]), minor=True)\n",
    "    ax.set_yticks(np.arange(extent[2], extent[3]), minor=True)\n",
    "\n",
    "    ax.grid(which='minor', color='w', linestyle='-', linewidth=2)\n",
    "    ax.tick_params(axis='x', which='minor', bottom='off')\n",
    "    plt.setp(ax.get_xticklabels(), rotation=30, ha='right')\n",
    "\n",
    "    ax.tick_params(axis='y', which='minor', left='off')\n",
    "    ax.set_xlabel('Time on {0} UTC'.format(time_grid[0].datetime.date()))\n",
    "    ax.set_title(name)\n",
    "    fig.subplots_adjust(left=0.35, right=0.9, top=0.9, bottom=0.1)\n",
    "    return fig, ax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "from astroplan.utils import time_grid_from_range\n",
    "\n",
    "# Define range of times to observe between\n",
    "start_time = Time('2017-11-14 22:00:01')\n",
    "end_time = Time('2017-11-15 04:00:01')\n",
    "time_resolution = 1 * u.hour\n",
    "\n",
    "# Create grid of times from ``start_time`` to ``end_time``\n",
    "# with resolution ``time_resolution``\n",
    "time_grid = time_grid_from_range([start_time, end_time],\n",
    "                                 time_resolution=time_resolution)\n",
    "\n",
    "for target in targets: \n",
    "    observability_grid = np.zeros((len(constraints), len(time_grid)))\n",
    "\n",
    "    for i, constraint in enumerate(constraints):\n",
    "        # Evaluate each constraint\n",
    "        observability_grid[i, :] = constraint(cape_town, target, times=time_grid)\n",
    "    \n",
    "    plot_constraints(observability_grid, constraints, time_grid, target.name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "In the above plot, grid squares that are purple represent times and constraints for which the target is observable, yellow is not observabable.\n",
    "\n",
    "Another way to visualize whether or not a target meets the constraints is with the `observability_table` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "from astroplan import observability_table\n",
    "observability_table(constraints, cape_town, targets, times=time_grid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "### Plots <a id='plots'></a>\n",
    "\n",
    "Let's track that target's motion through the sky for the next ten hours in a plot: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "from astroplan.plots import plot_sky\n",
    "\n",
    "# Plot at times: \n",
    "plot_times = Time.now() + np.linspace(0, 10, 10)*u.hour\n",
    "\n",
    "for target in targets:\n",
    "    plot_sky(target, cape_town, plot_times)\n",
    "plt.legend(loc=[1.2, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Perhaps you need to make a finder chart to help you find Proxima:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "from astroplan.plots import plot_finder_image\n",
    "\n",
    "plot_finder_image(targets[0], fov_radius=10*u.arcmin, \n",
    "                  reticle=True, survey='2MASS-K')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "source": [
    "<br /><br />\n",
    "***\n",
    "<br /><br />\n",
    "\n",
    "# Eclipsing objects <a id='eclipsing'></a>\n",
    "\n",
    "An ever-increasing number of eclipsing systems, such as transiting exoplanets and eclipsing binaries, are being discovered and followed-up. Here's how you can predict exoplanet transists or eclipses in binaries with `astroplan`.\n",
    "\n",
    "First you must define the ephemeris of your eclipsing system. Let's make one for the transiting exoplanet HAT-P-7 b:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# Measurements from Morris et al. 2013\n",
    "transit_epoch = Time(2454954.357463, format='jd', scale='tdb')\n",
    "orbital_period = 2.204737 * u.day \n",
    "transit_duration = 0.16430159 * u.day"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Now we'll create an `EclipsingSystem` object: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "from astroplan import EclipsingSystem\n",
    "\n",
    "hd189733b = EclipsingSystem(transit_epoch, orbital_period, \n",
    "                            duration=transit_duration, \n",
    "                            name='HD 189733 b')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "This object has handy convenience functions for computing the next primary or secondary eclipse: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "next_primary_eclipses = hd189733b.next_primary_eclipse_time(Time.now(), n_eclipses=3)\n",
    "\n",
    "print(\"Next three primary eclipses: {0}\".format(next_primary_eclipses.iso))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "We can combine these `EclipsingSystem` objects with the astroplan constraints framework to determine when transits are observable. Let's compute how many of the next 100 transits of HD 209458 b are observable from Apache Point Observatory:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "primary_eclipse_time = Time(2452826.628514, format='jd')\n",
    "orbital_period = 3.52474859 * u.day\n",
    "eclipse_duration = 0.1277 * u.day\n",
    "\n",
    "hd209458 = EclipsingSystem(primary_eclipse_time=primary_eclipse_time,\n",
    "                           orbital_period=orbital_period, duration=eclipse_duration,\n",
    "                           name='HD 209458 b')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "n_transits = 100  # This is the roughly number of transits per year\n",
    "obs_time = Time.now()\n",
    "midtransit_times = hd209458.next_primary_eclipse_time(obs_time, n_eclipses=n_transits)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "from astroplan import (PrimaryEclipseConstraint, AtNightConstraint, \n",
    "                       AltitudeConstraint, is_event_observable)\n",
    "import datetime as dt\n",
    "import astropy.units as u\n",
    "\n",
    "constraints = [AtNightConstraint.twilight_civil(),\n",
    "               AltitudeConstraint(min=30*u.deg)]\n",
    "\n",
    "apo = Observer.at_site('Apache Point Observatory')\n",
    "target = FixedTarget.from_name('HD 209458')\n",
    "\n",
    "observable = is_event_observable(constraints, apo, target, times=midtransit_times)\n",
    "\n",
    "observable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "Let's make a little table of the transits that can be observed from APO, with some basic observing info:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "observable_transits = midtransit_times[observable.ravel()]\n",
    "\n",
    "altaz = apo.altaz(observable_transits, target)\n",
    "alts = altaz.alt.degree\n",
    "times_utc = observable_transits.iso\n",
    "moon_separation = apo.moon_altaz(observable_transits).separation(altaz).degree\n",
    "moon_illumination = apo.moon_illumination(observable_transits)\n",
    "\n",
    "from astropy.table import Table\n",
    "\n",
    "Table(data=[times_utc, alts, moon_separation, moon_illumination], \n",
    "      names=['Mid-transit Time [UTC]', 'Altitude [deg]', \n",
    "             'Moon sep [deg]', 'Moon illum.'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.5.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}