NumericalPropagator.java

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

  18. import org.hipparchus.exception.LocalizedCoreFormats;
  19. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  20. import org.hipparchus.linear.DecompositionSolver;
  21. import org.hipparchus.linear.MatrixUtils;
  22. import org.hipparchus.linear.QRDecomposition;
  23. import org.hipparchus.linear.RealMatrix;
  24. import org.hipparchus.ode.ODEIntegrator;
  25. import org.hipparchus.util.FastMath;
  26. import org.hipparchus.util.Precision;
  27. import org.orekit.annotation.DefaultDataContext;
  28. import org.orekit.attitudes.Attitude;
  29. import org.orekit.attitudes.AttitudeProvider;
  30. import org.orekit.data.DataContext;
  31. import org.orekit.errors.OrekitException;
  32. import org.orekit.errors.OrekitIllegalArgumentException;
  33. import org.orekit.errors.OrekitInternalError;
  34. import org.orekit.errors.OrekitMessages;
  35. import org.orekit.forces.ForceModel;
  36. import org.orekit.forces.drag.AbstractDragForceModel;
  37. import org.orekit.forces.gravity.NewtonianAttraction;
  38. import org.orekit.forces.inertia.InertialForces;
  39. import org.orekit.forces.maneuvers.Maneuver;
  40. import org.orekit.forces.maneuvers.jacobians.Duration;
  41. import org.orekit.forces.maneuvers.jacobians.MassDepletionDelay;
  42. import org.orekit.forces.maneuvers.jacobians.MedianDate;
  43. import org.orekit.forces.maneuvers.jacobians.TriggerDate;
  44. import org.orekit.forces.maneuvers.trigger.ManeuverTriggerDetector;
  45. import org.orekit.forces.maneuvers.trigger.ResettableManeuverTriggers;
  46. import org.orekit.forces.radiation.RadiationForceModel;
  47. import org.orekit.frames.Frame;
  48. import org.orekit.orbits.Orbit;
  49. import org.orekit.orbits.OrbitType;
  50. import org.orekit.orbits.PositionAngleType;
  51. import org.orekit.propagation.AbstractMatricesHarvester;
  52. import org.orekit.propagation.AdditionalDataProvider;
  53. import org.orekit.propagation.CartesianToleranceProvider;
  54. import org.orekit.propagation.MatricesHarvester;
  55. import org.orekit.propagation.PropagationType;
  56. import org.orekit.propagation.Propagator;
  57. import org.orekit.propagation.SpacecraftState;
  58. import org.orekit.propagation.ToleranceProvider;
  59. import org.orekit.propagation.events.EventDetector;
  60. import org.orekit.propagation.events.ParameterDrivenDateIntervalDetector;
  61. import org.orekit.propagation.integration.AbstractIntegratedPropagator;
  62. import org.orekit.propagation.integration.AdditionalDerivativesProvider;
  63. import org.orekit.propagation.integration.StateMapper;
  64. import org.orekit.time.AbsoluteDate;
  65. import org.orekit.utils.AbsolutePVCoordinates;
  66. import org.orekit.utils.DoubleArrayDictionary;
  67. import org.orekit.utils.ParameterDriver;
  68. import org.orekit.utils.ParameterDriversList;
  69. import org.orekit.utils.ParameterDriversList.DelegatingDriver;
  70. import org.orekit.utils.ParameterObserver;
  71. import org.orekit.utils.TimeSpanMap;
  72. import org.orekit.utils.TimeSpanMap.Span;
  73. import org.orekit.utils.TimeStampedPVCoordinates;

  74. import java.util.ArrayList;
  75. import java.util.Arrays;
  76. import java.util.Collection;
  77. import java.util.Collections;
  78. import java.util.List;
  79. import java.util.Optional;
  80. import java.util.stream.Collectors;

  81. /** This class propagates {@link org.orekit.orbits.Orbit orbits} using
  82.  * numerical integration.
  83.  * <p>Numerical propagation is much more accurate than analytical propagation
  84.  * like for example {@link org.orekit.propagation.analytical.KeplerianPropagator
  85.  * Keplerian} or {@link org.orekit.propagation.analytical.EcksteinHechlerPropagator
  86.  * Eckstein-Hechler}, but requires a few more steps to set up to be used properly.
  87.  * Whereas analytical propagators are configured only thanks to their various
  88.  * constructors and can be used immediately after construction, numerical propagators
  89.  * configuration involve setting several parameters between construction time
  90.  * and propagation time.</p>
  91.  * <p>The configuration parameters that can be set are:</p>
  92.  * <ul>
  93.  *   <li>the initial spacecraft state ({@link #setInitialState(SpacecraftState)})</li>
  94.  *   <li>the central attraction coefficient ({@link #setMu(double)})</li>
  95.  *   <li>the various force models ({@link #addForceModel(ForceModel)},
  96.  *   {@link #removeForceModels()})</li>
  97.  *   <li>the {@link OrbitType type} of orbital parameters to be used for propagation
  98.  *   ({@link #setOrbitType(OrbitType)}),</li>
  99.  *   <li>the {@link PositionAngleType type} of position angle to be used in orbital parameters
  100.  *   to be used for propagation where it is relevant ({@link
  101.  *   #setPositionAngleType(PositionAngleType)}),</li>
  102.  *   <li>whether {@link MatricesHarvester state transition matrices and Jacobians matrices}
  103.  *   (with the option to include mass if a 7x7 initial matrix is passed) should be propagated along with orbital state
  104.  *   ({@link #setupMatricesComputation(String, RealMatrix, DoubleArrayDictionary)}),</li>
  105.  *   <li>whether {@link org.orekit.propagation.integration.AdditionalDerivativesProvider additional derivatives}
  106.  *   should be propagated along with orbital state ({@link
  107.  *   #addAdditionalDerivativesProvider(AdditionalDerivativesProvider)}),</li>
  108.  *   <li>the discrete events that should be triggered during propagation
  109.  *   ({@link #addEventDetector(EventDetector)},
  110.  *   {@link #clearEventsDetectors()})</li>
  111.  *   <li>the binding logic with the rest of the application ({@link #getMultiplexer()})</li>
  112.  * </ul>
  113.  * <p>From these configuration parameters, only the initial state is mandatory. The default
  114.  * propagation settings are in {@link OrbitType#EQUINOCTIAL equinoctial} parameters with
  115.  * {@link PositionAngleType#ECCENTRIC} longitude argument. If the central attraction coefficient
  116.  * is not explicitly specified, the one used to define the initial orbit will be used.
  117.  * However, specifying only the initial state and perhaps the central attraction coefficient
  118.  * would mean the propagator would use only Keplerian forces. In this case, the simpler {@link
  119.  * org.orekit.propagation.analytical.KeplerianPropagator KeplerianPropagator} class would
  120.  * perhaps be more effective.</p>
  121.  * <p>The underlying numerical integrator set up in the constructor may also have its own
  122.  * configuration parameters. Typical configuration parameters for adaptive stepsize integrators
  123.  * are the min, max and perhaps start step size as well as the absolute and/or relative errors
  124.  * thresholds.</p>
  125.  * <p>The state that is seen by the integrator is a simple seven elements double array.
  126.  * The six first elements are either:
  127.  * <ul>
  128.  *   <li>the {@link org.orekit.orbits.EquinoctialOrbit equinoctial orbit parameters} (a, e<sub>x</sub>,
  129.  *   e<sub>y</sub>, h<sub>x</sub>, h<sub>y</sub>, λ<sub>M</sub> or λ<sub>E</sub>
  130.  *   or λ<sub>v</sub>) in meters and radians,</li>
  131.  *   <li>the {@link org.orekit.orbits.KeplerianOrbit Keplerian orbit parameters} (a, e, i, ω, Ω,
  132.  *   M or E or v) in meters and radians,</li>
  133.  *   <li>the {@link org.orekit.orbits.CircularOrbit circular orbit parameters} (a, e<sub>x</sub>, e<sub>y</sub>, i,
  134.  *   Ω, α<sub>M</sub> or α<sub>E</sub> or α<sub>v</sub>) in meters
  135.  *   and radians,</li>
  136.  *   <li>the {@link org.orekit.orbits.CartesianOrbit Cartesian orbit parameters} (x, y, z, v<sub>x</sub>,
  137.  *   v<sub>y</sub>, v<sub>z</sub>) in meters and meters per seconds.
  138.  * </ul>
  139.  * <p> The last element is the mass in kilograms and changes only during thrusters firings
  140.  *
  141.  * <p>The following code snippet shows a typical setting for Low Earth Orbit propagation in
  142.  * equinoctial parameters and true longitude argument:</p>
  143.  * <pre>
  144.  * final double dP       = 0.001;
  145.  * final double minStep  = 0.001;
  146.  * final double maxStep  = 500;
  147.  * final double initStep = 60;
  148.  * final double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(dP).getTolerances(orbit, OrbitType.EQUINOCTIAL);
  149.  * AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tolerance[0], tolerance[1]);
  150.  * integrator.setInitialStepSize(initStep);
  151.  * propagator = new NumericalPropagator(integrator);
  152.  * </pre>
  153.  * <p>By default, at the end of the propagation, the propagator resets the initial state to the final state,
  154.  * thus allowing a new propagation to be started from there without recomputing the part already performed.
  155.  * This behaviour can be changed by calling {@link #setResetAtEnd(boolean)}.
  156.  * </p>
  157.  * <p>Beware the same instance cannot be used simultaneously by different threads, the class is <em>not</em>
  158.  * thread-safe.</p>
  159.  *
  160.  * @see SpacecraftState
  161.  * @see ForceModel
  162.  * @see org.orekit.propagation.sampling.OrekitStepHandler
  163.  * @see org.orekit.propagation.sampling.OrekitFixedStepHandler
  164.  * @see org.orekit.propagation.integration.IntegratedEphemeris
  165.  * @see TimeDerivativesEquations
  166.  *
  167.  * @author Mathieu Rom&eacute;ro
  168.  * @author Luc Maisonobe
  169.  * @author Guylaine Prat
  170.  * @author Fabien Maussion
  171.  * @author V&eacute;ronique Pommier-Maurussane
  172.  */
  173. public class NumericalPropagator extends AbstractIntegratedPropagator {

  174.     /** Default orbit type. */
  175.     public static final OrbitType DEFAULT_ORBIT_TYPE = OrbitType.EQUINOCTIAL;

  176.     /** Default position angle type. */
  177.     public static final PositionAngleType DEFAULT_POSITION_ANGLE_TYPE = PositionAngleType.ECCENTRIC;

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

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

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

  184.     /**
  185.      * boolean to know if a full attitude (with rates) is needed when computing derivatives for the ODE.
  186.      * since 12.1
  187.      */
  188.     private boolean needFullAttitudeForDerivatives = true;

  189.     /** Create a new instance of NumericalPropagator, based on orbit definition mu.
  190.      * After creation, the instance is empty, i.e. the attitude provider is set to an
  191.      * unspecified default law and there are no perturbing forces at all.
  192.      * This means that if {@link #addForceModel addForceModel} is not
  193.      * called after creation, the integrated orbit will follow a Keplerian
  194.      * evolution only. The defaults are {@link OrbitType#EQUINOCTIAL}
  195.      * for {@link #setOrbitType(OrbitType) propagation
  196.      * orbit type} and {@link PositionAngleType#ECCENTRIC} for {@link
  197.      * #setPositionAngleType(PositionAngleType) position angle type}.
  198.      *
  199.      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
  200.      *
  201.      * @param integrator numerical integrator to use for propagation.
  202.      * @see #NumericalPropagator(ODEIntegrator, AttitudeProvider)
  203.      */
  204.     @DefaultDataContext
  205.     public NumericalPropagator(final ODEIntegrator integrator) {
  206.         this(integrator,
  207.                 Propagator.getDefaultLaw(DataContext.getDefault().getFrames()));
  208.     }

  209.     /** Create a new instance of NumericalPropagator, based on orbit definition mu.
  210.      * After creation, the instance is empty, i.e. the attitude provider is set to an
  211.      * unspecified default law and there are no perturbing forces at all.
  212.      * This means that if {@link #addForceModel addForceModel} is not
  213.      * called after creation, the integrated orbit will follow a Keplerian
  214.      * evolution only. The defaults are {@link OrbitType#EQUINOCTIAL}
  215.      * for {@link #setOrbitType(OrbitType) propagation
  216.      * orbit type} and {@link PositionAngleType#ECCENTRIC} for {@link
  217.      * #setPositionAngleType(PositionAngleType) position angle type}.
  218.      * @param integrator numerical integrator to use for propagation.
  219.      * @param attitudeProvider the attitude law.
  220.      * @since 10.1
  221.      */
  222.     public NumericalPropagator(final ODEIntegrator integrator,
  223.                                final AttitudeProvider attitudeProvider) {
  224.         super(integrator, PropagationType.OSCULATING);
  225.         forceModels             = new ArrayList<>();
  226.         ignoreCentralAttraction = false;
  227.         initMapper();
  228.         setAttitudeProvider(attitudeProvider);
  229.         clearStepHandlers();
  230.         setOrbitType(DEFAULT_ORBIT_TYPE);
  231.         setPositionAngleType(DEFAULT_POSITION_ANGLE_TYPE);
  232.     }

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

  240.      /** Set the central attraction coefficient μ.
  241.       * <p>
  242.       * Setting the central attraction coefficient is
  243.       * equivalent to {@link #addForceModel(ForceModel) add}
  244.       * a {@link NewtonianAttraction} force model.
  245.       * * </p>
  246.       * @param mu central attraction coefficient (m³/s²)
  247.       * @see #addForceModel(ForceModel)
  248.       * @see #getAllForceModels()
  249.       */
  250.     @Override
  251.     public void setMu(final double mu) {
  252.         if (ignoreCentralAttraction) {
  253.             superSetMu(mu);
  254.         } else {
  255.             addForceModel(new NewtonianAttraction(mu));
  256.             superSetMu(mu);
  257.         }
  258.     }

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

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

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

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

  286.             try {
  287.                 // ensure we are notified of any mu change
  288.                 model.getParametersDrivers().get(0).addObserver(new ParameterObserver() {
  289.                     /** {@inheritDoc} */
  290.                     @Override
  291.                     public void valueSpanMapChanged(final TimeSpanMap<Double> previousValue, final ParameterDriver driver) {
  292.                         superSetMu(driver.getValue());
  293.                     }
  294.                     /** {@inheritDoc} */
  295.                     @Override
  296.                     public void valueChanged(final double previousValue, final ParameterDriver driver, final AbsoluteDate date) {
  297.                         superSetMu(driver.getValue());
  298.                     }
  299.                 });
  300.             } catch (OrekitException oe) {
  301.                 // this should never happen
  302.                 throw new OrekitInternalError(oe);
  303.             }

  304.             if (hasNewtonianAttraction()) {
  305.                 // there is already a central attraction model, replace it
  306.                 forceModels.set(forceModels.size() - 1, model);
  307.             } else {
  308.                 // there are no central attraction model yet, add it at the end of the list
  309.                 forceModels.add(model);
  310.             }
  311.         } else {
  312.             // we want to add a perturbing force model
  313.             if (hasNewtonianAttraction()) {
  314.                 // insert the new force model before Newtonian attraction,
  315.                 // which should always be the last one in the list
  316.                 forceModels.add(forceModels.size() - 1, model);
  317.             } else {
  318.                 // we only have perturbing force models up to now, just append at the end of the list
  319.                 forceModels.add(model);
  320.             }
  321.         }

  322.     }

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

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

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

  357.     /** Get propagation parameter type.
  358.      * @return orbit type used for propagation, null for
  359.      * propagating using {@link org.orekit.utils.AbsolutePVCoordinates} rather than {@link Orbit}
  360.      */
  361.     @Override
  362.     public OrbitType getOrbitType() {
  363.         return super.getOrbitType();
  364.     }

  365.     /** Set position angle type.
  366.      * <p>
  367.      * The position parameter type is meaningful only if {@link
  368.      * #getOrbitType() propagation orbit type}
  369.      * support it. As an example, it is not meaningful for propagation
  370.      * in {@link OrbitType#CARTESIAN Cartesian} parameters.
  371.      * </p>
  372.      * @param positionAngleType angle type to use for propagation
  373.      */
  374.     @Override
  375.     public void setPositionAngleType(final PositionAngleType positionAngleType) {
  376.         super.setPositionAngleType(positionAngleType);
  377.     }

  378.     /** Get propagation parameter type.
  379.      * @return angle type to use for propagation
  380.      */
  381.     @Override
  382.     public PositionAngleType getPositionAngleType() {
  383.         return super.getPositionAngleType();
  384.     }

  385.     /** Set the initial state.
  386.      * @param initialState initial state
  387.      */
  388.     public void setInitialState(final SpacecraftState initialState) {
  389.         resetInitialState(initialState);
  390.     }

  391.     /** {@inheritDoc} */
  392.     @Override
  393.     public void resetInitialState(final SpacecraftState state) {
  394.         super.resetInitialState(state);
  395.         if (!hasNewtonianAttraction()) {
  396.             // use the state to define central attraction
  397.             setMu(state.isOrbitDefined() ? state.getOrbit().getMu() : Double.NaN);
  398.         }
  399.         setStartDate(state.getDate());
  400.     }

  401.     /** Get the names of the parameters in the matrix returned by {@link MatricesHarvester#getParametersJacobian}.
  402.      * @return names of the parameters (i.e. columns) of the Jacobian matrix
  403.      */
  404.     List<String> getJacobiansColumnsNames() {
  405.         final List<String> columnsNames = new ArrayList<>();
  406.         for (final ForceModel forceModel : getAllForceModels()) {
  407.             for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
  408.                 if (driver.isSelected() && !columnsNames.contains(driver.getNamesSpanMap().getFirstSpan().getData())) {
  409.                     // As driver with same name should have same NamesSpanMap we only check if the first span is present,
  410.                     // if not we add all span names to columnsNames
  411.                     for (Span<String> span = driver.getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {
  412.                         columnsNames.add(span.getData());
  413.                     }
  414.                 }
  415.             }
  416.         }
  417.         Collections.sort(columnsNames);
  418.         return columnsNames;
  419.     }

  420.     /** {@inheritDoc}
  421.      * <p>
  422.      * Unlike other propagators, the numerical one can consider the mass as a state variable in the transition matrix.
  423.      * To do so, a 7x7 initial matrix is to be passed instead of 6x6.
  424.      * </p>
  425.      * */
  426.     @Override
  427.     protected AbstractMatricesHarvester createHarvester(final String stmName, final RealMatrix initialStm,
  428.                                                         final DoubleArrayDictionary initialJacobianColumns) {
  429.         return new NumericalPropagationHarvester(this, stmName, initialStm, initialJacobianColumns);
  430.     }

  431.     /** {@inheritDoc} */
  432.     @Override
  433.     public void clearMatricesComputation() {
  434.         final List<AdditionalDerivativesProvider> copiedDerivativesProviders = new ArrayList<>(getAdditionalDerivativesProviders());
  435.         copiedDerivativesProviders.stream().filter(AbstractStateTransitionMatrixGenerator.class::isInstance)
  436.                 .forEach(provider -> removeAdditionalDerivativesProvider(provider.getName()));
  437.         final List<AdditionalDataProvider<?>> copiedDataProviders = new ArrayList<>(getAdditionalDataProviders());
  438.         for (final AdditionalDataProvider<?> additionalDataProvider: copiedDataProviders) {
  439.             if (additionalDataProvider instanceof TriggerDate) {
  440.                 final TriggerDate triggerDate = (TriggerDate) additionalDataProvider;
  441.                 if (triggerDate.getMassDepletionDelay() != null) {
  442.                     removeAdditionalDerivativesProvider(triggerDate.getMassDepletionDelay().getName());
  443.                 }
  444.                 removeAdditionalDataProvider(additionalDataProvider.getName());
  445.             } else if (additionalDataProvider instanceof MedianDate || additionalDataProvider instanceof Duration) {
  446.                 removeAdditionalDataProvider(additionalDataProvider.getName());
  447.             }
  448.         }
  449.         super.clearMatricesComputation();
  450.     }

  451.     /** {@inheritDoc} */
  452.     @Override
  453.     protected void setUpStmAndJacobianGenerators() {

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

  456.             // set up the additional equations and additional state providers
  457.             final AbstractStateTransitionMatrixGenerator stmGenerator = setUpStmGenerator();
  458.             final List<String> triggersDates = setUpTriggerDatesJacobiansColumns(stmGenerator);
  459.             setUpRegularParametersJacobiansColumns(stmGenerator, triggersDates);

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

  463.         }

  464.     }

  465.     /** Set up the State Transition Matrix Generator.
  466.      * @return State Transition Matrix Generator
  467.      * @since 11.1
  468.      */
  469.     private AbstractStateTransitionMatrixGenerator setUpStmGenerator() {

  470.         final AbstractMatricesHarvester harvester = getHarvester();

  471.         // add the STM generator corresponding to the current settings, and setup state accordingly
  472.         AbstractStateTransitionMatrixGenerator stmGenerator = null;
  473.         for (final AdditionalDerivativesProvider equations : getAdditionalDerivativesProviders()) {
  474.             if (equations instanceof AbstractStateTransitionMatrixGenerator &&
  475.                 equations.getName().equals(harvester.getStmName())) {
  476.                 // the STM generator has already been set up in a previous propagation
  477.                 stmGenerator = (AbstractStateTransitionMatrixGenerator) equations;
  478.                 break;
  479.             }
  480.         }
  481.         if (stmGenerator == null) {
  482.             // this is the first time we need the STM generate, create it
  483.             if (harvester.getStateDimension() > 6) {
  484.                 stmGenerator = new ExtendedStateTransitionMatrixGenerator(harvester.getStmName(), getAllForceModels(),
  485.                         getAttitudeProvider());
  486.             } else {
  487.                 stmGenerator = new StateTransitionMatrixGenerator(harvester.getStmName(), getAllForceModels(),
  488.                         getAttitudeProvider());
  489.             }
  490.             addAdditionalDerivativesProvider(stmGenerator);
  491.         }

  492.         if (!getInitialIntegrationState().hasAdditionalData(harvester.getStmName())) {
  493.             // add the initial State Transition Matrix if it is not already there
  494.             // (perhaps due to a previous propagation)
  495.             setInitialState(stmGenerator.setInitialStateTransitionMatrix(getInitialState(),
  496.                                                                          harvester.getInitialStateTransitionMatrix(),
  497.                                                                          getOrbitType(),
  498.                                                                          getPositionAngleType()));
  499.         }

  500.         return stmGenerator;

  501.     }

  502.     /** Set up the Jacobians columns generator dedicated to trigger dates.
  503.      * @param stmGenerator State Transition Matrix generator
  504.      * @return names of the columns corresponding to trigger dates
  505.      * @since 13.1
  506.      */
  507.     private List<String> setUpTriggerDatesJacobiansColumns(final AbstractStateTransitionMatrixGenerator stmGenerator) {

  508.         final String stmName = stmGenerator.getName();
  509.         final boolean isMassInStm = stmGenerator instanceof ExtendedStateTransitionMatrixGenerator;
  510.         final List<String> names = new ArrayList<>();
  511.         for (final ForceModel forceModel : getAllForceModels()) {
  512.             if (forceModel instanceof Maneuver && ((Maneuver) forceModel).getManeuverTriggers() instanceof ResettableManeuverTriggers) {
  513.                 final Maneuver maneuver = (Maneuver) forceModel;
  514.                 final ResettableManeuverTriggers maneuverTriggers = (ResettableManeuverTriggers) maneuver.getManeuverTriggers();

  515.                 final Collection<EventDetector> selectedDetectors = maneuverTriggers.getEventDetectors().
  516.                         filter(ManeuverTriggerDetector.class::isInstance).
  517.                         map(triggerDetector -> ((ManeuverTriggerDetector<?>) triggerDetector).getDetector())
  518.                         .collect(Collectors.toList());
  519.                 for (final EventDetector detector: selectedDetectors) {
  520.                     if (detector instanceof ParameterDrivenDateIntervalDetector) {
  521.                         final ParameterDrivenDateIntervalDetector d = (ParameterDrivenDateIntervalDetector) detector;
  522.                         TriggerDate start;
  523.                         TriggerDate stop;

  524.                         if (d.getStartDriver().isSelected() || d.getMedianDriver().isSelected() || d.getDurationDriver().isSelected()) {
  525.                             // normally datedriver should have only 1 span but just in case the user defines several span, there will
  526.                             // be no problem here
  527.                             for (Span<String> span = d.getStartDriver().getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {
  528.                                 start = manageTriggerDate(stmName, maneuver, maneuverTriggers, span.getData(), true,
  529.                                         d.getThreshold(), isMassInStm);
  530.                                 names.add(start.getName());
  531.                                 start = null;
  532.                             }
  533.                         }
  534.                         if (d.getStopDriver().isSelected() || d.getMedianDriver().isSelected() || d.getDurationDriver().isSelected()) {
  535.                             // normally datedriver should have only 1 span but just in case the user defines several span, there will
  536.                             // be no problem here
  537.                             for (Span<String> span = d.getStopDriver().getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {
  538.                                 stop = manageTriggerDate(stmName, maneuver, maneuverTriggers, span.getData(), false,
  539.                                         d.getThreshold(), isMassInStm);
  540.                                 names.add(stop.getName());
  541.                                 stop = null;
  542.                             }
  543.                         }
  544.                         if (d.getMedianDriver().isSelected()) {
  545.                             // for first span
  546.                             Span<String> currentMedianNameSpan = d.getMedianDriver().getNamesSpanMap().getFirstSpan();
  547.                             MedianDate median =
  548.                                     manageMedianDate(d.getStartDriver().getNamesSpanMap().getFirstSpan().getData(),
  549.                                             d.getStopDriver().getNamesSpanMap().getFirstSpan().getData(), currentMedianNameSpan.getData());
  550.                             names.add(median.getName());
  551.                             // for all span
  552.                             // normally datedriver should have only 1 span but just in case the user defines several span, there will
  553.                             // be no problem here. /!\ medianDate driver, startDate driver and stopDate driver must have same span number
  554.                             for (int spanNumber = 1; spanNumber < d.getMedianDriver().getNamesSpanMap().getSpansNumber(); ++spanNumber) {
  555.                                 currentMedianNameSpan = d.getMedianDriver().getNamesSpanMap().getSpan(currentMedianNameSpan.getEnd());
  556.                                 median =
  557.                                         manageMedianDate(d.getStartDriver().getNamesSpanMap().getSpan(currentMedianNameSpan.getStart()).getData(),
  558.                                                 d.getStopDriver().getNamesSpanMap().getSpan(currentMedianNameSpan.getStart()).getData(),
  559.                                                 currentMedianNameSpan.getData());
  560.                                 names.add(median.getName());

  561.                             }

  562.                         }
  563.                         if (d.getDurationDriver().isSelected()) {
  564.                             // for first span
  565.                             Span<String> currentDurationNameSpan = d.getDurationDriver().getNamesSpanMap().getFirstSpan();
  566.                             Duration duration =
  567.                                     manageManeuverDuration(d.getStartDriver().getNamesSpanMap().getFirstSpan().getData(),
  568.                                             d.getStopDriver().getNamesSpanMap().getFirstSpan().getData(), currentDurationNameSpan.getData());
  569.                             names.add(duration.getName());
  570.                             // for all span
  571.                             for (int spanNumber = 1; spanNumber < d.getDurationDriver().getNamesSpanMap().getSpansNumber(); ++spanNumber) {
  572.                                 currentDurationNameSpan = d.getDurationDriver().getNamesSpanMap().getSpan(currentDurationNameSpan.getEnd());
  573.                                 duration =
  574.                                         manageManeuverDuration(d.getStartDriver().getNamesSpanMap().getSpan(currentDurationNameSpan.getStart()).getData(),
  575.                                                 d.getStopDriver().getNamesSpanMap().getSpan(currentDurationNameSpan.getStart()).getData(),
  576.                                                 currentDurationNameSpan.getData());
  577.                                 names.add(duration.getName());

  578.                             }
  579.                         }
  580.                     }
  581.                 }
  582.             }
  583.         }

  584.         return names;

  585.     }

  586.     /** Manage a maneuver trigger date.
  587.      * @param stmName name of the State Transition Matrix state
  588.      * @param maneuver maneuver force model
  589.      * @param mt trigger to which the driver is bound
  590.      * @param driverName name of the date driver
  591.      * @param start if true, the driver is a maneuver start
  592.      * @param threshold event detector threshold
  593.      * @param isMassInStm flag on presence on mass in STM
  594.      * @return generator for the date driver
  595.      * @since 13.1
  596.      */
  597.     private TriggerDate manageTriggerDate(final String stmName,
  598.                                           final Maneuver maneuver,
  599.                                           final ResettableManeuverTriggers mt,
  600.                                           final String driverName,
  601.                                           final boolean start,
  602.                                           final double threshold,
  603.                                           final boolean isMassInStm) {

  604.         TriggerDate triggerGenerator = null;

  605.         // check if we already have set up the provider
  606.         for (final AdditionalDataProvider<?> provider : getAdditionalDataProviders()) {
  607.             if (provider instanceof TriggerDate &&
  608.                 provider.getName().equals(driverName)) {
  609.                 // the Jacobian column generator has already been set up in a previous propagation
  610.                 triggerGenerator = (TriggerDate) provider;
  611.                 break;
  612.             }
  613.         }

  614.         if (triggerGenerator == null) {
  615.             // this is the first time we need the Jacobian column generator, create it
  616.             if (isMassInStm) {
  617.                 triggerGenerator = new TriggerDate(stmName, driverName, start, maneuver, threshold, true);
  618.             } else {
  619.                 final Optional<ForceModel> dragForce = getAllForceModels().stream().filter(AbstractDragForceModel.class::isInstance).findFirst();
  620.                 final Optional<ForceModel> srpForce = getAllForceModels().stream().filter(RadiationForceModel.class::isInstance).findFirst();
  621.                 final List<ForceModel> nonGravitationalForces = new ArrayList<>();
  622.                 dragForce.ifPresent(nonGravitationalForces::add);
  623.                 srpForce.ifPresent(nonGravitationalForces::add);
  624.                 triggerGenerator = new TriggerDate(stmName, driverName, start, maneuver, threshold, false,
  625.                         nonGravitationalForces.toArray(new ForceModel[0]));
  626.             }
  627.             mt.addResetter(triggerGenerator);
  628.             final MassDepletionDelay depletionDelay = triggerGenerator.getMassDepletionDelay();
  629.             if (depletionDelay != null) {
  630.                 addAdditionalDerivativesProvider(depletionDelay);
  631.             }
  632.             addAdditionalDataProvider(triggerGenerator);
  633.         }

  634.         if (!getInitialIntegrationState().hasAdditionalData(driverName)) {
  635.             // add the initial Jacobian column if it is not already there
  636.             // (perhaps due to a previous propagation)
  637.             final MassDepletionDelay depletionDelay = triggerGenerator.getMassDepletionDelay();
  638.             final double[] zeroes = new double[depletionDelay == null ? 7 : 6];
  639.             if (depletionDelay != null) {
  640.                 setInitialColumn(depletionDelay.getName(), zeroes);
  641.             }
  642.             setInitialColumn(driverName, getHarvester().getInitialJacobianColumn(driverName));
  643.         }

  644.         return triggerGenerator;

  645.     }

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

  654.         MedianDate medianGenerator = null;

  655.         // check if we already have set up the provider
  656.         for (final AdditionalDataProvider<?> provider : getAdditionalDataProviders()) {
  657.             if (provider instanceof MedianDate &&
  658.                 provider.getName().equals(medianName)) {
  659.                 // the Jacobian column generator has already been set up in a previous propagation
  660.                 medianGenerator = (MedianDate) provider;
  661.                 break;
  662.             }
  663.         }

  664.         if (medianGenerator == null) {
  665.             // this is the first time we need the Jacobian column generator, create it
  666.             medianGenerator = new MedianDate(startName, stopName, medianName);
  667.             addAdditionalDataProvider(medianGenerator);
  668.         }

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

  674.         return medianGenerator;

  675.     }

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

  684.         Duration durationGenerator = null;

  685.         // check if we already have set up the provider
  686.         for (final AdditionalDataProvider<?> provider : getAdditionalDataProviders()) {
  687.             if (provider instanceof Duration &&
  688.                 provider.getName().equals(durationName)) {
  689.                 // the Jacobian column generator has already been set up in a previous propagation
  690.                 durationGenerator = (Duration) provider;
  691.                 break;
  692.             }
  693.         }

  694.         if (durationGenerator == null) {
  695.             // this is the first time we need the Jacobian column generator, create it
  696.             durationGenerator = new Duration(startName, stopName, durationName);
  697.             addAdditionalDataProvider(durationGenerator);
  698.         }

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

  704.         return durationGenerator;

  705.     }

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

  713.         // first pass: gather all parameters (excluding trigger dates), binding similar names together
  714.         final ParameterDriversList selected = new ParameterDriversList();
  715.         for (final ForceModel forceModel : getAllForceModels()) {
  716.             for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
  717.                 if (!triggerDates.contains(driver.getNamesSpanMap().getFirstSpan().getData())) {
  718.                     // if the first span is not in triggerDates,
  719.                     // it means that the driver is not a trigger date and can be selected here
  720.                     selected.add(driver);
  721.                 }
  722.             }
  723.         }

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

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

  729.         // add the Jacobians column generators corresponding to parameters, and setup state accordingly
  730.         // a new column is needed for each value estimated so for each span of the parameterDriver
  731.         for (final DelegatingDriver driver : selected.getDrivers()) {

  732.             for (Span<String> currentNameSpan = driver.getNamesSpanMap().getFirstSpan(); currentNameSpan != null; currentNameSpan = currentNameSpan.next()) {

  733.                 IntegrableJacobianColumnGenerator generator = null;
  734.                 // check if we already have set up the providers
  735.                 for (final AdditionalDerivativesProvider provider : getAdditionalDerivativesProviders()) {
  736.                     if (provider instanceof IntegrableJacobianColumnGenerator &&
  737.                         provider.getName().equals(currentNameSpan.getData())) {
  738.                         // the Jacobian column generator has already been set up in a previous propagation
  739.                         generator = (IntegrableJacobianColumnGenerator) provider;
  740.                         break;
  741.                     }

  742.                 }

  743.                 if (generator == null) {
  744.                     // this is the first time we need the Jacobian column generator, create it
  745.                     final boolean isMassIncluded = stmGenerator.getStateDimension() == 7;
  746.                     generator = new IntegrableJacobianColumnGenerator(stmGenerator, currentNameSpan.getData(), isMassIncluded);
  747.                     addAdditionalDerivativesProvider(generator);
  748.                 }

  749.                 if (!getInitialIntegrationState().hasAdditionalData(currentNameSpan.getData())) {
  750.                     // add the initial Jacobian column if it is not already there
  751.                     // (perhaps due to a previous propagation)
  752.                     setInitialColumn(currentNameSpan.getData(), getHarvester().getInitialJacobianColumn(currentNameSpan.getData()));
  753.                 }

  754.             }

  755.         }

  756.     }

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

  767.         final SpacecraftState state = getInitialState();

  768.         final AbstractStateTransitionMatrixGenerator generator = (AbstractStateTransitionMatrixGenerator)
  769.                 getAdditionalDerivativesProviders().stream()
  770.                         .filter(AbstractStateTransitionMatrixGenerator.class::isInstance)
  771.                         .collect(Collectors.toList()).get(0);
  772.         final int expectedSize = generator.getStateDimension();
  773.         if (dYdQ.length != expectedSize) {
  774.             throw new OrekitException(LocalizedCoreFormats.DIMENSIONS_MISMATCH, dYdQ.length, expectedSize);
  775.         }

  776.         // convert to Cartesian Jacobian
  777.         final RealMatrix dYdC = MatrixUtils.createRealIdentityMatrix(expectedSize);
  778.         final double[][] jacobian = new double[6][6];
  779.         getOrbitType().convertType(state.getOrbit()).getJacobianWrtCartesian(getPositionAngleType(), jacobian);
  780.         dYdC.setSubMatrix(jacobian, 0, 0);
  781.         final DecompositionSolver solver = new QRDecomposition(dYdC, THRESHOLD).getSolver();
  782.         final double[] column = solver.solve(MatrixUtils.createRealVector(dYdQ)).toArray();

  783.         // set additional state
  784.         setInitialState(state.addAdditionalData(columnName, column));

  785.     }

  786.     /** {@inheritDoc} */
  787.     @Override
  788.     protected AttitudeProvider initializeAttitudeProviderForDerivatives() {
  789.         return needFullAttitudeForDerivatives ? getAttitudeProvider() : getFrozenAttitudeProvider();
  790.     }

  791.     /** {@inheritDoc} */
  792.     @Override
  793.     protected StateMapper createMapper(final AbsoluteDate referenceDate, final double mu,
  794.                                        final OrbitType orbitType, final PositionAngleType positionAngleType,
  795.                                        final AttitudeProvider attitudeProvider, final Frame frame) {
  796.         return new OsculatingMapper(referenceDate, mu, orbitType, positionAngleType, attitudeProvider, frame);
  797.     }

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

  800.         /** Simple constructor.
  801.          * <p>
  802.          * The position parameter type is meaningful only if {@link
  803.          * #getOrbitType() propagation orbit type}
  804.          * support it. As an example, it is not meaningful for propagation
  805.          * in {@link OrbitType#CARTESIAN Cartesian} parameters.
  806.          * </p>
  807.          * @param referenceDate reference date
  808.          * @param mu central attraction coefficient (m³/s²)
  809.          * @param orbitType orbit type to use for mapping (can be null for {@link AbsolutePVCoordinates})
  810.          * @param positionAngleType angle type to use for propagation
  811.          * @param attitudeProvider attitude provider
  812.          * @param frame inertial frame
  813.          */
  814.         OsculatingMapper(final AbsoluteDate referenceDate, final double mu,
  815.                          final OrbitType orbitType, final PositionAngleType positionAngleType,
  816.                          final AttitudeProvider attitudeProvider, final Frame frame) {
  817.             super(referenceDate, mu, orbitType, positionAngleType, attitudeProvider, frame);
  818.         }

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

  823.             final double mass = y[6];
  824.             if (mass <= 0.0) {
  825.                 throw new OrekitException(OrekitMessages.NOT_POSITIVE_SPACECRAFT_MASS, mass);
  826.             }

  827.             if (getOrbitType() == null) {
  828.                 // propagation uses absolute position-velocity-acceleration
  829.                 final Vector3D p = new Vector3D(y[0],    y[1],    y[2]);
  830.                 final Vector3D v = new Vector3D(y[3],    y[4],    y[5]);
  831.                 final Vector3D a;
  832.                 final AbsolutePVCoordinates absPva;
  833.                 if (yDot == null) {
  834.                     absPva = new AbsolutePVCoordinates(getFrame(), new TimeStampedPVCoordinates(date, p, v));
  835.                 } else {
  836.                     a = new Vector3D(yDot[3], yDot[4], yDot[5]);
  837.                     absPva = new AbsolutePVCoordinates(getFrame(), new TimeStampedPVCoordinates(date, p, v, a));
  838.                 }

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

  845.                 return new SpacecraftState(orbit, attitude).withMass(mass);
  846.             }

  847.         }

  848.         /** {@inheritDoc} */
  849.         public void mapStateToArray(final SpacecraftState state, final double[] y, final double[] yDot) {
  850.             if (getOrbitType() == null) {
  851.                 // propagation uses absolute position-velocity-acceleration
  852.                 final Vector3D p = state.getAbsPVA().getPosition();
  853.                 final Vector3D v = state.getAbsPVA().getVelocity();
  854.                 y[0] = p.getX();
  855.                 y[1] = p.getY();
  856.                 y[2] = p.getZ();
  857.                 y[3] = v.getX();
  858.                 y[4] = v.getY();
  859.                 y[5] = v.getZ();
  860.                 y[6] = state.getMass();
  861.             }
  862.             else {
  863.                 getOrbitType().mapOrbitToArray(state.getOrbit(), getPositionAngleType(), y, yDot);
  864.                 y[6] = state.getMass();
  865.             }
  866.         }

  867.     }

  868.     /** {@inheritDoc} */
  869.     protected MainStateEquations getMainStateEquations(final ODEIntegrator integrator) {
  870.         return new Main(integrator);
  871.     }

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

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

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

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

  880.         /** Flag keeping track whether Jacobian matrix needs to be recomputed or not. */
  881.         private boolean recomputingJacobian;

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

  886.             this.yDot     = new double[7];
  887.             this.jacobian = new double[6][6];
  888.             this.recomputingJacobian = true;

  889.             // feed internal event detectors
  890.             for (final ForceModel forceModel : forceModels) {
  891.                 forceModel.getEventDetectors().forEach(detector -> setUpEventDetector(integrator, detector));
  892.             }
  893.             getAttitudeProvider().getEventDetectors().forEach(detector -> setUpEventDetector(integrator, detector));

  894.             // default value for Jacobian is identity
  895.             for (int i = 0; i < jacobian.length; ++i) {
  896.                 Arrays.fill(jacobian[i], 0.0);
  897.                 jacobian[i][i] = 1.0;
  898.             }

  899.         }

  900.         /** {@inheritDoc} */
  901.         @Override
  902.         public void init(final SpacecraftState initialState, final AbsoluteDate target) {
  903.             needFullAttitudeForDerivatives = forceModels.stream().anyMatch(ForceModel::dependsOnAttitudeRate);

  904.             forceModels.forEach(fm -> fm.init(initialState, target));

  905.             final int numberOfForces = forceModels.size();
  906.             final OrbitType orbitType = getOrbitType();
  907.             if (orbitType != null && orbitType != OrbitType.CARTESIAN && numberOfForces > 0) {
  908.                 if (numberOfForces > 1) {
  909.                     recomputingJacobian = true;
  910.                 } else {
  911.                     recomputingJacobian = !(forceModels.get(0) instanceof NewtonianAttraction);
  912.                 }
  913.             } else {
  914.                 recomputingJacobian = false;
  915.             }
  916.         }

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

  920.             currentState = state;
  921.             Arrays.fill(yDot, 0.0);
  922.             if (recomputingJacobian) {
  923.                 // propagation uses Jacobian matrix of orbital parameters w.r.t. Cartesian ones
  924.                 currentState.getOrbit().getJacobianWrtCartesian(getPositionAngleType(), jacobian);
  925.             }

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

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


  940.             return yDot.clone();

  941.         }

  942.         /** {@inheritDoc} */
  943.         @Override
  944.         public void addKeplerContribution(final double mu) {
  945.             if (getOrbitType() == null) {
  946.                 // if mu is neither 0 nor NaN, we want to include Newtonian acceleration
  947.                 if (mu > 0) {
  948.                     // velocity derivative is Newtonian acceleration
  949.                     final Vector3D position = currentState.getPosition();
  950.                     final double r2         = position.getNormSq();
  951.                     final double coeff      = -mu / (r2 * FastMath.sqrt(r2));
  952.                     yDot[3] += coeff * position.getX();
  953.                     yDot[4] += coeff * position.getY();
  954.                     yDot[5] += coeff * position.getZ();
  955.                 }

  956.             } else {
  957.                 // propagation uses regular orbits
  958.                 currentState.getOrbit().addKeplerContribution(getPositionAngleType(), mu, yDot);
  959.             }
  960.         }

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

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

  976.     }

  977.     /** Estimate tolerance vectors for integrators when propagating in absolute position-velocity-acceleration.
  978.      * @param dP user specified position error
  979.      * @param absPva reference absolute position-velocity-acceleration
  980.      * @return a two rows array, row 0 being the absolute tolerance error and row 1
  981.      * being the relative tolerance error
  982.      * @deprecated since 13.0. Use {@link ToleranceProvider} for default and custom tolerances.
  983.      */
  984.     @Deprecated
  985.     public static double[][] tolerances(final double dP, final AbsolutePVCoordinates absPva) {
  986.         return ToleranceProvider.of(CartesianToleranceProvider.of(dP)).getTolerances(absPva);
  987.     }

  988.     /** Estimate tolerance vectors for integrators when propagating in orbits.
  989.      * <p>
  990.      * The errors are estimated from partial derivatives properties of orbits,
  991.      * starting from a scalar position error specified by the user.
  992.      * Considering the energy conservation equation V = sqrt(mu (2/r - 1/a)),
  993.      * we get at constant energy (i.e. on a Keplerian trajectory):
  994.      * <pre>
  995.      * V r² |dV| = mu |dr|
  996.      * </pre>
  997.      * <p> So we deduce a scalar velocity error consistent with the position error.
  998.      * From here, we apply orbits Jacobians matrices to get consistent errors
  999.      * on orbital parameters.
  1000.      *
  1001.      * <p>
  1002.      * The tolerances are only <em>orders of magnitude</em>, and integrator tolerances
  1003.      * are only local estimates, not global ones. So some care must be taken when using
  1004.      * these tolerances. Setting 1mm as a position error does NOT mean the tolerances
  1005.      * will guarantee a 1mm error position after several orbits integration.
  1006.      * </p>
  1007.      * @param dP user specified position error
  1008.      * @param orbit reference orbit
  1009.      * @param type propagation type for the meaning of the tolerance vectors elements
  1010.      * (it may be different from {@code orbit.getType()})
  1011.      * @return a two rows array, row 0 being the absolute tolerance error and row 1
  1012.      * being the relative tolerance error
  1013.      * @deprecated since 13.0. Use {@link ToleranceProvider} for default and custom tolerances.
  1014.      */
  1015.     @Deprecated
  1016.     public static double[][] tolerances(final double dP, final Orbit orbit, final OrbitType type) {
  1017.         return ToleranceProvider.getDefaultToleranceProvider(dP).getTolerances(orbit, type, PositionAngleType.TRUE);
  1018.     }

  1019.     /** Estimate tolerance vectors for integrators when propagating in orbits.
  1020.      * <p>
  1021.      * The errors are estimated from partial derivatives properties of orbits,
  1022.      * starting from scalar position and velocity errors specified by the user.
  1023.      * <p>
  1024.      * The tolerances are only <em>orders of magnitude</em>, and integrator tolerances
  1025.      * are only local estimates, not global ones. So some care must be taken when using
  1026.      * these tolerances. Setting 1mm as a position error does NOT mean the tolerances
  1027.      * will guarantee a 1mm error position after several orbits integration.
  1028.      * </p>
  1029.      * @param dP user specified position error
  1030.      * @param dV user specified velocity error
  1031.      * @param orbit reference orbit
  1032.      * @param type propagation type for the meaning of the tolerance vectors elements
  1033.      * (it may be different from {@code orbit.getType()})
  1034.      * @return a two rows array, row 0 being the absolute tolerance error and row 1
  1035.      * being the relative tolerance error
  1036.      * @since 10.3
  1037.      * @deprecated since 13.0. Use {@link ToleranceProvider} for default and custom tolerances.
  1038.      */
  1039.     @Deprecated
  1040.     public static double[][] tolerances(final double dP, final double dV,
  1041.                                         final Orbit orbit, final OrbitType type) {

  1042.         return ToleranceProvider.of(CartesianToleranceProvider.of(dP, dV, CartesianToleranceProvider.DEFAULT_ABSOLUTE_MASS_TOLERANCE))
  1043.                 .getTolerances(orbit, type, PositionAngleType.TRUE);
  1044.     }

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

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

  1049.             // inspect all force models to find InertialForces
  1050.             for (ForceModel force : forceModels) {
  1051.                 if (force instanceof InertialForces) {
  1052.                     return;
  1053.                 }
  1054.             }

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

  1057.         }

  1058.     }

  1059. }