When using Orekit for orbit prediction with an atmospheric drag model, the predicted orbit 10 days out shows a 6-second difference compared to STK

Sohnny · December 13, 2024, 8:36am

With atmospheric drag models disabled in both Orekit and STK, the prediction results are highly consistent. To clarify, for a 10-day orbit prediction, STK differs from the actual orbit by 1 second, while Orekit differs by 5 seconds, which translates to approximately 40 kilometers.The following are the configurations for Orekit and STK, respectively. I can ensure the input parameters are consistent:

The only difference I notice is that STK has an additional parameter for the area-to-mass ratio.

yzokras · December 13, 2024, 8:53am

Hi, can you provide a minimum working example code? In particular the propagator and the propagator parameters, as well as the mass.

Sohnny · December 13, 2024, 9:11am

    public static NumericalPropagator createNumericalPropagator(Orbit orbit, SpacecraftState spacecraftState) { // mass 85kg

            int degree = 32;
            int order = 32;
            
            double[][] tolerances = NumericalPropagator.tolerances(10.0, orbit, OrbitType.KEPLERIAN);

            NumericalPropagator propagator = new NumericalPropagator(new DormandPrince853Integrator(0.001, 60.0, tolerances[0], tolerances[1]));
            propagator.setInitialState(spacecraftState);
            propagator.setOrbitType(OrbitType.CARTESIAN);
            propagator.setPositionAngleType(PositionAngleType.MEAN);

            Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);

            NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(degree, order);

            ForceModel holmesFeatherstone = new HolmesFeatherstoneAttractionModel(itrf, provider);
            propagator.addForceModel(holmesFeatherstone);

            propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
            propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));


            OneAxisEllipsoid earth = ReferenceEllipsoid.getWgs84(itrf);

            CssiSpaceWeatherData inputParameters = new CssiSpaceWeatherData("SpaceWeather-All-v1.2.txt");

            Atmosphere atmosphere = new NRLMSISE00(inputParameters, CelestialBodyFactory.getSun(), earth);

            ForceModel dragModel = new DragForce(atmosphere, new IsotropicDrag(0.41396452, 4.115137));
            propagator.addForceModel(dragModel);


        return propagator;
    }

yzokras · December 13, 2024, 10:07am

Thanks. Are the force models in your Orekit propagator configured exactly like in STK? Did you set the spacecraft mass in the initial SpacecraftState? Otherwise by default the mass in Orekit is 1000kg IIRC.

Sohnny · December 20, 2024, 2:08am

As shown in my screenshots above, the atmospheric drag configuration parameters in Orekit are identical to those in STK, and both use the same SpaceWeather-All-v1.2.txt file. The mass is also the same at 85 kg (as noted in the code comments above). I am using version 12.0.1.

luc · December 20, 2024, 7:44am

I see two problems:
- tolerance NumericalPropagator.tolerances(10.0, orbit, OrbitType.KEPLERIAN); is huge, we recommend using 0.001m rather than 10m
- tolerance is computed using a state vector type set to OrbitType.KEPLERIAN but propagator is configured with OrbitType.CARTESIAN

Sohnny · December 23, 2024, 2:26am

I changed the tolerance settings as follows, but the results are essentially the same as before:

            double[][] tolerances = NumericalPropagator.tolerances(0.001, orbit, OrbitType.CARTESIAN);

            NumericalPropagator propagator = new NumericalPropagator(new DormandPrince853Integrator(0.001, 60.0, tolerances[0], tolerances[1]));
            propagator.setInitialState(spacecraftState);
            propagator.setOrbitType(OrbitType.CARTESIAN);
            propagator.setPositionAngleType(PositionAngleType.MEAN);

One point I mentioned earlier that I think is important is that the 5-second difference only appears when the atmospheric drag model is included. Without the drag model, the results match STK.

luc · December 23, 2024, 8:51am

Drag is one of the forces that are very sensitive to both initial conditions, modeling and integration. This is because density profile is exponential, hence whenever you have even the slightest difference in altitude between two propagation runs, it build up rather quickly.
In your case, you are propagating over a 10 days period, which is quite long when drag is involved and you still want accurate results. What its the altitude of your spacecraft?

Another point: are your Earth Orientation Parameters the same for Orekit and STK? If for example you don’t have EOP for this date in your orekit-data, then Orekit will just use the deterministic precession/nutation models, so you would end up with differences in position with respect to Earth (and hence to atmosphere) that could be around 15m due to polhody and up to a few hundred meters due to dUT1, so different density, then different acceleration then different propagation.

Sohnny · December 23, 2024, 9:13am

