Hello everyone,
I’ve started using Field propagation coupled with event detection with Orekit 12.0 and I’ve stumbled upon a couple of regressions. Basically, I’ve got a FieldNumericalPropagator with a custom ForceModel, whose intrinsic event detectors are instances of FieldDateDetector, let’s call the corresponding dates d_i.
I then propagate to several FieldAbsoludateDate equal to: d_i plus a Gradient variable with 0 constant term and non-zero derivatives. The original problem I came across was that with Orekit 12.0, the propagation would get stuck in a seemingly infinite loop at some point (not necessarily on the first date). This seems to happen in FieldDetectorBasedEventState from Hipparchus. I suspect it is due to one or several calls to isZero from CalculusFieldElement. However I couldn’t directly compare to version 2.3 as the class is new.
When I tried to extract a minimum reproduction example, I think I encountered another regression, but they might be linked.
The code is below, first for Orekit 12.0 and then for 11.3. The former leads to an MathRuntimeException while the latter runs fine.
@Test
void regressionDetection() {
final GradientField field = GradientField.getField(1);
final FieldAbsoluteDate<Gradient> fieldEpoch = FieldAbsoluteDate.getArbitraryEpoch(field);
final FieldPVCoordinates<Gradient> fieldPVCoordinates = new FieldPVCoordinates<>(
new FieldVector3D<>(field, new Vector3D(10000.e3, 0, 0)),
new FieldVector3D<>(field, new Vector3D(0, 7.5e3, 0)));
final FieldCartesianOrbit<Gradient> fieldOrbit = new FieldCartesianOrbit<>(fieldPVCoordinates,
FramesFactory.getGCRF(), fieldEpoch, field.getZero().add(Constants.EGM96_EARTH_MU));
final ClassicalRungeKuttaFieldIntegrator<Gradient> fieldIntegrator = new ClassicalRungeKuttaFieldIntegrator<>(field,
field.getZero().add(5.));
final FieldNumericalPropagator<Gradient> fieldNumericalPropagator = new FieldNumericalPropagator<>(
field, fieldIntegrator);
fieldNumericalPropagator.resetInitialState(new FieldSpacecraftState<>(fieldOrbit));
fieldNumericalPropagator.setResetAtEnd(true);
final AbsoluteDate epoch = fieldEpoch.toAbsoluteDate();
final TestForce testForce = new TestForce(epoch.shiftedBy(10.),
epoch.shiftedBy(100.), epoch.shiftedBy(1000.));
fieldNumericalPropagator.addForceModel(testForce);
for (AbsoluteDate date: testForce.dates) {
final FieldAbsoluteDate<Gradient> fieldDate = new FieldAbsoluteDate<>(field, date)
.shiftedBy(Gradient.variable(1, 0, 0.));
fieldNumericalPropagator.propagate(fieldDate);
}
}
private static class TestForce implements ForceModel {
private final AbsoluteDate[] dates;
TestForce (AbsoluteDate ...dates) {
this.dates = dates;
}
@Override
public boolean dependsOnPositionOnly() {
return true;
}
@Override
public Stream<EventDetector> getEventDetectors() {
return Arrays.stream(dates).map(DateDetector::new)
.map(d -> d.withHandler((a, b, c) -> Action.RESET_DERIVATIVES));
}
@Override
public <T extends CalculusFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventDetectors(Field<T> field) {
return Arrays.stream(dates).map(date -> new FieldDateDetector<>(field, new FieldAbsoluteDate<>(field, date)))
.map(d -> d.withHandler((a, b, c) -> Action.RESET_DERIVATIVES));
}
@Override
public Vector3D acceleration(SpacecraftState s, double[] parameters) {
return Vector3D.ZERO;
}
@Override
public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(FieldSpacecraftState<T> s, T[] parameters) {
return FieldVector3D.getZero(s.getDate().getField());
}
@Override
public List<ParameterDriver> getParametersDrivers() {
return new ArrayList<>();
}
}
and
@Test
void regressionDetection() {
final GradientField field = GradientField.getField(1);
final FieldAbsoluteDate<Gradient> fieldEpoch = FieldAbsoluteDate.getArbitraryEpoch(field);
final FieldPVCoordinates<Gradient> fieldPVCoordinates = new FieldPVCoordinates<>(
new FieldVector3D<>(field, new Vector3D(10000.e3, 0, 0)),
new FieldVector3D<>(field, new Vector3D(0, 7.5e3, 0)));
final FieldCartesianOrbit<Gradient> fieldOrbit = new FieldCartesianOrbit<>(fieldPVCoordinates,
FramesFactory.getGCRF(), fieldEpoch, field.getZero().add(Constants.EGM96_EARTH_MU));
final ClassicalRungeKuttaFieldIntegrator<Gradient> fieldIntegrator = new ClassicalRungeKuttaFieldIntegrator<>(field,
field.getZero().add(5.));
final FieldNumericalPropagator<Gradient> fieldNumericalPropagator = new FieldNumericalPropagator<>(
field, fieldIntegrator);
fieldNumericalPropagator.resetInitialState(new FieldSpacecraftState<>(fieldOrbit));
fieldNumericalPropagator.setResetAtEnd(true);
final AbsoluteDate epoch = fieldEpoch.toAbsoluteDate();
final TestForce testForce = new TestForce(epoch.shiftedBy(10.),
epoch.shiftedBy(100.), epoch.shiftedBy(1000.));
fieldNumericalPropagator.addForceModel(testForce);
for (AbsoluteDate date: testForce.dates) {
final FieldAbsoluteDate<Gradient> fieldDate = new FieldAbsoluteDate<>(field, date)
.shiftedBy(Gradient.variable(1, 0, 0.));
fieldNumericalPropagator.propagate(fieldDate);
}
}
private static class TestForce extends AbstractForceModel {
private final AbsoluteDate[] dates;
TestForce (AbsoluteDate ...dates) {
this.dates = dates;
}
@Override
public boolean dependsOnPositionOnly() {
return true;
}
@Override
public Stream<EventDetector> getEventsDetectors() {
return Arrays.stream(dates).map(DateDetector::new)
.map(d -> d.withHandler((a, b, c) -> Action.RESET_DERIVATIVES));
}
@Override
public <T extends CalculusFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventsDetectors(Field<T> field) {
return Arrays.stream(dates).map(date ->
new FieldDateDetector<>(field.getZero().add(AbstractDetector.DEFAULT_MAXCHECK),
field.getZero().add(AbstractDetector.DEFAULT_THRESHOLD), new FieldAbsoluteDate<>(field, date)))
.map(d -> d.withHandler((a, b, c) -> Action.RESET_DERIVATIVES));
}
@Override
public Vector3D acceleration(SpacecraftState s, double[] parameters) {
return Vector3D.ZERO;
}
@Override
public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(FieldSpacecraftState<T> s, T[] parameters) {
return FieldVector3D.getZero(s.getDate().getField());
}
@Override
public List<ParameterDriver> getParametersDrivers() {
return new ArrayList<>();
}
}
Cheers,
Romain.