Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Differences in ecliptic longitude between Skyfield and JPL's Horizons System for older dates #985

Open
ephilalethes opened this issue Jul 20, 2024 · 4 comments

Comments

@ephilalethes
Copy link

ephilalethes commented Jul 20, 2024

I use e.g. the following code to generate the ecliptic longitudes of Pluto and Juno seen from Earth on a given date and time in 1859:

from skyfield.api import load
from skyfield.data import mpc
from skyfield.constants import GM_SUN_Pitjeva_2005_km3_s2 as GM_SUN

ts = load.timescale()

planets = load('de440.bsp')

sun = planets['sun']
earth = planets['earth_barycenter']
pluto = planets['pluto_barycenter']

with load.open('skytest.txt') as f:
    minor_planets = mpc.load_mpcorb_dataframe(f)

minor_planets = minor_planets.set_index('designation', drop=False)
row = minor_planets.loc['(3) Juno']
juno = sun + mpc.mpcorb_orbit(row, ts, GM_SUN)

t = ts.utc(1859, 8, 23, 15, 20)
print("date and time:", t.utc_strftime())
ep = earth.at(t).observe(pluto).ecliptic_latlon('date')[1].degrees
print("ecliptic longitude of Pluto as seen from Earth, with precession:", ep)
ej = earth.at(t).observe(juno).ecliptic_latlon('date')[1].degrees
print("ecliptic longitude of Juno as seen from Earth, with precession:", ej)

Here the file skytest.txt contains just a single line for Juno:

00003    5.18  0.15 K243V  82.18651  247.77596  169.84069   12.98975  0.2561868  0.22595087   2.6696851  0 E2024-KE9  7606 117 1804-2024 0.63 M-v 3Ek MPCLINUX   4000      (3) Juno               20240529

The output is:

date and time: 1859-08-23 15:20:00 UTC
ecliptic longitude of Pluto as seen from Earth, with precession: 38.50194688084013
ecliptic longitude of Juno as seen from Earth, with precession: 217.3222419715898

However, when I use JPL's Horizons System to get the ecliptic longitudes at those exact moments, I get some differences; for Pluto:

1859-Aug-23 15:20      38.5040854

For Juno:

1859-Aug-23 15:20     216.5040346

As is evident, there's a difference of roughly 0.002 between the values for Pluto, and a whopping 0.8 between the values for Juno.

The Juno one I can only guess has to do with the lack of accuracy in the Kepler orbit method detailed in the documentation, and that JPL's Horizon System uses some other method that is more precise, but I'd be curious to know for sure.

The Pluto one is much stranger to me, because that's part of the ephemeris itself, and while the JPL value is using DE441, there shouldn't really be any difference between that and DE440. I also checked for a couple of the planets, and found a roughly 0.008 difference for Mercury and 0.004 difference for Jupiter.

Is there a clear answer for these inaccuracies? I guess what I'm wondering is if it's inherent to the precision of Skyfield itself, or whether there's something else strange going on. The 0.002 difference for Pluto above amounts to roughly 7.2 arcseconds, which while ultimately not very big still seems significantly larger than I'd expect from any imprecisions in bodies supported by the ephemerides themselves.

@brandon-rhodes
Copy link
Member

Take a look at the complete output from the HORIZONS system — or, if you like, cut and paste it here. My guess is that they are using UT1, based on the guessed historical orientation of the Earth during ancient times, whereas your Skyfield code is using UTC, which had no leap seconds in previous centuries and thus acted as a simple atomic clock.

@ephilalethes
Copy link
Author

ephilalethes commented Jul 21, 2024

Take a look at the complete output from the HORIZONS system — or, if you like, cut and paste it here. My guess is that they are using UT1, based on the guessed historical orientation of the Earth during ancient times, whereas your Skyfield code is using UTC, which had no leap seconds in previous centuries and thus acted as a simple atomic clock.

That appears to be entirely correct; the output does indeed say:

Times PRIOR to 1962 are UT1, a mean-solar time closely related to the
prior but now-deprecated GMT. Times AFTER 1962 are in UTC, the current
civil or "wall-clock" time-scale. UTC is kept within 0.9 seconds of UT1
using integer leap-seconds for 1972 and later years.

I'm not sure I entirely understand the implications of this though. The number of leap seconds in that time period doesn't seem to account for the differences, especially not the roughly 0.8 degree difference in the ecliptic longitude of Juno.

That aside, does this mean that there is any benefit to using UT1 instead of UTC for dates before 1962 when it comes to precision in general? I suppose the question is more about which to use if I see some old date (like the one in question) and I want to know what the ecliptic longitude actually was on that date, does UT1 coincide better with it?

As for this issue, I just tried running the above with t = ts.ut1(1859, 8, 23, 15, 20) instead, and the result appears to be almost exactly the same; the amount of seconds of difference between the two doesn't seem to impact the relevant precision at all:

date and time: 1859-08-23 15:20:00 UT1
ecliptic longitude of Pluto as seen from Earth, with precession: 38.50194886393843
ecliptic longitude of Juno as seen from Earth, with precession: 217.32215222137896

I also read that the MPC orbital elements being used should have precision to five digits, so I'm ultimately still confused not just about the smaller difference for Pluto (and for the main planets, like Mercury and Jupiter), but also the relatively large difference for Juno.

