Hi all! We came accross a failing of the root finding associated with event detection, when the following conditions are met:
-
The force field is non keplerian
-
The event we want to catch is very close to the current integration step: if the g function is going from positive to negative, then g is marginally negative at the current integration step, we are “right after the event”.
-
The detector is going through an event shifter with quite large time shift.
In the hipparchus routine DetectorBasedEventState#evaluateStep, a change of sign is detected and a root finding is triggered. Inside the root finder, the g function is again evaluated at the current integration step. However, the state used is not exactly the same as in when the change of sign is detected: the later is the soft current state sb, while the former is interpolator.getInterpolatedState(sb.getTime()). For some reason, these two states are not exactly equivalent: the primary state are equals to the last decimals, but not the primary derivatives. When these two states go through the event shifter with a quite large time shift, this non-keplerian difference yields variations on the 6-dimensional state, big enough to loose the change of sign, and an exception is raised.
Hereafter is a unit test to reproduce the bug:
@Test
public void reproduceBracketErrorBug() {
ClassicalRungeKuttaIntegratorBuilder integratorBuilder = new ClassicalRungeKuttaIntegratorBuilder(1.);
double positionScale = 10e3;
double a = 6879290;
double e = 0.0017;
double i = 1.7;
double pa = -5.3;
double raan = 0.42;
double anomaly = 3.04;
AbsoluteDate date = AbsoluteDate.ARBITRARY_EPOCH;
Frame frame = DataContext.getDefault().getFrames().getTOD(false);
KeplerianOrbit orbit = new KeplerianOrbit(a, e, i, pa, raan, anomaly, PositionAngleType.TRUE, frame, date, Constants.EGM96_EARTH_MU);
NumericalPropagatorBuilder numericalPropagatorBuilder = new NumericalPropagatorBuilder(orbit, integratorBuilder, PositionAngleType.TRUE, positionScale);
NumericalPropagator numericalPropagator = numericalPropagatorBuilder.buildPropagator();
HolmesFeatherstoneAttractionModel geoPotentialModel = new HolmesFeatherstoneAttractionModel(DataContext.getDefault().getFrames().getITRF(IERSConventions.IERS_2010, false),
DataContext.getDefault().getGravityFields().getNormalizedProvider(2, 0));
numericalPropagator.addForceModel(geoPotentialModel);
EventDetectionSettings detectionSettings = new EventDetectionSettings(AbstractDetector.DEFAULT_MAX_CHECK, AbstractDetector.DEFAULT_THRESHOLD,
AbstractDetector.DEFAULT_MAX_ITER);
EventOnIntegratorTimeStepDetector detector = new EventOnIntegratorTimeStepDetector(detectionSettings);
EventShifter eventShifter = new EventShifter(detector, false, 500, 500);
numericalPropagator.addEventDetector(eventShifter);
numericalPropagator.propagate(date.shiftedBy(10));
}
private class EventOnIntegratorTimeStepDetector extends AbstractDetector {
protected EventOnIntegratorTimeStepDetector(EventDetectionSettings detectionSettings) {
super(detectionSettings, new ContinueOnEvent());
}
@Override
public double g(SpacecraftState s) {
double exactY = -2391379.2800254025;
return s.getOrbit().getPosition().getY() - exactY + 1e-8;
}
@Override
protected AbstractDetector create(EventDetectionSettings detectionSettings, EventHandler newHandler) {
return new EventOnIntegratorTimeStepDetector(detectionSettings);
}
}
The expected stacktrace is
org.orekit.errors.OrekitException: org.orekit.propagation.integration.AbstractIntegratedPropagator$AdaptedEventDetector@563f38c4 failed to find root between 1 (g=-1.854245212844352E3) and 2 (g=1.0E-8)
Last iteration at 2 (g=-1.441338860988617E-7)
at org.orekit.errors.OrekitException.unwrap(OrekitException.java:154)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(AbstractIntegratedPropagator.java:572)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:475)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:435)
at org.orekit.propagation.numerical.NumericalPropagatorTest.reproduceBracketErrorBug(NumericalPropagatorTest.java:132)
at java.base/java.lang.reflect.Method.invoke(Method.java:569)
at java.base/java.util.ArrayList.forEach(ArrayList.java:1511)
at java.base/java.util.ArrayList.forEach(ArrayList.java:1511)
Caused by: org.hipparchus.exception.MathIllegalStateException: org.orekit.propagation.integration.AbstractIntegratedPropagator$AdaptedEventDetector@563f38c4 failed to find root between 1 (g=-1.854245212844352E3) and 2 (g=1.0E-8)
Last iteration at 2 (g=-1.441338860988617E-7)
at org.hipparchus.ode.events.DetectorBasedEventState.findRoot(DetectorBasedEventState.java:373)
at org.hipparchus.ode.events.DetectorBasedEventState.evaluateStep(DetectorBasedEventState.java:231)
at org.hipparchus.ode.AbstractIntegrator.lambda$acceptStep$4(AbstractIntegrator.java:328)
at java.base/java.util.ArrayList$ArrayListSpliterator.forEachRemaining(ArrayList.java:1625)
at java.base/java.util.stream.Streams$ConcatSpliterator.forEachRemaining(Streams.java:734)
at java.base/java.util.stream.ReferencePipeline$Head.forEach(ReferencePipeline.java:762)
at org.hipparchus.ode.AbstractIntegrator.acceptStep(AbstractIntegrator.java:328)
at org.hipparchus.ode.nonstiff.FixedStepRungeKuttaIntegrator.integrate(FixedStepRungeKuttaIntegrator.java:168)
at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(AbstractIntegratedPropagator.java:550)
... 6 more
Caused by: org.hipparchus.exception.MathIllegalArgumentException: interval does not bracket a root: f(1E0) = -1.854245212844352E3, f(2E0) = -144.13388609886168E-9
at org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver.doSolveInterval(BracketingNthOrderBrentSolver.java:209)
at org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver.solveInterval(BracketingNthOrderBrentSolver.java:427)
at org.hipparchus.analysis.solvers.BracketedUnivariateSolver.solveInterval(BracketedUnivariateSolver.java:126)
at org.hipparchus.ode.events.DetectorBasedEventState.findRoot(DetectorBasedEventState.java:364)
... 14 more