The spacecraft’s altitude is approximately 500 km, and it’s in a Sun-synchronous orbit. Your point about whether the Earth orientation parameters are the same between Orekit and STK is a good one – I’m not sure about that yet, and I’ll investigate it. Thank you for taking the time to analyze this issue.

Sohnny · December 24, 2024, 7:01am

I’m using Earth Orientation Parameters downloaded from here: Orekit / Orekit Data · GitLab. I haven’t added any extra data; I’m using the orekit-data directly after unpacking it. I noticed STK uses
EOP-v1.1.txt (213.9 KB)
. I don’t know how to ensure Orekit and STK use the same orientation file, as their formats seem different.

I tried deleting the “Earth-Orientation-Parameters” folder and found the difference with and without it is on the order of milliseconds. Given my current 5-second error, is this small enough to rule out the orientation file as the source of the error?

luc · December 24, 2024, 10:08am

The STK file you link to is also quite old, it contains only data up to september 2019.

The value you use depend on when you downloaded it, but I guess you have something far more recent than 2019 here. The orekit-data repository we maintain is just an example and updated only from time to time. It is the responsibility of users to refresh data regularly (their is an update.sh script that does this automatically, and people can just download the EOP files by themselves if they want).

I would not rule out EOP yet, as I fear STK EOP are too old and therefore inconsistent with orekit-data. Try updating the STK file and see if it changes something on their side.

Sohnny · December 25, 2024, 2:55am

I’m using actual ephemeris data from April 2024 for the comparison. The orekit-data was updated in December.
I updated the
EOP-v1.1.txt (239.9 KB)
file in STK to the latest version, and as I suspected, the difference is now only a few tens of milliseconds. It seems safe to rule out this factor.

baubin · January 6, 2025, 3:05am

@Sohnny I just thought of something. The NRLMSISE00 model has a bunch of switches, and turning them on/off can vary the drag result by km. I run NRLMSISE00 with switch 9 set to -1 and switch 1 set to 1:

NRLMSISE_CSSI_MODEL = NRLMSISE00(
    CSSI_WEATHER_MODEL,
    CelestialBodyFactory.getSun(),
    EARTH_BODY
).withSwitch(9, -1).withSwitch(1, 1)

Which as far as I can tell is the most high-fidelity way to run it. If you have different switch settings on in STK vs. Orekit you’d get very different results.

Sohnny · January 8, 2025, 7:17am

@baubin Thank you for your guidance. Your method did reduce the error, but in my actual tests, it only decreased by 100ms. In the example I provided above, there’s still a 5-6 second discrepancy compared to STK.

a_gior · January 17, 2025, 10:25am

Hello,
I am also trying to compare results and I was wondering the following:

which model is recommended between cssi and MSAFE? I have noticed that the latter, when set to AVERAGE, is similar (but not equal) to cssi. While when set to STRONG, provides much higher results than cssi.
is it possible to create a model with user-input constant values of Kp and F10.7?
is there a way to evaluate the Ap/Kp/F10.7 while using MSAFE? While with cssi I am able to retrieve this information, doing the same with MSAFE generates the following error.

strenghtLevel = MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE
supportedNames = MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES
inputParameters_Drag = MarshallSolarActivityFutureEstimation(supportedNames, strenghtLevel)
inputParameters_Drag.getDailyFlux(initialDate)
---------------------------------------------------------------------------------------------------------

inputParameters_Drag.getDailyFlux(initialDate)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
orekit.JavaError: <super: <class 'JavaError'>, <JavaError object>>
    Java stacktrace:
java.lang.NullPointerException
	at org.orekit.time.AbsoluteDate.durationFrom(AbsoluteDate.java:987)
	at org.orekit.models.earth.atmosphere.data.MarshallSolarActivityFutureEstimation.bracketDate(MarshallSolarActivityFutureEstimation.java:241)
	at org.orekit.models.earth.atmosphere.data.MarshallSolarActivityFutureEstimation.getMeanFlux(MarshallSolarActivityFutureEstimation.java:300)
	at org.orekit.models.earth.atmosphere.data.MarshallSolarActivityFutureEstimation.getDailyFlux(MarshallSolarActivityFutureEstimation.java:381)

baubin · January 17, 2025, 9:12pm

@Sohnny

I apologize if this is a dumb question, but did you actually set the mass in your NumericalPropagator. I only ask because I don’t see that value set in your code snippet. Your drag value 0.4139 m^2 is correct based on what you’ve told me, but you actually need to add

propagator.setMass(85) //mass in kg

To your code for the mass value in the propagator to be set correctly.