NumericalPropagator.java

/* Copyright 2002-2022 CS GROUP
 * Licensed to CS GROUP (CS) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * CS licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.orekit.propagation.numerical;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collections;
import java.util.List;

import org.hipparchus.exception.LocalizedCoreFormats;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.linear.MatrixUtils;
import org.hipparchus.linear.QRDecomposition;
import org.hipparchus.linear.RealMatrix;
import org.hipparchus.ode.ODEIntegrator;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.Precision;
import org.orekit.annotation.DefaultDataContext;
import org.orekit.attitudes.Attitude;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.data.DataContext;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitIllegalArgumentException;
import org.orekit.errors.OrekitInternalError;
import org.orekit.errors.OrekitMessages;
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;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngle;
import org.orekit.propagation.AbstractMatricesHarvester;
import org.orekit.propagation.AdditionalStateProvider;
import org.orekit.propagation.MatricesHarvester;
import org.orekit.propagation.PropagationType;
import org.orekit.propagation.Propagator;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.events.EventDetector;
import org.orekit.propagation.events.ParameterDrivenDateIntervalDetector;
import org.orekit.propagation.integration.AbstractIntegratedPropagator;
import org.orekit.propagation.integration.AdditionalDerivativesProvider;
import org.orekit.propagation.integration.StateMapper;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.AbsolutePVCoordinates;
import org.orekit.utils.DoubleArrayDictionary;
import org.orekit.utils.PVCoordinates;
import org.orekit.utils.ParameterDriver;
import org.orekit.utils.ParameterDriversList;
import org.orekit.utils.ParameterDriversList.DelegatingDriver;
import org.orekit.utils.ParameterObserver;
import org.orekit.utils.TimeStampedPVCoordinates;

/** This class propagates {@link org.orekit.orbits.Orbit orbits} using
 * numerical integration.
 * <p>Numerical propagation is much more accurate than analytical propagation
 * like for example {@link org.orekit.propagation.analytical.KeplerianPropagator
 * Keplerian} or {@link org.orekit.propagation.analytical.EcksteinHechlerPropagator
 * Eckstein-Hechler}, but requires a few more steps to set up to be used properly.
 * Whereas analytical propagators are configured only thanks to their various
 * constructors and can be used immediately after construction, numerical propagators
 * configuration involve setting several parameters between construction time
 * and propagation time.</p>
 * <p>The configuration parameters that can be set are:</p>
 * <ul>
 *   <li>the initial spacecraft state ({@link #setInitialState(SpacecraftState)})</li>
 *   <li>the central attraction coefficient ({@link #setMu(double)})</li>
 *   <li>the various force models ({@link #addForceModel(ForceModel)},
 *   {@link #removeForceModels()})</li>
 *   <li>the {@link OrbitType type} of orbital parameters to be used for propagation
 *   ({@link #setOrbitType(OrbitType)}),</li>
 *   <li>the {@link PositionAngle type} of position angle to be used in orbital parameters
 *   to be used for propagation where it is relevant ({@link
 *   #setPositionAngleType(PositionAngle)}),</li>
 *   <li>whether {@link MatricesHarvester state transition matrices and Jacobians matrices}
 *   should be propagated along with orbital state ({@link
 *   #setupMatricesComputation(String, RealMatrix, DoubleArrayDictionary)}),</li>
 *   <li>whether {@link org.orekit.propagation.integration.AdditionalDerivativesProvider additional derivatives}
 *   should be propagated along with orbital state ({@link
 *   #addAdditionalDerivativesProvider(AdditionalDerivativesProvider)}),</li>
 *   <li>the discrete events that should be triggered during propagation
 *   ({@link #addEventDetector(EventDetector)},
 *   {@link #clearEventsDetectors()})</li>
 *   <li>the binding logic with the rest of the application ({@link #getMultiplexer()})</li>
 * </ul>
 * <p>From these configuration parameters, only the initial state is mandatory. The default
 * propagation settings are in {@link OrbitType#EQUINOCTIAL equinoctial} parameters with
 * {@link PositionAngle#TRUE true} longitude argument. If the central attraction coefficient
 * is not explicitly specified, the one used to define the initial orbit will be used.
 * However, specifying only the initial state and perhaps the central attraction coefficient
 * would mean the propagator would use only Keplerian forces. In this case, the simpler {@link
 * org.orekit.propagation.analytical.KeplerianPropagator KeplerianPropagator} class would
 * perhaps be more effective.</p>
 * <p>The underlying numerical integrator set up in the constructor may also have its own
 * configuration parameters. Typical configuration parameters for adaptive stepsize integrators
 * are the min, max and perhaps start step size as well as the absolute and/or relative errors
 * thresholds.</p>
 * <p>The state that is seen by the integrator is a simple seven elements double array.
 * The six first elements are either:
 * <ul>
 *   <li>the {@link org.orekit.orbits.EquinoctialOrbit equinoctial orbit parameters} (a, e<sub>x</sub>,
 *   e<sub>y</sub>, h<sub>x</sub>, h<sub>y</sub>, λ<sub>M</sub> or λ<sub>E</sub>
 *   or λ<sub>v</sub>) in meters and radians,</li>
 *   <li>the {@link org.orekit.orbits.KeplerianOrbit Keplerian orbit parameters} (a, e, i, ω, Ω,
 *   M or E or v) in meters and radians,</li>
 *   <li>the {@link org.orekit.orbits.CircularOrbit circular orbit parameters} (a, e<sub>x</sub>, e<sub>y</sub>, i,
 *   Ω, α<sub>M</sub> or α<sub>E</sub> or α<sub>v</sub>) in meters
 *   and radians,</li>
 *   <li>the {@link org.orekit.orbits.CartesianOrbit Cartesian orbit parameters} (x, y, z, v<sub>x</sub>,
 *   v<sub>y</sub>, v<sub>z</sub>) in meters and meters per seconds.
 * </ul>
 * <p> The last element is the mass in kilograms and changes only during thrusters firings
 *
 * <p>The following code snippet shows a typical setting for Low Earth Orbit propagation in
 * equinoctial parameters and true longitude argument:</p>
 * <pre>
 * final double dP       = 0.001;
 * final double minStep  = 0.001;
 * final double maxStep  = 500;
 * final double initStep = 60;
 * final double[][] tolerance = NumericalPropagator.tolerances(dP, orbit, OrbitType.EQUINOCTIAL);
 * AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tolerance[0], tolerance[1]);
 * integrator.setInitialStepSize(initStep);
 * propagator = new NumericalPropagator(integrator);
 * </pre>
 * <p>By default, at the end of the propagation, the propagator resets the initial state to the final state,
 * thus allowing a new propagation to be started from there without recomputing the part already performed.
 * This behaviour can be chenged by calling {@link #setResetAtEnd(boolean)}.
 * </p>
 * <p>Beware the same instance cannot be used simultaneously by different threads, the class is <em>not</em>
 * thread-safe.</p>
 *
 * @see SpacecraftState
 * @see ForceModel
 * @see org.orekit.propagation.sampling.OrekitStepHandler
 * @see org.orekit.propagation.sampling.OrekitFixedStepHandler
 * @see org.orekit.propagation.integration.IntegratedEphemeris
 * @see TimeDerivativesEquations
 *
 * @author Mathieu Rom&eacute;ro
 * @author Luc Maisonobe
 * @author Guylaine Prat
 * @author Fabien Maussion
 * @author V&eacute;ronique Pommier-Maurussane
 */
