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