Doubts about DSST using (comparing with numerical integration and GMAT)

Orekit usage
1

Hello,
I have some problems with DSSTPropagator using, so I’ll be very grateful for your help.
My purpose is to propagate LEO satellite orbit for long periods, so I found recommendations about using DSSTPropagator for such calculations.
I used examples and tests and wrote a code, but the results of calculations are quite different from the same calculations, made with numerical propagator or in GMAT.
The differences are about x00 metres in semi-major axis and about 0.0x degrees in inclination for 10 - 30 days. I understand, that it’s incorrect to compare position-velocity components, so I compare orbit elements. I tried different combinations of “mean/osculating” tunings, but in all cases the difference seems to be too big.
At the beginning I used model with many forces (Earth gravity, Moon and Sun gravity, atmospheric drag), but I noticed, that the problem appears with the Earth gravity (the results for DSST and numerical integration are the same for spheric Earth gravity field, but there is a difference for 2 x 0 and further).
I examined all the tunings about constants, coordinate systems, etc., but I have no idea, how to explain such a difference. (I tried both position-velocity or keplerian orbit for initial SpacecraftState, nothing changed.)
Here is my code.

            File orekitData = new File("C:\\tempForOrekit");
            DataContext.getDefault().getDataProvidersManager().addProvider(new DirectoryCrawler(orekitData));

            AbsoluteDate initialDate = new AbsoluteDate(
                    2020, 3, 12, 4, 37, 0, TimeScalesFactory.getUTC());
            PVCoordinates pvCoordinates = new PVCoordinates(
                    new Vector3D(6823.698321e+3, -952.003960e+3, -13.187505e+3),
                    new Vector3D(-0.131194e+3, -0.976418e+3, 7.543914e+3));
            Frame inertialFrame = FramesFactory.getEME2000();
            KeplerianOrbit initialOrbit = new KeplerianOrbit(
                    6893.081283931438e+3,
                    0.001331122212474694,
                    Math.toRadians(97.44029236157961),
                    Math.toRadians(69.0070994141481),
                    Math.toRadians(352.0433593955477),
                    Math.toRadians(290.8823016332061),
                    PositionAngle.TRUE,
                    inertialFrame, initialDate, Constants.EGM96_EARTH_MU);

            final SpacecraftState initialState = new SpacecraftState(initialOrbit);

            // Initialize the DSST propagator
            final double minStep = initialState.getKeplerianPeriod();
            final double maxStep = 100. * minStep;
            final double[][] tol = DSSTPropagator.tolerances(1.0, initialState.getOrbit());
            final AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
            final DSSTPropagator dsstProp = new DSSTPropagator(integrator, PropagationType.OSCULATING);
            dsstProp.setInitialState(initialState, PropagationType.OSCULATING);

            // Central Body geopotential 2x0 (J2 only)
            final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(2, 0);
            final OneAxisEllipsoid earth = new OneAxisEllipsoid(provider.getAe(),
                    Constants.WGS84_EARTH_FLATTENING,
                    FramesFactory.getITRF(IERSConventions.IERS_2010, true));
            // Add the force models to the propagator
            dsstProp.addForceModel(new DSSTZonal(provider));
            dsstProp.addForceModel(new DSSTTesseral(earth.getBodyFrame(), Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider));

            // 10 days propagation
            final SpacecraftState finalState = dsstProp.propagate(initialDate.shiftedBy(10. * Constants.JULIAN_DAY));
            PVCoordinates pv = finalState.getPVCoordinates(inertialFrame);
            Orbit orbit = new KeplerianOrbit(pv, inertialFrame, finalState.getDate(), Constants.EGM96_EARTH_MU);
            double aM = orbit.getA();
            double e = orbit.getE();
            double iDegrees = Math.toDegrees(orbit.getI());

