PyEphem ephemeris differs from that returned by MPC

[blackhaz]

For some reason I am getting slightly different ephemerides from PyEphem, if compared to those returned by MPC. Example object: tracked_object = ephem.readdb("89958 2002 LY45,e,9.9092,188.1920,222.9011,1.641425,0.4686756,0.886287,250.1156,12/17.0/2020,2000,H17.00,0.15") Observatory code is U69 (long 240.5870 lat +37.070390, height 1376) PyEphem returns: 2021 06 26 082000 16 23 04.9 +51 05 36 Official MPC ephemeris: 2021 06 26 082000 16 22 36.3 +51 10 34

Any ideas where this discrepancy might be coming from? Coordinates used are the same in both cases. Many thanks! Max

PS: SkyField produces more accurate ephemeris.

[brandon-rhodes]

Thanks for the update, and for circling back on the issue so quickly!

[blackhaz]

Hi Brandon, Sorry, it appears I was too quick. SkyField also differs from MPC, gives me this: 2021 06 26 082000 16 22 24.7 +51 09 02. It's better but still far from the object. On my images from a 0.61 m telescope the object appears closer to the ephemeris from MPC. Any ideas where this discrepancy might be coming from?

[brandon-rhodes]

There are several possible sources of discrepancy. But to determine which one might be involved, I would need to see the Skyfield code you are using to produce those coordinates, and how you are using your telescope to produce coordinates?

[blackhaz]

Thank you for your time, Brandon. Here is the code used:

from skyfield.api import Star, load, wgs84
from skyfield.constants import GM_SUN_Pitjeva_2005_km3_s2 as GM_SUN
from skyfield.data import mpc
from astropy import units as u
from astropy.coordinates import SkyCoord

ts = load.timescale()
eph = load('de421.bsp')
sun = eph['sun']
earth = eph['earth']
with load.open('mpcorb.dat.target') as f:
    minor_planets = mpc.load_mpcorb_dataframe(f)
row = minor_planets.iloc[0]
tracked_object = sun + mpc.mpcorb_orbit(row, ts, GM_SUN)
observer = wgs84.latlon(240.5870, 37.070390, elevation_m=1376.643) # U69
t = ts.utc(2021, 6, 26, 0, 0, 0)
ra, dec, distance = (earth + observer).at(t).observe(tracked_object).radec()
print(ra,dec)

I am getting 16h 31m 05.19s +50deg 02' 25.0" as the output.

MPC returns via https://www.minorplanetcenter.net/iau/MPEph/MPEph.html:

89958              [H=17.00]
Date       UT      R.A. (J2000) Decl.    Delta     r     El.    Ph.   V      Sky Motion        Object    Sun   Moon                Uncertainty info
            h m s                                                            "/min    P.A.    Azi. Alt.  Alt.  Phase Dist. Alt.    3-sig/" P.A.
2021 06 26 000000 16 31 11.5 +50 04 01   0.196   1.077  102.8  67.0  16.0   12.55    310.5    225  +26   +37   0.98   085  -52       N/A   N/A / Map / Offsets

Note almost 2' off in declination.

[blackhaz]

Sorry, forgot the contents of mpcorb.dat.target which is the latest elements for 89958 from mpcorb.dat:

89958   17.00  0.15 K20CH 250.11558  222.90105  188.19200    9.90917  0.8862866  0.46867557   1.6414255  0 MPO594447   562  16 2000-2020 0.55 M-v 3Ek Pan        9803  (89958) 2002 LY45          20200516

[brandon-rhodes]

Thanks for sample code that I can run! (I've edited your comments to use triple-backticks around the code, so I can cut and paste it successfully). I don't see any big problems with your code at first glance — you're asking for UTC correctly, and are asking for J2000 coordinates not equinox-of-date.

Is there any chance that the MPC is perturbing the orbital parameters before applying them to the specific date you're asking about? Is there a checkbox in the interface you're using that would let you turn that off, if it's turned on? That can be one source of discrepancy.

[blackhaz]

Thank you for looking. In MPC's interface there is only an option to "generate perturbed ephemerides for unperturbed orbits" but they don't provide an option to generate unperturbed ephemerides.

[brandon-rhodes]

Another option to see if they are perturbing the elements might be for you to select the "Also display elements for epoch" field, enter the date for which you are generating the position, and compare the elements it returns to the elements from the latest mpcorb.dat file.

If there's no difference in the elements, then I'll have to search elsewhere for the discrepancy.

But if you do indeed see that the elements are different, then we will possibly have found the source of the difference!

[blackhaz]

Great idea. Here's what I've got from MPC: Elements at other epoch(s): Epoch (TT/JD) T (TT/JD) Peri. Node Incl. e q a Epoch (TT) 2459391.50000 2459434.9522715 222.90885 188.17803 9.90833 0.8862699 0.1867081 1.64168 2021/06/26.00

And here is the line in mpcorb.dat: 89958 17.00 0.15 K20CH 250.11558 222.90105 188.19200 9.90917 0.8862866 0.46867557 1.6414255 0 MPO594447 562 16 2000-2020 0.55 M-v 3Ek Pan 9803 (89958) 2002 LY45 20200516

Looks like they are slightly different indeed! I wonder, does that mean that the mpcorb.dat file is lagging?

[brandon-rhodes]

My guess is that the mpcorb.dat file contains the most recent "real" orbit, produced by researchers integrating real observations of the body. While the service is willing to produce synthetic perturbed orbits for the purpose of generating more recent positions, they're not willing to overwrite real observation-derived orbits in mpcorb.dat with speculative machine-generated perturbed orbits.

I wonder if their software for perturbation is available anywhere, or at least described in a publication?

[blackhaz]

I think we can close the case. I have used the latest orbital elements and SkyField puts it exactly where the asteroid is. I must apologize for raising a false alarm. Everything works perfectly.

[brandon-rhodes]

Thanks for trying out the latest orbital elements, and I'm glad the case is now closed! Let me know if you run across any guidance on how the MPC perturbs orbits; if we could include that information in Skyfield's documentation, we might avoid other folks running into the same issue you did.

[blackhaz]

Absolutely, will keep this in mind. Thank you for your time on this. I have tried searching for this, but found no information.