Fixed ephemeris generation with several derivatives providers.

Fixes #919

src/changes/changes.xml

@@ -28,6 +28,9 @@
       <action dev="louis" type="fix" issue="826">
         Fixed ephemeris generation using PropagatorParallelizer
       </action>
+      <action dev="luc" type="fix" issue="919">
+        Fixed ephemeris generation with several derivatives providers.
+      </action>
       <action dev="maxime" type="fix" issue="909">
         Fixed wrong implementation of NTW LOF frame.
       </action>

src/main/java/org/orekit/propagation/integration/AbstractIntegratedPropagator.java

@@ -1133,15 +1133,18 @@ public abstract class AbstractIntegratedPropagator extends AbstractPropagator {
             }
 
             // get the names of additional states managed by differential equations
-            final String[] names = new String[additionalDerivativesProviders.size()];
+            final String[] names      = new String[additionalDerivativesProviders.size()];
+            final int[]    dimensions = new int[additionalDerivativesProviders.size()];
             for (int i = 0; i < names.length; ++i) {
                 names[i] = additionalDerivativesProviders.get(i).getName();
+                dimensions[i] = additionalDerivativesProviders.get(i).getDimension();
             }
 
             // create the ephemeris
             ephemeris = new IntegratedEphemeris(startDate, minDate, maxDate,
                                                 stateMapper, propagationType, model,
-                                                unmanaged, getAdditionalStateProviders(), names);
+                                                unmanaged, getAdditionalStateProviders(),
+                                                names, dimensions);
 
         }
 

src/main/java/org/orekit/propagation/integration/IntegratedEphemeris.java

@@ -16,6 +16,7 @@
  */
 package org.orekit.propagation.integration;
 
+import java.util.Arrays;
 import java.util.List;
 import java.util.Map;
 
@@ -23,6 +24,7 @@ import org.hipparchus.ode.DenseOutputModel;
 import org.hipparchus.ode.ODEStateAndDerivative;
 import org.orekit.attitudes.AttitudeProvider;
 import org.orekit.errors.OrekitException;
+import org.orekit.errors.OrekitInternalError;
 import org.orekit.errors.OrekitMessages;
 import org.orekit.frames.Frame;
 import org.orekit.orbits.Orbit;
@@ -107,9 +109,9 @@ public class IntegratedEphemeris
      * @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 #IntegratedEphemeris(AbsoluteDate,
+     * @deprecated as of 11.1.2, replaced by {@link #IntegratedEphemeris(AbsoluteDate,
      * AbsoluteDate, AbsoluteDate, StateMapper, PropagationType, DenseOutputModel,
-     * DoubleArrayDictionary, List, String[])}
+     * DoubleArrayDictionary, List, String[], int[])}
      */
     @Deprecated
     public IntegratedEphemeris(final AbsoluteDate startDate,
@@ -134,7 +136,11 @@ public class IntegratedEphemeris
      * @param providers providers for pre-integrated states
      * @param equations names of additional equations
      * @since 11.1
+     * @deprecated as of 11.1.2, replaced by {@link #IntegratedEphemeris(AbsoluteDate,
+     * AbsoluteDate, AbsoluteDate, StateMapper, PropagationType, DenseOutputModel,
+     * DoubleArrayDictionary, List, String[], int[])}
      */
