Hello Orekit forum!

I am doing some tests regarding numerical propagation and I found out a strange behavior in the results. I noticed that the order I add the different force models to the numerical propagator affects the results of the propagation.

I would like to know if this strange phenomenon is normal or not, and (in case it is normal) if there is a systematic way/order to add the different force models to get the best result.

I add the following snippet to let you reproduce the results.

```
Frame j2000 = FramesFactory.getEME2000();
double mu = Constants.IERS2010_EARTH_MU;
TimeScale utc = TimeScalesFactory.getUTC();
IERSConventions iersConventions = IERSConventions.IERS_2010;
boolean simpleEOP = true;
Frame ecef = FramesFactory.getITRF(iersConventions, simpleEOP);
double m = 300.0;
AbsoluteDate date = new AbsoluteDate(2020, 05, 01, 12, 00, 00.000, utc);
double a = 10000000.0;
double e = 0.0;
double i = 45.0;
double raan = 50.0;
double pa = 60.0;
double ta = 0.0;
KeplerianOrbit orbit = new KeplerianOrbit(a, e, FastMath.toRadians(i), FastMath.toRadians(pa), FastMath.toRadians(raan), FastMath.toRadians(ta), PositionAngleType.TRUE, j2000, date, mu);
SpacecraftState state = new SpacecraftState(orbit, m);
int order = 50;
int degree= 50;
NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(degree, order);
ForceModel gravity = new HolmesFeatherstoneAttractionModel(ecef, provider);
ForceModel moon = new ThirdBodyAttraction(CelestialBodyFactory.getMoon());
ForceModel sun = new ThirdBodyAttraction(CelestialBodyFactory.getSun());
double dragArea = 2.0;
double srpArea = 3.0;
double cD = 2.5;
double cR = 1.4;
OneAxisEllipsoid oneAxisEllipsoid = new OneAxisEllipsoid(Constants.IERS2010_EARTH_EQUATORIAL_RADIUS, Constants.IERS2010_EARTH_FLATTENING, ecef);
RadiationSensitive radiationSensitive = new IsotropicRadiationSingleCoefficient(srpArea, cR);
ForceModel srp = new SolarRadiationPressure(CelestialBodyFactory.getSun(), oneAxisEllipsoid, radiationSensitive);
NRLMSISE00InputParameters parameters = new CssiSpaceWeatherData(CssiSpaceWeatherData.DEFAULT_SUPPORTED_NAMES);
Atmosphere atmosphere = new NRLMSISE00(parameters, CelestialBodyFactory.getSun(), oneAxisEllipsoid);
DragSensitive dragSensitive = new IsotropicDrag(dragArea, cD);
ForceModel drag = new DragForce(atmosphere, dragSensitive);
double minStep = 1e-4;
double maxStep = 1e2;
double positionTolerance = 1e-6;
OrbitType orbitType = OrbitType.CARTESIAN;
double [][] tolerances = NumericalPropagator.tolerances(positionTolerance, orbit, orbitType);
ODEIntegrator odeIntegrator = new DormandPrince853Integrator(minStep, maxStep, tolerances[0], tolerances[1]);
NumericalPropagator numericalPropagator = new NumericalPropagator(odeIntegrator);
numericalPropagator.setInitialState(state);
numericalPropagator.setOrbitType(orbitType);
// Set 1
numericalPropagator.addForceModel(gravity);
numericalPropagator.addForceModel(moon);
numericalPropagator.addForceModel(sun);
numericalPropagator.addForceModel(drag);
numericalPropagator.addForceModel(srp);
// Set 2
// numericalPropagator.addForceModel(moon);
// numericalPropagator.addForceModel(sun);
// numericalPropagator.addForceModel(gravity);
// numericalPropagator.addForceModel(drag);
// numericalPropagator.addForceModel(srp);
// Set 3
// numericalPropagator.addForceModel(moon);
// numericalPropagator.addForceModel(drag);
// numericalPropagator.addForceModel(sun);
// numericalPropagator.addForceModel(gravity);
// numericalPropagator.addForceModel(srp);
// Set 4
// numericalPropagator.addForceModel(moon);
// numericalPropagator.addForceModel(drag);
// numericalPropagator.addForceModel(srp);
// numericalPropagator.addForceModel(sun);
// numericalPropagator.addForceModel(gravity);
double tof = 15 * 24 * 3600.0;
SpacecraftState finalState = numericalPropagator.propagate(date.shiftedBy(tof));
```

I report the results in terms of Keplerian elements for the different sets of force models.

Set 1:

Epoch: 2020-05-16T12:00:00.000Z

a = 10004291.128031444

e = 0.000221236

i = 44.982032719

raan = 28.070068896

pa = 102.595153652

ta =160.101840116

Set 2:

Epoch: 2020-05-16T12:00:00.000Z

a = 10004291.128027441

e = 0.000221236

i = 44.982032719

raan = 28.070068896

pa = 102.595153705

ta = 160.101840136

Set 3:

Epoch: 2020-05-16T12:00:00.000Z

a = 10004291.128017593

e = 0.000221236

i = 44.982032719

raan = 28.070068896

pa = 102.595153832

ta = 160.101840189

Set 4:

Epoch: 2020-05-16T12:00:00.000Z

a = 10004291.127967408

e = 0.000221236

i = 44.982032719

raan = 28.070068896

pa = 102.595154509

ta = 160.101840458

I know that the orbit is circular and it is normal that the argument of perigee and true anomaly may differ, but I was not expecting differences in the semi-major axis. I add that these differences arise only when the gravity force model is sorted. If the gravity force model is the first force model to be added to the numerical propagator and the other force models are shuffled, all the orbital elements (also argument of perigee and true anomaly are exactly the same).

I also run an example considering an eccentricity equal to `e = 0.2`

and setting the `OrbitType.EQUINOCTIAL`

. The behaviour does not change for an elliptic orbit. I report the results for the different sets of force models.

Set 1:

Epoch: 2020-05-16T12:00:00.000Z

a = 10005876.549596444

e = 0.200138548

i = 44.982322448

raan = 26.201610444

pa = 85.303963999

ta =181.186988018

Set 2:

Epoch: 2020-05-16T12:00:00.000Z

a = 10005876.550832234

e = 0.200138548

i = 44.982322447

raan = 26.201610446

pa = 85.303963997

ta = 181.186985373

Set 3:

Epoch: 2020-05-16T12:00:00.000Z

a = 10005876.552108487

e = 0.200138548

i = 44.982322448

raan = 26.201610447

pa = 85.303964035

ta = 181.186981367

Set 4:

Epoch: 2020-05-16T12:00:00.000Z

a = 10005876.548659690

e = 0.200138548

i = 44.982322448

raan = 26.201610443

pa = 85.303963975

ta = 181.186990395

As you can see the order the force models are added to the numerical propagator is changing the semi-major axis result of some mm. If this is expected, what is the correct order to add the force models to generate the most accurate result?

I thank you all in advance for any help.

Regards,

Marco