Detector order modify propagation states behavior

Hello,

We are using numerical propagator with several node detector/handler pairs. We were facing a strange behavior and after investigations we found that swapping detectors definition resolved our problem.

We managed to show the problem with the code below. In our case, we define 2 detectors : one is counting number of orbit on ascending node, the other is logging each ascending/descending node.
The first test method defines a propagator with a node logger then a revolutions counter: it counts properly the number of orbit as seen in the console log and the passed assertions.
The other test method defines a propagator with a revolutions counter then a node logger : it goes twice through a node (seen at 2008-11-15T16:35:42.583 in the console log) and so the revolutions count is too high at the end.

////////////////////////////////////////////////////////////////////////////////
// Copyright Airbus Defence and Space (c), 2019
////////////////////////////////////////////////////////////////////////////////
package net.airbus.ds.orekit;

import org.hipparchus.ode.events.Action;
import org.junit.jupiter.api.Assertions;
import org.junit.jupiter.api.BeforeAll;
import org.junit.jupiter.api.BeforeEach;
import org.junit.jupiter.api.Test;
import org.orekit.data.DataProvidersManager;
import org.orekit.data.ZipJarCrawler;
import org.orekit.forces.gravity.potential.GravityFieldFactory;
import org.orekit.frames.FactoryManagedFrame;
import org.orekit.frames.FramesFactory;
import org.orekit.orbits.KeplerianOrbit;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngle;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.conversion.DormandPrince853IntegratorBuilder;
import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
import org.orekit.propagation.events.NodeDetector;
import org.orekit.propagation.events.handlers.EventHandler;
import org.orekit.propagation.integration.AdditionalEquations;
import org.orekit.propagation.numerical.NumericalPropagator;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScalesFactory;

import java.io.File;
import java.net.URISyntaxException;
import java.net.URL;
import java.net.URLClassLoader;
import java.util.Arrays;
import java.util.Optional;

/**
 * Bug raised : https://forum.orekit.org/t/detector-order-modify-propagation-states-behavior/626
 */
public class ResetStateOrderTest {

    private static final String STATE_NAME = "revolutions";
    private static final int STATE_INDEX = 0;

    private static final FactoryManagedFrame FRAME = FramesFactory.getTOD(false);
    public static final double MU = GravityFieldFactory.getNormalizedProvider(0, 0).getMu();

    private NumericalPropagatorBuilder builder;
    private Orbit orbit;

    public class RevolutionsCounterHandler implements EventHandler<NodeDetector> {

        private int count;

        public RevolutionsCounterHandler() {
            super();
        }

        @Override
        public Action eventOccurred(SpacecraftState s, NodeDetector detector, boolean increasing) {
            if (increasing) {
                if (detector.isForward()) {
                    count = (int) s.getAdditionalState(STATE_NAME)[STATE_INDEX] + 1;
                } else {
                    count = (int) s.getAdditionalState(STATE_NAME)[STATE_INDEX] - 1;
                }
                System.out.println(s.getDate() + "\t" + count);
                return Action.RESET_STATE;
            }
            return Action.CONTINUE;
        }

        @Override
        public SpacecraftState resetState(NodeDetector detector, SpacecraftState oldState) {
            return oldState.addAdditionalState(STATE_NAME, count);
        }
    }

    public class RevolutionsNbEquations implements AdditionalEquations {

        private final String name;

        public RevolutionsNbEquations(String name) {
            this.name = name;
        }

        @Override
        public String getName() {
            return this.name;
        }

        @Override
        public double[] computeDerivatives(SpacecraftState s, double[] pDot) {
            return null;
        }
    }

    private class LogHandler implements EventHandler<NodeDetector> {

        public LogHandler() {
        }

        @Override
        public Action eventOccurred(SpacecraftState s, NodeDetector detector, boolean increasing) {
            System.out.println(s.getDate() + "\t" + increasing);
            return Action.CONTINUE;
        }
    }

    @BeforeAll
    public static void loadOrekit() throws URISyntaxException {
        Optional<URL> orekitData = Arrays.stream(((URLClassLoader) ClassLoader.getSystemClassLoader()).getURLs()).filter(url -> url.getFile().contains("orekit-data")).findFirst();
        DataProvidersManager.getInstance().addProvider(new ZipJarCrawler(new File(orekitData.get().getPath())));
    }

