Event States not Equivalent to Propagated States

Orekit usage
1

I’m running into a strange issue while using event detectors. I’ve created an example that is similar enough to my use case. I’m using a DateDetector with a RecordAndContinue handler to store states at a certain time step. But the states that are stored do not match what I get if I just propagate to that epochs directly.

In my example, I first propagate to each test epoch and calculate the range between the states using two identical propagators. I get back 0 as expected. Next, I add DateDetectors with a RecordAndContinue handler to one of the propagators. I propagate to the stop time and record all the events. For each event, I then propagate to the event epoch with a second propagator that is identical to the first, but it does not have event detectors. When I calculate the range between the two states (the event state and the propagated state), I get a non-zero result.

@Log4j2
public class OrekitTest {

  @Test
  public void testOrekit() {

    // Initialize Orekit.
    File orekitData = new File("PATH_TO_OREKIT_DATA");
    DataProvidersManager manager = DataContext.getDefault().getDataProvidersManager();
    manager.addProvider(new DirectoryCrawler(orekitData));

    // Establish time range.
    TimeScale utc = TimeScalesFactory.getUTC();
    AbsoluteDate initialEpoch = new AbsoluteDate("2023-09-01T00:00:00.000", utc);
    AbsoluteDate finalEpoch = new AbsoluteDate("2023-09-03T00:00:00.000", utc);

    // Set initial state.
    TimeStampedPVCoordinates pv = new TimeStampedPVCoordinates(initialEpoch,
        new Vector3D(3.202221011875539e+07, 2.743293506737402e+07, -7.418022649217934e+04),
        new Vector3D(-2.000287709983509e+03, 2.334929394216369e+03, 4.490568450445024e-00));
    KeplerianOrbit orbit = new KeplerianOrbit(pv, FramesFactory.getGCRF(),
        Constants.IERS2010_EARTH_MU);
    SpacecraftState initialState = new SpacecraftState(orbit);

    // Create propagators.
    Propagator prop = createPropagator(initialState);
    Propagator propNoDetectors = createPropagator(initialState);

    // Get the epochs to propagate to.
    double step = 6 * 3600.0; // 6 hours.
    List<AbsoluteDate> epochs = new ArrayList<>();
    for (AbsoluteDate epoch = initialEpoch; epoch.compareTo(finalEpoch.shiftedBy(-step)) < 0; ) {
      epoch = epoch.shiftedBy(step);
      epochs.add(epoch);
    }

    // Propagate to dates and get range.
    for (AbsoluteDate epoch : epochs) {
      SpacecraftState state1 = prop.propagate(epoch);
      SpacecraftState state2 = propNoDetectors.propagate(epoch);
      double range = calculateRange(state1, state2);
      log.error("Range: {} m", range);
    }

    // Re-create the propagators.
    prop = createPropagator(initialState);
    propNoDetectors = createPropagator(initialState);

    // Add the date detectors.
    RecordAndContinue<EventDetector> rac = new RecordAndContinue<>();
    for (AbsoluteDate epoch : epochs) {
      prop.addEventDetector(new DateDetector(epoch).withHandler(rac));
    }

    // Propagate (and record events).
    prop.propagate(finalEpoch);

    // Compare states at event times.
    for (int i = 0; i < rac.getEvents().size(); ++i) {

      Event<?> event = rac.getEvents().get(i);
      AbsoluteDate epoch = epochs.get(i);
      double epochDifference = event.getState().getDate().durationFrom(epoch);
      log.error("Epoch Difference: {} sec", epochDifference);

      AbsoluteDate eventEpoch = event.getState().getDate();
      SpacecraftState state = propNoDetectors.propagate(eventEpoch);
      double range = calculateRange(event.getState(), state);
      log.error("Range: {} m", range);
    }

  }

  private double calculateRange(SpacecraftState first, SpacecraftState second) {
    return first.getPVCoordinates().getPosition().subtract(second.getPVCoordinates().getPosition())
        .getNorm();
  }

  private Propagator createPropagator(SpacecraftState initialState) {

    // Create the integrator.
    double minStep = initialState.getKeplerianPeriod() / 100.0;
    double maxStep = initialState.getKeplerianPeriod();
    double[][] tol = DSSTPropagator.tolerances(10.0, initialState.getOrbit());
    AdaptiveStepsizeIntegrator integrator =
        new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);

    // Initialize the propagator.
    NumericalPropagator propagator = new NumericalPropagator(integrator);
    propagator.setInitialState(initialState);

