pyephem conversion of (Alt, Az) to (Ra, Dec) and back not internally consistent

Question

I am finding that when I convert an (Alt, Az) to an (Ra, Dec) and then back with PyEphem, I don't get what I started with. Below is a simple example.

import ephem
print ephem.__version__
# '3.7.3.4'

gbt = ephem.Observer()
gbt.long = '-79:50:23.4'
gbt.lat = '38:25:59.23'
gbt.pressure = 0 # no refraction correction.
gbt.epoch = ephem.J2000
# Set the date to the epoch so there is nothing changing.
gbt.date = '2000/01/01 12:00:00'

# Should get the north pole right?
ra, dec = gbt.radec_of(0, '38:25:59.23')
# Not the north pole... error might be abberation.
print dec
# 89:59:30.5

# Now check internal consistancy by reversing the calculation.
pole = ephem.FixedBody()
pole._ra = ra
pole._dec = dec
pole._epoch = ephem.J2000
pole.compute(gbt)
# Should get what I started with right?
alt = pole.alt
# Not what I started with... error unknown.
print alt
# 38:26:26.7

As noted in the comments, not getting exactly the north pole might just be stellar aberration, although the 30" is more than Wikipedia stated maximum effect of 20".

The fact that I don't get the same thing when I do the reverse calculation is truly puzzling me. Any suggestions?

Answer 1

The result you are getting is off because of both aberration and nutation. If you were to compile PyEphem yourself and comment out lines lines 271 and 272 of circum.h then you would find that you get exactly the result that you expect — with those edits, the code would look like:

    /* correct EOD equatoreal for nutation/aberation to form apparent 
     * geocentric
     */
    /* nut_eq(mjed, &ra, &dec); */
    /* ab_eq(mjed, lsn, &ra, &dec); */
    op->s_gaera = ra;
    op->s_gaedec = dec;

When you ask PyEphem to "go backwards" from an observed RA and dec to the sky-position behind them, it merely reverses refraction (which you have already turned off) and precession to generate its answer.

Why does it stop there? Why does it not try to reverse nutation and aberration? (Besides the practical reason: that those quantities are estimated through expensive polynomials that cannot be easily reversed!)

The reason why it does not attempt to reverse-compensate for nutation and aberration is that it does not know the range to the object that sits at the RA and dec you are asking about. If you are asking about that RA and dec because you saw a satellite passing overhead, for example, then aberration would be irrelevant — earth satellites travel in the same relativistic frame as the Earth does — and nutation would also be irrelevant, since you would not be interested in where the "ideal" pole of earth points — you would be interested in where the pole was pointing on the particular night that you looked up and observed the satellite passing overhead.

So without knowing whether the object you saw is an earth satellite, or the Moon, or a planet out in the different relativistic frame of the solar system, or something even farther away, the "libastro" library does the simplest thing possible and stops there, rather than mucking up the answer with reverse-effects that might not even apply to your situation. I will keep this in mind, though, as I approach PyEphem's next version, and think about whether multiple radec_of() techniques might not be appropriate.