+    @Deprecated
     public IntegratedEphemeris(final AbsoluteDate startDate,
                                final AbsoluteDate minDate, final AbsoluteDate maxDate,
                                final StateMapper mapper, final PropagationType type,
@@ -142,6 +148,31 @@ public class IntegratedEphemeris
                                final DoubleArrayDictionary unmanaged,
                                final List<AdditionalStateProvider> 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 providers for pre-integrated states
+     * @param equations names of additional equations
+     * @param dimensions dimensions of additional equations
+     * @since 11.1.2
+     */
+    public IntegratedEphemeris(final AbsoluteDate startDate,
+                               final AbsoluteDate minDate, final AbsoluteDate maxDate,
+                               final StateMapper mapper, final PropagationType type,
+                               final DenseOutputModel model,
+                               final DoubleArrayDictionary unmanaged,
+                               final List<AdditionalStateProvider> providers,
+                               final String[] equations, final int[] dimensions) {
 
         super(mapper.getAttitudeProvider());
 
@@ -160,11 +191,37 @@ public class IntegratedEphemeris
 
         // set up providers to map the final elements of the model array to additional states
         for (int i = 0; i < equations.length; ++i) {
-            addAdditionalStateProvider(new LocalGenerator(equations[i], i));
+            addAdditionalStateProvider(new LocalGenerator(equations[i], i, dimensions[i]));
         }
 
     }
 
+    /** Compute remaining dimensions for additional equations.
+     * @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.1.2 this method is temporary and should be removed
+     * when the calling constructors are removed
+     * @since 11.1.2
+     */
+    @Deprecated
+    private static int[] remainingDimensions(final DenseOutputModel model,
+                                             final DoubleArrayDictionary unmanaged,
+                                             final List<AdditionalStateProvider> providers,
+                                             final String[] equations) {
+        final ODEStateAndDerivative 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
@@ -270,13 +327,18 @@ public class IntegratedEphemeris
         /** 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) {
-            this.name  = name;
-            this.index = index;
+        LocalGenerator(final String name, final int index, final int dimension) {
+            this.name      = name;
+            this.index     = index;
+            this.dimension = dimension;
         }
 
         /** {@inheritDoc} */
@@ -288,7 +350,8 @@ public class IntegratedEphemeris
         public double[] getAdditionalState(final SpacecraftState state) {
 
             // extract the part of the interpolated array corresponding to the additional state
-            return getInterpolatedState(state.getDate()).getSecondaryState(index + 1);
+            final double[] combined = getInterpolatedState(state.getDate()).getSecondaryState(1);
+            return Arrays.copyOfRange(combined, index, index + dimension);
 
         }
 

src/test/java/org/orekit/propagation/integration/IntegratedEphemerisTest.java

@@ -16,11 +16,19 @@
  */
 package org.orekit.propagation.integration;
 
+import java.util.Collections;
+
 import org.hipparchus.geometry.euclidean.threed.Vector3D;
 import org.hipparchus.linear.MatrixUtils;
 import org.hipparchus.linear.RealMatrix;
+import org.hipparchus.ode.DenseOutputModel;
+import org.hipparchus.ode.ExpandableODE;
+import org.hipparchus.ode.ODEState;
+import org.hipparchus.ode.OrdinaryDifferentialEquation;
+import org.hipparchus.ode.SecondaryODE;
 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
+import org.hipparchus.ode.nonstiff.EulerIntegrator;
 import org.junit.Assert;
 import org.junit.Before;
 import org.junit.Test;
@@ -28,15 +36,18 @@ import org.orekit.Utils;
 import org.orekit.attitudes.CelestialBodyPointed;
 import org.orekit.attitudes.InertialProvider;
 import org.orekit.bodies.CelestialBodyFactory;
+import org.orekit.errors.OrekitInternalError;
 import org.orekit.forces.gravity.potential.GravityFieldFactory;
 import org.orekit.forces.gravity.potential.ICGEMFormatReader;
 import org.orekit.frames.FramesFactory;
 import org.orekit.orbits.EquinoctialOrbit;
 import org.orekit.orbits.Orbit;
 import org.orekit.orbits.OrbitType;
+import org.orekit.orbits.PositionAngle;
 import org.orekit.propagation.BoundedPropagator;
 import org.orekit.propagation.EphemerisGenerator;
 import org.orekit.propagation.MatricesHarvester;
+import org.orekit.propagation.PropagationType;
 import org.orekit.propagation.SpacecraftState;
 import org.orekit.propagation.analytical.KeplerianPropagator;
 import org.orekit.propagation.numerical.NumericalPropagator;
@@ -155,6 +166,54 @@ public class IntegratedEphemerisTest {
 
     }
 
+    @Deprecated
+    @Test
+    public void testDeprecated() {
+
+        EulerIntegrator integ = new EulerIntegrator(1.0);
+        DenseOutputModel dom = new DenseOutputModel();
+        integ.addStepHandler(dom);
+        ExpandableODE eode = new ExpandableODE(new OrdinaryDifferentialEquation() {
+            public int getDimension() { return 1; }
+            public double[] computeDerivatives(double t, double[] y) { return y; }
+        });
+        eode.addSecondaryEquations(new SecondaryODE() {
+            public int getDimension() { return 1; }
+            public double[] computeDerivatives(double t, double[] primary,
+                                               double[] primaryDot, double[] secondary) { return secondary; }
+        });
+        integ.integrate(eode, new ODEState(0.0, new double[1], new double[1][1]), 1.0);
+
+        StateMapper mapper = new StateMapper(AbsoluteDate.ARBITRARY_EPOCH, Constants.EIGEN5C_EARTH_MU,
+                                             OrbitType.CARTESIAN, PositionAngle.TRUE,
+                                             new InertialProvider(FramesFactory.getEME2000()),
+                                             FramesFactory.getEME2000()) {
+            public void mapStateToArray(SpacecraftState state, double[] y, double[] yDot) {}
+            public SpacecraftState mapArrayToState(AbsoluteDate date, double[] y, double[] yDot, PropagationType type) {
+                return null;
+            }
+        };
+
+        try {
+            new IntegratedEphemeris(AbsoluteDate.ARBITRARY_EPOCH,
+                                    AbsoluteDate.PAST_INFINITY, AbsoluteDate.FUTURE_INFINITY,
+                                    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
+        }
+
+        IntegratedEphemeris ie = new IntegratedEphemeris(AbsoluteDate.ARBITRARY_EPOCH,
+                                                         AbsoluteDate.PAST_INFINITY, AbsoluteDate.FUTURE_INFINITY,
+                                                         mapper, PropagationType.OSCULATING,
+                                                         dom, Collections.emptyMap(), Collections.emptyList(),
+                                                         new String[] { "equation-1" });
+        Assert.assertNotNull(ie);
+
+    }
+
     @Before
     public void setUp() {
 

src/test/java/org/orekit/propagation/numerical/NumericalPropagatorTest.java

@@ -25,6 +25,7 @@ import java.util.List;
 import java.util.concurrent.ExecutionException;
 import java.util.concurrent.atomic.AtomicInteger;
 import java.util.function.Consumer;
+import java.util.stream.IntStream;
 import java.util.stream.Stream;
 
 import org.hamcrest.CoreMatchers;
@@ -81,6 +82,8 @@ import org.orekit.propagation.BoundedPropagator;
 import org.orekit.propagation.EphemerisGenerator;
 import org.orekit.propagation.FieldSpacecraftState;
 import org.orekit.propagation.SpacecraftState;
+import org.orekit.propagation.conversion.DormandPrince853IntegratorBuilder;
+import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
 import org.orekit.propagation.events.AbstractDetector;
 import org.orekit.propagation.events.ApsideDetector;
 import org.orekit.propagation.events.DateDetector;
@@ -1523,6 +1526,71 @@ public class NumericalPropagatorTest {
 
     }
 
+    @Test
+    public void testAdditionalDerivatives() {
+
+        final Frame eme2000 = FramesFactory.getEME2000();
+        final AbsoluteDate date = new AbsoluteDate(new DateComponents(2008, 6, 23),
+                                                   new TimeComponents(14, 0, 0),
+                                                   TimeScalesFactory.getUTC());
+        final Orbit initialOrbit = new KeplerianOrbit(8000000.0, 0.01, 0.87, 2.44, 0.21, -1.05, PositionAngle.MEAN,
+                                           eme2000,
+                                           date, Constants.EIGEN5C_EARTH_MU);
+        NumericalPropagatorBuilder builder = new NumericalPropagatorBuilder(initialOrbit, new DormandPrince853IntegratorBuilder(0.02, 0.2, 1.), PositionAngle.TRUE, 10);
+        NumericalPropagator propagator = builder.buildPropagator(builder.getSelectedNormalizedParameters());
+
+        IntStream.
+        range(0, 2).
+        mapToObj(i -> new EmptyDerivativeProvider("test_provider_" + i, new double[] { 10 * i })).
+        forEach(provider -> addDerivativeProvider(propagator, provider));
+
+        EphemerisGenerator generator = propagator.getEphemerisGenerator();
+        propagator.propagate(initialOrbit.getDate().shiftedBy(600));
+        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(0.0,  finalState.getAdditionalState("test_provider_0")[0], 1.0e-15);
+        Assert.assertEquals(1,    finalState.getAdditionalState("test_provider_1").length);
+        Assert.assertEquals(10.0, finalState.getAdditionalState("test_provider_1")[0], 1.0e-15);
+    }
+
+    private void addDerivativeProvider(NumericalPropagator propagator, EmptyDerivativeProvider provider) {
+        SpacecraftState initialState = propagator.getInitialState();
+        propagator.setInitialState(initialState.addAdditionalState(provider.getName(), provider.getInitialState()));
+        propagator.addAdditionalDerivativesProvider(provider);
+    }
+
+    private static class EmptyDerivativeProvider implements AdditionalDerivativesProvider {
+
+        private final String name;
+        private final double[] state;
+
+        public EmptyDerivativeProvider(String name, double[] state) {
+            this.name = name;
+            this.state = state;
+        }
+
+        @Override
+        public double[] derivatives(SpacecraftState s) {
+            return new double[getDimension()];
+        }
+
+        @Override
+        public String getName() {
+            return name;
+        }
+
+        @Override
+        public int getDimension() {
+            return state.length;
+        }
+
+        public double[] getInitialState() {
+            return state;
+        }
+    }
+
     /** Record the dates treated by the handler.
      *  If they are out of an interval defined by a start and final date.
      */