We are in the process of upgrading an application from Orekit 12 to 13 and noticed a behavioral change involving extrapolation. It seems that forward and backward extrapolation are no longer symmetric: backward extrapolation behaves as before, while forward extrapolation shows differences on the order of ~6 meters.
It is not clear to us whether this indicates:
- an issue or regression,
- an expected behavioral change in Orekit 13, or
- a problem with our test assumptions or configuration.
We would appreciate any guidance. We have consolidated the relevant functionality into the unit test below. The test passes with Orekit 12 but fails with Orekit 13.
import static org.junit.jupiter.api.Assertions.assertEquals;
import static org.junit.jupiter.api.Assertions.assertThrows;
import java.io.StringReader;
import java.util.List;
import java.util.Map;
import org.junit.jupiter.api.Test;
import org.orekit.attitudes.FrameAlignedProvider;
import org.orekit.data.DataContext;
import org.orekit.data.DataSource;
import org.orekit.files.stk.STKEphemerisFile;
import org.orekit.files.stk.STKEphemerisFileParser;
import org.orekit.files.stk.STKEphemerisFile.STKCoordinateSystem;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.orbits.Orbit;
import org.orekit.propagation.BoundedPropagator;
import org.orekit.propagation.Propagator;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.analytical.AbstractAnalyticalPropagator;
import org.orekit.propagation.analytical.AggregateBoundedPropagator;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScalesFactory;
import org.orekit.time.UTCScale;
class ExtrapolationTest {
@Test
void testExtrapolation() {
final double maximumExtrapolation = 60.0;
// create propagator (mirrors original method)
final BoundedPropagator propagator = createStkEphemerisFilePropagator(maximumExtrapolation);
final SpacecraftState initialState = propagator.getInitialState();
final UTCScale utc = TimeScalesFactory.getUTC();
// propagate to nominal final state
final SpacecraftState finalState = propagator.propagate(
new AbsoluteDate(2007, 1, 12, 0, 5, 0., utc));
// left extrapolation (one minute before start)
final SpacecraftState leftExtrapolatedState = propagator.propagate(
new AbsoluteDate(2007, 1, 11, 23, 59, 0., utc));
final SpacecraftState expectedLeft = initialState.shiftedBy(-maximumExtrapolation);
assertEquals(expectedLeft.getPVCoordinates().getPosition(),
leftExtrapolatedState.getPVCoordinates().getPosition());
assertEquals(expectedLeft.getPVCoordinates().getVelocity(),
leftExtrapolatedState.getPVCoordinates().getVelocity());
// right extrapolation (one minute after final nominal state)
final SpacecraftState rightExtrapolatedState = propagator.propagate(
new AbsoluteDate(2007, 1, 12, 0, 6, 0., utc));
final SpacecraftState expectedRight = finalState.shiftedBy(+maximumExtrapolation);
assertEquals(expectedRight.getPVCoordinates().getPosition().getX(),
rightExtrapolatedState.getPVCoordinates().getPosition().getX(), 1e-9);
assertEquals(expectedRight.getPVCoordinates().getPosition().getY(),
rightExtrapolatedState.getPVCoordinates().getPosition().getY(), 1e-9);
assertEquals(expectedRight.getPVCoordinates().getPosition().getZ(),
rightExtrapolatedState.getPVCoordinates().getPosition().getZ(), 1e-9);
assertEquals(expectedRight.getPVCoordinates().getVelocity().getX(),
rightExtrapolatedState.getPVCoordinates().getVelocity().getX(), 1e-12);
assertEquals(expectedRight.getPVCoordinates().getVelocity().getY(),
rightExtrapolatedState.getPVCoordinates().getVelocity().getY(), 2e-12);
assertEquals(expectedRight.getPVCoordinates().getVelocity().getZ(),
rightExtrapolatedState.getPVCoordinates().getVelocity().getZ(), 1e-12);
// outside extrapolation bounds should throw
assertThrows(IllegalArgumentException.class,
() -> propagator.propagate(
new AbsoluteDate(2007, 1, 11, 23, 58, 59.999999999, utc)));
assertThrows(IllegalArgumentException.class,
() -> propagator.propagate(
new AbsoluteDate(2007, 1, 12, 0, 6, 0.000000001, utc)));
}
/** Encapsulated helper method for building the aggregate propagator with optional extrapolation */
private static BoundedPropagator createStkEphemerisFilePropagator(final double maximumExtrapolation) {
final UTCScale utc = TimeScalesFactory.getUTC();
final double mu = 3.986004415e14;
// parse ephemeris from in-memory string
final DataSource source = new DataSource("stk.ephemeris", () -> new StringReader(stkEphemeris()));
final Map<STKCoordinateSystem, Frame> frameMapping = Map.of(
STKCoordinateSystem.ICRF, FramesFactory.getGCRF(),
STKCoordinateSystem.J2000, FramesFactory.getEME2000());
final STKEphemerisFileParser parser = new STKEphemerisFileParser("99999", mu, utc, frameMapping);
final STKEphemerisFile file = parser.parse(source);
final BoundedPropagator basePropagator =
file.getSatellites().values().iterator().next().getPropagator();
// add optional left/right extrapolation
final SpacecraftState initialState = basePropagator.getInitialState();
final SpacecraftState finalState = basePropagator.propagate(basePropagator.getMaxDate());
final BoundedPropagator left =
new ShiftingBoundedAnalyticalPropagator(initialState, -maximumExtrapolation, utc);
final BoundedPropagator right =
new ShiftingBoundedAnalyticalPropagator(finalState, +maximumExtrapolation, utc);
return new AggregateBoundedPropagator(List.of(left, basePropagator, right));
}
/** Simple shifting propagator for testing extrapolation */
private static final class ShiftingBoundedAnalyticalPropagator
extends AbstractAnalyticalPropagator implements BoundedPropagator {
private final SpacecraftState reference;
private final AbsoluteDate minDate;
private final AbsoluteDate maxDate;
private final UTCScale utc;
ShiftingBoundedAnalyticalPropagator(final SpacecraftState reference,
final double maxShift,
final UTCScale utc) {
super(FrameAlignedProvider.of(reference.getFrame()));
this.reference = reference;
this.utc = utc;
resetInitialState(reference);
if (maxShift > 0) {
minDate = reference.getDate();
maxDate = minDate.shiftedBy(maxShift);
} else {
maxDate = reference.getDate();
minDate = maxDate.shiftedBy(maxShift);
}
}
@Override
public Orbit propagateOrbit(final AbsoluteDate date) {
if (date.isBefore(minDate) || date.isAfter(maxDate)) {
throw new IllegalArgumentException(
String.format("Date must be in [%s, %s], was %s",
minDate.toStringRfc3339(utc),
maxDate.toStringRfc3339(utc),
date.toStringRfc3339(utc)));
}
return reference.getOrbit().shiftedBy(date.durationFrom(reference.getDate()));
}
@Override
protected double getMass(final AbsoluteDate date) {
return reference.getMass();
}
@Override
protected void resetIntermediateState(final SpacecraftState state, final boolean forward) {
throw new UnsupportedOperationException();
}
@Override
public AbsoluteDate getMinDate() { return minDate; }
@Override
public AbsoluteDate getMaxDate() { return maxDate; }
}
/** Minimal STK ephemeris for test */
private static String stkEphemeris() {
return "stk.v.12.0\n" +
"BEGIN Ephemeris\n" +
" NumberOfEphemerisPoints 6\n" +
" ScenarioEpoch 12 Jan 2007 00:00:00.000\n" +
" InterpolationMethod Lagrange\n" +
" InterpolationSamplesM1 5\n" +
" CentralBody Earth\n" +
" CoordinateSystem J2000\n" +
" EphemerisTimePosVel\n" +
" 0.0 -7.3217622956859889e+06 -6.5112200072246802e+05 -2.5954523517972045e+06 -1.3783361548820392e+03 -4.8106444659146173e+03 5.1079817147894801e+03\n" +
" 60.0 -7.3933276740881586e+06 -9.3862795014265948e+05 -2.2851879696136559e+06 -1.0065654101328107e+03 -4.7704797808443709e+03 5.2315577564422465e+03\n" +
" 120.0 -7.4424870313313799e+06 -1.2232893062177158e+06 -1.9679831047277695e+06 -6.3165697134411255e+02 -4.7158467735648492e+03 5.3392717215232315e+03\n" +
" 180.0 -7.4690862155061280e+06 -1.5042425362945164e+06 -1.6447997891084312e+06 -2.5474943725825887e+02 -4.6469049447876278e+03 5.4307880975566077e+03\n" +
" 240.0 -7.4730396315300921e+06 -1.7806349998268955e+06 -1.3166187108281432e+06 1.2301114478231821e+02 -4.5638578109296532e+03 5.5058207040365623e+03\n" +
" 300.0 -7.4543305805486357e+06 -2.0516275743448907e+06 -9.8443622230934782e+05 5.0047487187925583e+02 -4.4669523821145613e+03 5.5641337301821031e+03\n" +
"END Ephemeris\n";
}
}