Skip to content
GitLab
Commit edefe97c authored by Luc Maisonobe's avatar Luc Maisonobe
Browse files

Started tests.

parent 235f59e0
Hide whitespace changes
Inline Side-by-side
src/main/java/org/orekit/forces/maneuvers/trigger/DateBasedManeuverTriggers.java
View file @ edefe97c
......@@ -116,6 +116,8 @@ public class DateBasedManeuverTriggers extends IntervalEventTrigger<ParameterDri
detector.getStopDriver().getBaseDate());
fd.getStartDriver().setName(detector.getStartDriver().getName());
fd.getStopDriver().setName(detector.getStopDriver().getName());
fd.getMedianDriver().setName(detector.getMedianDriver().getName());
fd.getDurationDriver().setName(detector.getDurationDriver().getName());
@SuppressWarnings("unchecked")
final FieldAbstractDetector<D, S> converted = (FieldAbstractDetector<D, S>) fd;
......@@ -127,7 +129,9 @@ public class DateBasedManeuverTriggers extends IntervalEventTrigger<ParameterDri
@Override
public List<ParameterDriver> getParametersDrivers() {
return Arrays.asList(getFiringIntervalDetector().getStartDriver(),
getFiringIntervalDetector().getStopDriver());
getFiringIntervalDetector().getStopDriver(),
getFiringIntervalDetector().getMedianDriver(),
getFiringIntervalDetector().getDurationDriver());
}
}
src/main/java/org/orekit/propagation/events/FieldParameterDrivenDateIntervalDetector.java
View file @ edefe97c
......@@ -16,16 +16,32 @@
*/
package org.orekit.propagation.events;
import java.util.List;
import java.util.stream.Collectors;
import org.hipparchus.CalculusFieldElement;
import org.hipparchus.Field;
import org.hipparchus.util.FastMath;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;
import org.orekit.propagation.FieldSpacecraftState;
import org.orekit.propagation.events.handlers.FieldEventHandler;
import org.orekit.propagation.events.handlers.FieldStopOnEvent;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.DateDriver;
import org.orekit.utils.ParameterDriver;
import org.orekit.utils.ParameterObserver;
/** Detector for date intervals that may be offset thanks to parameter drivers.
* <p>
* Two dual views can be used for date intervals: either start date/stop date or
* median date/duration. {@link #getStartDriver() start}/{@link #getStopDriver() stop}
* drivers and {@link #getMedianDriver() median}/{@link #getDurationDriver() duration}
* drivers work in pair. Both drivers in one pair can be selected and their changes will
* be propagated to the other pair, but attempting to select drivers in both
* pairs at the same time will trigger an exception. Changing the value of a driver
* that is not selected should be avoided as it leads to inconsistencies between the pairs.
* </p>
* @see org.orekit.propagation.Propagator#addEventDetector(EventDetector)
* @param <T> type of the field elements
* @author Luc Maisonobe
......@@ -39,12 +55,37 @@ public class FieldParameterDrivenDateIntervalDetector<T extends CalculusFieldEle
/** Default suffix for stop driver. */
public static final String STOP_SUFFIX = "_STOP";
/** Default suffix for median driver. */
public static final String MEDIAN_SUFFIX = "_MEDIAN";
/** Default suffix for duration driver. */
public static final String DURATION_SUFFIX = "_DURATION";
/** Reference interval start driver. */
private DateDriver start;
/** Reference interval stop driver. */
private DateDriver stop;
/** Median date driver. */
private DateDriver median;
/** Duration driver. */
private ParameterDriver duration;
/** Build a new instance.
* @param field field to which the elements belong
* @param prefix prefix to use for parameter drivers names
* @param refMedian reference interval median date
* @param refDuration reference duration
*/
public FieldParameterDrivenDateIntervalDetector(final Field<T> field, final String prefix,
final AbsoluteDate refMedian, final double refDuration) {
this(field, prefix,
refMedian.shiftedBy(-0.5 * refDuration),
refMedian.shiftedBy(+0.5 * refDuration));
}
/** Build a new instance.
* @param field field to which the elements belong
* @param prefix prefix to use for parameter drivers names
......@@ -58,7 +99,9 @@ public class FieldParameterDrivenDateIntervalDetector<T extends CalculusFieldEle
DEFAULT_MAX_ITER,
new FieldStopOnEvent<FieldParameterDrivenDateIntervalDetector<T>, T>(),
new DateDriver(refStart, prefix + START_SUFFIX, true),
new DateDriver(refStop, prefix + STOP_SUFFIX, false));
new DateDriver(refStop, prefix + STOP_SUFFIX, false),
new DateDriver(refStart.shiftedBy(0.5 * refStop.durationFrom(refStart)), prefix + MEDIAN_SUFFIX, true),
new ParameterDriver(prefix + DURATION_SUFFIX, refStop.durationFrom(refStart), 1.0, 0.0, Double.POSITIVE_INFINITY));
}
/** Private constructor with full parameters.
......@@ -73,13 +116,43 @@ public class FieldParameterDrivenDateIntervalDetector<T extends CalculusFieldEle
* @param handler event handler to call at event occurrences
* @param start reference interval start driver
* @param stop reference interval stop driver
* @param median median date driver
* @param duration duration driver
*/
private FieldParameterDrivenDateIntervalDetector(final T maxCheck, final T threshold, final int maxIter,
final FieldEventHandler<? super FieldParameterDrivenDateIntervalDetector<T>, T> handler,
final DateDriver start, final DateDriver stop) {
final DateDriver start, final DateDriver stop,
final DateDriver median, final ParameterDriver duration) {
super(maxCheck, threshold, maxIter, handler);
this.start = start;
this.stop = stop;
this.start = start;
this.stop = stop;
this.median = median;
this.duration = duration;
// set up delegation between drivers
replaceBindingObserver(start, new StartObserver());
replaceBindingObserver(stop, new StopObserver());
replaceBindingObserver(median, new MedianObserver());
replaceBindingObserver(duration, new DurationObserver());
}
/** Replace binding observers.
* @param driver driver for whose binding observers should be replaced
* @param bindingObserver new binding observer
*/
private void replaceBindingObserver(final ParameterDriver driver, final BindingObserver bindingObserver) {
// remove the previous binding observers
final List<ParameterObserver> original = driver.
getObservers().
stream().
filter(observer -> observer instanceof FieldParameterDrivenDateIntervalDetector.BindingObserver).
collect(Collectors.toList());
original.forEach(observer -> driver.removeObserver(observer));
driver.addObserver(bindingObserver);
}
/** {@inheritDoc} */
......@@ -87,10 +160,15 @@ public class FieldParameterDrivenDateIntervalDetector<T extends CalculusFieldEle
protected FieldParameterDrivenDateIntervalDetector<T> create(final T newMaxCheck, final T newThreshold, final int newMaxIter,
final FieldEventHandler<? super FieldParameterDrivenDateIntervalDetector<T>, T> newHandler) {
return new FieldParameterDrivenDateIntervalDetector<>(newMaxCheck, newThreshold, newMaxIter, newHandler,
start, stop);
start, stop, median, duration);
}
/** Get the driver for start date.
* <p>
* Note that the start date is automatically adjusted if either
* {@link #getMedianDriver() median date} or {@link #getDurationDriver() duration}
* are {@link ParameterDriver#isSelected() selected} and changed.
* </p>
* @return driver for start date
*/
public DateDriver getStartDriver() {
......@@ -98,12 +176,41 @@ public class FieldParameterDrivenDateIntervalDetector<T extends CalculusFieldEle
}
/** Get the driver for stop date.
* <p>
* Note that the stop date is automatically adjusted if either
* {@link #getMedianDriver() median date} or {@link #getDurationDriver() duration}
* are {@link ParameterDriver#isSelected() selected} changed.
* </p>
* @return driver for stop date
*/
public DateDriver getStopDriver() {
return stop;
}
/** Get the driver for median date.
* <p>
* Note that the median date is automatically adjusted if either
* {@link #getStartDriver()} start date or {@link #getStopDriver() stop date}
* are {@link ParameterDriver#isSelected() selected} changed.
* </p>
* @return driver for median date
*/
public DateDriver getMedianDriver() {
return median;
}
/** Get the driver for duration.
* <p>
* Note that the duration is automatically adjusted if either
* {@link #getStartDriver()} start date or {@link #getStopDriver() stop date}
* are {@link ParameterDriver#isSelected() selected} changed.
* </p>
* @return driver for duration
*/
public ParameterDriver getDurationDriver() {
return duration;
}
/** Compute the value of the switching function.
* <p>
* The function is positive for dates within the interval defined
......@@ -121,4 +228,73 @@ public class FieldParameterDrivenDateIntervalDetector<T extends CalculusFieldEle
s.getDate().durationFrom(stop.getDate()).negate());
}
/** Base observer. */
private abstract class BindingObserver implements ParameterObserver {
/** {@inheritDoc} */
@Override
public void valueChanged(final double previousValue, final ParameterDriver driver) {
if (driver.isSelected()) {
setDelta(driver.getValue() - previousValue);
}
}
/** {@inheritDoc} */
@Override
public void selectionChanged(final boolean previousSelection, final ParameterDriver driver) {
if ((start.isSelected() || stop.isSelected()) &&
(median.isSelected() || duration.isSelected())) {
throw new OrekitException(OrekitMessages.INCONSISTENT_SELECTION,
start.getName(), stop.getName(),
median.getName(), duration.getName());
}
}
/** Change a value.
* @param delta change of value
*/
protected abstract void setDelta(double delta);
}
/** Observer for start date. */
private class StartObserver extends BindingObserver {
/** {@inheritDoc} */
@Override
protected void setDelta(final double delta) {
median.setValue(median.getValue() + 0.5 * delta);
duration.setValue(duration.getValue() - delta);
}
}
/** Observer for stop date. */
private class StopObserver extends BindingObserver {
/** {@inheritDoc} */
@Override
protected void setDelta(final double delta) {
median.setValue(median.getValue() + 0.5 * delta);
duration.setValue(duration.getValue() + delta);
}
}
/** Observer for median date. */
private class MedianObserver extends BindingObserver {
/** {@inheritDoc} */
@Override
protected void setDelta(final double delta) {
start.setValue(start.getValue() + delta);
stop.setValue(stop.getValue() + delta);
}
}
/** Observer for duration. */
private class DurationObserver extends BindingObserver {
/** {@inheritDoc} */
@Override
protected void setDelta(final double delta) {
start.setValue(start.getValue() - 0.5 * delta);
stop.setValue(stop.getValue() + 0.5 * delta);
}
}
}
src/main/java/org/orekit/propagation/events/ParameterDrivenDateIntervalDetector.java
View file @ edefe97c
......@@ -140,7 +140,7 @@ public class ParameterDrivenDateIntervalDetector extends AbstractDetector<Parame
final List<ParameterObserver> original = driver.
getObservers().
stream().
filter(observer -> observer instanceof BindingObserver).
filter(observer -> observer instanceof ParameterDrivenDateIntervalDetector.BindingObserver).
collect(Collectors.toList());
original.forEach(observer -> driver.removeObserver(observer));
......
src/test/java/org/orekit/propagation/numerical/TriggersDerivativesTest.java
View file @ edefe97c
......@@ -17,9 +17,14 @@
package org.orekit.propagation.numerical;
import java.io.File;
import java.io.IOException;
import java.io.PrintStream;
import java.nio.charset.StandardCharsets;
import java.util.ArrayList;
import java.util.Collections;
import java.util.List;
import java.util.Locale;
import org.hipparchus.analysis.differentiation.UnivariateDerivative1;
import org.hipparchus.geometry.euclidean.threed.Rotation;
......@@ -73,18 +78,34 @@ public class TriggersDerivativesTest {
final double duration = 200.0;
// first propagation, covering the maneuver
DateBasedManeuverTriggers triggers1 = new DateBasedManeuverTriggers("MAN_0", firing, duration);
final NumericalPropagator propagator1 = buildPropagator(OrbitType.EQUINOCTIAL, PositionAngle.TRUE, 20,
firing, duration, true, 0);
setUpManeuverJacobianComputation(true, propagator1);
firing, duration, triggers1);
propagator1.
getAllForceModels().
forEach(fm -> fm.
getParametersDrivers().
stream().
filter(d -> d.getName().equals("MAN_0_START") ||
d.getName().equals(NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT)).
forEach(d -> d.setSelected(true)));
final MatricesHarvester harvester1 = propagator1.setupMatricesComputation("stm", null, null);
final SpacecraftState state1 = propagator1.propagate(firing.shiftedBy(2 * duration));
final RealMatrix stm1 = harvester1.getStateTransitionMatrix(state1);
final RealMatrix jacobian1 = harvester1.getParametersJacobian(state1);
// second propagation, interrupted during maneuver
final NumericalPropagator propagator2 = buildPropagator(OrbitType.EQUINOCTIAL, PositionAngle.TRUE, 20,
firing, duration, true, 0);
setUpManeuverJacobianComputation(true, propagator2);
DateBasedManeuverTriggers triggers2 = new DateBasedManeuverTriggers("MAN_0", firing, duration);
final NumericalPropagator propagator2 = buildPropagator(OrbitType.EQUINOCTIAL, PositionAngle.TRUE, 20,
firing, duration, triggers2);
propagator2.
getAllForceModels().
forEach(fm -> fm.
getParametersDrivers().
stream().
filter(d -> d.getName().equals("MAN_0_START") ||
d.getName().equals(NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT)).
forEach(d -> d.setSelected(true)));
// some additional providers for test coverage
final StateTransitionMatrixGenerator dummyStmGenerator =
......@@ -128,29 +149,89 @@ public class TriggersDerivativesTest {
@Test
public void testDerivativeWrtStartTimeCartesian() {
doTestDerivativeWrtTime(true, OrbitType.CARTESIAN,
doTestDerivativeWrtStartStopTime(true, OrbitType.CARTESIAN,
0.022, 0.012, 0.012, 0.013, 0.012, 0.021);
}
@Test
public void testDerivativeWrtStartTimeKeplerian() {
doTestDerivativeWrtTime(true, OrbitType.KEPLERIAN,
doTestDerivativeWrtStartStopTime(true, OrbitType.KEPLERIAN,
0.012, 0.011, 0.011, 0.012, 0.011, 0.017);
}
@Test
public void testDerivativeWrtStopTimeCartesian() {
doTestDerivativeWrtTime(false, OrbitType.CARTESIAN,
doTestDerivativeWrtStartStopTime(false, OrbitType.CARTESIAN,
0.00033, 0.00045, 0.00040, 0.00022, 0.00020, 0.00010);
}
@Test
public void testDerivativeWrtStopTimeKeplerian() {
doTestDerivativeWrtTime(false, OrbitType.KEPLERIAN,
doTestDerivativeWrtStartStopTime(false, OrbitType.KEPLERIAN,
0.0011, 0.00020, 0.00002, 0.00082, 0.00008, 0.00019);
}
private void doTestDerivativeWrtTime(final boolean start, final OrbitType orbitType, final double...tolerance) {
@Test
public void testDerivativeWrtMedianKeplerian() {
doTestDerivativeWrtMedianTime(OrbitType.KEPLERIAN,
0.0011, 0.00020, 0.00002, 0.00082, 0.00008, 0.00019);
}
private void doTestDerivativeWrtStartStopTime(final boolean start, final OrbitType orbitType, final double...tolerance) {
final AbsoluteDate firing = new AbsoluteDate(new DateComponents(2004, 1, 2),
new TimeComponents(4, 15, 34.080),
TimeScalesFactory.getUTC());
final List<Propagator> propagators = new ArrayList<>();
// propagators will be combined using finite differences to compute derivatives
final PositionAngle positionAngle = PositionAngle.TRUE;
final int degree = 20;
final double duration = 200.0;
final double h = 1.0;
final double samplingtep = 240.0;
for (int k = -4; k <= 4; ++k) {
final DateBasedManeuverTriggers trigger = start ?
new DateBasedManeuverTriggers("MAN_0", firing.shiftedBy(k * h), duration - k * h) :
new DateBasedManeuverTriggers("MAN_0", firing, duration + k * h);
propagators.add(buildPropagator(orbitType, positionAngle, degree, firing, duration, trigger));
}
// the central propagator (k = 4) will compute derivatives autonomously using State and TriggersDerivatives
final NumericalPropagator autonomous = (NumericalPropagator) propagators.get(4);
final MatricesHarvester harvester = autonomous.setupMatricesComputation("stm", null, null);
autonomous.
getAllForceModels().
forEach(fm -> fm.
getParametersDrivers().
stream().
filter(d -> d.getName().equals(start ? "MAN_0_START" : "MAN_0_STOP") ||
d.getName().equals(NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT)).
forEach(d -> d.setSelected(true)));
DerivativesSampler sampler = new DerivativesSampler(harvester, orbitType, positionAngle,
firing, duration, h, samplingtep);
new PropagatorsParallelizer(propagators, sampler).
propagate(firing.shiftedBy(-30 * samplingtep), firing.shiftedBy(duration + 300 * samplingtep));
double[] maxRelativeError = new double[tolerance.length];
for (final Entry entry : sampler.sample) {
for (int i = 0; i < entry.finiteDifferences.length; ++i) {
double f = entry.finiteDifferences[i].getFirstDerivative();
double c = entry.closedForm[i].getFirstDerivative();
maxRelativeError[i] = FastMath.max(maxRelativeError[i], FastMath.abs(f - c) / FastMath.max(1.0e-10, FastMath.abs(f)));
}
}
for (int i = 0; i < tolerance.length; ++i) {
Assert.assertEquals(0.0, maxRelativeError[i], tolerance[i]);
}
}
private void doTestDerivativeWrtMedianTime(final OrbitType orbitType, final double...tolerance) {
final AbsoluteDate firing = new AbsoluteDate(new DateComponents(2004, 1, 2),
new TimeComponents(4, 15, 34.080),
......@@ -165,13 +246,22 @@ public class TriggersDerivativesTest {
final double h = 1.0;
final double samplingtep = 240.0;
for (int k = -4; k <= 4; ++k) {
propagators.add(buildPropagator(orbitType, positionAngle, degree, firing, duration, start, k * h));
final DateBasedManeuverTriggers triggers =
new DateBasedManeuverTriggers("MAN_0", firing.shiftedBy(k * h), duration);
propagators.add(buildPropagator(orbitType, positionAngle, degree, firing, duration, triggers));
}
// the central propagator (k = 4) will compute derivatives autonomously using State and TriggersDerivatives
final NumericalPropagator autonomous = (NumericalPropagator) propagators.get(4);
final MatricesHarvester harvester = autonomous.setupMatricesComputation("stm", null, null);
setUpManeuverJacobianComputation(start, autonomous);
autonomous.
getAllForceModels().
forEach(fm -> fm.
getParametersDrivers().
stream().
filter(d -> d.getName().equals("MAN_0_MEDIAN") ||
d.getName().equals(NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT)).
forEach(d -> d.setSelected(true)));
DerivativesSampler sampler = new DerivativesSampler(harvester, orbitType, positionAngle,
firing, duration, h, samplingtep);
......@@ -188,36 +278,91 @@ public class TriggersDerivativesTest {
}
}
analyzeSample(sampler, orbitType, firing, degree, true, null);
for (int i = 0; i < tolerance.length; ++i) {
Assert.assertEquals(0.0, maxRelativeError[i], tolerance[i]);
}
}
private void setUpManeuverJacobianComputation(final boolean start, final NumericalPropagator propagator) {
propagator.
getAllForceModels().
forEach(fm -> fm.
getParametersDrivers().
stream().
filter(d -> d.getName().equals(start ? "MAN_0_START" : "MAN_0_STOP") ||
d.getName().equals(NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT)).
forEach(d -> d.setSelected(true)));
private void analyzeSample(final DerivativesSampler sampler, final OrbitType orbitType, final AbsoluteDate firing,
final int degree, final boolean plot, final String outputDir) {
if (!plot) {
return;
}
final ProcessBuilder pb = new ProcessBuilder("gnuplot").
redirectOutput(ProcessBuilder.Redirect.INHERIT).
redirectError(ProcessBuilder.Redirect.INHERIT);
pb.environment().remove("XDG_SESSION_TYPE");
Process gnuplot;
try {
gnuplot = pb.start();
try (PrintStream out = new PrintStream(gnuplot.getOutputStream(), false, StandardCharsets.UTF_8.name())) {
final String fileName;
if (outputDir == null) {
fileName = null;
out.format(Locale.US, "set terminal qt size %d, %d title 'complex plotter'%n", 1000, 1000);
} else {
fileName = new File(outputDir,
"triggers-partials-" + orbitType + "-degree-" + degree + ".png").
getAbsolutePath();
out.format(Locale.US, "set terminal pngcairo size %d, %d%n", 1000, 1000);
out.format(Locale.US, "set output '%s'%n", fileName);
}
out.format(Locale.US, "set offset graph 0.05, 0.05, 0.05, 0.05%n");
out.format(Locale.US, "set view map scale 1%n");
out.format(Locale.US, "set xlabel 't - t_{start}'%n");
if (orbitType == OrbitType.CARTESIAN) {
out.format(Locale.US, "set ylabel 'd\\{X,Y,Z\\}/dt_{start} (m)'%n");
out.format(Locale.US, "set key bottom left%n");
} else {
out.format(Locale.US, "set ylabel 'da/dt_{start} (m)'%n");
out.format(Locale.US, "set key top left%n");
}
out.format(Locale.US, "set title 'derivatives of %s state, gravity field degree %d'%n", orbitType, degree);
out.format(Locale.US, "$data <<EOD%n");
for (final Entry entry : sampler.sample) {
out.format(Locale.US, "%.6f", entry.date.durationFrom(firing));
for (int i = 0; i < entry.finiteDifferences.length; ++i) {
out.format(Locale.US, " %.9f %.9f",
entry.finiteDifferences[i].getFirstDerivative(),
entry.closedForm[i].getFirstDerivative());
}
out.format(Locale.US, "%n");
}
out.format(Locale.US, "EOD%n");
if (orbitType == OrbitType.CARTESIAN) {
out.print("plot $data using 1:2 with lines title 'dX/dt_{start} finite differences', \\\n" +
" '' using 1:3 with points title 'dX/dt_{start} closed form based on STM', \\\n" +
" '' using 1:4 with lines title 'dY/dt_{start} finite differences', \\\n" +
" '' using 1:5 with points title 'dY/dt_{start} closed form based on STM', \\\n" +
" '' using 1:6 with lines title 'dZ/dt_{start} finite differences', \\\n" +
" '' using 1:7 with points title 'dZ/dt_{start} closed form based on STM'\n");
} else {
out.print("plot $data using 1:2 with lines title 'da/dt_{start} finite differences', \\\n" +
" '' using 1:3 with points title 'da/dt_{start} closed form based on STM'\n");
}
if (fileName == null) {
out.format(Locale.US, "pause mouse close%n");
} else {
System.out.format(Locale.US, "output written to %s%n", fileName);
}
}
} catch (IOException ioe) {
Assert.fail(ioe.getLocalizedMessage());
}
}
private NumericalPropagator buildPropagator(final OrbitType orbitType, final PositionAngle positionAngle,
final int degree, final AbsoluteDate firing, final double duration,
final boolean start, final double shift) {
final DateBasedManeuverTriggers triggers) {
final double delta = FastMath.toRadians(-7.4978);
final double alpha = FastMath.toRadians(351);
AttitudeProvider attitudeProvider = new InertialProvider(new Rotation(new Vector3D(alpha, delta), Vector3D.PLUS_I));
final DateBasedManeuverTriggers triggers =
start ?
new DateBasedManeuverTriggers("MAN_0", firing