From 1f41a049803653062e22df2a4a60c830b248e5b3 Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Wed, 11 May 2022 10:57:40 +0200 Subject: [PATCH 1/2] Fixed additional states in integrated ephemerides. Fixes #925 --- .../AbstractIntegratedPropagator.java | 3 +- .../FieldAbstractIntegratedPropagator.java | 29 +++++- .../integration/FieldIntegratedEphemeris.java | 76 ++++++++++++++-- .../integration/IntegratedEphemeris.java | 4 +- .../FieldIntegratedEphemerisTest.java | 70 +++++++++++++++ .../FieldNumericalPropagatorTest.java | 88 +++++++++++++++++++ .../numerical/NumericalPropagatorTest.java | 8 +- 7 files changed, 266 insertions(+), 12 deletions(-) diff --git a/src/main/java/org/orekit/propagation/integration/AbstractIntegratedPropagator.java b/src/main/java/org/orekit/propagation/integration/AbstractIntegratedPropagator.java index 4ec638af6..ce60408a8 100644 --- a/src/main/java/org/orekit/propagation/integration/AbstractIntegratedPropagator.java +++ b/src/main/java/org/orekit/propagation/integration/AbstractIntegratedPropagator.java @@ -1051,7 +1051,7 @@ public abstract class AbstractIntegratedPropagator extends AbstractPropagator { /** Generated ephemeris. */ private BoundedPropagator ephemeris; - /** Variable used to store the last interpolator handled by the object.*/ + /** Last interpolator handled by the object.*/ private ODEStateInterpolator lastInterpolator; /** Set the end date. @@ -1107,6 +1107,7 @@ public abstract class AbstractIntegratedPropagator extends AbstractPropagator { // Update the model's finalTime with the last interpolator. model.finish(lastInterpolator.getCurrentState()); + // set up the boundary dates final double tI = model.getInitialTime(); final double tF = model.getFinalTime(); diff --git a/src/main/java/org/orekit/propagation/integration/FieldAbstractIntegratedPropagator.java b/src/main/java/org/orekit/propagation/integration/FieldAbstractIntegratedPropagator.java index daff0b661..66faf46c2 100644 --- a/src/main/java/org/orekit/propagation/integration/FieldAbstractIntegratedPropagator.java +++ b/src/main/java/org/orekit/propagation/integration/FieldAbstractIntegratedPropagator.java @@ -1022,6 +1022,9 @@ public abstract class FieldAbstractIntegratedPropagator ephemeris; + /** Last interpolator handled by the object.*/ + private FieldODEStateInterpolator lastInterpolator; + /** Set the end date. * @param endDate end date */ @@ -1038,11 +1041,15 @@ public abstract class FieldAbstractIntegratedPropagator getGeneratedEphemeris() { + // Each time we try to get the ephemeris, rebuild it using the last data. + buildEphemeris(); return ephemeris; } @@ -1050,11 +1057,26 @@ public abstract class FieldAbstractIntegratedPropagator interpolator) { model.handleStep(interpolator); + lastInterpolator = interpolator; } /** {@inheritDoc} */ @Override public void finish(final FieldODEStateAndDerivative finalState) { + buildEphemeris(); + } + + /** Method used to produce ephemeris at a given time. + * Can be used at multiple times, updating the ephemeris to + * its last state. + */ + private void buildEphemeris() { + // buildEphemeris was built in order to allow access to what was previously the finish method. + // This now allows to call it through getGeneratedEphemeris, therefore through an external call, + // which was not previously the case. + + // Update the model's finalTime with the last interpolator. + model.finish(lastInterpolator.getCurrentState()); // set up the boundary dates final T tI = model.getInitialTime(); @@ -1085,15 +1107,18 @@ public abstract class FieldAbstractIntegratedPropagator(startDate, minDate, maxDate, stateMapper, propagationType, model, - unmanaged, getAdditionalStateProviders(), names); + unmanaged, getAdditionalStateProviders(), + names, dimensions); } diff --git a/src/main/java/org/orekit/propagation/integration/FieldIntegratedEphemeris.java b/src/main/java/org/orekit/propagation/integration/FieldIntegratedEphemeris.java index e29dae4fd..1cec33028 100644 --- a/src/main/java/org/orekit/propagation/integration/FieldIntegratedEphemeris.java +++ b/src/main/java/org/orekit/propagation/integration/FieldIntegratedEphemeris.java @@ -16,6 +16,7 @@ */ package org.orekit.propagation.integration; +import java.util.Arrays; import java.util.Collections; import java.util.List; import java.util.Map; @@ -24,6 +25,7 @@ import org.hipparchus.CalculusFieldElement; import org.hipparchus.ode.FieldDenseOutputModel; import org.hipparchus.ode.FieldODEStateAndDerivative; import org.orekit.errors.OrekitException; +import org.orekit.errors.OrekitInternalError; import org.orekit.errors.OrekitMessages; import org.orekit.frames.Frame; import org.orekit.orbits.FieldOrbit; @@ -109,9 +111,9 @@ public class FieldIntegratedEphemeris > * @param unmanaged unmanaged additional states that must be simply copied * @param providers providers for pre-integrated states * @param equations names of additional equations - * @deprecated as of 11.1, replaced by {@link #FieldIntegratedEphemeris(FieldAbsoluteDate, + * @deprecated as of 11.2, replaced by {@link #FieldIntegratedEphemeris(FieldAbsoluteDate, * FieldAbsoluteDate, FieldAbsoluteDate, FieldStateMapper, PropagationType, - * FieldDenseOutputModel, FieldArrayDictionary, List, String[])} + * FieldDenseOutputModel, FieldArrayDictionary, List, String[], int[])} */ @Deprecated public FieldIntegratedEphemeris(final FieldAbsoluteDate startDate, @@ -137,7 +139,11 @@ public class FieldIntegratedEphemeris > * @param providers generators for pre-integrated states * @param equations names of additional equations * @since 11.1 + * @deprecated as of 11.2, replaced by {@link #FieldIntegratedEphemeris(FieldAbsoluteDate, + * FieldAbsoluteDate, FieldAbsoluteDate, FieldStateMapper, PropagationType, + * FieldDenseOutputModel, FieldArrayDictionary, List, String[], int[]) */ + @Deprecated public FieldIntegratedEphemeris(final FieldAbsoluteDate startDate, final FieldAbsoluteDate minDate, final FieldAbsoluteDate maxDate, final FieldStateMapper mapper, final PropagationType type, @@ -145,6 +151,31 @@ public class FieldIntegratedEphemeris > final FieldArrayDictionary unmanaged, final List> providers, final String[] equations) { + this(startDate, minDate, maxDate, mapper, type, model, + unmanaged, providers, equations, + remainingDimensions(model, unmanaged, providers, equations)); + } + + /** Creates a new instance of IntegratedEphemeris. + * @param startDate Start date of the integration (can be minDate or maxDate) + * @param minDate first date of the range + * @param maxDate last date of the range + * @param mapper mapper between raw double components and spacecraft state + * @param type type of orbit to output (mean or osculating) + * @param model underlying raw mathematical model + * @param unmanaged unmanaged additional states that must be simply copied + * @param providers generators for pre-integrated states + * @param equations names of additional equations + * @param dimensions dimensions of additional equations + * @since 11.2 + */ + public FieldIntegratedEphemeris(final FieldAbsoluteDate startDate, + final FieldAbsoluteDate minDate, final FieldAbsoluteDate maxDate, + final FieldStateMapper mapper, final PropagationType type, + final FieldDenseOutputModel model, + final FieldArrayDictionary unmanaged, + final List> providers, + final String[] equations, final int[] dimensions) { super(startDate.getField(), mapper.getAttitudeProvider()); @@ -162,12 +193,41 @@ public class FieldIntegratedEphemeris > } // set up providers to map the final elements of the model array to additional states + int index = 0; for (int i = 0; i < equations.length; ++i) { - addAdditionalStateProvider(new LocalGenerator(equations[i], i)); + addAdditionalStateProvider(new LocalGenerator(equations[i], index, dimensions[i])); + index += dimensions[i]; } } + /** Compute remaining dimensions for additional equations. + * @param tupe of the field elements + * @param model underlying raw mathematical model + * @param unmanaged unmanaged additional states that must be simply copied + * @param providers providers for pre-integrated states + * @param equations names of additional equations + * @return dimensions of additional equations + * @deprecated as of 11.2 this method is temporary and should be removed + * when the calling constructors are removed + * @since 11.2 + */ + @Deprecated + private static > int[] remainingDimensions(final FieldDenseOutputModel model, + final FieldArrayDictionary unmanaged, + final List> providers, + final String[] equations) { + final FieldODEStateAndDerivative osd = model.getInterpolatedState(model.getInitialTime()); + if (equations.length != osd.getNumberOfSecondaryStates()) { + throw new OrekitInternalError(null); + } + final int[] dimensions = new int[equations.length]; + for (int i = 0; i < dimensions.length; ++i) { + dimensions[i] = osd.getSecondaryStateDimension(i + 1); + } + return dimensions; + } + /** Interpolate the model at some date. * @param date desired interpolation date * @return state interpolated at date @@ -271,13 +331,18 @@ public class FieldIntegratedEphemeris > /** Index of the additional state. */ private final int index; + /** Dimension of the additional state. */ + private final int dimension; + /** Simple constructor. * @param name name of the additional state * @param index index of the additional state + * @param dimension dimension of the additional state */ - LocalGenerator(final String name, final int index) { + LocalGenerator(final String name, final int index, final int dimension) { this.name = name; this.index = index; + this.dimension = dimension; } /** {@inheritDoc} */ @@ -291,7 +356,8 @@ public class FieldIntegratedEphemeris > public T[] getAdditionalState(final FieldSpacecraftState state) { // extract the part of the interpolated array corresponding to the additional state - return getInterpolatedState(state.getDate()).getSecondaryState(index + 1); + final T[] combined = getInterpolatedState(state.getDate()).getSecondaryState(1); + return Arrays.copyOfRange(combined, index, index + dimension); } diff --git a/src/main/java/org/orekit/propagation/integration/IntegratedEphemeris.java b/src/main/java/org/orekit/propagation/integration/IntegratedEphemeris.java index 6c3785dfc..1483a96e3 100644 --- a/src/main/java/org/orekit/propagation/integration/IntegratedEphemeris.java +++ b/src/main/java/org/orekit/propagation/integration/IntegratedEphemeris.java @@ -190,8 +190,10 @@ public class IntegratedEphemeris } // set up providers to map the final elements of the model array to additional states + int index = 0; for (int i = 0; i < equations.length; ++i) { - addAdditionalStateProvider(new LocalGenerator(equations[i], i, dimensions[i])); + addAdditionalStateProvider(new LocalGenerator(equations[i], index, dimensions[i])); + index += dimensions[i]; } } diff --git a/src/test/java/org/orekit/propagation/integration/FieldIntegratedEphemerisTest.java b/src/test/java/org/orekit/propagation/integration/FieldIntegratedEphemerisTest.java index 059f32b08..d34221da9 100644 --- a/src/test/java/org/orekit/propagation/integration/FieldIntegratedEphemerisTest.java +++ b/src/test/java/org/orekit/propagation/integration/FieldIntegratedEphemerisTest.java @@ -18,21 +18,30 @@ package org.orekit.propagation.integration; import java.lang.reflect.InvocationTargetException; import java.lang.reflect.Method; +import java.util.Collections; import java.util.List; import org.hipparchus.CalculusFieldElement; import org.hipparchus.Field; import org.hipparchus.geometry.euclidean.threed.FieldVector3D; +import org.hipparchus.ode.FieldDenseOutputModel; +import org.hipparchus.ode.FieldExpandableODE; +import org.hipparchus.ode.FieldODEState; +import org.hipparchus.ode.FieldOrdinaryDifferentialEquation; +import org.hipparchus.ode.FieldSecondaryODE; import org.hipparchus.ode.events.Action; import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator; import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator; +import org.hipparchus.ode.nonstiff.EulerFieldIntegrator; import org.hipparchus.util.Decimal64Field; import org.hipparchus.util.MathArrays; import org.junit.Assert; import org.junit.Before; import org.junit.Test; import org.orekit.Utils; +import org.orekit.attitudes.InertialProvider; import org.orekit.errors.OrekitException; +import org.orekit.errors.OrekitInternalError; import org.orekit.errors.OrekitMessages; import org.orekit.forces.gravity.potential.GravityFieldFactory; import org.orekit.forces.gravity.potential.ICGEMFormatReader; @@ -40,10 +49,12 @@ import org.orekit.frames.FramesFactory; import org.orekit.orbits.FieldEquinoctialOrbit; import org.orekit.orbits.FieldOrbit; import org.orekit.orbits.OrbitType; +import org.orekit.orbits.PositionAngle; import org.orekit.propagation.FieldAdditionalStateProvider; import org.orekit.propagation.FieldBoundedPropagator; import org.orekit.propagation.FieldEphemerisGenerator; import org.orekit.propagation.FieldSpacecraftState; +import org.orekit.propagation.PropagationType; import org.orekit.propagation.analytical.FieldKeplerianPropagator; import org.orekit.propagation.events.FieldDateDetector; import org.orekit.propagation.numerical.FieldNumericalPropagator; @@ -74,6 +85,12 @@ public class FieldIntegratedEphemerisTest { doTestNoReset(Decimal64Field.getInstance()); } + @Deprecated + @Test + public void testDeprecated() { + doTestDeprecated(Decimal64Field.getInstance()); + } + private > void doTestNormalKeplerIntegration(Field field) { FieldOrbit initialOrbit = createOrbit(field); FieldNumericalPropagator numericalPropagator = createPropagator(field); @@ -215,6 +232,59 @@ public class FieldIntegratedEphemerisTest { } + @Deprecated + private > void doTestDeprecated(Field field) { + + EulerFieldIntegrator integ = new EulerFieldIntegrator<>(field, field.getZero().newInstance(1.0)); + FieldDenseOutputModel dom = new FieldDenseOutputModel<>(); + integ.addStepHandler(dom); + FieldExpandableODE eode = new FieldExpandableODE(new FieldOrdinaryDifferentialEquation() { + public int getDimension() { return 1; } + public T[] computeDerivatives(T t, T[] y) { return y; } + }); + eode.addSecondaryEquations(new FieldSecondaryODE() { + public int getDimension() { return 1; } + public T[] computeDerivatives(T t, T[] primary, T[] primaryDot, T[] secondary) { return secondary; } + }); + integ.integrate(eode, + new FieldODEState<>(field.getZero().newInstance(0.0), + MathArrays.buildArray(field, 1), + MathArrays.buildArray(field, 1, 1)), + field.getZero().newInstance(1.0)); + + FieldStateMapper mapper = new FieldStateMapper(FieldAbsoluteDate.getArbitraryEpoch(field), + field.getZero().newInstance(Constants.EIGEN5C_EARTH_MU), + OrbitType.CARTESIAN, PositionAngle.TRUE, + new InertialProvider(FramesFactory.getEME2000()), + FramesFactory.getEME2000()) { + public void mapStateToArray(FieldSpacecraftState state, T[] y, T[] yDot) {} + public FieldSpacecraftState mapArrayToState(FieldAbsoluteDate date, T[] y, T[] yDot, PropagationType type) { + return null; + } + }; + + try { + new FieldIntegratedEphemeris(FieldAbsoluteDate.getArbitraryEpoch(field), + FieldAbsoluteDate.getPastInfinity(field), + FieldAbsoluteDate.getFutureInfinity(field), + mapper, PropagationType.OSCULATING, + dom, Collections.emptyMap(), Collections.emptyList(), + new String[] { "equation-1", "equation-2" }); + Assert.fail("an exception should have been thrown"); + } catch (OrekitInternalError oie) { + // expected as only one equation could be handled properly by this deprecated constructor + } + + FieldIntegratedEphemeris ie = new FieldIntegratedEphemeris(FieldAbsoluteDate.getArbitraryEpoch(field), + FieldAbsoluteDate.getPastInfinity(field), + FieldAbsoluteDate.getFutureInfinity(field), + mapper, PropagationType.OSCULATING, + dom, Collections.emptyMap(), Collections.emptyList(), + new String[] { "equation-1" }); + Assert.assertNotNull(ie); + + } + @Before public void setUp() { mu = 3.9860047e14; diff --git a/src/test/java/org/orekit/propagation/numerical/FieldNumericalPropagatorTest.java b/src/test/java/org/orekit/propagation/numerical/FieldNumericalPropagatorTest.java index c6dfb81b4..54a27db5d 100644 --- a/src/test/java/org/orekit/propagation/numerical/FieldNumericalPropagatorTest.java +++ b/src/test/java/org/orekit/propagation/numerical/FieldNumericalPropagatorTest.java @@ -17,6 +17,7 @@ package org.orekit.propagation.numerical; import java.lang.reflect.Array; +import java.util.stream.IntStream; import org.hamcrest.MatcherAssert; import org.hipparchus.CalculusFieldElement; @@ -77,8 +78,10 @@ import org.orekit.propagation.integration.FieldAbstractIntegratedPropagator; import org.orekit.propagation.integration.FieldAdditionalDerivativesProvider; import org.orekit.propagation.sampling.FieldOrekitStepHandler; import org.orekit.propagation.sampling.FieldOrekitStepInterpolator; +import org.orekit.time.DateComponents; import org.orekit.time.FieldAbsoluteDate; import org.orekit.time.FieldTimeStamped; +import org.orekit.time.TimeComponents; import org.orekit.time.TimeScale; import org.orekit.time.TimeScalesFactory; import org.orekit.utils.Constants; @@ -1633,6 +1636,91 @@ public class FieldNumericalPropagatorTest { 1.14, 9.1, 140.3, 1066.7, 3306.9); } + @Test + public void testAdditionalDerivatives() { + doTestAdditionalDerivatives(Decimal64Field.getInstance()); + } + + private > void doTestAdditionalDerivatives(final Field field) { + + final Frame eme2000 = FramesFactory.getEME2000(); + final FieldAbsoluteDate date = new FieldAbsoluteDate<>(field, + new DateComponents(2008, 6, 23), + new TimeComponents(14, 0, 0), + TimeScalesFactory.getUTC()); + final FieldOrbit initialOrbit = new FieldKeplerianOrbit<>(field.getZero().newInstance(8000000.0), + field.getZero().newInstance(0.01), + field.getZero().newInstance(0.87), + field.getZero().newInstance(2.44), + field.getZero().newInstance(0.21), + field.getZero().newInstance(-1.05), + PositionAngle.MEAN, + eme2000, + date, field.getZero().newInstance(Constants.EIGEN5C_EARTH_MU)); + FieldNumericalPropagator propagator = new FieldNumericalPropagator<>(field, + new DormandPrince853FieldIntegrator<>(field, 0.02, 0.2, 1.0, 1.0e-3)); + propagator.resetInitialState(new FieldSpacecraftState<>(initialOrbit)); + + IntStream. + range(0, 2). + mapToObj(i -> new EmptyDerivativeProvider<>(field, "test_provider_" + i, new double[] { 10 * i, 20 * i })). + forEach(provider -> addDerivativeProvider(propagator, provider)); + + FieldEphemerisGenerator generator = propagator.getEphemerisGenerator(); + propagator.propagate(initialOrbit.getDate().shiftedBy(600)); + FieldBoundedPropagator ephemeris = generator.getGeneratedEphemeris(); + final FieldSpacecraftState finalState = ephemeris.propagate(initialOrbit.getDate().shiftedBy(300)); + Assert.assertEquals(2, finalState.getAdditionalStatesValues().size()); + Assert.assertEquals(2, finalState.getAdditionalState("test_provider_0").length); + Assert.assertEquals(0.0, finalState.getAdditionalState("test_provider_0")[0].getReal(), 1.0e-15); + Assert.assertEquals(0.0, finalState.getAdditionalState("test_provider_0")[1].getReal(), 1.0e-15); + Assert.assertEquals(2, finalState.getAdditionalState("test_provider_1").length); + Assert.assertEquals(10.0, finalState.getAdditionalState("test_provider_1")[0].getReal(), 1.0e-15); + Assert.assertEquals(20.0, finalState.getAdditionalState("test_provider_1")[1].getReal(), 1.0e-15); + } + + private > void addDerivativeProvider(FieldNumericalPropagator propagator, + EmptyDerivativeProvider provider) { + FieldSpacecraftState initialState = propagator.getInitialState(); + propagator.setInitialState(initialState.addAdditionalState(provider.getName(), provider.getInitialState())); + propagator.addAdditionalDerivativesProvider(provider); + } + + private static class EmptyDerivativeProvider> implements FieldAdditionalDerivativesProvider { + + private final Field field; + private final String name; + private final T[] state; + + public EmptyDerivativeProvider(Field field, String name, double[] state) { + this.field = field; + this.name = name; + this.state = MathArrays.buildArray(field, state.length); + for (int i = 0; i < state.length; ++i) { + this.state[i] = field.getZero().newInstance(state[i]); + } + } + + @Override + public T[] derivatives(FieldSpacecraftState s) { + return MathArrays.buildArray(field, getDimension()); + } + + @Override + public String getName() { + return name; + } + + @Override + public int getDimension() { + return state.length; + } + + public T[] getInitialState() { + return state; + } + } + private static > void doTestShift(final FieldCartesianOrbit orbit, final OrbitType orbitType, final PositionAngle angleType, final boolean withDerivatives, final double error60s, final double error120s, diff --git a/src/test/java/org/orekit/propagation/numerical/NumericalPropagatorTest.java b/src/test/java/org/orekit/propagation/numerical/NumericalPropagatorTest.java index 6634fa6e2..e0982face 100644 --- a/src/test/java/org/orekit/propagation/numerical/NumericalPropagatorTest.java +++ b/src/test/java/org/orekit/propagation/numerical/NumericalPropagatorTest.java @@ -1541,7 +1541,7 @@ public class NumericalPropagatorTest { IntStream. range(0, 2). - mapToObj(i -> new EmptyDerivativeProvider("test_provider_" + i, new double[] { 10 * i })). + mapToObj(i -> new EmptyDerivativeProvider("test_provider_" + i, new double[] { 10 * i, 20 * i })). forEach(provider -> addDerivativeProvider(propagator, provider)); EphemerisGenerator generator = propagator.getEphemerisGenerator(); @@ -1549,10 +1549,12 @@ public class NumericalPropagatorTest { BoundedPropagator ephemeris = generator.getGeneratedEphemeris(); final SpacecraftState finalState = ephemeris.propagate(initialOrbit.getDate().shiftedBy(300)); Assert.assertEquals(2, finalState.getAdditionalStatesValues().size()); - Assert.assertEquals(1, finalState.getAdditionalState("test_provider_0").length); + Assert.assertEquals(2, finalState.getAdditionalState("test_provider_0").length); Assert.assertEquals(0.0, finalState.getAdditionalState("test_provider_0")[0], 1.0e-15); - Assert.assertEquals(1, finalState.getAdditionalState("test_provider_1").length); + Assert.assertEquals(0.0, finalState.getAdditionalState("test_provider_0")[1], 1.0e-15); + Assert.assertEquals(2, finalState.getAdditionalState("test_provider_1").length); Assert.assertEquals(10.0, finalState.getAdditionalState("test_provider_1")[0], 1.0e-15); + Assert.assertEquals(20.0, finalState.getAdditionalState("test_provider_1")[1], 1.0e-15); } private void addDerivativeProvider(NumericalPropagator propagator, EmptyDerivativeProvider provider) { -- GitLab From e2ac8a5f38538d9dacb08b81f194a9ffd7c75bdb Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Wed, 11 May 2022 17:03:13 +0200 Subject: [PATCH 2/2] Don't loose additional derivatives when generating ephemeris. Fixes #927 --- .../integration/FieldIntegratedEphemeris.java | 74 ++++++++---------- .../integration/IntegratedEphemeris.java | 72 ++++++++--------- .../FieldIntegratedEphemerisTest.java | 77 +++++++++++++++++++ .../integration/IntegratedEphemerisTest.java | 64 +++++++++++++++ 4 files changed, 205 insertions(+), 82 deletions(-) diff --git a/src/main/java/org/orekit/propagation/integration/FieldIntegratedEphemeris.java b/src/main/java/org/orekit/propagation/integration/FieldIntegratedEphemeris.java index 1cec33028..aa5924ac5 100644 --- a/src/main/java/org/orekit/propagation/integration/FieldIntegratedEphemeris.java +++ b/src/main/java/org/orekit/propagation/integration/FieldIntegratedEphemeris.java @@ -101,6 +101,16 @@ public class FieldIntegratedEphemeris > /** Unmanaged additional states that must be simply copied. */ private final FieldArrayDictionary unmanaged; + /** Names of additional equations. + * @since 11.2 + */ + private final String[] equations; + + /** Dimensions of additional equations. + * @since 11.2 + */ + private final int[] dimensions; + /** Creates a new instance of IntegratedEphemeris. * @param startDate Start date of the integration (can be minDate or maxDate) * @param minDate first date of the range @@ -192,12 +202,8 @@ public class FieldIntegratedEphemeris > addAdditionalStateProvider(provider); } - // set up providers to map the final elements of the model array to additional states - int index = 0; - for (int i = 0; i < equations.length; ++i) { - addAdditionalStateProvider(new LocalGenerator(equations[i], index, dimensions[i])); - index += dimensions[i]; - } + this.equations = equations.clone(); + this.dimensions = dimensions.clone(); } @@ -322,44 +328,28 @@ public class FieldIntegratedEphemeris > return updateAdditionalStates(basicPropagate(getMinDate())); } - /** Local generator for additional state data. */ - private class LocalGenerator implements FieldAdditionalStateProvider { - - /** Name of the additional state. */ - private final String name; - - /** Index of the additional state. */ - private final int index; - - /** Dimension of the additional state. */ - private final int dimension; - - /** Simple constructor. - * @param name name of the additional state - * @param index index of the additional state - * @param dimension dimension of the additional state - */ - LocalGenerator(final String name, final int index, final int dimension) { - this.name = name; - this.index = index; - this.dimension = dimension; - } - - /** {@inheritDoc} */ - @Override - public String getName() { - return name; + /** {@inheritDoc} */ + @Override + protected FieldSpacecraftState updateAdditionalStates(final FieldSpacecraftState original) { + + FieldSpacecraftState updated = super.updateAdditionalStates(original); + + if (equations.length > 0) { + final FieldODEStateAndDerivative osd = getInterpolatedState(updated.getDate()); + final T[] combinedState = osd.getSecondaryState(1); + final T[] combinedDerivative = osd.getSecondaryDerivative(1); + int index = 0; + for (int i = 0; i < equations.length; ++i) { + final T[] state = Arrays.copyOfRange(combinedState, index, index + dimensions[i]); + final T[] derivative = Arrays.copyOfRange(combinedDerivative, index, index + dimensions[i]); + updated = updated. + addAdditionalState(equations[i], state). + addAdditionalStateDerivative(equations[i], derivative); + index += dimensions[i]; + } } - /** {@inheritDoc} */ - @Override - public T[] getAdditionalState(final FieldSpacecraftState state) { - - // extract the part of the interpolated array corresponding to the additional state - final T[] combined = getInterpolatedState(state.getDate()).getSecondaryState(1); - return Arrays.copyOfRange(combined, index, index + dimension); - - } + return updated; } diff --git a/src/main/java/org/orekit/propagation/integration/IntegratedEphemeris.java b/src/main/java/org/orekit/propagation/integration/IntegratedEphemeris.java index 1483a96e3..825ea90e9 100644 --- a/src/main/java/org/orekit/propagation/integration/IntegratedEphemeris.java +++ b/src/main/java/org/orekit/propagation/integration/IntegratedEphemeris.java @@ -99,6 +99,16 @@ public class IntegratedEphemeris /** Unmanaged additional states that must be simply copied. */ private final DoubleArrayDictionary unmanaged; + /** Names of additional equations. + * @since 11.2 + */ + private final String[] equations; + + /** Dimensions of additional equations. + * @since 11.2 + */ + private final int[] dimensions; + /** Creates a new instance of IntegratedEphemeris. * @param startDate Start date of the integration (can be minDate or maxDate) * @param minDate first date of the range @@ -189,12 +199,8 @@ public class IntegratedEphemeris addAdditionalStateProvider(provider); } - // set up providers to map the final elements of the model array to additional states - int index = 0; - for (int i = 0; i < equations.length; ++i) { - addAdditionalStateProvider(new LocalGenerator(equations[i], index, dimensions[i])); - index += dimensions[i]; - } + this.equations = equations.clone(); + this.dimensions = dimensions.clone(); } @@ -320,42 +326,28 @@ public class IntegratedEphemeris return updateAdditionalStates(basicPropagate(getMinDate())); } - /** Local generator for additional state data. */ - private class LocalGenerator implements AdditionalStateProvider { - - /** Name of the additional state. */ - private final String name; - - /** Index of the additional state. */ - private final int index; - - /** Dimension of the additional state. */ - private final int dimension; - - /** Simple constructor. - * @param name name of the additional state - * @param index index of the additional state - * @param dimension dimension of the additional state - */ - LocalGenerator(final String name, final int index, final int dimension) { - this.name = name; - this.index = index; - this.dimension = dimension; - } - - /** {@inheritDoc} */ - public String getName() { - return name; + /** {@inheritDoc} */ + @Override + protected SpacecraftState updateAdditionalStates(final SpacecraftState original) { + + SpacecraftState updated = super.updateAdditionalStates(original); + + if (equations.length > 0) { + final ODEStateAndDerivative osd = getInterpolatedState(updated.getDate()); + final double[] combinedState = osd.getSecondaryState(1); + final double[] combinedDerivative = osd.getSecondaryDerivative(1); + int index = 0; + for (int i = 0; i < equations.length; ++i) { + final double[] state = Arrays.copyOfRange(combinedState, index, index + dimensions[i]); + final double[] derivative = Arrays.copyOfRange(combinedDerivative, index, index + dimensions[i]); + updated = updated. + addAdditionalState(equations[i], state). + addAdditionalStateDerivative(equations[i], derivative); + index += dimensions[i]; + } } - /** {@inheritDoc} */ - public double[] getAdditionalState(final SpacecraftState state) { - - // extract the part of the interpolated array corresponding to the additional state - final double[] combined = getInterpolatedState(state.getDate()).getSecondaryState(1); - return Arrays.copyOfRange(combined, index, index + dimension); - - } + return updated; } diff --git a/src/test/java/org/orekit/propagation/integration/FieldIntegratedEphemerisTest.java b/src/test/java/org/orekit/propagation/integration/FieldIntegratedEphemerisTest.java index d34221da9..9e0c187f8 100644 --- a/src/test/java/org/orekit/propagation/integration/FieldIntegratedEphemerisTest.java +++ b/src/test/java/org/orekit/propagation/integration/FieldIntegratedEphemerisTest.java @@ -91,6 +91,11 @@ public class FieldIntegratedEphemerisTest { doTestDeprecated(Decimal64Field.getInstance()); } + @Test + public void testAdditionalDerivatives() { + doTestAdditionalDerivatives(Decimal64Field.getInstance()); + } + private > void doTestNormalKeplerIntegration(Field field) { FieldOrbit initialOrbit = createOrbit(field); FieldNumericalPropagator numericalPropagator = createPropagator(field); @@ -285,6 +290,56 @@ public class FieldIntegratedEphemerisTest { } + private > void doTestAdditionalDerivatives(final Field field) { + + final FieldOrbit initialOrbit = createOrbit(field); + FieldAbsoluteDate finalDate = initialOrbit.getDate().shiftedBy(10.0); + double[][] tolerances = FieldNumericalPropagator.tolerances(field.getZero().newInstance(1.0e-3), + initialOrbit, OrbitType.CARTESIAN); + DormandPrince853FieldIntegrator integrator = + new DormandPrince853FieldIntegrator<>(field, 1.0e-6, 10.0, tolerances[0], tolerances[1]); + integrator.setInitialStepSize(1.0e-3); + FieldNumericalPropagator propagator = new FieldNumericalPropagator<>(field, integrator); + final DerivativesProvider provider1 = new DerivativesProvider<>(field, "provider-1", 3); + propagator.addAdditionalDerivativesProvider(provider1); + final DerivativesProvider provider2 = new DerivativesProvider<>(field, "provider-2", 1); + propagator.addAdditionalDerivativesProvider(provider2); + final FieldEphemerisGenerator generator = propagator.getEphemerisGenerator(); + propagator.setInitialState(new FieldSpacecraftState<>(initialOrbit). + addAdditionalState(provider1.getName(), MathArrays.buildArray(field, provider1.getDimension())). + addAdditionalState(provider2.getName(), MathArrays.buildArray(field, provider2.getDimension()))); + propagator.propagate(finalDate); + FieldBoundedPropagator ephemeris = generator.getGeneratedEphemeris(); + + for (double dt = 0; dt < ephemeris.getMaxDate().durationFrom(ephemeris.getMinDate()).getReal(); dt += 0.1) { + FieldSpacecraftState state = ephemeris.propagate(ephemeris.getMinDate().shiftedBy(dt)); + checkState(dt, state, provider1); + checkState(dt, state, provider2); + } + + } + + private > void checkState(final double dt, final FieldSpacecraftState state, + final DerivativesProvider provider) { + + Assert.assertTrue(state.hasAdditionalState(provider.getName())); + Assert.assertEquals(provider.getDimension(), state.getAdditionalState(provider.getName()).length); + for (int i = 0; i < provider.getDimension(); ++i) { + Assert.assertEquals(i * dt, + state.getAdditionalState(provider.getName())[i].getReal(), + 4.0e-15 * i * dt); + } + + Assert.assertTrue(state.hasAdditionalStateDerivative(provider.getName())); + Assert.assertEquals(provider.getDimension(), state.getAdditionalStateDerivative(provider.getName()).length); + for (int i = 0; i < provider.getDimension(); ++i) { + Assert.assertEquals(i, + state.getAdditionalStateDerivative(provider.getName())[i].getReal(), + 2.0e-14 * i); + } + + } + @Before public void setUp() { mu = 3.9860047e14; @@ -319,4 +374,26 @@ public class FieldIntegratedEphemerisTest { } double mu; + + private static class DerivativesProvider> implements FieldAdditionalDerivativesProvider { + private final String name; + private final T[] derivatives; + DerivativesProvider(final Field field, final String name, final int dimension) { + this.name = name; + this.derivatives = MathArrays.buildArray(field, dimension); + for (int i = 0; i < dimension; ++i) { + derivatives[i] = field.getZero().newInstance(i); + } + } + public String getName() { + return name; + } + public int getDimension() { + return derivatives.length; + } + public T[] derivatives(final FieldSpacecraftState s) { + return derivatives; + } + } + } diff --git a/src/test/java/org/orekit/propagation/integration/IntegratedEphemerisTest.java b/src/test/java/org/orekit/propagation/integration/IntegratedEphemerisTest.java index ace8600f3..7ada76b83 100644 --- a/src/test/java/org/orekit/propagation/integration/IntegratedEphemerisTest.java +++ b/src/test/java/org/orekit/propagation/integration/IntegratedEphemerisTest.java @@ -214,6 +214,49 @@ public class IntegratedEphemerisTest { } + @Test + public void testAdditionalDerivatives() { + + AbsoluteDate finalDate = initialOrbit.getDate().shiftedBy(10.0); + double[][] tolerances = NumericalPropagator.tolerances(1.0e-3, initialOrbit, OrbitType.CARTESIAN); + DormandPrince853Integrator integrator = new DormandPrince853Integrator(1.0e-6, 10.0, tolerances[0], tolerances[1]); + integrator.setInitialStepSize(1.0e-3); + NumericalPropagator propagator = new NumericalPropagator(integrator); + final DerivativesProvider provider1 = new DerivativesProvider("provider-1", 3); + propagator.addAdditionalDerivativesProvider(provider1); + final DerivativesProvider provider2 = new DerivativesProvider("provider-2", 1); + propagator.addAdditionalDerivativesProvider(provider2); + final EphemerisGenerator generator = propagator.getEphemerisGenerator(); + propagator.setInitialState(new SpacecraftState(initialOrbit). + addAdditionalState(provider1.getName(), new double[provider1.getDimension()]). + addAdditionalState(provider2.getName(), new double[provider2.getDimension()])); + propagator.propagate(finalDate); + BoundedPropagator ephemeris = generator.getGeneratedEphemeris(); + + for (double dt = 0; dt < ephemeris.getMaxDate().durationFrom(ephemeris.getMinDate()); dt += 0.1) { + SpacecraftState state = ephemeris.propagate(ephemeris.getMinDate().shiftedBy(dt)); + checkState(dt, state, provider1); + checkState(dt, state, provider2); + } + + } + + private void checkState(final double dt, final SpacecraftState state, final DerivativesProvider provider) { + + Assert.assertTrue(state.hasAdditionalState(provider.getName())); + Assert.assertEquals(provider.getDimension(), state.getAdditionalState(provider.getName()).length); + for (int i = 0; i < provider.getDimension(); ++i) { + Assert.assertEquals(i * dt, state.getAdditionalState(provider.getName())[i], 4.0e-15 * i * dt); + } + + Assert.assertTrue(state.hasAdditionalStateDerivative(provider.getName())); + Assert.assertEquals(provider.getDimension(), state.getAdditionalStateDerivative(provider.getName()).length); + for (int i = 0; i < provider.getDimension(); ++i) { + Assert.assertEquals(i, state.getAdditionalStateDerivative(provider.getName())[i], 2.0e-14 * i); + } + + } + @Before public void setUp() { @@ -247,4 +290,25 @@ public class IntegratedEphemerisTest { private Orbit initialOrbit; private NumericalPropagator numericalPropagator; + private static class DerivativesProvider implements AdditionalDerivativesProvider { + private final String name; + private final double[] derivatives; + DerivativesProvider(final String name, final int dimension) { + this.name = name; + this.derivatives = new double[dimension]; + for (int i = 0; i < dimension; ++i) { + derivatives[i] = i; + } + } + public String getName() { + return name; + } + public int getDimension() { + return derivatives.length; + } + public double[] derivatives(final SpacecraftState s) { + return derivatives; + } + } + } -- GitLab