brandon-rhodes
commented
3 months ago Thank you for providing a working notebook! I was pleasantly surprised to find a "play" button I could press and re-run the example right in my browser, which made diagnosis very quick. I hadn't used colab.research.google.com before.
The problem is that you are using two different coordinate systems. By default Skyfield uses the modern ICRS = J2000 coordinate system. But the old PyEphem library uses, for apparent coordinates, the dynamical precessing reference system of the Earth's mean equator and pole. Try telling Skyfield to apply precession, like this:
ra, dec, distance = apparent.radec(epoch='date')Does that fix the inconsistency in your copy of the notebook as well?
alexanderflint
While working on a Nautical Almanac calculator, I've noticed that the RA and Decl calculations for any body (Sun, planet, star) are different when using Skyfield vs when using pyephem.
It looks like pyephem yields the correct result when compared to GHA and Decl in a printed Nautical Almanac [ and also the tables from nauticalalmanac.com, but that's not surprising because they use pyephem for their calculations. ] While the differences aren't massive, they are on the minutes scale, so using the Skyfield version would put each body's geographic position off by a couple nautical miles, and the resulting fix from several bodies potentially off by more than that.
For example, here are Skyfield's results for Mars on 1 Jan 2024 00:00:00h GMT:
Sight: Mars Sight Datetime: 2024/1/1 0:0:0 RA = 17h 46m 45.64s Dec = -23° 57.1m
And here are pyephem's results for Mars at the same datetime:
Sight: Mars Sight Datetime: 2024/1/1 0:0:0 RA = 17h 48m 13.1s Decl = -23° 57.7m
Using the pyephem version to convert RA to GHA, you get:
Sight: Mars Sight Datetime: 2024/1/1 0:0:0 GHA: 193° 5.8, Decl: -23° 57.7
And this exactly matches the printed Nautical Almanac, so pyephem appears to be correct.
The Skyfield version of the calculation uses the usual:
eph = load('de421.bsp') earth, mars = eph[3], eph[4] sight_datetime = f'{sight_year}/{sight_month}/{sight_day} {sight_hour}:{sight_minute}:{sight_second}' barycentric = earth.at(t) astrometric = barycentric.observe(mars) apparent = astrometric.apparent() ra, dec, distance = apparent.radec()
Colab Jupyter notebook with example code: https://colab.research.google.com/drive/1SPg3wAb_0JSG6jwim_rILLlpi5DUZ2Sy?usp=sharing