public class NumericalPropagator extends AbstractIntegratedPropagator {

    /** Space dimension. */
    private static final int SPACE_DIMENSION = 3;

    /** State dimension. */
    private static final int STATE_DIMENSION = 2 * SPACE_DIMENSION;

    /** Threshold for matrix solving. */
    private static final double THRESHOLD = Precision.SAFE_MIN;

    /** Force models used during the extrapolation of the orbit. */
    private final List<ForceModel> forceModels;

    /** boolean to ignore or not the creation of a NewtonianAttraction. */
    private boolean ignoreCentralAttraction;

    /** Create a new instance of NumericalPropagator, based on orbit definition mu.
     * After creation, the instance is empty, i.e. the attitude provider is set to an
     * unspecified default law and there are no perturbing forces at all.
     * This means that if {@link #addForceModel addForceModel} is not
     * called after creation, the integrated orbit will follow a Keplerian
     * evolution only. The defaults are {@link OrbitType#EQUINOCTIAL}
     * for {@link #setOrbitType(OrbitType) propagation
     * orbit type} and {@link PositionAngle#TRUE} for {@link
     * #setPositionAngleType(PositionAngle) position angle type}.
     *
     * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
     *
     * @param integrator numerical integrator to use for propagation.
     * @see #NumericalPropagator(ODEIntegrator, AttitudeProvider)
     */
    @DefaultDataContext
    public NumericalPropagator(final ODEIntegrator integrator) {
        this(integrator,
                Propagator.getDefaultLaw(DataContext.getDefault().getFrames()));
    }

    /** Create a new instance of NumericalPropagator, based on orbit definition mu.
     * After creation, the instance is empty, i.e. the attitude provider is set to an
     * unspecified default law and there are no perturbing forces at all.
     * This means that if {@link #addForceModel addForceModel} is not
     * called after creation, the integrated orbit will follow a Keplerian
     * evolution only. The defaults are {@link OrbitType#EQUINOCTIAL}
     * for {@link #setOrbitType(OrbitType) propagation
     * orbit type} and {@link PositionAngle#TRUE} for {@link
     * #setPositionAngleType(PositionAngle) position angle type}.
     * @param integrator numerical integrator to use for propagation.
     * @param attitudeProvider the attitude law.
     * @since 10.1
     */
    public NumericalPropagator(final ODEIntegrator integrator,
                               final AttitudeProvider attitudeProvider) {
        super(integrator, PropagationType.MEAN);
        forceModels             = new ArrayList<ForceModel>();
        ignoreCentralAttraction = false;
        initMapper();
        setAttitudeProvider(attitudeProvider);
        clearStepHandlers();
        setOrbitType(OrbitType.EQUINOCTIAL);
        setPositionAngleType(PositionAngle.TRUE);
    }

