SGP4: STK vs Orekit

Hello,

I just finished running a quick comparison between STK’s and Orekit’s SGP4. I have attached the documents showing everything. I basically used a TLE from the ISS and propagated for one day. The CSV attached shows a comparison between the cartesian position and also for the keplerian elements. They both give practically the same results, my only concern is that, even though in terms of percentage we are talking about less than 1% difference, in absolute terms, for the X, Y, and Z positions the absolute difference is about 20km on average. Is this still acceptable and what are the reasons that could lead to this?

Thank you!

Iliass

ps: all in ICRF

STK vs Orekit.xlsx (1.3 MB)
TLE.tce (140 Bytes)

Looking at the smooth sinusoidal errors at orbital period, and assuming you use the exact same TLE in both programs, I would think there is a discrepency in frames. This discrepency is really large.
Are you sure about the frames conversions you use in both Orekit and STK? Note that TLE are given in a very strange and ill-defined frame called TEME, which is neither EME200, nor GCRF, nor MOD, nor TOD. It is almost inertial like the other ones, but specific to TLE.
Another point is that converting to Keplerian elements just don"t work. TLE are mean elements and point-to-point conversion is not meaningful.

2 Likes

If you use the same TLE in both programs, no it is not acceptable. TLE have very crude accuracy, and we often see errors that count in hundreds of kilometers int TLE, so you should never expect the truth to be extracted from TLE. However, the model is widespread and all programs lying around use the same, so they should all give the same (wrong) position if given the same TLE. Orekit handling of TLE is well validated, STK handling of TLE is well validated, and the same is true for probably all other programs. The only reasons I see to explain such huge errors are either not using the same TLE or some frames conversions problems, or using Keplerian elements as an intermediate step which is a clear no-go.

Hello Luc and thank you for your answer. I have attached a test case file testCase.zip (21.2 KB).

I use the exact same TLE for both STK and Orekit. I have not made any previous conversions to Keplerian elements, I use directly the TLE with the TLE Propagator available in Orekit. With regards to the frames, I specify the satellite one to be GCRF, does not that mean the propagator will do the conversion itself when giving the results?
In STK I simply generated a report selecting ICRF (which I believe is equivalent to GCRF in this case) for the cartesian and keplerian coordinates.

Thank you.

Iliass

In TestCase.propagate, at lines 61 to 66, replace

    org.orekit.propagation.Propagator propagator = TLEPropagator.selectExtrapolator(
        new TLE(input.getTle1(), input.getTle2()),
        attitudeProvider,
        input.getMass(),
        satelliteFrame
    );

by

    org.orekit.propagation.Propagator propagator = TLEPropagator.selectExtrapolator(
        new TLE(input.getTle1(), input.getTle2()),
        attitudeProvider,
        input.getMass(),
       FramesFactory.getTEME());

Then in EphemerisHandler.getCartesianParameter, replace

    Vector3D orekitPosition = state.getPVCoordinates().getPosition();
    Coordinate position = new Coordinate(orekitPosition.getX(), orekitPosition.getY(), orekitPosition.getZ());
    Vector3D orekitVelocity = state.getPVCoordinates().getVelocity();
    Coordinate velocity = new Coordinate(orekitVelocity.getX(), orekitVelocity.getY(), orekitVelocity.getZ());
    Vector3D orekitAcceleration = state.getPVCoordinates().getAcceleration();
    Coordinate acceleration = new Coordinate(orekitAcceleration.getX(), orekitAcceleration.getY(), orekitAcceleration.getZ());

    return new CartesianParameter(position, velocity, acceleration);

by

   PVCoordinates pva = state.getPVCoordinates(satelliteFrame);
  return new CartesianParameter(pva.getPosition(),
                                pva.getVelocity(),
                                pva.getAcceleration());

This is because the TLE propagator must be in TEME frame, and then you should convert its output into another frame as desired.

There are also other uses of Cartesian coordinates in you code, at least in EphemerisHandler.getOther and EphemerisHandler.getSubSatellitePoint where you should also perform the frame transform as the state will be in TEME but you want it in satelliteFrame.

1 Like

Yes, that makes sense. I have applied the changes and it works. I am getting submeter differences on the cartesian position now. What about the keplerian elements though?

TLE correspond to a mean orbit, this implies that short periods have been removed from them. A consequence is that the Keplerian elements that you get from a single point may not be significant and may be inconsistent. The only generally accepted way to convert TLE to anything else is to perform fitting over a long range that covers at least one orbit, generally more.
TLE are a strange beast, apart from everything else in space flight dynamics, and compatible only with itself. It is so widely used because data are publicly available in this format, not because it has any technical advantage. The model is poor, the frame is ill-defined, the data has unknown accuracy (but generally very crude accuracy), but it is used, widely used, much too often in my opinion.

1 Like

Yes, I understand the TLE is a mean orbit and its different limitations, but here we are talking about the ephemeris generated after the propagation, if the conversion is done from the cartesian elements, which should be possible, as these are consistent the Keplerian elements should also be and if for the cartesian elements we are getting the same as STK so it should be for the keplerian, don’t you think? As in, I am not trying to generate a keplerian orbit from the TLE, but rather simply expressing the position as the six classic elements instead of the cartesian ones. I hope this makes sense.

Thank you as always Luc :smiley:

Yes, you could use this just as a change of variable and not as something with physical meaning. In this case, you could compare it to STK but have to check the central attraction coefficient, which is specific for the SGP models (µ = 3.986008e+14 m³/s² if the comment in TLEConstants.java is correct).
But don’t try to do any computation with that, including a small shiftedBy(dt) with a small dt, I would not trust the result.