Inconsistent TLE results with Astropy/SGP4

[clarktristan]

Why are the TEME vectors off with this code? They are off by tens of kms. The julian dates are slightly off, and i can't seem to get them to align

from skyfield.api import EarthSatellite, load
from astropy.time import Time
from sgp4.api import Satrec

tle1 = '1 25544U 98067A   19343.69339541  .00001764  00000-0  38792-4 0  9991'
tle2 = '2 25544  51.6439 211.2001 0007417  17.6667  85.6398 15.50103472202482'

satellite1 = Satrec.twoline2rv(tle1, tle2)
t1 = Time(satellite1.jdsatepoch, format='jd')
error_code, teme_p, teme_v = satellite1.sgp4(t1.jd1, t1.jd2)  # in km and km/s
print(t1, teme_p, teme_v)

satellite2 = EarthSatellite(tle1, tle2)
t2 = load.timescale().from_astropy(t1)
v = satellite2.at(t2).velocity.km_per_s
teme = satellite2.at(t2).position.km
print(t2, teme, v)

[brandon-rhodes]

1. You cannot use satellite1.jdsatepoch by itself. The sgp4 documentation explains that it is merely the whole-day part of the epoch; you need to add in satellite1.jdsatepochF which is the corresponding fraction, expressing (as a fraction of a day) how many hours and minutes and seconds the epoch has.
2. The .position.km is a GCRS position, not a TEME position, so I wouldn't expect them to match closely. Let me know where I might be able to more thoroughly document that fact so that users are more likely to find it when reading. Thanks!

[brandon-rhodes]

Good morning! I wanted to follow up and see if there were any steps for Skyfield to take here — like updating or improving documentation. If not, then hopefully your code is working better. I'll plan on closing this issue late in the week.

[brandon-rhodes]

As we haven't heard back in more than a week, I hope that means that your issue is resolved. If not, simply comment here with a further update on where things stand, and we can continue the discussion.

[russdot]

Hello! Related to the original post, I'm having some trouble constructing skyfield.timelib.Time object using jdsatepoch and jdsatepochF components.

What is the proper way to do that? I'm still a bit new to the different timescales, and it's not quite clear how I can accomplish it. (or, equivalently with sgp4.api.jday or skyfield.timelib.julian_date presumably).

I notice in the original question the author is using astropy.time.Time, which appears to work just fine.

from datetime import datetime, timezone
from astropy.time import Time
from skyfield.timelib import julian_date
from skyfield.api import load

ts = load.timescale()

def test_jday_to_astropy_utc():

    date_jd = ts.from_datetime(datetime(2023, 7, 19, 14, 14, 17, tzinfo=timezone.utc))

    jd = julian_date(2023, 7, 19, 14, 14, 17)
    astropy_t = Time(jd, format='jd')
    assert ts.from_astropy(astropy_t).utc_iso() == date_jd.utc_iso() # Succeeds

    assert ts.tt_jd(jd).utc_iso() == date_jd.utc_iso() # Fails with a difference of 69 seconds

Am I missing something?

Thanks for this fantastic package!