Addendum:

From looking more at the Horizons output, it occurred to me that they were inherently using apparent positions rather than astrometric ones for the ecliptic longitude; using apparent() instead I get:

date and time: 1859-08-23 15:20:00 UT1
ecliptic longitude of Pluto as seen from Earth, with precession: 38.5041606254349
ecliptic longitude of Juno as seen from Earth, with precession: 217.31988130429593

Here at least Pluto has been brought into agreement to 3 decimal places, which I can be happy with (although I'm curious why it's different beyond that, trying the same with UTC didn't really make a significant difference, the two are still roughly 0.00008 off, but that doesn't matter as much to me, I guess it could have to do with differences in assumptions about atmospheric refraction when it's that small).

Juno on the other hand still appears to be more or less just as wrong (around 0.8 degrees off), so presumably there's something about the way Horizons computes its position relative to Skyfield; perhaps you or someone else might know.

Another thing I just found out is that Stellarium has the ecliptic longitude at around 217.63 degrees at the same date and time, which is closer to the Skyfield value, but confusingly off by around 0.3 degrees in the other direction instead.

Looking further into it I checked for the other three largest small bodies, Ceres, Pallas, and Vesta, and they seem to be even further off. Here are the ecliptic longitudes given by Horizons:

Ceres
1859-Aug-23 15:20     347.2725205

Pallas
1859-Aug-23 15:20     314.5696718

Vesta
1859-Aug-23 15:20      19.2065536

And here are the ones found with Skyfield, for each combination of astrometric/apparent and with/without precession:

date and time: 1859-08-23 15:20:00 UT1

ASTROMETRIC

with precession
ecliptic longitude of Ceres as seen from Earth: 344.73079032198683
ecliptic longitude of Pallas as seen from Earth: 309.4286171228486
ecliptic longitude of Juno as seen from Earth: 217.32215222137896
ecliptic longitude of Vesta as seen from Earth: 8.732412577670233

without precession
ecliptic longitude of Ceres as seen from Earth: 346.69256032731795
ecliptic longitude of Pallas as seen from Earth: 311.3776741611754
ecliptic longitude of Juno as seen from Earth: 219.28144731078157
ecliptic longitude of Vesta as seen from Earth: 10.692461485856374

APPARENT

with precession
ecliptic longitude of Ceres as seen from Earth: 344.73646204842953
ecliptic longitude of Pallas as seen from Earth: 309.4346864486547
ecliptic longitude of Juno as seen from Earth: 217.31988130429593
ecliptic longitude of Vesta as seen from Earth: 8.736923786463954

without precession
ecliptic longitude of Ceres as seen from Earth: 346.69823200260623
ecliptic longitude of Pallas as seen from Earth: 311.3837423935764
ecliptic longitude of Juno as seen from Earth: 219.27917673002062
ecliptic longitude of Vesta as seen from Earth: 10.696972416002792

These are many degrees off even in the best cases (Vesta being ~10° off), and strangely it seems like the correct combination of accounting for precession and using the apparent position (which is what Horizons uses) is even wronger than other combinations in some cases.

I also checked e.g. Vesta on a contemporary date, and in that case the Horizons value agrees almost perfectly with the Skyfield value for apparent position and accounting for precession; Horizons:

2024-Jul-20 12:00     133.5711239

Skyfield:

ecliptic longitude of Vesta as seen from Earth: 133.571096822202

So it seems to be a problem that creeps in over longer spans of time, but I can't really figure out what is causing these huge discrepancies. I'm of course assuming Horizons is correct and Skyfield is doing something inaccurate, but I guess I don't know, maybe Skyfield is correct and Horizons isn't, heh.

@ephilalethes
Copy link
Author

@brandon-rhodes:

From what I gather JPL uses something called a "small-body perturber" (SB441-N16); do you think that could that account for the large differences observed over time? On the Wikipedia page for the JPL Small-Body Database it also says:

Most objects such as asteroids get a two-body solution (Sun+object) recomputed twice a year. Comets generally have their two-body orbits computed at a time near the perihelion passage (closest approach to the Sun) as to have the two-body orbit more reasonably accurate for both before and after perihelion. For most asteroids, the epoch used to define an orbit is updated twice a year. Orbital uncertainties in the JPL Small-Body Database are listed at the 1-sigma level.

@brandon-rhodes
Copy link
Member

That aside, does this mean that there is any benefit to using UT1 instead of UTC for dates before 1962 when it comes to precision in general? I suppose the question is more about which to use if I see some old date (like the one in question) and I want to know what the ecliptic longitude actually was on that date, does UT1 coincide better with it?

If you are interested in older human experience and human records, you will want to use UT1, so that your "days" correspond to the sunrise/sunset days actually experienced by earlier peoples. They didn't know the Earth's rotation was slowing, or that someday we'd invent atomic clocks, so UTC times would be meaningless to them and not correspond to their own idea of a "day".

From what I gather JPL uses something called a "small-body perturber" (SB441-N16); do you think that could that account for the large differences observed over time?

Yes, that could make an important difference. Try using HORIZONS with the first setting as "Ephemeris Type: Osculating Orbital Elements" — I think that will tell you what instantaneous orbital elements they're using for Juno at a distant date in the past? Then you could try giving those same elements to Skyfield.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

2 participants