Trying to match EGM2008 with STK

Orekit usage
1

Hello,

I am trying to use a NumericPropagator with the EGM2008 gravity model for a project that I am working on. Unfortunately, when I check the results against STK, they are way off after 15 days. When I visualize them in CZML, what I see is pretty good agreement for a day or so, and then the STK propagator appears to pull ahead of the Orekit propagator. The error by the end of day 15 is entirely in the in-track direction, and the orbital planes appear to remain co-planar. The STK propagator is roughly 2.5 minutes (in time) ahead at this point.

To initialize the EGM2008 model, downloaded the EGM2008 file (record 104) from the ICGEM website here:

I placed the EGM2008.gfc file into my orekit-data/Potential directory and renamed the eigen-6s.gfc file to eigen-6s.old.

My code to create the gravity model looks like this:

    public static ForceModel createGravityModel(int degree, int order, Frame centralBodyFrame) {

        NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(
                degree,
                order);
        ForceModel holmesFeatherstone = new HolmesFeatherstoneAttractionModel(
                centralBodyFrame,
                provider);

        return holmesFeatherstone;
    }

Note that I am using degree 70, order 70 for testing.

I have tried various centralBodyFrames [IERS_2010, IERS_2003, IERS_1996] with ‘simpleEOP’ set to [true,false]. Just for fun, I tried GTOD with applyEOPCorr set to false as well.

My propagator is initialized with this code

public static NumericalPropagator createOrekitPropagator(
            List<ForceModel> forceModelList,
            Orbit orbit) {
        final double minStep = 0.001;
        final double maxStep = 1000;
        
        // error control parameters (absolute and relative)
        final double positionError = 10.0;
        ODEIntegrator integrator = new DormandPrince853IntegratorBuilder(
                minStep,
                maxStep,
                positionError).buildIntegrator(
                orbit,
                OrbitType.CARTESIAN);

        // set up space dynamics propagator
        NumericalPropagator propagator = new NumericalPropagator(integrator);

        for (ForceModel forceModel : forceModelList) {
            propagator.addForceModel(forceModel);
        }

        return propagator;
    }

I’m relatively confident that I have my initial orbit setup right as STK and the Orekit propagator appear in the same location when viewing the cesium.

With that being said, I’m at a bit of a loss here. Does anyone have any suggestions for next steps? If I can get things to line up, Orekit will do pretty much exactly what I need it to do.

FWIW, this is a Sun Synch / ~550km orbit.

Thanks!

Screenshot of Cesium Viewer showing the difference between Orekit and STK for one centralBodyFrame (though all show the same behavior)

Screen Shot 2021-06-29 at 12.27.02 AM
Screen Shot 2021-06-29 at 12.27.02 AM1319×820 202 KB

2

Hi @fbl100, welcome onboard

I would guess two things are wrong with the adaptive step size settings.
You could first try to use a lower value for positionError; 10m is good for quick analysis and short term, but insufficient for either accurate propagation or long term. I would suggest something much lower, say 0.1m or 0.01m.
Then, it appears you set up the underlying mathematical integrator tolerances notifying it that it should spread the scalar position (i.e. the 10m distance above) to a state vector that will at the end represent Cartesian elements, but you don’t configure the Numerical propagator orbit type and you end up with the default OrbitType.EQUINOCTIAL. The state vector really used for integration is therefore not consistent with how the tolerance vector has been set in the mathematical integrator. The mathematical integrator has no way to understand what the state vector represents, for it is is only an array of double, this is why you have to pass it the correct information when it computes the tolerance vector from the scalar position tolerance. As you specified OrbitType.CARTESIAN, basically the three first elements of the tolerance vectors are set to 10m (and the three following ones probably set to something of the order of magnitude of 0.01m/s) as they are assumed to be x, y, z, vx, vy, vz. But as the NumericalPropagator will really use a EquinoctialOrbit, the second and third components will really be an eccentricity vector, not a Y and Z Cartesian components, and a tolerance of 10.0 on eccentricity is… huge.

If you want to be sure you remain consistent between the integrator settings and the numerical propagator setting, you should call setOrbitType on the NumericalPropagator after it has been built and before it is started, making sure you use the same type as you used for integrator configuration.

3

@luc Thank you for the quick reply.

Changing the tolerance and the orbitType on the propagator definitely helped. After further discussion with my teammate (who is generating the STK data), we realized that we were using different integrators. He was using RK4. When I changed to RK4 with the same step size as him, we’re now 1500m apart after 15 days. That’s not fantastic, but I think it’ll work for now.

Can you think of any other settings that we may need to check? We did double check that we’re using the same value for Mu.

4

You should take care of the body equatorial radius in the gravity field. In fact, it has not really any geometrical meaning, it is rather a scaling coefficient, which is often confusing. If for example you use a model with ae= 6378136.3 m (which happens to be Constants.EGM96_EARTH_EQUATORIAL_RADIUS) but you initialized the value using Constants.WGS84_EARTH_EQUATORIAL_RADIUS which is 63781367m, then you get inconsistencies. This does not mean that you select one of these two values and use it everywhere, it means that you should use the equatorial radius that is consistent with the gravity model for computing gravity acceleration, even if you use another equatorial radius for models corresponding to the geometry of the surface of the Earth : these are two different concepts and they can use different values, and it is confusing. So check this radius too making sure you compare data with the same meaning (either acceleration or geometry) in both programs.

You may look at your other force models if any.
If for example you have drag there are many things that can be different, but I guess you have no drag, otherwise you would have a much larger differences after 15 days.
If you have 3rd body attraction, check the ephemerides of the body (Moon and Sun).

5

Thanks again for the explanation. It looks like my gravity model is correctly reading the 6378136.3 for EGM2008 (it’s also listed in the gfc file) and I’m not initializing it anywhere else. In the debugger, I see ae set to 6378136.3 on the ConstantSphericalHarmonics provider. I feel like 1500m is a little much over 15 days, but getting things to line up with STK is always a bit of a challenge.

When it comes to atmosphere, I was unable to get the NRLMSISE00 model to compare within reason. Orekit wants Ap, but STK wants Kp, but I don’t think that was the difference. Similarly, the HarrisPriester model diverged rather quickly. I suspect that it’s a setup problem, but for the time being, I was able to wrap my implementation of Jacchia Roberts that has been validated against STK and provide that Atmosphere instead which produced reasonable results.

Thanks for all the help.

6

Could also check you’re using the same EOP and leap second file. That makes a small difference as to the direction of gravity.

7

In my experience, using DormandPrince853Integrator you need a positionError of 1e-4 m to achieve sub-meter level precision (<1m) when propagating a high-drag orbit for one day.
But I recommend testing different positionError to check the precision.