    /** Set the flag to ignore or not the creation of a {@link NewtonianAttraction}.
     * @param ignoreCentralAttraction if true, {@link NewtonianAttraction} is <em>not</em>
     * added automatically if missing
     */
    public void setIgnoreCentralAttraction(final boolean ignoreCentralAttraction) {
        this.ignoreCentralAttraction = ignoreCentralAttraction;
    }

     /** Set the central attraction coefficient μ.
      * <p>
      * Setting the central attraction coefficient is
      * equivalent to {@link #addForceModel(ForceModel) add}
      * a {@link NewtonianAttraction} force model.
      * </p>
     * @param mu central attraction coefficient (m³/s²)
     * @see #addForceModel(ForceModel)
     * @see #getAllForceModels()
     */
    public void setMu(final double mu) {
        if (ignoreCentralAttraction) {
            superSetMu(mu);
        } else {
            addForceModel(new NewtonianAttraction(mu));
        }
    }

    /** Set the central attraction coefficient μ only in upper class.
     * @param mu central attraction coefficient (m³/s²)
     */
    private void superSetMu(final double mu) {
        super.setMu(mu);
    }

    /** Check if Newtonian attraction force model is available.
     * <p>
     * Newtonian attraction is always the last force model in the list.
     * </p>
     * @return true if Newtonian attraction force model is available
     */
    private boolean hasNewtonianAttraction() {
        final int last = forceModels.size() - 1;
        return last >= 0 && forceModels.get(last) instanceof NewtonianAttraction;
    }

    /** Add a force model.
     * <p>If this method is not called at all, the integrated orbit will follow
     * a Keplerian evolution only.</p>
     * @param model {@link ForceModel} to add (it can be either a perturbing force
     * model or an instance of {@link NewtonianAttraction})
     * @see #removeForceModels()
     * @see #setMu(double)
     */
    public void addForceModel(final ForceModel model) {

        if (model instanceof NewtonianAttraction) {
            // we want to add the central attraction force model

            try {
                // ensure we are notified of any mu change
                model.getParametersDrivers().get(0).addObserver(new ParameterObserver() {
                    /** {@inheritDoc} */
                    @Override
                    public void valueChanged(final double previousValue, final ParameterDriver driver) {
                        superSetMu(driver.getValue());
                    }
                });
            } catch (OrekitException oe) {
                // this should never happen
                throw new OrekitInternalError(oe);
            }

            if (hasNewtonianAttraction()) {
                // there is already a central attraction model, replace it
                forceModels.set(forceModels.size() - 1, model);
            } else {
                // there are no central attraction model yet, add it at the end of the list
                forceModels.add(model);
            }
        } else {
            // we want to add a perturbing force model
            if (hasNewtonianAttraction()) {
                // insert the new force model before Newtonian attraction,
                // which should always be the last one in the list
                forceModels.add(forceModels.size() - 1, model);
            } else {
                // we only have perturbing force models up to now, just append at the end of the list
                forceModels.add(model);
            }
        }

    }

    /** Remove all force models (except central attraction).
     * <p>Once all perturbing forces have been removed (and as long as no new force
     * model is added), the integrated orbit will follow a Keplerian evolution
     * only.</p>
     * @see #addForceModel(ForceModel)
     */
    public void removeForceModels() {
        final int last = forceModels.size() - 1;
        if (hasNewtonianAttraction()) {
            // preserve the Newtonian attraction model at the end
            final ForceModel newton = forceModels.get(last);
            forceModels.clear();
            forceModels.add(newton);
        } else {
            forceModels.clear();
        }
    }

    /** Get all the force models, perturbing forces and Newtonian attraction included.
     * @return list of perturbing force models, with Newtonian attraction being the
     * last one
     * @see #addForceModel(ForceModel)
     * @see #setMu(double)
     */
    public List<ForceModel> getAllForceModels() {
        return Collections.unmodifiableList(forceModels);
    }

    /** Set propagation orbit type.
     * @param orbitType orbit type to use for propagation, null for
     * propagating using {@link org.orekit.utils.AbsolutePVCoordinates} rather than {@link Orbit}
     */
    public void setOrbitType(final OrbitType orbitType) {
        super.setOrbitType(orbitType);
    }

    /** Get propagation parameter type.
     * @return orbit type used for propagation, null for
     * propagating using {@link org.orekit.utils.AbsolutePVCoordinates} rather than {@link Orbit}
     */
    public OrbitType getOrbitType() {
        return super.getOrbitType();
    }

    /** Set position angle type.
     * <p>
     * The position parameter type is meaningful only if {@link
     * #getOrbitType() propagation orbit type}
     * support it. As an example, it is not meaningful for propagation
     * in {@link OrbitType#CARTESIAN Cartesian} parameters.
     * </p>
     * @param positionAngleType angle type to use for propagation
     */
    public void setPositionAngleType(final PositionAngle positionAngleType) {
        super.setPositionAngleType(positionAngleType);
    }

    /** Get propagation parameter type.
     * @return angle type to use for propagation
     */
    public PositionAngle getPositionAngleType() {
        return super.getPositionAngleType();
    }

    /** Set the initial state.
     * @param initialState initial state
     */
    public void setInitialState(final SpacecraftState initialState) {
        resetInitialState(initialState);
    }

