Asteroids.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import ephem\n",
    "import gzip\n",
    "\n",
    "tau = pi * 2.0\n",
    "white = (1.0, 1.0, 1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "f = gzip.open('data/ELEMENTS.NUMBR.gz')\n",
    "lines = iter(f)\n",
    "header = next(lines)\n",
    "dashes = next(lines)\n",
    "\n",
    "bodies = []\n",
    "\n",
    "for line in f:\n",
    "    fields = line[24:].split()\n",
    "    b = ephem.EllipticalBody()\n",
    "    b._epoch = int(fields[0])     # \"Epoch\"\n",
    "    b._a = float(fields[1])       # \"a\"\n",
    "    b._e = float(fields[2])       # \"e\"\n",
    "    b._inc = float(fields[3])     # \"i\"\n",
    "    b._om = float(fields[4])      # \"w\"\n",
    "    b._Om = float(fields[5])      # \"Node\"\n",
    "    b._M = float(fields[6])       # \"M\"\n",
    "    b._epoch_M = b._epoch  # complete guess\n",
    "    b._g = float(fields[7])\n",
    "    bodies.append(b)\n",
    "    \n",
    "print 'Loaded', len(bodies), 'asteroids'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How did the number of asteroids known to our species become so large?\n",
    "Because since the first LINEAR system started operating in 1998,\n",
    "asteroid discovery has been an automated process\n",
    "that now involves dozens of automated telescopes across the world.\n",
    "\n",
    "Here is a video animation of our asteroid discoveries so far:\n",
    "\n",
    "http://www.youtube.com/watch?v=S_d-gs0WoUw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def compute_asteroid_positions(mjd):\n",
    "    for b in bodies:\n",
    "        b.compute(mjd)\n",
    "    theta = np.array([ b.hlon for b in bodies ])\n",
    "    r = np.array([ b.sun_distance for b in bodies ])\n",
    "    return theta, r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mjd = ephem.Date('2013/1/11 15:45')\n",
    "theta, r = compute_asteroid_positions(mjd)\n",
    "\n",
    "print 'First asteroid is at angle', theta[0]\n",
    "print 'and at distance', r[0], 'AU from the Sun'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def generate_planet_positions(cls, mjd0, mjd1):\n",
    "    planet = cls()\n",
    "    dates = np.arange(mjd0, mjd1, 1.0)\n",
    "    for date in dates:\n",
    "        planet.compute(date)\n",
    "        yield (planet.hlon, planet.sun_distance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "subplot = figure().add_subplot(111, polar=True)\n",
    "\n",
    "subplot.scatter(theta, r, s=0.5)\n",
    "\n",
    "for cls in ephem.Jupiter, ephem.Mars, ephem.Saturn:\n",
    "    g = generate_planet_positions(cls, mjd - 30, mjd + 30)\n",
    "    theta_r = np.array(list(g)).T\n",
    "    subplot.plot(theta_r[0], theta_r[1], linewidth=8)\n",
    "\n",
    "subplot.axis([None, None, 0.0, 12.0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x =  r * cos(theta)\n",
    "y = -r * sin(theta)\n",
    "\n",
    "heatmap, xedges, yedges = np.histogram2d(\n",
    "    y, x, bins=200, range=[[-7, 7], [-7, 7]])\n",
    "\n",
    "extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]\n",
    "imshow(heatmap, extent=extent)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "imshow(log(heatmap + 1.0), extent=extent)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n, bins, patches = hist(r, bins=60, range=[0.0, 12.0])\n",
    "\n",
    "for cls in ephem.Mars, ephem.Jupiter, ephem.Saturn:\n",
    "    planet = cls()\n",
    "    planet.compute(mjd)\n",
    "    annotate(\n",
    "        planet.name,\n",
    "        xy=(planet.sun_distance, 30000),\n",
    "        xytext=(planet.sun_distance, 50000),\n",
    "        arrowprops=dict(facecolor='black', shrink=0.2),\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "jupiter = ephem.Jupiter()\n",
    "jupiter.compute(mjd)\n",
    "\n",
    "print 'Jupiter heliocentric longitude =', jupiter.hlon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n, bins, patches = hist(theta, bins=48)\n",
    "jpos = (jupiter.hlon, 1500.0)\n",
    "annotate('J', xy=jpos, fontsize=24, color=white)\n",
    "axis([0.0, tau, None, None])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n, bins, patches = hist(theta[r < 3.0], bins=48)\n",
    "jpos = (jupiter.hlon, 850.0)\n",
    "annotate('J', xy=jpos, fontsize=24, color=white)\n",
    "axis([0.0, tau, None, None])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "n, bins, patches = hist(theta[r > 4.5], bins=48)\n",
    "jpos = (jupiter.hlon, 200.0)\n",
    "annotate('J', xy=jpos, fontsize=24)\n",
    "axis([0.0, tau, None, None])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}