Fixed event bracketing problem induced by numerical noise.

Fixes #921 Merge branch 'issue-919' into develop

Fixed event bracketing problem induced by numerical noise.
Fixes #921 Merge branch 'issue-919' into develop
d8d2ee6d · Luc Maisonobe · 41551176 · d8d2ee6d · d8d2ee6d · d8d2ee6d
Commit d8d2ee6d authored Apr 22, 2022 by Luc Maisonobe
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -21,6 +21,9 @@
  </properties>
  <body>
    <release version="11.2" date="TBD" description="TBD">
+      <action dev="luc" type="fix" issue="921">
+        Fixed event bracketing problem induced by numerical noise at end of search interval.
+      </action>
      <action dev="luc" type="fix" issue="919">
        Fixed ephemeris generation with several derivatives providers.
      </action>

--- a/src/main/java/org/orekit/propagation/events/EventState.java
+++ b/src/main/java/org/orekit/propagation/events/EventState.java
@@ -312,11 +312,25 @@ public class EventState<T extends EventDetector> {
                // both non-zero, the usual case, use a root finder.
                // time zero for evaluating the function f. Needs to be final
                final AbsoluteDate fT0 = loopT;
+                final double tbDouble = tb.durationFrom(fT0);
+                final double middle = 0.5 * tbDouble;
                final UnivariateFunction f = dt -> {
-                    return g(interpolator.getInterpolatedState(fT0.shiftedBy(dt)));
+                    // use either fT0 or tb as the base time for shifts
+                    // in order to ensure we reproduce exactly those times
+                    // using only one reference time like fT0 would imply
+                    // to use ft0.shiftedBy(tbDouble), which may be different
+                    // from tb due to numerical noise (see issue 921)
+                    final AbsoluteDate t;
+                    if (forward == dt <= middle) {
+                        // use start of interval as reference
+                        t = fT0.shiftedBy(dt);
+                    } else {
+                        // use end of interval as reference
+                        t = tb.shiftedBy(dt - tbDouble);
+                    }
+                    return g(interpolator.getInterpolatedState(t));
                };
                // tb as a double for use in f
-                final double tbDouble = tb.durationFrom(fT0);
                if (forward) {
                    try {
                        final Interval interval =

--- a/src/test/java/org/orekit/propagation/events/EventDetectorTest.java
+++ b/src/test/java/org/orekit/propagation/events/EventDetectorTest.java
@@ -354,6 +354,27 @@ public class EventDetectorTest {

    }

+    @Test
+    public void testNumericalNoiseAtIntervalEnd() {
+
+        Frame eme2000 = FramesFactory.getEME2000();
+        TimeScale utc = TimeScalesFactory.getUTC();
+        final AbsoluteDate initialDate   = new AbsoluteDate(2011, 5, 11, utc);
+        final Orbit orbit = new EquinoctialOrbit(new PVCoordinates(new Vector3D(4008462.4706055815, -3155502.5373837613, -5044275.9880020910),
+                                                                   new Vector3D(-5012.9298276860990, 1920.3567095973078, -5172.7403501801580)),
+                                                 eme2000, initialDate, Constants.WGS84_EARTH_MU);
+        Propagator propagator = new KeplerianPropagator(orbit);
+        double base  = 3600.0;
+        double noise = FastMath.scalb(base, -60);
+        // introduce some numerical noise by using two separate shifts
+        AbsoluteDate finalTime = initialDate.shiftedBy(base).shiftedBy(2 * noise);
+        AbsoluteDate eventTime = finalTime.shiftedBy(-noise);
+        propagator.addEventDetector(new DateDetector(eventTime).withMaxCheck(base).withThreshold(noise / 2));
+        SpacecraftState finalState = propagator.propagate(finalTime);
+        Assert.assertEquals(0.0, finalState.getDate().durationFrom(eventTime), noise);
+
+    }
+
    @Test
    public void testForwardAnalytical() {
        doTestScheduling(0.0, 1.0, 21, this::buildAnalytical);