//            reference data from GMAT:
//            Cartesian State                       Keplerian State
//            ---------------------------           --------------------------------
//            X  =   2617.7180296979 km             SMA  =   6877.0840271856 km
//            Y  =   916.41629724835 km             ECC  =   0.0022472704215
//            Z  =  -6308.9965773698 km             INC  =   97.468200487389 deg
//            VX =   7.0209016138680 km/sec         RAAN =   1.9452039595408 deg
//            VY =  -0.1417284428202 km/sec         AOP  =   93.814590070776 deg
//            VZ =   2.8985288247191 km/sec         TA   =   198.77517269129 deg
//            MA   =   198.85818878177 deg
//            EA   =   198.81665869991 deg

            assertEquals(6877.0840271856e+3, aM, 1); // got 6877039.514159619
            assertEquals(0.0022472704215, e, 1e-6); // got 0.002250298057808621
            assertEquals(97.468200487389, iDegrees, 1e-6); // got 97.44902204199919

Thank you in advance for your answers.

2

Hi @IrinaPonomareva

In your code, you seem to retrieve just the position and velocity from the initial state and rebuild an orbit, setting the central attraction coefficient to a constant. It would probably be better to just get the orbit using finalState.getOrbit(). The orbit type used by DSST is equinoctial, but all orbits implement getA(), getE() and getI() so you don’t have to convert it.

I have doubts about the central attraction coefficient. You use Constants.EGM96_EARTH_MU, you should perhaps retrieve the one from the gravity field, using provider.getMu(). I don’t know if it matches the one from the gravity field but it would be safer.

I don’t see any obvious error any error in the code. You may try to reduce tolerance when you call DSSTPropagator.tolerances to something like 0.01 but I doubt it will change a lot the results.

Do you see any pattern in the error (some short period terms for example) if you plot the difference as a time series instead of just one point?
Are all three propagators (DSST, numerical and GMAT) different or do at least two of them agree with each other?

3

Hi, @luc ,
thank you for your quick answer.
I followed your recommendations and modified my code, but there is no positive effect.
Answering your questions:

  1. Numerical integration and GMAT give nearly the same results (for 10 days propagation):
----------
Difference GMAT <-> DSST in aM is 44.51298189163208
Difference  <-> DSST in e is 3.027633090022311E-6
Difference  <-> DSST in iDegrees is 0.019178445388746468
----------
Difference GMAT <-> num in aM is 0.7631760407239199
Difference  <-> num in e is 2.688015611877932E-6
Difference  <-> num in iDegrees is 1.1065084581218798E-4
----------
  1. I compared numerical integration and DSST for time series, there is obvious systematic effect in inclination difference:

изображение
изображение2745×1795 236 KB

