Comparison between DSST and numerical propagator

Hello everyone,

I’m comparing the results obtained from a numerical propagation to those obtained from a DSST (mean) propagation. However, I found some results which I’m not sure how to interpret: you can see that mean SMA and eccentricity obtained from DSST are somewhat capturing the behaviour of the osculating elements. However I don’t understand what’s happening with the inclination.
numerical vs dsst

The x-axis is the time expressed in seconds, with a total propagation time of 1 year, and the orbit is a 550km-altitude Sun-Synchronous orbit. I don’t know if this can be causing any issue but for DSSTPropagator I’m using the same integrator that I use for NumericalPropagator (with maxStep = 1e+4). Also, the force models used in both propagators are:

  • HolmesFeatherston(4,4) for numerical, Zonal(4,4) and Tesseral(4,4) for DSST
  • Sun and Moon third body
  • Drag
  • SRP

Do you think that there is something wrong? Also, is it correct that the mean eccentricity is pretty far from the osculating one for the first part of the graph? If I were to propagate for a much shorter period of time (say, 5 days), I would get this
numerical vs dsst short

Just by looking at this I would assume that there is something wrong with the eccentricity too. Is this normal?

Hi Emiliano,

I’m not sure this behavior is unreasonable, but I’m no semi-analytical expert, so let’s go through some checks about the API.

What do you mean exactly by the “same integrator”?
Is it the same scheme e.g. Dormand-Prince 8 or also the same min/max steps? 10000s is definitely too much for a NumericalPropagator.
On another note, do you pass the same initial SpacecraftState to the DSST as a MEAN or OSCULATING condition?
Finally, what happens if you output the osculating orbits with the DSST (with corresponding PropagationType at constructor level)?
Do you use the EphemerisGenerator by the way?

Best,
Romain.

Hello Romain,

The integrator is the same both in terms of scheme (DoPri853) and min/max steps. I know that 1e+4 is a bit large for a NumericalPropagator, but I initially tested it with 1e+2 (again for both) and the results are the same, so I think the integrator is able to select a suitable step for the NumericalPropagator, far from the max step allowed. The large maxStep is however very useful to DSSTPropagator to run significantly faster then the test with a maxStep of 1e+2.
The initial state is given to DSST as MEAN and the DSSTPropagator is set to PropagationType.MEAN, while to the NumericalPropagator I give the corresponding osculating state (computed via the computeOsculatingState of the DSSTPropagator).
I tried using the PropagationType.OSCULATING condition, these are the results
numerical vs dsst osc

Again, the inclination is capturing the short-term and secular variation, but not the long-term oscillations. I can provide a snippet of my code if you think that can be helpful.
Finally, yes, I’m using an EphemerisGenerator for each propagator.

Best regards,
Emiliano

Hello,

I also observed a small difference between numerical and DSST propagation on GEO orbit some time ago. After investigations, for my case, it seems linked to the moon perturbation. Without the moon, the two propagations are really close. But when I add the moon, I see a small difference that is more clearly visible on the evolution of the terrestrial longitude.

Nominal force models:
36x36 (not needed for a GEO orbit but set at this value for the test)
sun
moon
SRP

GEO propagation on 14 days with moon, there is a slight divergence on the terrestrial mean longitude.
with “mean” for DSST propagation type


with “osculating” for DSST propagation type

same GEO propagation on 14 days without moon. Here, the numerical and DSST propagations are really close.


The initial state of the numerical propagator is obtained with DSSTPropagator.computeOsculatingState on the mean state used to propagate with DSST.

Christophe

PS: on the last plot, it’s not the “true longitude” but the terrestrial longitude.

Hi guys,

@Emiliano you need to pass the initial state to DSST as OSCULATING. That should explain your difference on the short term
Moreover, I would suggest using different integrators between numerical and semi analytical. The latter could be given a fixed step scheme like RK4 with large steps.

@Christophe
Again I’m no semi analytical expert, but I would assume that the divergence on the anomaly (or its equivalent) to be bound to appear over time. Otherwise these models would have replaced full numerical ones no?
Maybe the Moon exacerbates the problem because the short periods terms are less accurate?

Best,
Romain.

Hello,

Why do I need to give the initial state as OSCULATING?
What I do in my code is:
→ I define an initial mean orbit
→ I construct the two propagators with the DSSTPropagator set to PropagationType.MEAN
→ I convert the initial mean orbit to an osculating one and give it to the NumericalPropagator
→ I feed the mean orbit to the DSSTPropagator as a MEAN state

