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

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");

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,
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

// 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
``````

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?

Hi, @luc ,
I followed your recommendations and modified my code, but there is no positive effect.

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:

My modified code is below:

``````            File orekitData = new File("C:\\tempForOrekit");

// 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,
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

// Initialize numerical propagator
// By example from https://www.orekit.org/site-orekit-tutorials-11.1/tutorials/propagation.html
// 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);
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.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("----------");
``````

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.

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

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.

Thank you!