Orekit v12 KeplerianPropagator propagateOrbit while-loop?

[Ryan]

I’ve long wondered about the loop in KeplerianPropagator::propagateOrbit and, in particular, the comment -

// we use a loop here to compensate for very small date shifts error
// that occur with long propagation time

- that the loop is required due to loss of precision with long propagation durations, but I don’t think that this is accurate. (This is present in both v12 and v13, at least.)

I’m surprised this loop exists at all, as it should be the responsibility of the shiftedBy impl.s to ensure that the shift worked as intended. And I don’t think I’ve seen equivalent loops in any other propagation.

I digress: black-box testing and debugging (see below for black-box testing/debugging details) suggests that, in fact, this is/seems to be due to the AbsoluteDate(final AbsoluteDate since, final double elapsedDuration) constructor (v12.2).

Moreover, inspection of the Orbit::shiftedBy implementations - checked KeplerianOrbit, CartesianOrbit, CircularOrbit, and EquinoctialOrbit - reveals no inappropriate date manipulation, so I think it can only be due to the aforementioned constructor.

I see that in v13 AbsoluteDate is now backed by TimeOffset which strictly uses longs to keep track of everything (seconds and attoseconds, IIRC) - which is great to see - so I wonder whether this issue is still present.

For v <= 12, I’m wondering whether the AbsoluteDate constructor should be updated with a loop like we see in propagateOrbit.

Expand for full AbsoluteDate black-box testing procedure

Using

        final TimeScale scale1 = ...;
        final TimeScale scale2 = ...;
        
        final int y = 2000;
        final AbsoluteDate t1 = new AbsoluteDate(new DateComponents(y, 1, 1),
                                                 TimeComponents.H12,
                                                 scale1); // J2000 defn if scale1 == TT
        
        System.out.println("start");
        for (int dy = 0; dy < 100; dy++) {
            for (int m = 1; m <= 12; m++) {
                final AbsoluteDate t2 = new AbsoluteDate(new DateComponents(2000 + dy, m, 1),
                                                         TimeComponents.H12,
                                                         scale2);
                
                final AbsoluteDate shifted = t1.shiftedBy(t2.durationFrom(t1));
                final double       error   = t2.durationFrom(shifted);
                
                if (error != 0.0) {
                    System.out.format("[%s, %s], [%3d, %3d]: % e%n", scale1, scale2, dy, m, error);
                }
            }
        }
        
        System.out.println("end");

with two TimeScales that have a constant offset - or using two TimeScales that are equal - error == 0 for each iteration.

However, using a pair like [TT, UTC]:

final TimeScale scale1 = TimeScalesFactory.getTT();
final TimeScale scale2 = TimeScalesFactory.getUTC();

We see the first non-zero error in month 2, and it’s non-zero for each moth from then on:

start
[TT, UTC], [  0,   2]:  1.080309e-10
[TT, UTC], [  0,   3]: -3.576304e-10
...
[TT, UTC], [  4,   5]:  1.454353e-08
...
[TT, UTC], [  8,   8]: -1.525879e-08
...
end

Of course my first thought is that this is really due to something w.r.t. leap seconds, and testing with TT and TAI (constant offset and zero offset) more-or-less confirmed this, as there is never a non-zero difference.

Adding a few extra pieces to the if (error != 0.0) block for easier debugging:

final var dt = t2.durationFrom(t1);
final var sd = t1.shiftedBy(dt);

t2.durationFrom(sd);

and stepping in, we arrive at the aforementioned AbsoluteDate constructor, where the value of regularOffset always matches - in magnitude but not necessarily sign - the error in the output.

[luc]

You are right.
This loop was introduced because with the primitive double-based dates, we sometimes very slightly missed the target date, i.e. propagator.propagate(date).getDate() did not return exactly the input date.
I agree that with the new long-based dates, this loop is probably not needed anymore.

I think there will still be small inaccuracies if the date is computed as date0.shiftedBy(dt) where dt is a primitive double instead of a TimeOffset, but at least the inaccuracy will occur before we propagate, hence the input date and output date should be the identical (i.e. exhibit the same inaccuracy).

Could you check if removing the loop still allows all unit tests to pass? If so, then you could open an issue to have the loop removed.

[Ryan]

Sure thing! I’m working on getting v13 set-up locally and I’ll let you know.

[Ryan]

O.k., so, in KeplerianPropagator using just durationFrom doesn’t cut it, but I’ve got it working by using accurateDurationFrom:

    public Orbit propagateOrbit(final AbsoluteDate date) {
        final Orbit      orbit = states.get(date).getOrbit();
        final TimeOffset dt    = date.accurateDurationFrom(orbit.getDate());
        return orbit.shiftedBy(dt);
    }

However, I’ve not had any luck with FieldKeplerianPropagator. It seems there isn’t a Field equivalent to accurateDurationFrom, and durationFrom results in the same problems as in the normal KeplerianPropagator when using the base durationFrom; i.e., some failed tests due to dates not matching internally, and an infinite loop in testVariableStep (see below for a screenshot).

Any suggestions for the Field side of things?

[luc]

I don’t see any obvious solution yet.
At least the primitive double version you propose should work.

[Ryan]

Would it problematic to keep the loop in the Field version and not in the original?

I have a few minutes to keep looking this afternoon (US, EST). I’d like to figure out why this isn’t an issue everywhere else.

E.g., why should we only have this loop here?

Moreover, why is it that, when the loop is removed, errors are thrown due to dates not being equal, such as in check calls within the larger propagation structure.

[luc]

@Ryan

Would it problematic to keep the loop in the Field version and not in the original?

No, I think it would be fair.

@Ryan

E.g., why should we only have this loop here?

We only add this kind of loop when we encounter a problem. I don’t remember exactly when it happened in this case, but for sure it happened and it was a problem, even if we missed the date by only one femtosecond : a test on the returned state would simply say “date is strictly earlier than target”, despite the error was really really small.

@Ryan

why is it that, when the loop is removed, errors are thrown due to dates not being equal

If users specify a date and there is no event to stop the propagation earlier, then this date should really be reached. doing differently would put the burden on the caller that could not use anymore logic with tests like “propagate until date t is reached” but would need more complex logic like “propagate until date t is almost reached with a tolerance of delta_t”.