shashwatak / satellite-js

Modular set of functions for SGP4 and SDP4 propagation of TLEs.
MIT License
902 stars 145 forks source link

Possible issue with accuracy #110

Closed letmaik closed 1 year ago

letmaik commented 1 year ago

I noticed that this library produces slightly different results than Python's skyfield library. Is this expected? Both libraries seem to have test vectors.

satellite-js: XYZ: [1366.441281990948, -5843.920852076301, 3180.698197781816 ] Longitude: : 58.053917147708106 Latitude: 28.072343404662533 Elevation: 418.87478629257475

skyfield: XYZ: [ 1343.67765613 -5850.69932306 3177.91747594] Longitude: 58.053774471070135 Latitude: 28.072163998036647 Elevation: 418.8672200797935

Minimal repro code:

<head>
<script src="https://unpkg.com/satellite.js@4.1.4/dist/satellite.min.js"></script>
</head>

<body>

<script>

function computeLocation(satrec, time) {
    const gmst = satellite.gstime(time);
    const eci = satellite.propagate(satrec, time);
    const gdPos = satellite.eciToGeodetic(eci.position, gmst);
    return {
        eci: eci.position,
        lat: satellite.radiansToDegrees(gdPos.latitude),
        lon: satellite.radiansToDegrees(gdPos.longitude),
        alt: gdPos.height
    }
}

const tle = [
    '1 25544U 98067A   22256.53615712  .00008923  00000+0  16307-3 0  9998',
    '2 25544  51.6425 258.8876 0002207 229.9964 274.2805 15.50186801358886']
const satData = satellite.twoline2satrec(tle[0], tle[1])
const pos = computeLocation(satData, new Date("2022-09-13T15:30:00Z"))
document.body.innerText = JSON.stringify(pos, null, 2)

</script>
</body>

And as reference, skyfield:

import datetime
import skyfield.sgp4lib as sgp4lib
import skyfield.api as skyfield

sat = sgp4lib.EarthSatellite(
    '1 25544U 98067A   22256.53615712  .00008923  00000+0  16307-3 0  9998',
    '2 25544  51.6425 258.8876 0002207 229.9964 274.2805 15.50186801358886')
t = sat.ts.from_datetime(datetime.datetime(2022,9,13,15,30,0, tzinfo=datetime.timezone.utc))
sat_pos = sat.at(t)
print(sat_pos.xyz.km)
sat_geo = skyfield.wgs84.geographic_position_of(sat_pos)
print(f"Epoch: {sat.epoch.utc_iso()}")
print(f"Longitude: {sat_geo.longitude.degrees}")
print(f"Latitude: {sat_geo.latitude.degrees}")
print(f"Elevation: {sat_geo.elevation.km}")
thkruz commented 1 year ago

Take a look at https://github.com/shashwatak/satellite-js/pull/108, but be sure to use WGS 72 in skyfield or you will surely get different answers. Without reviewing their code, it could be the use of 1st/2nd derivatives of mean motion causing differences or any of the bugs I address in that PR.

letmaik commented 1 year ago

I tried with your branch and get this: "x": 1366.422830298306, "y": -5843.92852074945, "z": 3180.6759955027815

"lat": 28.0721645587034, "lon": 58.053728950245045, "alt": 418.86722017424836

Interestingly, the XYZ values are still quite different from skyfield but the geographic coordinates now are nearly identical.

letmaik commented 1 year ago

FWIW, the XYZ values of skyfield are GCRS coordinates. Does satellite-js use something slightly different? Maybe I shouldn't look too closely at those and be happy with matching geographic coordinates.

thkruz commented 1 year ago

@letmaik I would recommend looking at https://space-track.org for the official propagator used by USSF. I have tested 20,000 or so TLEs against mine and got answers within 1mm of USSF for the ECI coordinates. I would treat their tool as the truth and then work backwards to determine which library has the issue.

If you can't figure it out I will take a look when I am done traveling in a week or so.

Depending on your application it may not matter much. If you are just trying to visualize stuff for fun it is a trivial amount of error between those three answers.

letmaik commented 1 year ago

@thkruz OK let's close this. I'll consider opening an issue on skyfield's repo, but for my use cases it's good enough. The issue only becomes apparent when you combine coordinates computed by both libraries in a single visualization, but I guess that's not a common use case.