pyephem, libnova, stellarium, JPL Horizons disagree on moon RA/DEC?

Question

MINOR EDIT: I say below that JPL's Horizons library is not open source. Actually, it is, and it's available here: http://naif.jpl.nasa.gov/naif/tutorials.html

At 2013-01-01 00:00:00 UTC at 0 degrees north latitude, 0 degrees east latitude, sea level elevation, what is the J2000 epoch right ascension and declination of the moon?

Sadly, different libraries give slightly different answers. Converted to degrees, the summarized results (RA first):

Stellarium: 141.9408333000, 9.8899166666 [precision: .0004166640, .0000277777] 
Pyephem: 142.1278749990,  9.8274722221 [precision .0000416655, .0000277777] 
Libnova: 141.320712606865, 9.76909442356909 [precision unknown] 
Horizons: 141.9455833320, 9.8878888888 [precision: .0000416655, .0000277777] 

My question: why? Notes:

• I realize these differences are small, but:

  • I use pyephem and libnova to calculate sun/moon rise/set, and these times can be very sensitive to position at higher latitudes (eg, midnight sun).

  • I can understand JPL's Horizons library not being open source, but the other three are. Shouldn't someone work out the differences in these libraries and merge them? This is my main complaint. Do the stellarium/pyephem/libnova library authors have a fundamental difference in how to make these calculations, or do they just need to merge their code?

• I also realize there might be other reasons the calculations are different, and would appreciate any help in rectifying these possible errors:

  • Pyephem and Libnova may be using the epoch of the date instead of J2000

  • The moon is close enough that observer location can affect its RA/DEC (parallax effect).

  • I'm using Perl's Astro::Nova and Python's pyephem, not the original C implementations of these libraries. However, if these differences are caused by using Perl/Python, that is important in my opinion.

• My code (w/ raw results):

  • First, Perl and Astro::Nova:

#!/bin/perl

# RA/DEC of moon at 0N 0E at 0000 UTC 01 Jan 2013 
use Astro::Nova; 
# 1356998400 == 01 Jan 2013 0000 UTC 
$jd = Astro::Nova::get_julian_from_timet(1356998400); 
$coords = Astro::Nova::get_lunar_equ_coords($jd); 
print join(",",($coords->get_ra(), $coords->get_dec())),"\n"; 

RESULT: 141.320712606865,9.76909442356909 

- Second, Python and pyephem:

#!/usr/local/bin/python 

# RA/DEC of moon at 0N 0E at 0000 UTC 01 Jan 2013 
import ephem; e = ephem.Observer(); e.date = '2013/01/01 00:00:00'; 
moon = ephem.Moon(); moon.compute(e); print moon.ra, moon.dec 

RESULT: 9:28:30.69 9:49:38.9

- The stellarium result (snapshot): 

- The JPL Horizons result (snapshot): 

[JPL Horizons requires POST data (not really, but pretend), so I couldn't post a URL].

• I haven't linked them (lazy), but I believe there are many unanswered questions on stackoverflow that effectively reduce to this question (inconsistency of precision astronomical libraries), including some of my own questions.

• I'm playing w this stuff at: https://github.com/barrycarter/bcapps/tree/master/ASTRO

Answer 1

I have no idea what Stellarium is doing, but I think I know about the other three. You are correct that only Horizons is using J2000 instead of the epoch-of-date for this apparent, locale-specific observation. You can bring it into close agreement with PyEphem by clicking "change" next to the "Table Settings" and switching from "1. Astrometric RA & DEC" to "2. Apparent RA & DEC."

The difference with Libnova is a bit trickier, but my late-night guess is that Libnova uses UT instead of Ephemeris Time, and so to make PyEphem give the same answer you have to convert from one time to the other:

import ephem
moon, e = ephem.Moon(), ephem.Observer()
e.date = '2013/01/01 00:00:00'
e.date -= ephem.delta_t() * ephem.second
moon.compute(e)
print moon.a_ra / ephem.degree, moon.a_dec / ephem.degree

This outputs:

141.320681918 9.77023197401

Which is, at least, much closer than before. Note that you might also want to do this in your PyEphem code if you want it to ignore refraction like you have asked Horizons to; though for this particular observation I am not seeing it make any difference:

e.pressure = 0

Any residual difference is probably (but not definitely; there could be other sources of error that are not occurring to me right now) due to the different programs using different formulae to predict where the planets will be. PyEphem uses the old but popular VSOP87. Horizons uses the much more recent — and exact — DE405 and DE406, as stated in its output. I do not know what models of the solar system the other products use.