SGP4 propagation and reference frame issue

[axelsauvage]

Hello,

I tested the Lox project to compute GEO satellite position from TLE using the SGP4 propagator.
I found incorrect results in my application when computing satellite position in ICRF reference frame.

Here is a basic program to compute such position using Lox:

use lox_orbits::propagators::sgp4::{Elements, Sgp4};
use lox_time::{time_scales::Tai, TimeBuilder};

use lox_orbits::propagators::Propagator;

fn main() {
    let tle = Elements::from_tle(
        Some("INTELSAT 36 (IS-36)  ".to_string()),
        "1 41747U 16053A   25026.69560333 -.00000008  00000+0  00000+0 0  9995".as_bytes(),
        "2 41747   0.0093  36.1594 0001240 298.1404 110.8362  1.00272270 30870".as_bytes(),
    )
    .unwrap();
    let sgp4 = Sgp4::new(tle).unwrap();

    let time = TimeBuilder::new(Tai).with_ymd(2025, 1, 27).with_hms(0, 0, 0.0).build().unwrap();
    let state  = sgp4.propagate(time).unwrap();

    println!("x = {} km", state.position()[0]);
    println!("y = {} km", state.position()[1]);
    println!("z = {} km", state.position()[2]);
}

Output:

x = -40755.39647027197 km
y = -10823.118977946837 km
z = 12.227111947919987 km

I use Python lib called Skyfield to check the position, here is a simple program for test:

from skyfield.api import load
from skyfield.iokit import parse_tle_file
from skyfield.api import wgs84
from skyfield.framelib import itrs, ICRS

ts = load.timescale()

tle = [
    b"INTELSAT 36 (IS-36)  ",
    b"1 41747U 16053A   25026.69560333 -.00000008  00000+0  00000+0 0  9995",
    b"2 41747   0.0093  36.1594 0001240 298.1404 110.8362  1.00272270 30870"
]
satellite = list(parse_tle_file(tle, ts))[0]
print(satellite)

t = ts.tai(2025, 1, 27, 0, 0, 0)
x, y, z = satellite.at(t).frame_xyz(ICRS).km

print("icrs")
print('  x = {:.3f} km'.format(x))
print('  y = {:.3f} km'.format(y))
print('  z = {:.3f} km'.format(z))

Output:

INTELSAT 36 (IS-36) catalog #41747 epoch 2025-01-26 16:41:40 UTC
icrs
  x = -40815.298 km
  y = -10594.440 km
  z = 112.144 km

After checking the internal SGP4 propagator, I found that it output position in TEME frame (related to SGP4 algorithm), and not exactly the ICRF frame, which could expain the difference in position in the Lox and Skyfield lib.

In the Lox lib, the output from SGP4 is stored in ICRF frame:

impl Propagator<Tai, Earth, Icrf> for Sgp4 {
    type Error = Sgp4Error;

    fn propagate(&self, time: Time<Tai>) -> Result<State<Tai, Earth, Icrf>, Self::Error> {
        let dt = (time - self.time).to_decimal_seconds() / SECONDS_PER_MINUTE;
        let prediction = self.constants.propagate(MinutesSinceEpoch(dt))?;
        let position = DVec3::from_array(prediction.position);
        let velocity = DVec3::from_array(prediction.velocity);
        Ok(State::new(time, position, velocity, Earth, Icrf))
    }
}

But the SGP4 prediction stores position in the TEME frame, which is different frame:

/// Predicted satellite position and velocity after SGP4 propagation
///
/// The position and velocity are given in the True Equator, Mean Equinox (TEME) of epoch reference frame.
#[derive(Debug, Clone)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct Prediction {
    /// The three position components (x, y, z) in km
    pub position: [f64; 3],

    /// The three velocity components (x, y, z) in km.s⁻¹
    pub velocity: [f64; 3],
}

Have you used SGP4 satellite propagation in this project yet ? Am I doing something incorrect ?
Any help on this point would be helpful, thank you.

Cargo.toml:

[package]
name = "lox-sgp4"
version = "0.1.0"
edition = "2021"

[dependencies]
lox-bodies = "0.1.0-alpha.5"
lox-orbits = "0.1.0-alpha.8"
lox-time = "0.1.0-alpha.3"

[helgee]

Thanks for raising this issue and looking into this. Your analysis seems to be correct and we are missing the reference frame transformation.