Commit d64ebe44 authored by Luc Maisonobe's avatar Luc Maisonobe
Browse files

Moved around classes for maneuvers jacobians.

parent 92ce876f
......@@ -14,7 +14,10 @@
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.propagation;
package org.orekit.forces.maneuvers.jacobians;
import org.orekit.propagation.AdditionalStateProvider;
import org.orekit.propagation.SpacecraftState;
/** Generator for one column of a Jacobian matrix for special case of maneuver duration.
* <p>
......@@ -23,10 +26,10 @@ package org.orekit.propagation;
* </p>
* @author Luc Maisonobe
* @since 11.1
* @see MedianDateJacobianColumnGenerator
* @see TriggerDateJacobianColumnGenerator
* @see MedianDate
* @see TriggerDate
*/
public class DurationJacobianColumnGenerator implements AdditionalStateProvider {
public class Duration implements AdditionalStateProvider {
/** Name of the parameter corresponding to the start date. */
private final String startName;
......@@ -42,7 +45,7 @@ public class DurationJacobianColumnGenerator implements AdditionalStateProvider
* @param stopName name of the parameter corresponding to the stop date
* @param columnName name of the parameter corresponding to the column
*/
public DurationJacobianColumnGenerator(final String startName, final String stopName, final String columnName) {
public Duration(final String startName, final String stopName, final String columnName) {
this.startName = startName;
this.stopName = stopName;
this.columnName = columnName;
......
......@@ -14,7 +14,10 @@
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.propagation;
package org.orekit.forces.maneuvers.jacobians;
import org.orekit.propagation.AdditionalStateProvider;
import org.orekit.propagation.SpacecraftState;
/** Generator for one column of a Jacobian matrix for special case of maneuver median date.
* <p>
......@@ -23,10 +26,10 @@ package org.orekit.propagation;
* </p>
* @author Luc Maisonobe
* @since 11.1
* @see DurationJacobianColumnGenerator
* @see TriggerDateJacobianColumnGenerator
* @see Duration
* @see TriggerDate
*/
public class MedianDateJacobianColumnGenerator implements AdditionalStateProvider {
public class MedianDate implements AdditionalStateProvider {
/** Name of the parameter corresponding to the start date. */
private final String startName;
......@@ -42,7 +45,7 @@ public class MedianDateJacobianColumnGenerator implements AdditionalStateProvide
* @param stopName name of the parameter corresponding to the stop date
* @param columnName name of the parameter corresponding to the column
*/
public MedianDateJacobianColumnGenerator(final String startName, final String stopName, final String columnName) {
public MedianDate(final String startName, final String stopName, final String columnName) {
this.startName = startName;
this.stopName = stopName;
this.columnName = columnName;
......
......@@ -14,7 +14,7 @@
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.propagation;
package org.orekit.forces.maneuvers.jacobians;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.linear.MatrixUtils;
......@@ -23,6 +23,8 @@ import org.hipparchus.linear.RealMatrix;
import org.hipparchus.linear.RealVector;
import org.orekit.forces.maneuvers.Maneuver;
import org.orekit.forces.maneuvers.trigger.ManeuverTriggersResetter;
import org.orekit.propagation.AdditionalStateProvider;
import org.orekit.propagation.SpacecraftState;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.TimeSpanMap;
......@@ -78,10 +80,10 @@ import org.orekit.utils.TimeSpanMap;
* </p>
* @author Luc Maisonobe
* @since 11.1
* @see MedianDateJacobianColumnGenerator
* @see DurationJacobianColumnGenerator
* @see MedianDate
* @see Duration
*/
public class TriggerDateJacobianColumnGenerator
public class TriggerDate
implements AdditionalStateProvider, ManeuverTriggersResetter {
/** Dimension of the state. */
......@@ -97,10 +99,13 @@ public class TriggerDateJacobianColumnGenerator
private final String stmName;
/** Name of the parameter corresponding to the column. */
private final String columnName;
private final String triggerName;
/** Mass depletion effect. */
private final MassDepletionDelay massDepletionDelay;
/** Start/stop management flag. */
private boolean manageStart;
private final boolean manageStart;
/** Maneuver force model. */
private final Maneuver maneuver;
......@@ -119,28 +124,28 @@ public class TriggerDateJacobianColumnGenerator
/** Simple constructor.
* @param stmName name of State Transition Matrix state
* @param columnName name of the parameter corresponding to the column
* @param triggerName name of the parameter corresponding to the trigger date column
* @param manageStart if true, we compute derivatives with respect to maneuver start
* @param maneuver maneuver force model
* @param threshold event detector threshold
*/
public TriggerDateJacobianColumnGenerator(final String stmName, final String columnName,
final boolean manageStart, final Maneuver maneuver,
final double threshold) {
this.stmName = stmName;
this.columnName = columnName;
this.manageStart = manageStart;
this.maneuver = maneuver;
this.threshold = threshold;
this.contribution = null;
this.trigger = null;
this.forward = true;
public TriggerDate(final String stmName, final String triggerName, final boolean manageStart,
final Maneuver maneuver, final double threshold) {
this.stmName = stmName;
this.triggerName = triggerName;
this.massDepletionDelay = new MassDepletionDelay(triggerName, manageStart, maneuver);
this.manageStart = manageStart;
this.maneuver = maneuver;
this.threshold = threshold;
this.contribution = null;
this.trigger = null;
this.forward = true;
}
/** {@inheritDoc} */
@Override
public String getName() {
return columnName;
return triggerName;
}
/** {@inheritDoc}
......@@ -150,7 +155,14 @@ public class TriggerDateJacobianColumnGenerator
*/
@Override
public boolean yield(final SpacecraftState state) {
return !state.hasAdditionalState(stmName);
return !(state.hasAdditionalState(stmName) && state.hasAdditionalState(massDepletionDelay.getName()));
}
/** Get the mass depletion effect processor.
* @return mass depletion effect processor
*/
public MassDepletionDelay getMassDepletionDelay() {
return massDepletionDelay;
}
/** {@inheritDoc} */
......@@ -174,34 +186,29 @@ public class TriggerDateJacobianColumnGenerator
/** {@inheritDoc} */
@Override
public double[] getAdditionalState(final SpacecraftState state) {
// we check contribution rather than triggered because this method
// is called after maneuverTriggered and before resetState,
// when preparing the old state to be reset
final double[] c = contribution == null ? null : contribution.get(state.getDate());
if (c == null) {
// no thrust, no effect
return ZERO;
} else {
// part of the effect due to the acceleration performed at trigger time
// primary effect: full maneuver contribution at (delayed) trigger date
final double[] effect = getStm(state).operate(c);
if (trigger != null) {
// part of the effect due to mass change influence on current acceleration
final double[] parameters = maneuver.getParameters();
final Vector3D acceleration = maneuver.acceleration(state, parameters);
final double flowRate = maneuver.getPropulsionModel().getMassDerivatives(state, parameters);
final double ratio = state.getDate().durationFrom(trigger) * flowRate / state.getMass();
effect[3] += ratio * acceleration.getX();
effect[4] += ratio * acceleration.getY();
effect[5] += ratio * acceleration.getZ();
// secondary effect: maneuver change throughout thrust as mass depletion is delayed
final double[] secondary = state.getAdditionalState(massDepletionDelay.getName());
// sum up both effects
for (int i = 0; i < effect.length; ++i) {
effect[i] += secondary[i];
}
return effect;
}
}
/** {@inheritDoc}*/
......@@ -224,7 +231,7 @@ public class TriggerDateJacobianColumnGenerator
final Vector3D acceleration = maneuver.acceleration(stateWhenFiring, maneuver.getParameters());
// initialize derivatives computation
final double sign = (forward ^ manageStart) ? +1 : -1;
final double sign = (forward == manageStart) ? -1 : +1;
final RealVector rhs = MatrixUtils.createRealVector(STATE_DIMENSION);
rhs.setEntry(3, sign * acceleration.getX());
rhs.setEntry(4, sign * acceleration.getY());
......
......@@ -41,6 +41,9 @@ import org.orekit.forces.ForceModel;
import org.orekit.forces.gravity.NewtonianAttraction;
import org.orekit.forces.inertia.InertialForces;
import org.orekit.forces.maneuvers.Maneuver;
import org.orekit.forces.maneuvers.jacobians.Duration;
import org.orekit.forces.maneuvers.jacobians.MedianDate;
import org.orekit.forces.maneuvers.jacobians.TriggerDate;
import org.orekit.forces.maneuvers.trigger.AbstractManeuverTriggers;
import org.orekit.forces.maneuvers.trigger.ManeuverTriggers;
import org.orekit.frames.Frame;
......@@ -49,13 +52,10 @@ import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngle;
import org.orekit.propagation.AbstractMatricesHarvester;
import org.orekit.propagation.AdditionalStateProvider;
import org.orekit.propagation.DurationJacobianColumnGenerator;
import org.orekit.propagation.MatricesHarvester;
import org.orekit.propagation.MedianDateJacobianColumnGenerator;
import org.orekit.propagation.PropagationType;
import org.orekit.propagation.Propagator;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.TriggerDateJacobianColumnGenerator;
import org.orekit.propagation.events.EventDetector;
import org.orekit.propagation.events.ParameterDrivenDateIntervalDetector;
import org.orekit.propagation.integration.AbstractIntegratedPropagator;
......@@ -503,22 +503,22 @@ public class NumericalPropagator extends AbstractIntegratedPropagator {
map (d -> (ParameterDrivenDateIntervalDetector) d).
forEach(d -> {
if (d.getStartDriver().isSelected() || d.getMedianDriver().isSelected() || d.getDurationDriver().isSelected()) {
final TriggerDateJacobianColumnGenerator start =
final TriggerDate start =
manageTriggerDate(stmName, maneuver, amt, d.getStartDriver().getName(), true, d.getThreshold());
names.add(start.getName());
}
if (d.getStopDriver().isSelected() || d.getMedianDriver().isSelected() || d.getDurationDriver().isSelected()) {
final TriggerDateJacobianColumnGenerator stop =
final TriggerDate stop =
manageTriggerDate(stmName, maneuver, amt, d.getStopDriver().getName(), false, d.getThreshold());
names.add(stop.getName());
}
if (d.getMedianDriver().isSelected()) {
final MedianDateJacobianColumnGenerator median =
final MedianDate median =
manageMedianDate(d.getStartDriver().getName(), d.getStopDriver().getName(), d.getMedianDriver().getName());
names.add(median.getName());
}
if (d.getDurationDriver().isSelected()) {
final DurationJacobianColumnGenerator duration =
final Duration duration =
manageManeuverDuration(d.getStartDriver().getName(), d.getStopDriver().getName(), d.getDurationDriver().getName());
names.add(duration.getName());
}
......@@ -542,36 +542,37 @@ public class NumericalPropagator extends AbstractIntegratedPropagator {
* @return generator for the date driver
* @since 11.1
*/
private TriggerDateJacobianColumnGenerator manageTriggerDate(final String stmName,
final Maneuver maneuver,
final AbstractManeuverTriggers amt,
final String driverName,
final boolean start,
final double threshold) {
private TriggerDate manageTriggerDate(final String stmName,
final Maneuver maneuver,
final AbstractManeuverTriggers amt,
final String driverName,
final boolean start,
final double threshold) {
TriggerDateJacobianColumnGenerator triggerGenerator = null;
TriggerDate triggerGenerator = null;
// check if we already have set up the provider
for (final AdditionalStateProvider provider : getAdditionalStateProviders()) {
if (provider instanceof TriggerDateJacobianColumnGenerator &&
if (provider instanceof TriggerDate &&
provider.getName().equals(driverName)) {
// the Jacobian column generator has already been set up in a previous propagation
triggerGenerator = (TriggerDateJacobianColumnGenerator) provider;
triggerGenerator = (TriggerDate) provider;
break;
}
}
if (triggerGenerator == null) {
// this is the first time we need the Jacobian column generator, create it
triggerGenerator = new TriggerDateJacobianColumnGenerator(stmName, driverName,
start, maneuver, threshold);
triggerGenerator = new TriggerDate(stmName, driverName, start, maneuver, threshold);
amt.addResetter(triggerGenerator);
addAdditionalDerivativesProvider(triggerGenerator.getMassDepletionDelay());
addAdditionalStateProvider(triggerGenerator);
}
if (!getInitialIntegrationState().hasAdditionalState(driverName)) {
// add the initial Jacobian column if it is not already there
// (perhaps due to a previous propagation)
setInitialColumn(triggerGenerator.getMassDepletionDelay().getName(), new double[6]);
setInitialColumn(driverName, getHarvester().getInitialJacobianColumn(driverName));
}
......@@ -586,23 +587,23 @@ public class NumericalPropagator extends AbstractIntegratedPropagator {
* @return generator for the median driver
* @since 11.1
*/
private MedianDateJacobianColumnGenerator manageMedianDate(final String startName, final String stopName, final String medianName) {
private MedianDate manageMedianDate(final String startName, final String stopName, final String medianName) {
MedianDateJacobianColumnGenerator medianGenerator = null;
MedianDate medianGenerator = null;
// check if we already have set up the provider
for (final AdditionalStateProvider provider : getAdditionalStateProviders()) {
if (provider instanceof MedianDateJacobianColumnGenerator &&
if (provider instanceof MedianDate &&
provider.getName().equals(medianName)) {
// the Jacobian column generator has already been set up in a previous propagation
medianGenerator = (MedianDateJacobianColumnGenerator) provider;
medianGenerator = (MedianDate) provider;
break;
}
}
if (medianGenerator == null) {
// this is the first time we need the Jacobian column generator, create it
medianGenerator = new MedianDateJacobianColumnGenerator(startName, stopName, medianName);
medianGenerator = new MedianDate(startName, stopName, medianName);
addAdditionalStateProvider(medianGenerator);
}
......@@ -623,23 +624,23 @@ public class NumericalPropagator extends AbstractIntegratedPropagator {
* @return generator for the median driver
* @since 11.1
*/
private DurationJacobianColumnGenerator manageManeuverDuration(final String startName, final String stopName, final String durationName) {
private Duration manageManeuverDuration(final String startName, final String stopName, final String durationName) {
DurationJacobianColumnGenerator durationGenerator = null;
Duration durationGenerator = null;
// check if we already have set up the provider
for (final AdditionalStateProvider provider : getAdditionalStateProviders()) {
if (provider instanceof DurationJacobianColumnGenerator &&
if (provider instanceof Duration &&
provider.getName().equals(durationName)) {
// the Jacobian column generator has already been set up in a previous propagation
durationGenerator = (DurationJacobianColumnGenerator) provider;
durationGenerator = (Duration) provider;
break;
}
}
if (durationGenerator == null) {
// this is the first time we need the Jacobian column generator, create it
durationGenerator = new DurationJacobianColumnGenerator(startName, stopName, durationName);
durationGenerator = new Duration(startName, stopName, durationName);
addAdditionalStateProvider(durationGenerator);
}
......
......@@ -14,7 +14,7 @@
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.propagation.numerical;
package org.orekit.forces.maneuvers.jacobians;
import java.io.File;
......@@ -22,9 +22,10 @@ import java.io.IOException;
import java.io.PrintStream;
import java.nio.charset.StandardCharsets;
import java.util.ArrayList;
import java.util.Collections;
import java.util.Arrays;
import java.util.List;
import java.util.Locale;
import java.util.Optional;
import org.hipparchus.analysis.differentiation.UnivariateDerivative1;
import org.hipparchus.geometry.euclidean.threed.Rotation;
......@@ -58,8 +59,7 @@ import org.orekit.propagation.MatricesHarvester;
import org.orekit.propagation.Propagator;
import org.orekit.propagation.PropagatorsParallelizer;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.TriggerDateJacobianColumnGenerator;
import org.orekit.propagation.integration.AdditionalDerivativesProvider;
import org.orekit.propagation.numerical.NumericalPropagator;
import org.orekit.propagation.sampling.MultiSatStepHandler;
import org.orekit.propagation.sampling.OrekitStepInterpolator;
import org.orekit.time.AbsoluteDate;
......@@ -67,181 +67,105 @@ import org.orekit.time.DateComponents;
import org.orekit.time.TimeComponents;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.Constants;
import org.orekit.utils.DoubleArrayDictionary;
import org.orekit.utils.IERSConventions;
public class TriggersDerivativesTest {
@Test
public void testInterrupt() {
final AbsoluteDate firing = new AbsoluteDate(new DateComponents(2004, 1, 2),
new TimeComponents(4, 15, 34.080),
TimeScalesFactory.getUTC());
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, 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
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 =
new StateTransitionMatrixGenerator("dummy-1",
Collections.emptyList(),
propagator2.getAttitudeProvider());
propagator2.addAdditionalDerivativesProvider(dummyStmGenerator);
propagator2.setInitialState(propagator2.getInitialState().addAdditionalState(dummyStmGenerator.getName(), new double[36]));
propagator2.addAdditionalDerivativesProvider(new IntegrableJacobianColumnGenerator(dummyStmGenerator, "dummy-2"));
propagator2.setInitialState(propagator2.getInitialState().addAdditionalState("dummy-2", new double[6]));
propagator2.addAdditionalDerivativesProvider(new AdditionalDerivativesProvider() {
public String getName() { return "dummy-3"; }
public int getDimension() { return 1; }
public double[] derivatives(SpacecraftState s) { return new double[1]; }
});
propagator2.setInitialState(propagator2.getInitialState().addAdditionalState("dummy-3", new double[1]));
propagator2.addAdditionalStateProvider(new TriggerDateJacobianColumnGenerator(dummyStmGenerator.getName(), "dummy-4", true, null, 1.0e-6));
propagator2.addAdditionalStateProvider(new AdditionalStateProvider() {
public String getName() { return "dummy-5"; }
public double[] getAdditionalState(SpacecraftState s) { return new double[1]; }
});
final MatricesHarvester harvester2 = propagator2.setupMatricesComputation("stm", null, null);
final SpacecraftState intermediate = propagator2.propagate(firing.shiftedBy(0.5 * duration));
final RealMatrix stmI = harvester2.getStateTransitionMatrix(intermediate);
final RealMatrix jacobianI = harvester2.getParametersJacobian(intermediate);
// intermediate state has really different matrices, they are still building up
Assert.assertEquals(0.1253, stmI.subtract(stm1).getNorm1() / stm1.getNorm1(), 1.0e-4);
Assert.assertEquals(0.0172, jacobianI.subtract(jacobian1).getNorm1() / jacobian1.getNorm1(), 1.0e-4);
// restarting propagation where we left it
final SpacecraftState state2 = propagator2.propagate(firing.shiftedBy(2 * duration));
final RealMatrix stm2 = harvester2.getStateTransitionMatrix(state2);
final RealMatrix jacobian2 = harvester2.getParametersJacobian(state2);
// after completing the two-stage propagation, we get the same matrices
Assert.assertEquals(0.0, stm2.subtract(stm1).getNorm1(), 1.0e-13 * stm1.getNorm1());
Assert.assertEquals(0.0, jacobian2.subtract(jacobian1).getNorm1(), 2.0e-14 * jacobian1.getNorm1());
}
@Test
public void testDerivativeWrtStartTimeCartesianForward() {
doTestDerivativeWrtStartStopTime(true, OrbitType.CARTESIAN, true,
0.022, 0.012, 0.012, 0.013, 0.012, 0.021);
4.3e-6, 2.8e-6, 4.3e-6, 1.5e-5, 8.1e-6, 1.5e-5);
}
@Test
public void testDerivativeWrtStartTimeKeplerianForward() {
doTestDerivativeWrtStartStopTime(true, OrbitType.KEPLERIAN, true,
0.012, 0.011, 0.011, 0.012, 0.011, 0.017);
1.5e-5, 1.5e-5, 1.5e-5, 6.8e-6, 1.1e-4, 3.2e-5);
}
@Test
public void testDerivativeWrtStartTimeCartesianBackward() {
doTestDerivativeWrtStartStopTime(true, OrbitType.CARTESIAN, false,
0.022, 0.012, 0.012, 0.013, 0.012, 0.021);
9.7e-5, 9.7e-5, 9.7e-5, 4.7e-8, 4.7e-8, 4.7e-8);
}
@Test
public void testDerivativeWrtStartTimeKeplerianBackward() {
doTestDerivativeWrtStartStopTime(true, OrbitType.KEPLERIAN, false,
0.012, 0.011, 0.011, 0.012, 0.011, 0.017);
6.5e-8, 6.5e-8, 6.5e-8, 1.1e-6, 1.8e-6, 6.2e-7);
}
@Test
public void testDerivativeWrtStopTimeCartesianForward() {
doTestDerivativeWrtStartStopTime(false, OrbitType.CARTESIAN, true,
0.00033, 0.00045, 0.00040, 0.00022, 0.00020, 0.00010);
5.2e-10, 6.0e-8, 6.6e-11, 6.8e-12, 3.9e-11, 7.5e-12);
}
@Test
public void testDerivativeWrtStopTimeKeplerianForward() {
doTestDerivativeWrtStartStopTime(false, OrbitType.KEPLERIAN, true,
0.0011, 0.00020, 0.00002, 0.00082, 0.00008, 0.00019);
1.8e-11, 9.7e-12, 2.9e-12, 3.0e-9, 2.6e-9, 1.8e-9);
}
@Test
public void testDerivativeWrtStopTimeCartesianBackward() {
doTestDerivativeWrtStartStopTime(false, OrbitType.CARTESIAN, false,
0.00033, 0.00045, 0.00040, 0.00022, 0.00020, 0.00010);
9.7e-5, 9.7e-5, 9.7e-5, 1.5e-5, 2.1e-5, 1.5e-5);
}
@Test
public void testDerivativeWrtStopTimeKeplerianBackward() {
doTestDerivativeWrtStartStopTime(false, OrbitType.KEPLERIAN, false,
0.0011, 0.00020, 0.00002, 0.00082, 0.00008, 0.00019);
1.5e-5, 1.5e-5, 1.5e-5, 3.8e-6, 1.2e-4, 3.0e-4);
}
@Test
public void testDerivativeWrtMedianCartesianForward() {
doTestDerivativeWrtMedianDuration(true, OrbitType.CARTESIAN, true,
0.0011, 0.00020, 0.00002, 0.00082, 0.00008, 0.00019);
1.8e-5, 1.2e-5, 1.8e-5, 2.0e-2, 1.7e-2, 2.0e-2);
}