    /** {@inheritDoc} */
    public void resetInitialState(final SpacecraftState state) {
        super.resetInitialState(state);
        if (!hasNewtonianAttraction()) {
            // use the state to define central attraction
            setMu(state.getMu());
        }
        setStartDate(state.getDate());
    }

    /** Get the names of the parameters in the matrix returned by {@link MatricesHarvester#getParametersJacobian}.
     * @return names of the parameters (i.e. columns) of the Jacobian matrix
     */
    List<String> getJacobiansColumnsNames() {
        final List<String> columnsNames = new ArrayList<>();
        for (final ForceModel forceModel : getAllForceModels()) {
            for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
                if (driver.isSelected() && !columnsNames.contains(driver.getName())) {
                    columnsNames.add(driver.getName());
                }
            }
        }
        Collections.sort(columnsNames);
        return columnsNames;
    }

    /** {@inheritDoc} */
    @Override
    protected AbstractMatricesHarvester createHarvester(final String stmName, final RealMatrix initialStm,
                                                        final DoubleArrayDictionary initialJacobianColumns) {
        return new NumericalPropagationHarvester(this, stmName, initialStm, initialJacobianColumns);
    }

    /** {@inheritDoc} */
    @Override
    protected void setUpStmAndJacobianGenerators() {

        final AbstractMatricesHarvester harvester = getHarvester();
        if (harvester != null) {

            // set up the additional equations and additional state providers
            final StateTransitionMatrixGenerator stmGenerator = setUpStmGenerator();
            final List<String> triggersDates = setUpTriggerDatesJacobiansColumns(stmGenerator.getName());
            setUpRegularParametersJacobiansColumns(stmGenerator, triggersDates);

            // as we are now starting the propagation, everything is configured
            // we can freeze the names in the harvester
            harvester.freezeColumnsNames();

        }

    }

    /** Set up the State Transition Matrix Generator.
     * @return State Transition Matrix Generator
     * @since 11.1
     */
    private StateTransitionMatrixGenerator setUpStmGenerator() {

        final AbstractMatricesHarvester harvester = getHarvester();

        // add the STM generator corresponding to the current settings, and setup state accordingly
        StateTransitionMatrixGenerator stmGenerator = null;
        for (final AdditionalDerivativesProvider equations : getAdditionalDerivativesProviders()) {
            if (equations instanceof StateTransitionMatrixGenerator &&
                equations.getName().equals(harvester.getStmName())) {
                // the STM generator has already been set up in a previous propagation
                stmGenerator = (StateTransitionMatrixGenerator) equations;
                break;
            }
        }
        if (stmGenerator == null) {
            // this is the first time we need the STM generate, create it
            stmGenerator = new StateTransitionMatrixGenerator(harvester.getStmName(), getAllForceModels(), getAttitudeProvider());
            addAdditionalDerivativesProvider(stmGenerator);
        }

        if (!getInitialIntegrationState().hasAdditionalState(harvester.getStmName())) {
            // add the initial State Transition Matrix if it is not already there
            // (perhaps due to a previous propagation)
            setInitialState(stmGenerator.setInitialStateTransitionMatrix(getInitialState(),
                                                                         harvester.getInitialStateTransitionMatrix(),
                                                                         getOrbitType(),
                                                                         getPositionAngleType()));
        }

        return stmGenerator;

    }

    /** Set up the Jacobians columns generator dedicated to trigger dates.
     * @param stmName name of the State Transition Matrix state
     * @return names of the columns corresponding to trigger dates
     * @since 11.1
     */
    private List<String> setUpTriggerDatesJacobiansColumns(final String stmName) {

        final List<String> names = new ArrayList<>();
        for (final ForceModel forceModel : getAllForceModels()) {
            if (forceModel instanceof Maneuver) {
                final Maneuver         maneuver         = (Maneuver) forceModel;
                final ManeuverTriggers maneuverTriggers = maneuver.getManeuverTriggers();
                if (maneuverTriggers instanceof AbstractManeuverTriggers) {

                    // FIXME: when issue https://gitlab.orekit.org/orekit/orekit/-/issues/854 is solved
                    // the previous if statement and the following cast should be removed as the following
                    // code should really be done for all ManeuverTriggers and not only AbstractManeuverTriggers
                    final AbstractManeuverTriggers amt = (AbstractManeuverTriggers) maneuverTriggers;

                    amt.getEventsDetectors().
                        filter(d -> d instanceof ParameterDrivenDateIntervalDetector).
                        map (d -> (ParameterDrivenDateIntervalDetector) d).
                        forEach(d -> {
                            if (d.getStartDriver().isSelected() || d.getMedianDriver().isSelected() || d.getDurationDriver().isSelected()) {
                                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 TriggerDate stop =
                                                manageTriggerDate(stmName, maneuver, amt, d.getStopDriver().getName(),  false, d.getThreshold());
                                names.add(stop.getName());
                            }
                            if (d.getMedianDriver().isSelected()) {
                                final MedianDate median =
                                                manageMedianDate(d.getStartDriver().getName(), d.getStopDriver().getName(), d.getMedianDriver().getName());
                                names.add(median.getName());
                            }
                            if (d.getDurationDriver().isSelected()) {
                                final Duration duration =
                                                manageManeuverDuration(d.getStartDriver().getName(), d.getStopDriver().getName(), d.getDurationDriver().getName());
                                names.add(duration.getName());
                            }
                        });

                }
            }
        }

        return names;

    }

    /** Manage a maneuver trigger date.
     * @param stmName name of the State Transition Matrix state
     * @param maneuver maneuver force model
     * @param amt trigger to which the driver is bound
     * @param driverName name of the date driver
     * @param start if true, the driver is a maneuver start
     * @param threshold event detector threshold
     * @return generator for the date driver
     * @since 11.1
     */
    private TriggerDate manageTriggerDate(final String stmName,
                                          final Maneuver maneuver,
                                          final AbstractManeuverTriggers amt,
                                          final String driverName,
                                          final boolean start,
                                          final double threshold) {

        TriggerDate triggerGenerator = null;

        // check if we already have set up the provider
        for (final AdditionalStateProvider provider : getAdditionalStateProviders()) {
            if (provider instanceof TriggerDate &&
                provider.getName().equals(driverName)) {
                // the Jacobian column generator has already been set up in a previous propagation
                triggerGenerator = (TriggerDate) provider;
                break;
            }
        }

        if (triggerGenerator == null) {
            // this is the first time we need the Jacobian column generator, create it
            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));
        }

        return triggerGenerator;

    }

    /** Manage a maneuver median date.
     * @param startName name of the start driver
     * @param stopName name of the stop driver
     * @param medianName name of the median driver
     * @return generator for the median driver
     * @since 11.1
     */
    private MedianDate manageMedianDate(final String startName, final String stopName, final String medianName) {

        MedianDate medianGenerator = null;

        // check if we already have set up the provider
        for (final AdditionalStateProvider provider : getAdditionalStateProviders()) {
            if (provider instanceof MedianDate &&
                provider.getName().equals(medianName)) {
                // the Jacobian column generator has already been set up in a previous propagation
                medianGenerator = (MedianDate) provider;
                break;
            }
        }

        if (medianGenerator == null) {
            // this is the first time we need the Jacobian column generator, create it
            medianGenerator = new MedianDate(startName, stopName, medianName);
            addAdditionalStateProvider(medianGenerator);
        }

        if (!getInitialIntegrationState().hasAdditionalState(medianName)) {
            // add the initial Jacobian column if it is not already there
            // (perhaps due to a previous propagation)
            setInitialColumn(medianName, getHarvester().getInitialJacobianColumn(medianName));
        }

        return medianGenerator;

    }

    /** Manage a maneuver duration.
     * @param startName name of the start driver
     * @param stopName name of the stop driver
     * @param durationName name of the duration driver
     * @return generator for the median driver
     * @since 11.1
     */
    private Duration manageManeuverDuration(final String startName, final String stopName, final String durationName) {

        Duration durationGenerator = null;

        // check if we already have set up the provider
        for (final AdditionalStateProvider provider : getAdditionalStateProviders()) {
            if (provider instanceof Duration &&
                provider.getName().equals(durationName)) {
                // the Jacobian column generator has already been set up in a previous propagation
                durationGenerator = (Duration) provider;
                break;
            }
        }

        if (durationGenerator == null) {
            // this is the first time we need the Jacobian column generator, create it
            durationGenerator = new Duration(startName, stopName, durationName);
            addAdditionalStateProvider(durationGenerator);
        }

        if (!getInitialIntegrationState().hasAdditionalState(durationName)) {
            // add the initial Jacobian column if it is not already there
            // (perhaps due to a previous propagation)
            setInitialColumn(durationName, getHarvester().getInitialJacobianColumn(durationName));
        }

        return durationGenerator;

    }

    /** Set up the Jacobians columns generator for regular parameters.
     * @param stmGenerator generator for the State Transition Matrix
     * @param triggerDates names of the columns already managed as trigger dates
     * @since 11.1
     */
    private void setUpRegularParametersJacobiansColumns(final StateTransitionMatrixGenerator stmGenerator,
                                                        final List<String> triggerDates) {

        // first pass: gather all parameters (excluding trigger dates), binding similar names together
        final ParameterDriversList selected = new ParameterDriversList();
        for (final ForceModel forceModel : getAllForceModels()) {
            for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
                if (!triggerDates.contains(driver.getName())) {
                    selected.add(driver);
                }
            }
        }

        // second pass: now that shared parameter names are bound together,
        // their selections status have been synchronized, we can filter them
        selected.filter(true);

        // third pass: sort parameters lexicographically
        selected.sort();

        // add the Jacobians column generators corresponding to parameters, and setup state accordingly
        for (final DelegatingDriver driver : selected.getDrivers()) {

            IntegrableJacobianColumnGenerator generator = null;

            // check if we already have set up the providers
            for (final AdditionalDerivativesProvider provider : getAdditionalDerivativesProviders()) {
                if (provider instanceof IntegrableJacobianColumnGenerator &&
                    provider.getName().equals(driver.getName())) {
                    // the Jacobian column generator has already been set up in a previous propagation
                    generator = (IntegrableJacobianColumnGenerator) provider;
                    break;
                }
            }

            if (generator == null) {
                // this is the first time we need the Jacobian column generator, create it
                generator = new IntegrableJacobianColumnGenerator(stmGenerator, driver.getName());
                addAdditionalDerivativesProvider(generator);
            }

            if (!getInitialIntegrationState().hasAdditionalState(driver.getName())) {
                // add the initial Jacobian column if it is not already there
                // (perhaps due to a previous propagation)
                setInitialColumn(driver.getName(), getHarvester().getInitialJacobianColumn(driver.getName()));
            }

        }

    }

    /** Add the initial value of the column to the initial state.
     * <p>
     * The initial state must already contain the Cartesian State Transition Matrix.
     * </p>
     * @param columnName name of the column
     * @param dYdQ column of the Jacobian ∂Y/∂qₘ with respect to propagation type,
     * if null (which is the most frequent case), assumed to be 0
     * @since 11.1
     */
    private void setInitialColumn(final String columnName, final double[] dYdQ) {

        final SpacecraftState state = getInitialState();

        if (dYdQ.length != STATE_DIMENSION) {
            throw new OrekitException(LocalizedCoreFormats.DIMENSIONS_MISMATCH,
                                      dYdQ.length, STATE_DIMENSION);
        }

        // convert to Cartesian Jacobian
        final double[][] dYdC = new double[STATE_DIMENSION][STATE_DIMENSION];
        getOrbitType().convertType(state.getOrbit()).getJacobianWrtCartesian(getPositionAngleType(), dYdC);
        final double[] column = new QRDecomposition(MatrixUtils.createRealMatrix(dYdC), THRESHOLD).
                        getSolver().
                        solve(MatrixUtils.createRealVector(dYdQ)).
                        toArray();

        // set additional state
        setInitialState(state.addAdditionalState(columnName, column));

    }

    /** {@inheritDoc} */
    @Override
    public TimeStampedPVCoordinates getPVCoordinates(final AbsoluteDate date, final Frame frame) {
        return propagate(date).getPVCoordinates(frame);
    }

    /** {@inheritDoc} */
    @Override
    protected StateMapper createMapper(final AbsoluteDate referenceDate, final double mu,
                                       final OrbitType orbitType, final PositionAngle positionAngleType,
                                       final AttitudeProvider attitudeProvider, final Frame frame) {
        return new OsculatingMapper(referenceDate, mu, orbitType, positionAngleType, attitudeProvider, frame);
    }

    /** Internal mapper using directly osculating parameters. */
    private static class OsculatingMapper extends StateMapper {

        /** Simple constructor.
         * <p>
         * The position parameter type is meaningful only if {@link
         * #getOrbitType() propagation orbit type}
         * support it. As an example, it is not meaningful for propagation
         * in {@link OrbitType#CARTESIAN Cartesian} parameters.
         * </p>
         * @param referenceDate reference date
         * @param mu central attraction coefficient (m³/s²)
         * @param orbitType orbit type to use for mapping (can be null for {@link AbsolutePVCoordinates})
         * @param positionAngleType angle type to use for propagation
         * @param attitudeProvider attitude provider
         * @param frame inertial frame
         */
        OsculatingMapper(final AbsoluteDate referenceDate, final double mu,
                         final OrbitType orbitType, final PositionAngle positionAngleType,
                         final AttitudeProvider attitudeProvider, final Frame frame) {
            super(referenceDate, mu, orbitType, positionAngleType, attitudeProvider, frame);
        }

        /** {@inheritDoc} */
        public SpacecraftState mapArrayToState(final AbsoluteDate date, final double[] y, final double[] yDot,
                                               final PropagationType type) {
            // the parameter type is ignored for the Numerical Propagator

            final double mass = y[6];
            if (mass <= 0.0) {
                throw new OrekitException(OrekitMessages.SPACECRAFT_MASS_BECOMES_NEGATIVE, mass);
            }

            if (getOrbitType() == null) {
                // propagation uses absolute position-velocity-acceleration
                final Vector3D p = new Vector3D(y[0],    y[1],    y[2]);
                final Vector3D v = new Vector3D(y[3],    y[4],    y[5]);
                final Vector3D a;
                final AbsolutePVCoordinates absPva;
                if (yDot == null) {
                    absPva = new AbsolutePVCoordinates(getFrame(), new TimeStampedPVCoordinates(date, p, v));
                } else {
                    a = new Vector3D(yDot[3], yDot[4], yDot[5]);
                    absPva = new AbsolutePVCoordinates(getFrame(), new TimeStampedPVCoordinates(date, p, v, a));
                }

                final Attitude attitude = getAttitudeProvider().getAttitude(absPva, date, getFrame());
                return new SpacecraftState(absPva, attitude, mass);
            } else {
                // propagation uses regular orbits
                final Orbit orbit       = getOrbitType().mapArrayToOrbit(y, yDot, getPositionAngleType(), date, getMu(), getFrame());
                final Attitude attitude = getAttitudeProvider().getAttitude(orbit, date, getFrame());

                return new SpacecraftState(orbit, attitude, mass);
            }

        }

        /** {@inheritDoc} */
        public void mapStateToArray(final SpacecraftState state, final double[] y, final double[] yDot) {
            if (getOrbitType() == null) {
                // propagation uses absolute position-velocity-acceleration
                final Vector3D p = state.getAbsPVA().getPosition();
                final Vector3D v = state.getAbsPVA().getVelocity();
                y[0] = p.getX();
                y[1] = p.getY();
                y[2] = p.getZ();
                y[3] = v.getX();
                y[4] = v.getY();
                y[5] = v.getZ();
                y[6] = state.getMass();
            }
            else {
                getOrbitType().mapOrbitToArray(state.getOrbit(), getPositionAngleType(), y, yDot);
                y[6] = state.getMass();
            }
        }

    }

    /** {@inheritDoc} */
    protected MainStateEquations getMainStateEquations(final ODEIntegrator integrator) {
        return new Main(integrator);
    }

    /** Internal class for osculating parameters integration. */
    private class Main implements MainStateEquations, TimeDerivativesEquations {

        /** Derivatives array. */
        private final double[] yDot;

        /** Current state. */
        private SpacecraftState currentState;

        /** Jacobian of the orbital parameters with respect to the Cartesian parameters. */
        private double[][] jacobian;

        /** Simple constructor.
         * @param integrator numerical integrator to use for propagation.
         */
        Main(final ODEIntegrator integrator) {

            this.yDot     = new double[7];
            this.jacobian = new double[6][6];

            for (final ForceModel forceModel : forceModels) {
                forceModel.getEventsDetectors().forEach(detector -> setUpEventDetector(integrator, detector));
            }

            if (getOrbitType() == null) {
                // propagation uses absolute position-velocity-acceleration
                // we can set Jacobian once and for all
                for (int i = 0; i < jacobian.length; ++i) {
                    Arrays.fill(jacobian[i], 0.0);
                    jacobian[i][i] = 1.0;
                }
            }

        }

        /** {@inheritDoc} */
        @Override
        public void init(final SpacecraftState initialState, final AbsoluteDate target) {
            forceModels.forEach(fm -> fm.init(initialState, target));
        }

        /** {@inheritDoc} */
        @Override
        public double[] computeDerivatives(final SpacecraftState state) {

            currentState = state;
            Arrays.fill(yDot, 0.0);
            if (getOrbitType() != null) {
                // propagation uses regular orbits
                currentState.getOrbit().getJacobianWrtCartesian(getPositionAngleType(), jacobian);
            }

            // compute the contributions of all perturbing forces,
            // using the Kepler contribution at the end since
            // NewtonianAttraction is always the last instance in the list
            for (final ForceModel forceModel : forceModels) {
                forceModel.addContribution(state, this);
            }

            if (getOrbitType() == null) {
                // position derivative is velocity, and was not added above in the force models
                // (it is added when orbit type is non-null because NewtonianAttraction considers it)
                final Vector3D velocity = currentState.getPVCoordinates().getVelocity();
                yDot[0] += velocity.getX();
                yDot[1] += velocity.getY();
                yDot[2] += velocity.getZ();
            }


            return yDot.clone();

        }

        /** {@inheritDoc} */
        @Override
        public void addKeplerContribution(final double mu) {
            if (getOrbitType() == null) {

                // if mu is neither 0 nor NaN, we want to include Newtonian acceleration
                if (mu > 0) {
                    // velocity derivative is Newtonian acceleration
                    final Vector3D position = currentState.getPVCoordinates().getPosition();
                    final double r2         = position.getNormSq();
                    final double coeff      = -mu / (r2 * FastMath.sqrt(r2));
                    yDot[3] += coeff * position.getX();
                    yDot[4] += coeff * position.getY();
                    yDot[5] += coeff * position.getZ();
                }

            } else {
                // propagation uses regular orbits
                currentState.getOrbit().addKeplerContribution(getPositionAngleType(), mu, yDot);
            }
        }

        /** {@inheritDoc} */
        public void addNonKeplerianAcceleration(final Vector3D gamma) {
            for (int i = 0; i < 6; ++i) {
                final double[] jRow = jacobian[i];
                yDot[i] += jRow[3] * gamma.getX() + jRow[4] * gamma.getY() + jRow[5] * gamma.getZ();
            }
        }

        /** {@inheritDoc} */
        @Override
        public void addMassDerivative(final double q) {
            if (q > 0) {
                throw new OrekitIllegalArgumentException(OrekitMessages.POSITIVE_FLOW_RATE, q);
            }
            yDot[6] += q;
        }

    }

    /** Estimate tolerance vectors for integrators when propagating in absolute position-velocity-acceleration.
     * @param dP user specified position error
     * @param absPva reference absolute position-velocity-acceleration
     * @return a two rows array, row 0 being the absolute tolerance error and row 1
     * being the relative tolerance error
     * @see NumericalPropagator#tolerances(double, Orbit, OrbitType)
     */
    public static double[][] tolerances(final double dP, final AbsolutePVCoordinates absPva) {

        final double relative = dP / absPva.getPosition().getNorm();
        final double dV = relative * absPva.getVelocity().getNorm();

        final double[] absTol = new double[7];
        final double[] relTol = new double[7];

        absTol[0] = dP;
        absTol[1] = dP;
        absTol[2] = dP;
        absTol[3] = dV;
        absTol[4] = dV;
        absTol[5] = dV;

        // we set the mass tolerance arbitrarily to 1.0e-6 kg, as mass evolves linearly
        // with trust, this often has no influence at all on propagation
        absTol[6] = 1.0e-6;

        Arrays.fill(relTol, relative);

        return new double[][] {
            absTol, relTol
        };

    }

    /** Estimate tolerance vectors for integrators when propagating in orbits.
     * <p>
     * The errors are estimated from partial derivatives properties of orbits,
     * starting from a scalar position error specified by the user.
     * Considering the energy conservation equation V = sqrt(mu (2/r - 1/a)),
     * we get at constant energy (i.e. on a Keplerian trajectory):
     * <pre>
     * V r² |dV| = mu |dr|
     * </pre>
     * <p> So we deduce a scalar velocity error consistent with the position error.
     * From here, we apply orbits Jacobians matrices to get consistent errors
     * on orbital parameters.
     *
     * <p>
     * The tolerances are only <em>orders of magnitude</em>, and integrator tolerances
     * are only local estimates, not global ones. So some care must be taken when using
     * these tolerances. Setting 1mm as a position error does NOT mean the tolerances
     * will guarantee a 1mm error position after several orbits integration.
     * </p>
     * @param dP user specified position error
     * @param orbit reference orbit
     * @param type propagation type for the meaning of the tolerance vectors elements
     * (it may be different from {@code orbit.getType()})
     * @return a two rows array, row 0 being the absolute tolerance error and row 1
     * being the relative tolerance error
     */
    public static double[][] tolerances(final double dP, final Orbit orbit, final OrbitType type) {

        // estimate the scalar velocity error
        final PVCoordinates pv = orbit.getPVCoordinates();
        final double r2 = pv.getPosition().getNormSq();
        final double v  = pv.getVelocity().getNorm();
        final double dV = orbit.getMu() * dP / (v * r2);

        return tolerances(dP, dV, orbit, type);

    }

    /** Estimate tolerance vectors for integrators when propagating in orbits.
     * <p>
     * The errors are estimated from partial derivatives properties of orbits,
     * starting from scalar position and velocity errors specified by the user.
     * <p>
     * The tolerances are only <em>orders of magnitude</em>, and integrator tolerances
     * are only local estimates, not global ones. So some care must be taken when using
     * these tolerances. Setting 1mm as a position error does NOT mean the tolerances
     * will guarantee a 1mm error position after several orbits integration.
     * </p>
     * @param dP user specified position error
     * @param dV user specified velocity error
     * @param orbit reference orbit
     * @param type propagation type for the meaning of the tolerance vectors elements
     * (it may be different from {@code orbit.getType()})
     * @return a two rows array, row 0 being the absolute tolerance error and row 1
     * being the relative tolerance error
     * @since 10.3
     */
    public static double[][] tolerances(final double dP, final double dV,
                                        final Orbit orbit, final OrbitType type) {

        final double[] absTol = new double[7];
        final double[] relTol = new double[7];

        // we set the mass tolerance arbitrarily to 1.0e-6 kg, as mass evolves linearly
        // with trust, this often has no influence at all on propagation
        absTol[6] = 1.0e-6;

        if (type == OrbitType.CARTESIAN) {
            absTol[0] = dP;
            absTol[1] = dP;
            absTol[2] = dP;
            absTol[3] = dV;
            absTol[4] = dV;
            absTol[5] = dV;
        } else {

            // convert the orbit to the desired type
            final double[][] jacobian = new double[6][6];
            final Orbit converted = type.convertType(orbit);
            converted.getJacobianWrtCartesian(PositionAngle.TRUE, jacobian);

            for (int i = 0; i < 6; ++i) {
                final double[] row = jacobian[i];
                absTol[i] = FastMath.abs(row[0]) * dP +
                            FastMath.abs(row[1]) * dP +
                            FastMath.abs(row[2]) * dP +
                            FastMath.abs(row[3]) * dV +
                            FastMath.abs(row[4]) * dV +
                            FastMath.abs(row[5]) * dV;
                if (Double.isNaN(absTol[i])) {
                    throw new OrekitException(OrekitMessages.SINGULAR_JACOBIAN_FOR_ORBIT_TYPE, type);
                }
            }

        }

        Arrays.fill(relTol, dP / FastMath.sqrt(orbit.getPVCoordinates().getPosition().getNormSq()));

        return new double[][] {
            absTol, relTol
        };

    }

    /** {@inheritDoc} */
    @Override
    protected void beforeIntegration(final SpacecraftState initialState, final AbsoluteDate tEnd) {

        if (!getFrame().isPseudoInertial()) {

            // inspect all force models to find InertialForces
            for (ForceModel force : forceModels) {
                if (force instanceof InertialForces) {
                    return;
                }
            }

            // throw exception if no inertial forces found
            throw new OrekitException(OrekitMessages.INERTIAL_FORCE_MODEL_MISSING, getFrame().getName());

        }

    }

}