{
 "metadata": {
  "name": "Asteroids"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import ephem\n",
      "import gzip\n",
      "\n",
      "tau = pi * 2.0\n",
      "white = (1.0, 1.0, 1.0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "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'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "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",
     "collapsed": false,
     "input": [
      "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"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "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'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "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)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "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])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "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)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "imshow(log(heatmap + 1.0), extent=extent)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "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",
      "        )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "jupiter = ephem.Jupiter()\n",
      "jupiter.compute(mjd)\n",
      "\n",
      "print 'Jupiter heliocentric longitude =', jupiter.hlon"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "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])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "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])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "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])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}