    @BeforeEach
    public void initializeBuilder() throws URISyntaxException {
        DormandPrince853IntegratorBuilder integratorBuilder = new DormandPrince853IntegratorBuilder(0.001, 300, 0.002);
        orbit = OrbitType.EQUINOCTIAL.convertType(new KeplerianOrbit(7000000, 0, Math.toRadians(98), 0, 0, 0, PositionAngle.MEAN, FRAME, new AbsoluteDate(2008, 11, 15, 8, 30, 0, TimeScalesFactory.getUTC()), MU));
        builder = new NumericalPropagatorBuilder(orbit, integratorBuilder, PositionAngle.TRUE, 10000);
        builder.setMass(3500);
    }

    @Test
    public void testNodesBeforeCounting() {
        NumericalPropagator propagator = builder.buildPropagator(builder.getSelectedNormalizedParameters());
        addNodeDetector(propagator);
        addRevolutionsCounter(propagator);
        SpacecraftState state = propagator.propagate(orbit.getDate().shiftedBy(3600 * 12));
        Assertions.assertEquals(7, (int) state.getAdditionalState(STATE_NAME)[STATE_INDEX]);
    }

    @Test
    public void testCountingBeforeNodes() throws Exception {
        NumericalPropagator propagator = builder.buildPropagator(builder.getSelectedNormalizedParameters());
        addRevolutionsCounter(propagator);
        addNodeDetector(propagator);
        SpacecraftState state = propagator.propagate(orbit.getDate().shiftedBy(3600 * 12));

        Assertions.assertEquals(7, (int) state.getAdditionalState(STATE_NAME)[STATE_INDEX]);
    }


    private void addNodeDetector(NumericalPropagator propagator) {
        propagator.addEventDetector(new NodeDetector(orbit, FRAME).withHandler(new LogHandler()));
    }

    private void addRevolutionsCounter(NumericalPropagator propagator) {
        propagator.addAdditionalEquations(new RevolutionsNbEquations(STATE_NAME));
        propagator.addEventDetector(new NodeDetector(orbit, FRAME).withHandler(new RevolutionsCounterHandler()));
        propagator.setInitialState(propagator.getInitialState().addAdditionalState(STATE_NAME, 0));
    }
}

We also replaced RESET_STATE action by a CONTINUE action and it does not go twice through the same node.

May we have misused detectors ?

Best regards,

Anne-Laure

The quick solution is don’t use RESET_STATE.

If you want to see what is happening add an OrekitStepHandler and print the satellite’s z coordinate around the event. Since your tolerance is about 1e-10 you’ll want to print out states within +/- 2e-10 of the detected event time with a step size of about 1e-11.

What I expect you will find is that the states from before and after the RESET_STATE are different for times after the event time. This is because the RESET_STATE action forces the propagator to recompute the trajectory after the event time. Resetting the state introduces some error because the interpolation error at the event time is now included in all future states. For maneuvers this interpolation error is usually less than the uncertainty in the actual delta-v applied.

I was not able to run your example due to some missing classes, but here is my guess as to what happened in more detail. An event is detected and the event time is determined to be t_e which is less than t_0, the time of the actual root with the constraint that |t_e - t_0| < tol, where tol is the tolerance set for the event detector (1e-10 in this case). The event detector will not look for another zero crossing until t_e + tol, which is after t_0. When you reset state the numerical propagator computes a new trajectory with a slightly different orbital state due to the interpolation error. I’m guessing this new trajectory slightly changes the time of the zero such that it is greater than t_e + tol at which point the event detector is allowed to find another root and it does.

Thank you Evan for your fast reply.

I will keep the solution by setting a larger tolerance as we need to edit states.
I reedit the first code with no custom class (sorry!) and here the dependencies, just in case :

<dependency>
        <groupId>org.orekit</groupId>
        <artifactId>orekit</artifactId>
        <version>10.0</version>
    </dependency>    
    <dependency>
        <groupId>org.junit.jupiter</groupId>
        <artifactId>junit-jupiter-engine</artifactId>
        <version>5.4.2</version>
        <scope>test</scope>
    </dependency>
    <dependency>
        <groupId>org.junit.jupiter</groupId>
        <artifactId>junit-jupiter-api</artifactId>
        <version>5.4.2</version>
        <scope>test</scope>
    </dependency>
    <dependency>
        <groupId>org.junit.jupiter</groupId>
        <artifactId>junit-jupiter</artifactId>
        <version>5.4.2</version>
        <scope>test</scope>
    </dependency>

Best regards,

Anne-Laure