Axis-Wobble.ipynb

{
 "metadata": {
  "name": ""
 },
 "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 pylab as pl\n",
      "import numpy as np\n",
      "import math\n",
      "\n",
      "tau = 2.0 * math.pi\n",
      "halfpi = 0.5 * math.pi\n",
      "\n",
      "degree = tau / 360.0\n",
      "minute = degree / 60.\n",
      "second = minute / 60."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import ephem"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "pole = ephem.FixedBody()\n",
      "pole._ra = '0'    # hours\n",
      "pole._dec = '90'  # degrees\n",
      "\n",
      "def pole_positions(dates):\n",
      "    theta = np.zeros(len(dates))\n",
      "    r = np.zeros(len(dates))\n",
      "    for i, d in enumerate(dates):\n",
      "        pole.compute(d)\n",
      "        theta[i] = pole.ra\n",
      "        r[i] = (halfpi - pole.dec) / degree\n",
      "    return theta, r"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "start_date = ephem.Date((-162, 1, 1))\n",
      "end_date = ephem.Date((2200, 1, 1))\n",
      "\n",
      "print float(start_date), float(end_date)\n",
      "\n",
      "dates = np.linspace(start_date, end_date, 365.0)\n",
      "theta, r = pole_positions(dates)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot = pl.figure().add_subplot(111, polar=True)\n",
      "plot.plot(theta, r)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# When will it be closest?\n",
      "# TODO: this makes no sense; re-do this as Polaris later?\n",
      "\n",
      "print min(theta)\n",
      "print np.where(theta==min(theta))\n",
      "print dates[theta==min(theta)]\n",
      "\n",
      "closest_dates = dates[theta==min(theta)]\n",
      "print ephem.Date(closest_dates[0])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "start_date = ephem.Date((1950, 1, 1))\n",
      "end_date = ephem.Date((2050, 1, 1))\n",
      "\n",
      "dates = range(int(start_date), int(end_date))\n",
      "theta, r = pole_positions(dates)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot = pl.figure().add_subplot(111, polar=True)\n",
      "plot.plot(theta, r)\n",
      "plot.axis([None, None, 0.0, 0.20])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot = pl.figure().add_subplot(111, polar=True)\n",
      "plot.plot(theta, r)\n",
      "plot.axis([None, None, 0.0, 0.02])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot = pl.figure().add_subplot(111)\n",
      "plot.plot(theta, r)\n",
      "plot.axis([89.0 * degree, 91.0 * degree,\n",
      "           0.0, 0.1])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}