    // Add J22 perturbations.
    Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
    HolmesFeatherstoneAttractionModel gravityModel = new HolmesFeatherstoneAttractionModel(itrf,
        GravityFieldFactory.getNormalizedProvider(2, 2));

    // Add the force models to the propagator
    propagator.addForceModel(gravityModel);

    return propagator;
  }

}

It seems to have something to do with the event detection and handling itself. Any idea what’s going on here? Maybe I’m getting the event state incorrectly? Thanks in advance!

2

Hello there,

Well the event detection affect the stepsize of the propagator, so you’re getting different integration noise with and without them. How big are we talking here with your non-zero difference?
Side note: why do you use the DSST routines for the tolerances? The numerical propagator has its own.

Best,
Romain.

3

I’m getting differences on the order of tens of meters. I thought that was a little bit too large for just some numerical precision issues or integrator noise. Here’s what I get when running that sample:

Range at 2023-09-01T06:00:00.000Z: 1.2061311055629265 m
Range at 2023-09-01T12:00:00.000Z: 19.34622091652894 m
Range at 2023-09-01T18:00:00.000Z: 25.024743341962665 m
Range at 2023-09-02T00:00:00.000Z: 21.396167066155222 m
Range at 2023-09-02T06:00:00.000Z: 99.42941722453584 m
Range at 2023-09-02T12:00:00.000Z: 77.07795291428513 m
Range at 2023-09-02T18:00:00.000Z: 78.80801529463407 m

And no reason for the DSST tolerances. Just an oversight by me when I was making this post :slight_smile:

Edit:
When I changed the tolerances initialization to this:

double[][] tol = NumericalPropagator.tolerances(10.0, initialState.getOrbit(),
        OrbitType.CARTESIAN);

I actually got larger range differences:

Range at 2023-09-01T06:00:00.000Z: 92.06407346257473 m
Range at 2023-09-01T12:00:00.000Z: 389.23255346679235 m
Range at 2023-09-01T18:00:00.000Z: 456.2105038438005 m
Range at 2023-09-02T00:00:00.000Z: 137.64949024330465 m
Range at 2023-09-02T06:00:00.000Z: 2006.97414159891 m
Range at 2023-09-02T12:00:00.000Z: 2079.406474987599 m
Range at 2023-09-02T18:00:00.000Z: 1829.9862702700873 m
4

The larger difference with the new tolerances is because you computed them for Cartesian coordinates while your integration variables are Keplerian (inferred from the orbit you passed to the SpacecraftState initializing your propagator).
Edit: also, your maximum stepsize is really big for a full numerical model.

Romain.

5

@Serrof Thanks. I made those changes, but my primary concern is still why there is a difference at all. The epochs in the event handler are correct, but the states are significantly off. Even with different step sizes between the detector/non-detector propagators, I wouldn’t expect to see such large differences in states.

6

It seems to me that in one case you are performing several propagations with intermediate end dates set to the epochs and in the other case you are performing one propagation with one end date.
This, associated with the very large tolerance you put on the propagators may generate quite large noise.

The position tolerance value in the calls to {DSST|Numerical}Propagator.tolerances is a local error, not a global error and it is an estimation made by the underlying integrator. This local error is therefore an approximation valid only at the end of each step. When several steps are performed, the error will often grow up.

Another factor explaining large errors is linked to events. When an event is triggered, it generally occurs within one step, not at the end. The integrator underlying models are optimized to reduce error at end of step only (this is done by solving the order conditions), not in the middle of the step. This implies that even when using an order 8 integrator that will estimated a local error at end of step of the form k_1 \times h^8, the state vector error when event is triggered may be something like k_2\times h^7 or worse. If the propagator is stopped at this date and then restarted, the intermediate point will be quite bad and errors will grow a lot.

This is the reason why, when you have discontinuities in the dynamics that are handled by events, you should set a tolerance that is really really small. In a recent study I made when dealing with eclipses, I had to use tolerances of the order of magnitude of 10^{-6} m or even smaller in order to get really fine results in some comparison (I was dealing with orbits at a few centimeters level). In your case, 10m on the tolerance, plus events that interrupt the propagation, plus restarts from approximated states all have nasty effects on the accuracy and cannot be compared to a single run propagation.

7

@luc Thanks for the explanation. That makes sense to me. Tightening the tolerances in my propagator settings makes my states much less noisy. And, for my simple case, it didn’t seem to have much of an impact on runtime. :slight_smile:

Thanks @luc and @Serrof for the help!