I tried giving DSSTPropagator the same initial osculating orbit I give to NumericalPropagator and propagating with the MEAN condition, but the results are the same as if I give the initial mean orbit as MEAN…

Best regards,
Emiliano

My bad, I thought you were starting from an osculating orbit.

@Christophe I suspect the difference is caused by Week Time Dependent terms which are not taken into account in Orekit’s DSST Third body model.

There is an old issue opened to add those terms : Add weak-time-dependent (WTD) terms in the DSST lunar solar short period motion. (#576) · Issues · Orekit / Orekit · GitLab

You can find more theoretical details in the two following references:

  1. https://www.researchgate.net/publication/349311171_Towards_Accurate_Orbit_Determination_using_Semi-analytical_Satellite_Theory (Chapter “OREKIT DSST ORBIT DETERMINATION RESULTS”)
  2. https://elib.dlr.de/140658/1/Thesis_Setty-2020.pdf (Page 59 + Comparison between Fig. 4.3 without WTD and Fig. 4.7 with WTD)

To summarize, they are short period terms that are not present in Orekit DSST third body model and that improve considerably the consistency between DSST and numerical propagation models.

One of my “Orekit dream” is to find time to implement these terms :slight_smile:

@Emiliano for a LEO satellite you can use a fixed step of about 10 orbital periods

Thanks to both @Serrof and @bcazabonne for your suggestions on the integrator to use.

However I still have no clue as to why my DSSTPropagator isn’t correctly representing the evolution of the mean inclination. By playing a little with the force models, I found out that my problem is caused by the geopotential model. Indeed, if I account for all other perturbations except for the spherical expansion of the gravity field, the mean inclination correctly follows the osculating one.

Additionally, what’s even weirder to me is the fact that if I use computeMeanState of my DSSTPropagator on the osculating state that I obtain from NumericalPropagator, the mean inclination is correctly following the osculating inclination. This makes me think that the way I construct the force models is correct. This happens with as many force models as I want. No matter how many forces I account for, the state returned by the computeMeanState is always accurate. So basically I think that my DSSTPropagator is correctly converting between osculating and mean states, but for some reason is not correctly propagating

Perhaps you could check the body frame?
I think DSST uses a simplified model for Earth frame, with a separate Frame and centralBodyRotationRate at tesseral terms construction, whereas numerical propagator takes all small variations of the Earth frame into account (precession, nutation tides, …). If the frames do not match, you may experience drift.

One way to check this would be to use a very simplified Earth frame also for numerical propagator (perhaps using CelestialBodyFactory.getEarth().getBodyOrientedFrame() or something similar).

I tried using the frame you suggested but I still get the same results.

Here is my code, maybe that way it can be easier to locate a possible error:

n, m = 4, 4;
gravityProvider = GravityFieldFactory.getConstantNormalizedProvider(n,m);
nonSphericalGravity = HolmesFeatherstoneAttractionModel(CelestialBodyFactory.getEarth().getBodyOrientedFrame(), gravityProvider);        
forces = [nonSphericalGravity];

gravityProviderDsst = GravityFieldFactory.getConstantUnnormalizedProvider(n,m);
zonalDsst = DSSTZonal(gravityProviderDsst);
tesseralDsst = DSSTTesseral(CelestialBodyFactory.getEarth().getBodyOrientedFrame(),
                            Constants.WGS84_EARTH_ANGULAR_VELOCITY,
                            gravityProviderDsst);
forcesDsst = [zonalDsst, tesseralDsst];

minStep = 1e-2;
maxStep = 1e+3;
absTol = 1e-9;
relTol = 1e-9;
integrator = DormandPrince853Integrator(minStep, maxStep, absTol, relTol);
integrator.setInitialStepSize(60.0);

propagator = NumericalPropagator(integrator);
propagator.setResetAtEnd(False);
propagator.setAttitudeProvider(earthCenterAttitudeLaw);
for force in forces:
    propagator.addForceModel(force);
ephGen = propagator.getEphemerisGenerator();

initOrbitMean = KeplerianOrbit(a0, e0, i0, pa0, raan0, f0, PositionAngle.TRUE,
                               FramesFactory.getEME2000(), initDate, Constants.WGS84_EARTH_MU);
initStateMean = SpacecraftState(initOrbitMean, mass);

integratorDsst = ClassicalRungeKuttaIntegrator(10*initOrbitMean.getKeplerianPeriod());  
propagatorDsst = DSSTPropagator(integratorDsst, PropagationType.MEAN);
propagatorDsst.setResetAtEnd(False);
for force in forcesDsst:
    propagatorDsst.addForceModel(force);
ephGenDsst = propagatorDsst.getEphemerisGenerator();
initStateOsc = propagatorDsst.computeOsculatingState(initStateMean, earthCenterAttitudeLaw, propagatorDsst.getAllForceModels());

propagator.setInitialState(initStateOsc);
propagatorDsst.setInitialState(initStateMean, PropagationType.MEAN);

nDays = 14;
tSpan = nDays*day;
finalDate = initDate.shiftedBy(tSpan);

propagator.propagate(finalDate);
print("NumericalPropagator done");
propagatorDsst.propagate(finalDate);
print("DSSTPropagator done");

timeVec = linspace(0, tSpan, 30*nDays);
aOsc, eOsc, iOsc = [], [], [];
aMean, eMean, iMean = [], [], [];
for time in timeVec:
    dateNow = initDate.shiftedBy(float(time));
    pvOsc = ephGen.getGeneratedEphemeris().getPVCoordinates(dateNow, FramesFactory.getEME2000());
    pvMean = ephGenDsst.getGeneratedEphemeris().getPVCoordinates(dateNow, FramesFactory.getEME2000());
    
    orbitOsc = KeplerianOrbit(pvOsc, FramesFactory.getEME2000(), Constants.WGS84_EARTH_MU);
    orbitMean = KeplerianOrbit(pvMean, FramesFactory.getEME2000(), Constants.WGS84_EARTH_MU);
    
    aOsc.append(orbitOsc.getA());
    eOsc.append(orbitOsc.getE());
    iOsc.append(orbitOsc.getI());
    aMean.append(orbitMean.getA());
    eMean.append(orbitMean.getE());
    iMean.append(orbitMean.getI());

Hello,

why do you use a simple RK4 with a fixed step of 10 revolutions (ClassicalRungeKuttaIntegrator) for DSST and not DormandPrince853Integrator ?

Thanks Bryan,

You’re right, this is a good candidate for the differences I observed. I wasn’t aware that the model for the moon was not “full” in DSST.

Thanks!

Looking at this reference, it is stated that:

It is worth noting that the step size in the case of the slow-dynamics component integration assumes that the mean element motion frequencies due to the zonal harmonics in the geopotential are slow relative to the frequencies due to the lunar perturbations. The integration step for a LEO orbit is usually around half a day, which corresponds to several revolutions. However, as the orbital period increases, the step size remains the same because the lunar-solar frequencies do not change, although the number of steps per revolution may be different due to the variation in the revolution duration.

DoPri853 is probably just overkill and would be unnecessarily slower.

I don’t agree. You can use Dormand-Prince 8-5(3) with settings that correspond to the dynamics of the process. If the process is slow, you will simply end up with step sizes that can reach several days, there is nothing wrong with that. As the integrator performs 17 evaluations of the ODE at each step, if you have a step size of one day you just evaluate the ODE about once every 5000 seconds.

With semi-analytical theory, the numerical part managed by the integrator is not the one that is time consuming. The overhead is spread over everything else. Even just writing the output on a file may consume more cycles than integrating the ODE.

Interesting. So the recommendation of using a fixed-step integration scheme is not really something that must be respected at all costs?

I don’t think so.
If you are dealing with stable orbits, variable stepsize integrator would in fact end up selecting one size and stick to it. So you can just select whatever you want for integrating, and variable step size integrator may select something that could be counter-intuitive (several days).
If on the other hand you are computing orbits that change on the long term, for example to evaluate life time due to reentry with propagation that cover several years, then you will notice that the step size decreases as dynamics become more noisy before reentry.

Thanks for the clarification!

Anyway, about the issue with the mean vs osculating inclination, is there any gross error in the code that I posted before?

I suspect the problem is that you convert mean elements to Cartesian. Could you just use getOrbit on your propagated states (that might prevent the use of the ephemeris generator I have not checked)?
And preferably make comparisons in coordinates of the semi analytical theory so here equinoctial.

Romain.