My modified code is below:

            File orekitData = new File("C:\\tempForOrekit");
            DataContext.getDefault().getDataProvidersManager().addProvider(new DirectoryCrawler(orekitData));

            // Central Body geopotential 2x0 (J2 only)
            final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(2, 0);
            final OneAxisEllipsoid earth = new OneAxisEllipsoid(provider.getAe(),
                    Constants.WGS84_EARTH_FLATTENING,
                    FramesFactory.getITRF(IERSConventions.IERS_2010, true));

            AbsoluteDate initialDate = new AbsoluteDate(
                    2020, 3, 12, 4, 37, 0, TimeScalesFactory.getUTC());
            PVCoordinates pvCoordinates = new PVCoordinates(
                    new Vector3D(6823.698321e+3, -952.003960e+3, -13.187505e+3),
                    new Vector3D(-0.131194e+3, -0.976418e+3, 7.543914e+3));
            Frame inertialFrame = FramesFactory.getEME2000();
            KeplerianOrbit initialOrbit = new KeplerianOrbit(
                    6893.081283931438e+3,
                    0.001331122212474694,
                    Math.toRadians(97.44029236157961),
                    Math.toRadians(69.0070994141481),
                    Math.toRadians(352.0433593955477),
                    Math.toRadians(290.8823016332061),
                    PositionAngle.TRUE,
                    inertialFrame, initialDate, provider.getMu());

            final SpacecraftState initialState = new SpacecraftState(initialOrbit);

            // Initialize the DSST propagator
            final double minStepDSST = initialState.getKeplerianPeriod();
            final double maxStepDSST = 100. * minStepDSST;
            final double[][] tol = DSSTPropagator.tolerances(1.0, initialState.getOrbit());
            final AdaptiveStepsizeIntegrator integratorDSST = new DormandPrince853Integrator(minStepDSST, maxStepDSST, tol[0], tol[1]);
            final DSSTPropagator dsstProp = new DSSTPropagator(integratorDSST, PropagationType.OSCULATING);
            dsstProp.setInitialState(initialState, PropagationType.OSCULATING);
            // Add the force models to the DSST propagator
            dsstProp.addForceModel(new DSSTZonal(provider));
            dsstProp.addForceModel(new DSSTTesseral(earth.getBodyFrame(), Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider));

            // Initialize numerical propagator
            // By example from https://www.orekit.org/site-orekit-tutorials-11.1/tutorials/propagation.html
            // Adaptive step integrator
            // with a minimum step of 0.001 and a maximum step of 1000
            double minStepNum = 0.001;
            double maxStepNum = 1000.0;
            double positionTolerance = 0.01;
            OrbitType propagationType = OrbitType.CARTESIAN;
            double[][] tolerances =
                    NumericalPropagator.tolerances(positionTolerance, initialOrbit, propagationType);
            AdaptiveStepsizeIntegrator integratorNum =
                    new DormandPrince853Integrator(minStepNum, maxStepNum, tolerances[0], tolerances[1]);
            NumericalPropagator numProp = new NumericalPropagator(integratorNum);
            numProp.setOrbitType(propagationType);
            NormalizedSphericalHarmonicsProvider provider2 =
                    GravityFieldFactory.getNormalizedProvider(2, 0);
            ForceModel holmesFeatherstone =
                    new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010,
                            true),
                            provider2);
            numProp.addForceModel(holmesFeatherstone);
            numProp.setInitialState(initialState);

            //30 days propagation with 1 hour step - comparing of DDST and numerical integration
            AbsoluteDate intermediateDate = initialDate;
            AbsoluteDate finalDate = initialDate.shiftedBy(30. * Constants.JULIAN_DAY);
            while (intermediateDate.isBefore(finalDate)) {
                //intermediate result for DSST
                SpacecraftState finalStateDSST = dsstProp.propagate(intermediateDate);
                Orbit orbitDSST = finalStateDSST.getOrbit();
                double aMDSST = orbitDSST.getA();
                double eDSST = orbitDSST.getE();
                double iDegreesDSST = Math.toDegrees(orbitDSST.getI());
                //intermediate result for numerical integration
                SpacecraftState finalStateNum = numProp.propagate(intermediateDate);
                Orbit orbitNum = finalStateNum.getOrbit();
                double aMNum = orbitNum.getA();
                double eNum = orbitNum.getE();
                double iDegreesNum = Math.toDegrees(orbitNum.getI());
                System.out.println(finalStateDSST.getDate().toString() + "; "
                        + (aMDSST - aMNum) + "; "
                        + (eDSST - eNum) + "; "
                        + (iDegreesDSST - iDegreesNum) + ";");
                intermediateDate = intermediateDate.shiftedBy(1.0 * Constants.JULIAN_DAY / 24.0);
            }

            // 10 days propagation - comparing with GMAT for final date
            final SpacecraftState finalStateDSST = dsstProp.propagate(initialDate.shiftedBy(10. * Constants.JULIAN_DAY));
            Orbit orbitDSST = finalStateDSST.getOrbit();
            double aMDSST = orbitDSST.getA();
            double eDSST = orbitDSST.getE();
            double iDegreesDSST = Math.toDegrees(orbitDSST.getI());

