Hello everyone,
I came accross a potential bug in the DSST propagator in osculating elements when used with at least one detector. It leads to hyperbolic orbits during the propagation.
To reproduce the bug, we added two apside detectors with two different thresholds:
@Test
public void dsstInOsculatingElementsYieldsHyperbolicOrbits() {
Utils.setDataRoot("regular-data:potential/grgs-format");
GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
// initial state
final AbsoluteDate orbitEpoch = new AbsoluteDate(2023, 2, 18, TimeScalesFactory.getUTC());
final Frame inertial = FramesFactory.getCIRF(IERSConventions.IERS_2010, true);
final KeplerianOrbit orbit = new KeplerianOrbit(42166000.0, 0.00028, FastMath.toRadians(0.05), FastMath.toRadians(66.0),
FastMath.toRadians(270.0), FastMath.toRadians(11.94), PositionAngleType.MEAN,
inertial, orbitEpoch, Constants.WGS84_EARTH_MU);
final EquinoctialOrbit equinoctial = new EquinoctialOrbit(orbit);
// create propagator
final double[][] tol = DSSTPropagator.tolerances(0.001, equinoctial);
final AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(3600.0, 86400.0, tol[0], tol[1]);
final DSSTPropagator propagator = new DSSTPropagator(integrator, PropagationType.OSCULATING);
// add force models
final Frame ecefFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
final UnnormalizedSphericalHarmonicsProvider gravityProvider = GravityFieldFactory.getUnnormalizedProvider(8, 8);
propagator.addForceModel(new DSSTZonal(gravityProvider));
propagator.addForceModel(new DSSTTesseral(ecefFrame, Constants.WGS84_EARTH_ANGULAR_VELOCITY, gravityProvider));
propagator.addForceModel(new DSSTThirdBody(CelestialBodyFactory.getSun(), gravityProvider.getMu()));
propagator.addForceModel(new DSSTThirdBody(CelestialBodyFactory.getMoon(), gravityProvider.getMu()));
// add two dummy detectors with different thresholds
propagator.addEventDetector(new ApsideDetector(1e-3, equinoctial).withHandler((spacecraftState, detector, increasing) -> {
return Action.CONTINUE;
}));
propagator.addEventDetector(new ApsideDetector(1e-4, equinoctial).withHandler((spacecraftState, detector, increasing) -> {
return Action.CONTINUE;
}));
// propagation leads to hyperbolic orbits
propagator.setInitialState(new SpacecraftState(equinoctial, 6000.0), PropagationType.MEAN);
propagator.propagate(orbitEpoch.shiftedBy(20.0 * Constants.JULIAN_DAY));
}
This test is an extension of testIssue1029() of src/test/java/org/orekit/propagation/semianalytical/dsst/DSSTPropagatorTest.java · develop · Orekit / Orekit · GitLab with two detectors.
Note that the bug also randomly appears in other contexts, with only one detector.
The root cause seems to be a bad interpolation during the computation of the short periodic contribution of either zonal or tesseral geopotential terms, when evaluated during the “mean to osculating” conversion of the state in
org.orekit.propagation.semianalytical.dsst.DSSTPropagator.MeanPlusShortPeriodicMapper#mapArrayToState.
I can see that, at this level of the code, the short periodic terms are computed using slots that have been previously updated using org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel#updateShortPeriodTerms(double, org.orekit.propagation.SpacecraftState…) in org.orekit.propagation.semianalytical.dsst.DSSTPropagator.ShortPeriodicsHandler#handleStep.
I wonder why there is no immediate call to updateShortPeriodTerms in DSSTPropagator.MeanPlusShortPeriodicMapper#mapArrayToState right before the conversion to make sure that the slots are correctly updated, as it is done in org.orekit.propagation.semianalytical.dsst.DSSTPropagator#computeOsculatingState ?
Thanks in advance!
Anne-Laure, on behafl of Bastien Le Bihan (he failed to create its own account due to too restrictive company mails filter)
on the issue description.