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