//            Cartesian State                       Keplerian State
//            ---------------------------           --------------------------------
//            X  =   2617.7180296979 km             SMA  =   6877.0840271856 km
//            Y  =   916.41629724835 km             ECC  =   0.0022472704215
//            Z  =  -6308.9965773698 km             INC  =   97.468200487389 deg
//            VX =   7.0209016138680 km/sec         RAAN =   1.9452039595408 deg
//            VY =  -0.1417284428202 km/sec         AOP  =   93.814590070776 deg
//            VZ =   2.8985288247191 km/sec         TA   =   198.77517269129 deg
//            MA   =   198.85818878177 deg
//            EA   =   198.81665869991 deg

            System.out.println("----------");
            System.out.println("Difference GMAT <-> DSST in aM is " + Math.abs(aMDSST - 6877.0840271856e+3));
            System.out.println("Difference  <-> DSST in e is " + Math.abs(eDSST - 0.0022472704215));
            System.out.println("Difference  <-> DSST in iDegrees is " + Math.abs(iDegreesDSST - 97.468200487389));

            final SpacecraftState finalStateNum = numProp.propagate(initialDate.shiftedBy(10. * Constants.JULIAN_DAY));
            Orbit orbitNum = finalStateNum.getOrbit();
            double aMNum = orbitNum.getA();
            double eNum = orbitNum.getE();
            double iDegreesNum = Math.toDegrees(orbitNum.getI());

            System.out.println("----------");
            System.out.println("Difference GMAT <-> num in aM is " + Math.abs(aMNum - 6877.0840271856e+3));
            System.out.println("Difference  <-> num in e is " + Math.abs(eNum - 0.0022472704215));
            System.out.println("Difference  <-> num in iDegrees is " + Math.abs(iDegreesNum - 97.468200487389));
            System.out.println("----------");

Thank you in advance for your answers.

4

I don’t understand what happens.
I guess you could remove the DSSTTesseral force model as this test case only uses J₂.
I first suspected the gravity fields models were different between GMAT and Orekit, but as they are consistent when numerical propagation is used, this rules out this hypothesis.
There are some J₂-related effects we don’t have implemented yet in Orekit: the J₂ squared terms and the J₂ secular m-daily coupling terms. Maybe this is what you observe.

See for example the PDF file attached at the end of the following message from Paul Cefola in the old mailing lists : thread 2016-04-15, it references these missing terms. These terms are still missing now.

5

Hi Irina and Luc,

I confirm that the differences between GMAT and Orekit-DSST for your example are caused by the J2 squared terms that are not considered in Orekit-DSST. According to [1](Figures 1 and 2) for an inclined near circular orbit, like the one in your example, the error is about 50 meters in position.

Adding J2 squared terms in Orekit-DSST is something we want to do since years. Unfortunately, we have not yet found the time to do this development.
In [1], Folcik and Cefola propose an implementation based on a Gaussian quadrature. It is a very accurate implementation, but Gaussian quadrature are usually time consuming. Indeed, in DSST theory, atmospheric drag and solar radiation pressure effects are also based on a Gaussian quadrature and adding these perturbations increase a lot the computation time.
Last year, San Juan et al. [2] propose a closed-form implementation of the J2 squared effect in DSST allowing to not use the Gaussian quadrature to evaluate the effect. This has the advantage of being very fast.

In my opinion, the best implementation to perform in Orekit is the one in Reference [2]. Could you open an issue in order to add these terms? Could you put the link of this forum’s discussion in the issue description?

Thank you and best regards,
Bryan

[1] : Z. Folcik and P. J. Cefola, A general solution to the second order J2 contribution in a mean element semianalytical satellite theory,” Proceedings Advanced Maui Optical and Space Surveillance Technologies Conference, AMOS 2012, Maui, HI, USA, Maui Economic Development Board, September 2012 (385.1 KB)

[2] : San-Juan J. F., Lopez R., Perez I., Folcik Z., and Cefola P. J., Including the close-form J 2 2 in DSST, AAS Paper 21-233, 31st AAS/AIAA Space Flight Dynamics Meeting, Charlotte, NC, January 31 – February 4, 2021. (431.7 KB)

6

Hi @bcazabonne , thank you for your answer. Issue #931 (Differences between DSST and GMAT / numerical integrations (J2 squared terms not implemented) (#931) · Issues · Orekit / Orekit · GitLab) is created.

7

Thank you!