Epoch is half-day too early for EarthSatellites built from_satrec

[chylu]

First of all, thanks Brandon for making such a thorough, powerful astronomy library. It's simplified my job tremendously.

Recently I've had the need to create EarthSatellites from my own orbital elements. Per the documentation, I've made a Satrec object and used sgp4init, then used EarthSatellite.from_satrec. To test this method, I tried using this method on an existing satellite with a TLE to compare results.

I found that the model calculates position and velocity correctly, but the purported epoch differs by exactly a half-day. I believe the issue is probably in the calculation of epoch in the EarthSatellite.from_satrec function.

Here is some sample code you can run to reproduce this issue:

from sgp4.api import Satrec, WGS72, days2mdhms
from sgp4.ext import jday

from skyfield.api import EarthSatellite, load

import numpy as np

TIMESCALE = load.timescale()
DEGREES_PER_RADIANS = 360 / (2 * np.pi)

name = "FORMOSAT 7-3"
tle1 = "1 44343U 19036E   20286.64549342  .00000777  00000-0  27989-4 0 9992"
tle2 = "2 44343  24.0014  31.6920 0010169 255.8372 104.1002 15.07436211 69802"

satnum = 44343
ndot = .00000777
bstar = 27989e-4

# manually convert the epoch string (wish this were an official skyfield function)
string_epoch = "20286.64549342"
two_digit_year = int(string_epoch[0:2])
fractional_day = float(string_epoch[2:])
if two_digit_year < 57:
    year = 2000 + two_digit_year
else:
    year = 1900 + two_digit_year

month, day, hour, minute, second = days2mdhms(two_digit_year, fractional_day)
julian_day = jday(year, month, day, hour, minute, second)

# orbital elements from the TLE
inclination = 24.0014 / DEGREES_PER_RADIANS
raan = 31.6920 / DEGREES_PER_RADIANS
eccentricity = .0010169
argument_of_perigee = 255.8372 / DEGREES_PER_RADIANS
mean_anomaly = 104.1002 / DEGREES_PER_RADIANS
mean_motion = 15.07436211

# convert things into the format expected by Satrec
epoch = julian_day - 2433281.5
mean_motion_radians_per_minute = mean_motion * (2 * np.pi) / (24 * 60)

satrec = Satrec()
satrec.sgp4init(
    WGS72,  # gravity model
    'i',  # 'a' = old AFSPC mode, 'i' = improved mode
    satnum,  # satnum: Satellite number
    epoch,  # epoch: days since 1949 December 31 00:00 UT
    bstar,  # bstar: drag coefficient (/earth radii)
    ndot,  # ndot: ballistic coefficient (revs/day)
    0.0,  # nddot: second derivative of mean motion (revs/day^3)
    eccentricity,  # ecco: eccentricity
    argument_of_perigee,  # argpo: argument of perigee (radians)
    inclination,  # inclo: inclination (radians)
    mean_anomaly,  # mo: mean anomaly (radians)
    mean_motion_radians_per_minute,  # no_kozai: mean motion (radians/minute)
    raan,  # nodeo: right ascension of ascending node (radians)
)

satrec_satellite = EarthSatellite.from_satrec(satrec, TIMESCALE)

tle_satellite = EarthSatellite(tle1, tle2, name, TIMESCALE)

print(tle_satellite.epoch)  # <Time tt=2459135.1462941607>
print(satrec_satellite.epoch)  # <Time tt=2459134.6462941607> - a half-day early!

[chylu]

As an aside, when I calculate epoch manually, I'm consistently getting results about 30 seconds removed from the epoch of EarthSatellites generated from TLEs. I don't think it's of grave importance, but I'm not sure why this happens.

[brandon-rhodes]

First of all, thanks Brandon for making such a thorough, powerful astronomy library. It's simplified my job tremendously.

Thank you for the kind words! That's very encouraging, and I'm happy to hear that Skyfield has made your job easier (rather than more complicated!).

I looked, and it does look like Skyfield was off by a half-day in its new from_satrec() method when building the .epoch parameter for users. If you would like to try out the fix without waiting for the next version of Skyfield to be released, you can try this:

pip install https://github.com/skyfielders/python-skyfield/archive/master.zip

Let me know if it works better, and thanks for letting me know about the discrepancy! Your example code was excellent; I was able to paste it and run it immediately.