Utils-HORIZONS-data.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "How to fetch current physical data for the planets from the HORIZONS system.\n",
      "\n",
      "This data is then cut-and-pasted into Skyfield."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from telnetlib import Telnet\n",
      "\n",
      "def fetch_planet_data():\n",
      "    t = Telnet('horizons.jpl.nasa.gov', 6775)\n",
      "    for n in [10] + [100*i + 99 for i in range(1, 10)]:\n",
      "        t.read_until('Horizons>')\n",
      "        t.write('%s\\r\\n' % n)\n",
      "        yield t.read_until('<cr>:')\n",
      "        t.write('\\r\\n')\n",
      "    t.write('q\\r\\n')\n",
      "    print t.read_all()\n",
      "\n",
      "texts = list(fetch_planet_data())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " \r\n",
        " \r",
        "\r\n",
        "Horizons> q\r\n",
        "                                                                          \r",
        "\r\n",
        "     ___    _____     ___                                                 \r",
        "\r\n",
        "    /_ /|  /____/ \\  /_ /|       Horizons On-line Ephemeris System v3.85  \r",
        "\r\n",
        "    | | | |  __ \\ /| | | |       Solar System Dynamics Group              \r",
        "\r\n",
        " ___| | | | |__) |/  | | |__     Jet Propulsion Laboratory                \r",
        "\r\n",
        "/___| | | |  ___/    | |/__ /|   Pasadena, CA, USA                        \r",
        "\r\n",
        "|_____|/  |_|/       |_____|/                                             \r",
        "\r\n",
        " \r",
        "\r\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "There!\n",
      "\n",
      "With that simple routine complete,\n",
      "we have completed the data fetch,\n",
      "closed the open Telnet connection to the JPL,\n",
      "and can now process the data in peace without bothering them with a connection dangling open."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Here is what the block for the Sun looks like.\n",
      "\n",
      "print texts[0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 10\r\n",
        "\u001b[?1h\u001b=\r",
        "*******************************************************************************\r\n",
        " Revised : Jul 31, 2013                  Sun                                 10\r\n",
        "\r\n",
        " PHYSICAL PROPERTIES (revised Jan 16, 2014):\r\n",
        "  GM (10^11 km^3/s^2)   = 1.3271244004193938  Mass (10^30 kg)   ~ 1.988544\r\n",
        "  Radius (photosphere)  = 6.963(10^5) km  Angular diam at 1 AU  = 1919.3\"\r\n",
        "  Solar Radius (IAU)    = 6.955(10^5) km  Mean density          = 1.408 g/cm^3\r\n",
        "  Surface gravity       =  274.0 m/s^2    Moment of inertia     = 0.059\r\n",
        "  Escape velocity       =  617.7 km/s     Adopted sidereal per  = 25.38 d\r\n",
        "  Pole (RA,DEC in deg.) =  286.13,63.87   Obliquity to ecliptic = 7 deg 15'        \r\n",
        "  Solar constant (1 AU) = 1367.6 W/m^2    Solar lumin.(erg/s)   =  3.846(10^33)\r\n",
        "  Mass-energy conv rate = 4.3(10^12 gm/s) Effective temp (K)    =  5778\r\n",
        "  Surf. temp (photosphr)= 6600 K (bottom) Surf. temp (photosphr)=  4400 K (top)\r\n",
        "  Photospheric depth    = ~400 km         Chromospheric depth   = ~2500 km\r\n",
        "  Sunspot cycle         = 11.4 yr         Cycle 22 sunspot min. =  1991 A.D.\r\n",
        "\r\n",
        "  Motn. rel to nrby strs= apex : RA=271 deg; DEC=+30 deg\r\n",
        "                          speed: 19.4 km/s = 0.0112 AU/day\r\n",
        "  Motn. rel to 2.73K BB = apex : l=264.7+-0.8; b=48.2+-0.5\r\n",
        "                          speed: 369 +-11 km/s\r\n",
        "*******************************************************************************\r\n",
        "\r",
        "\u001b[K\u001b[?1l\u001b> Select ... [E]phemeris, [F]tp, [M]ail, [R]edisplay, ?, <cr>:\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# And one of the blocks for a planet.\n",
      "\n",
      "print texts[1]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 199\r\n",
        "\u001b[?1h\u001b=\r",
        "*******************************************************************************\r\n",
        " Revised: Jul 31, 2013              Mercury                             199 / 1\r\n",
        "\r\n",
        " GEOPHYSICAL DATA (updated 2008-Feb-07):\r\n",
        "  Mean radius (km)      =  2440(+-1)     Density (g cm^-3)     =  5.427\r\n",
        "  Mass (10^23 kg )      =     3.302      Flattening, f         =  \r\n",
        "  Volume (x10^10 km^3)  =     6.085      Semi-major axis       = \r\n",
        "  Sidereal rot. period  =    58.6462 d   Rot. Rate (x10^5 s)   =  0.124001\r\n",
        "  Mean solar day        =   175.9421 d   Polar gravity ms^-2   =\r\n",
        "  Mom. of Inertia       =     0.33       Equ. gravity  ms^-2   =  3.701     \r\n",
        "  Core radius (km)      = ~1600          Potential Love # k2   =\r\n",
        "\r\n",
        "  GM (km^3 s^-2)        = 22032.09       Equatorial Radius, Re =    2440 km\r\n",
        "  GM 1-sigma (km^3 s^-2)=   +-0.91       Mass ratio (sun/plnt) = 6023600\r\n",
        " \r\n",
        "  Atmos. pressure (bar) =                Max. angular diam.    = 11.0\"   \r\n",
        "  Mean Temperature (K)  =                Visual mag. V(1,0)    = -0.42 \r\n",
        "  Geometric albedo      =   0.106        Obliquity to orbit[1] = 2.11' +/- 0.1'  \r\n",
        "  Sidereal orb. per.    =   0.2408467 y  Mean Orbit vel.  km/s = 47.362 \r\n",
        "  Sidereal orb. per.    =  87.969257  d  Escape vel. km/s      =  4.435\r\n",
        "  Hill's sphere rad. Rp =  94.4          Planetary Solar Const = 9936.9 (Wm^2)\r\n",
        "\r\n",
        "[1] Margot et al., Science 316, 2007\r\n",
        "*******************************************************************************\r\n",
        "\r",
        "\u001b[K\u001b[?1l\u001b> Select ... [E]phemeris, [F]tp, [M]ail, [R]edisplay, ?, <cr>:\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def parse(text):\n",
      "    for line in text.splitlines():\n",
      "        if line.strip().startswith('Revised'):\n",
      "            name = line[28:60].strip()\n",
      "            yield 'name', name\n",
      "            continue\n",
      "        halfway = line.find(' ', 40)\n",
      "        halves = [line] if (halfway == -1) else [line[:halfway], line[halfway:]]\n",
      "        for half in halves:\n",
      "            if '\\x1b' in half:\n",
      "                continue\n",
      "            if '=' in half:\n",
      "                name, value = half.rsplit('=', 1)\n",
      "            elif '~' in half:\n",
      "                name, value = half.rsplit('~', 1)\n",
      "            else:\n",
      "                continue\n",
      "            yield name.strip(), value.strip()\n",
      "\n",
      "data = [list(parse(text)) for text in texts]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "masses = []\n",
      "\n",
      "for items in data:\n",
      "    for k, v in items:\n",
      "        if k == 'name':\n",
      "            name = v\n",
      "            continue\n",
      "        if 'Mass' not in k:\n",
      "            continue\n",
      "        if '10^' not in k:\n",
      "            continue\n",
      "        exponent = k.split('^')[1].split()[0]\n",
      "        mantissa = v.split('+-')[0]\n",
      "        mass = float(mantissa) * 10.0 ** float(exponent)\n",
      "        masses.append((name, mass))\n",
      "\n",
      "print len(masses)\n",
      "masses"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "[('Sun', 1.9885440000000002e+30),\n",
        " ('Mercury', 3.302e+23),\n",
        " ('Venus', 4.8685e+24),\n",
        " ('Earth', 5.97219e+24),\n",
        " ('Mars', 6.4185e+23),\n",
        " ('Jupiter', 1.89813e+27),\n",
        " ('Saturn', 5.68319e+26),\n",
        " ('Uranus', 8.681029999999999e+25),\n",
        " ('Neptune', 1.0241e+26),\n",
        " ('134340 Pluto', 1.3069999999999998e+22)]"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "radii = []\n",
      "\n",
      "for items in data:\n",
      "    for k, v in items:\n",
      "        if k == 'name':\n",
      "            name = v\n",
      "            continue\n",
      "        if 'adius' not in k:\n",
      "            continue\n",
      "        if 'mean' not in k.lower() and 'Pluto' not in k and 'Solar' not in k:\n",
      "            continue\n",
      "        print k, v\n",
      "        radius_km = float(v.split()[0].split('+-')[0].split('(')[0])\n",
      "        if '10^' in v:\n",
      "            radius_km *= 10 ** int(v.split('10^')[1].split(')')[0])\n",
      "        radii.append((name, radius_km))\n",
      "\n",
      "print len(radii)\n",
      "radii"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Solar Radius (IAU) 6.955(10^5) km\n",
        "Mean radius (km) 2440(+-1)\n",
        "Mean radius (km) 6051.8(4+-1)\n",
        "Mean radius, km 6371.01+-0.01\n",
        "Mean radius (km) 3389.9(2+-4)\n",
        "Volumetric mean radius 69911+-6 km\n",
        "Volumetric mean radius 58232+-6 km\n",
        "Volumetric mean radius 25362+-12 km\n",
        "Volumetric mean radius 24624+-21 km\n",
        "Radius of Pluto, Rp 1195 km\n",
        "10\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "[('Sun', 695500.0),\n",
        " ('Mercury', 2440.0),\n",
        " ('Venus', 6051.8),\n",
        " ('Earth', 6371.01),\n",
        " ('Mars', 3389.9),\n",
        " ('Jupiter', 69911.0),\n",
        " ('Saturn', 58232.0),\n",
        " ('Uranus', 25362.0),\n",
        " ('Neptune', 24624.0),\n",
        " ('134340 Pluto', 1195.0)]"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    }
   ],
   "metadata": {}
  }
 ]
}