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  
19  import org.hipparchus.exception.LocalizedCoreFormats;
20  import org.hipparchus.geometry.euclidean.threed.Vector3D;
21  import org.hipparchus.linear.DecompositionSolver;
22  import org.hipparchus.linear.MatrixUtils;
23  import org.hipparchus.linear.QRDecomposition;
24  import org.hipparchus.linear.RealMatrix;
25  import org.hipparchus.ode.ODEIntegrator;
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.OrekitInternalError;
33  import org.orekit.errors.OrekitMessages;
34  import org.orekit.forces.ForceModel;
35  import org.orekit.forces.drag.AbstractDragForceModel;
36  import org.orekit.forces.gravity.NewtonianAttraction;
37  import org.orekit.forces.inertia.InertialForces;
38  import org.orekit.forces.maneuvers.Maneuver;
39  import org.orekit.forces.maneuvers.jacobians.Duration;
40  import org.orekit.forces.maneuvers.jacobians.MassDepletionDelay;
41  import org.orekit.forces.maneuvers.jacobians.MedianDate;
42  import org.orekit.forces.maneuvers.jacobians.TriggerDate;
43  import org.orekit.forces.maneuvers.trigger.ManeuverTriggerDetector;
44  import org.orekit.forces.maneuvers.trigger.ResettableManeuverTriggers;
45  import org.orekit.forces.radiation.RadiationForceModel;
46  import org.orekit.frames.Frame;
47  import org.orekit.orbits.Orbit;
48  import org.orekit.orbits.OrbitType;
49  import org.orekit.orbits.PositionAngleType;
50  import org.orekit.propagation.AbstractMatricesHarvester;
51  import org.orekit.propagation.AdditionalDataProvider;
52  import org.orekit.propagation.CartesianToleranceProvider;
53  import org.orekit.propagation.MatricesHarvester;
54  import org.orekit.propagation.PropagationType;
55  import org.orekit.propagation.Propagator;
56  import org.orekit.propagation.SpacecraftState;
57  import org.orekit.propagation.ToleranceProvider;
58  import org.orekit.propagation.events.DetectorModifier;
59  import org.orekit.propagation.events.EventDetector;
60  import org.orekit.propagation.events.ParameterDrivenDateIntervalDetector;
61  import org.orekit.propagation.events.handlers.EventHandler;
62  import org.orekit.propagation.integration.AbstractIntegratedPropagator;
63  import org.orekit.propagation.integration.AdditionalDerivativesProvider;
64  import org.orekit.propagation.integration.StateMapper;
65  import org.orekit.time.AbsoluteDate;
66  import org.orekit.utils.AbsolutePVCoordinates;
67  import org.orekit.utils.DoubleArrayDictionary;
68  import org.orekit.utils.ParameterDriver;
69  import org.orekit.utils.ParameterDriversList;
70  import org.orekit.utils.ParameterDriversList.DelegatingDriver;
71  import org.orekit.utils.ParameterObserver;
72  import org.orekit.utils.TimeSpanMap;
73  import org.orekit.utils.TimeSpanMap.Span;
74  import org.orekit.utils.TimeStampedPVCoordinates;
75  
76  import java.util.ArrayList;
77  import java.util.Collection;
78  import java.util.Collections;
79  import java.util.List;
80  import java.util.Optional;
81  import java.util.stream.Collectors;
82  
83  /** This class propagates {@link org.orekit.orbits.Orbit orbits} using
84   * numerical integration.
85   * <p>Numerical propagation is much more accurate than analytical propagation
86   * like for example {@link org.orekit.propagation.analytical.KeplerianPropagator
87   * Keplerian} or {@link org.orekit.propagation.analytical.EcksteinHechlerPropagator
88   * Eckstein-Hechler}, but requires a few more steps to set up to be used properly.
89   * Whereas analytical propagators are configured only thanks to their various
90   * constructors and can be used immediately after construction, numerical propagators
91   * configuration involve setting several parameters between construction time
92   * and propagation time.</p>
93   * <p>The configuration parameters that can be set are:</p>
94   * <ul>
95   *   <li>the initial spacecraft state ({@link #setInitialState(SpacecraftState)})</li>
96   *   <li>the central attraction coefficient ({@link #setMu(double)})</li>
97   *   <li>the various force models ({@link #addForceModel(ForceModel)},
98   *   {@link #removeForceModels()})</li>
99   *   <li>the {@link OrbitType type} of orbital parameters to be used for propagation
100  *   ({@link #setOrbitType(OrbitType)}),</li>
101  *   <li>the {@link PositionAngleType type} of position angle to be used in orbital parameters
102  *   to be used for propagation where it is relevant ({@link
103  *   #setPositionAngleType(PositionAngleType)}),</li>
104  *   <li>whether {@link MatricesHarvester state transition matrices and Jacobians matrices}
105  *   (with the option to include mass if a 7x7 initial matrix is passed) should be propagated along with orbital state
106  *   ({@link #setupMatricesComputation(String, RealMatrix, DoubleArrayDictionary)}),</li>
107  *   <li>whether {@link org.orekit.propagation.integration.AdditionalDerivativesProvider additional derivatives}
108  *   should be propagated along with orbital state ({@link
109  *   #addAdditionalDerivativesProvider(AdditionalDerivativesProvider)}),</li>
110  *   <li>the discrete events that should be triggered during propagation
111  *   ({@link #addEventDetector(EventDetector)},
112  *   {@link #clearEventsDetectors()})</li>
113  *   <li>the binding logic with the rest of the application ({@link #getMultiplexer()})</li>
114  * </ul>
115  * <p>From these configuration parameters, only the initial state is mandatory. The default
116  * propagation settings are in {@link OrbitType#EQUINOCTIAL equinoctial} parameters with
117  * {@link PositionAngleType#ECCENTRIC} longitude argument. If the central attraction coefficient
118  * is not explicitly specified, the one used to define the initial orbit will be used.
119  * However, specifying only the initial state and perhaps the central attraction coefficient
120  * would mean the propagator would use only Keplerian forces. In this case, the simpler {@link
121  * org.orekit.propagation.analytical.KeplerianPropagator KeplerianPropagator} class would
122  * perhaps be more effective.</p>
123  * <p>The underlying numerical integrator set up in the constructor may also have its own
124  * configuration parameters. Typical configuration parameters for adaptive stepsize integrators
125  * are the min, max and perhaps start step size as well as the absolute and/or relative errors
126  * thresholds.</p>
127  * <p>The state that is seen by the integrator is a simple seven elements double array.
128  * The six first elements are either:
129  * <ul>
130  *   <li>the {@link org.orekit.orbits.EquinoctialOrbit equinoctial orbit parameters} (a, e<sub>x</sub>,
131  *   e<sub>y</sub>, h<sub>x</sub>, h<sub>y</sub>, λ<sub>M</sub> or λ<sub>E</sub>
132  *   or λ<sub>v</sub>) in meters and radians,</li>
133  *   <li>the {@link org.orekit.orbits.KeplerianOrbit Keplerian orbit parameters} (a, e, i, ω, Ω,
134  *   M or E or v) in meters and radians,</li>
135  *   <li>the {@link org.orekit.orbits.CircularOrbit circular orbit parameters} (a, e<sub>x</sub>, e<sub>y</sub>, i,
136  *   Ω, α<sub>M</sub> or α<sub>E</sub> or α<sub>v</sub>) in meters
137  *   and radians,</li>
138  *   <li>the {@link org.orekit.orbits.CartesianOrbit Cartesian orbit parameters} (x, y, z, v<sub>x</sub>,
139  *   v<sub>y</sub>, v<sub>z</sub>) in meters and meters per seconds.
140  * </ul>
141  * <p> The last element is the mass in kilograms and changes only during thrusters firings
142  *
143  * <p>The following code snippet shows a typical setting for Low Earth Orbit propagation in
144  * equinoctial parameters and true longitude argument:</p>
145  * <pre>
146  * final double dP       = 0.001;
147  * final double minStep  = 0.001;
148  * final double maxStep  = 500;
149  * final double initStep = 60;
150  * final double[][] tolerance = ToleranceProvider.getDefaultToleranceProvider(dP).getTolerances(orbit, OrbitType.EQUINOCTIAL);
151  * AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tolerance[0], tolerance[1]);
152  * integrator.setInitialStepSize(initStep);
153  * propagator = new NumericalPropagator(integrator);
154  * </pre>
155  * <p>By default, at the end of the propagation, the propagator resets the initial state to the final state,
156  * thus allowing a new propagation to be started from there without recomputing the part already performed.
157  * This behaviour can be changed by calling {@link #setResetAtEnd(boolean)}.
158  * </p>
159  * <p>Beware the same instance cannot be used simultaneously by different threads, the class is <em>not</em>
160  * thread-safe.</p>
161  *
162  * @see SpacecraftState
163  * @see ForceModel
164  * @see org.orekit.propagation.sampling.OrekitStepHandler
165  * @see org.orekit.propagation.sampling.OrekitFixedStepHandler
166  * @see org.orekit.propagation.integration.IntegratedEphemeris
167  * @see TimeDerivativesEquations
168  *
169  * @author Mathieu Rom&eacute;ro
170  * @author Luc Maisonobe
171  * @author Guylaine Prat
172  * @author Fabien Maussion
173  * @author V&eacute;ronique Pommier-Maurussane
174  */
175 public class NumericalPropagator extends AbstractIntegratedPropagator {
176 
177     /** Default orbit type. */
178     public static final OrbitType DEFAULT_ORBIT_TYPE = OrbitType.EQUINOCTIAL;
179 
180     /** Default position angle type. */
181     public static final PositionAngleType DEFAULT_POSITION_ANGLE_TYPE = PositionAngleType.ECCENTRIC;
182 
183     /** Threshold for matrix solving. */
184     private static final double THRESHOLD = Precision.SAFE_MIN;
185 
186     /** Force models used during the extrapolation of the orbit. */
187     private final List<ForceModel> forceModels;
188 
189     /** boolean to ignore or not the creation of a NewtonianAttraction. */
190     private boolean ignoreCentralAttraction;
191 
192     /**
193      * boolean to know if a full attitude (with rates) is needed when computing derivatives for the ODE.
194      * since 12.1
195      */
196     private boolean needFullAttitudeForDerivatives = true;
197 
198     /** Create a new instance of NumericalPropagator, based on orbit definition mu.
199      * After creation, the instance is empty, i.e. the attitude provider is set to an
200      * unspecified default law and there are no perturbing forces at all.
201      * This means that if {@link #addForceModel addForceModel} is not
202      * called after creation, the integrated orbit will follow a Keplerian
203      * evolution only. The defaults are {@link OrbitType#EQUINOCTIAL}
204      * for {@link #setOrbitType(OrbitType) propagation
205      * orbit type} and {@link PositionAngleType#ECCENTRIC} for {@link
206      * #setPositionAngleType(PositionAngleType) position angle type}.
207      *
208      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
209      *
210      * @param integrator numerical integrator to use for propagation.
211      * @see #NumericalPropagator(ODEIntegrator, AttitudeProvider)
212      */
213     @DefaultDataContext
214     public NumericalPropagator(final ODEIntegrator integrator) {
215         this(integrator,
216                 Propagator.getDefaultLaw(DataContext.getDefault().getFrames()));
217     }
218 
219     /** Create a new instance of NumericalPropagator, based on orbit definition mu.
220      * After creation, the instance is empty, i.e. the attitude provider is set to an
221      * unspecified default law and there are no perturbing forces at all.
222      * This means that if {@link #addForceModel addForceModel} is not
223      * called after creation, the integrated orbit will follow a Keplerian
224      * evolution only. The defaults are {@link OrbitType#EQUINOCTIAL}
225      * for {@link #setOrbitType(OrbitType) propagation
226      * orbit type} and {@link PositionAngleType#ECCENTRIC} for {@link
227      * #setPositionAngleType(PositionAngleType) position angle type}.
228      * @param integrator numerical integrator to use for propagation.
229      * @param attitudeProvider the attitude law.
230      * @since 10.1
231      */
232     public NumericalPropagator(final ODEIntegrator integrator,
233                                final AttitudeProvider attitudeProvider) {
234         super(integrator, PropagationType.OSCULATING);
235         forceModels             = new ArrayList<>();
236         ignoreCentralAttraction = false;
237         initMapper();
238         setAttitudeProvider(attitudeProvider);
239         clearStepHandlers();
240         setOrbitType(DEFAULT_ORBIT_TYPE);
241         setPositionAngleType(DEFAULT_POSITION_ANGLE_TYPE);
242     }
243 
244     /** Set the flag to ignore or not the creation of a {@link NewtonianAttraction}.
245      * @param ignoreCentralAttraction if true, {@link NewtonianAttraction} is <em>not</em>
246      * added automatically if missing
247      */
248     public void setIgnoreCentralAttraction(final boolean ignoreCentralAttraction) {
249         this.ignoreCentralAttraction = ignoreCentralAttraction;
250     }
251 
252      /** Set the central attraction coefficient μ.
253       * <p>
254       * Setting the central attraction coefficient is
255       * equivalent to {@link #addForceModel(ForceModel) add}
256       * a {@link NewtonianAttraction} force model.
257       * * </p>
258       * @param mu central attraction coefficient (m³/s²)
259       * @see #addForceModel(ForceModel)
260       * @see #getAllForceModels()
261       */
262     @Override
263     public void setMu(final double mu) {
264         if (ignoreCentralAttraction) {
265             superSetMu(mu);
266         } else {
267             addForceModel(new NewtonianAttraction(mu));
268             superSetMu(mu);
269         }
270     }
271 
272     /** Set the central attraction coefficient μ only in upper class.
273      * @param mu central attraction coefficient (m³/s²)
274      */
275     private void superSetMu(final double mu) {
276         super.setMu(mu);
277     }
278 
279     /** Check if Newtonian attraction force model is available.
280      * <p>
281      * Newtonian attraction is always the last force model in the list.
282      * </p>
283      * @return true if Newtonian attraction force model is available
284      */
285     private boolean hasNewtonianAttraction() {
286         final int last = forceModels.size() - 1;
287         return last >= 0 && forceModels.get(last) instanceof NewtonianAttraction;
288     }
289 
290     /** Add a force model.
291      * <p>If this method is not called at all, the integrated orbit will follow
292      * a Keplerian evolution only.</p>
293      * @param model {@link ForceModel} to add (it can be either a perturbing force
294      * model or an instance of {@link NewtonianAttraction})
295      * @see #removeForceModels()
296      * @see #setMu(double)
297      */
298     public void addForceModel(final ForceModel model) {
299 
300         if (model instanceof NewtonianAttraction) {
301             // we want to add the central attraction force model
302 
303             try {
304                 // ensure we are notified of any mu change
305                 model.getParametersDrivers().get(0).addObserver(new ParameterObserver() {
306                     /** {@inheritDoc} */
307                     @Override
308                     public void valueSpanMapChanged(final TimeSpanMap<Double> previousValue, final ParameterDriver driver) {
309                         superSetMu(driver.getValue());
310                     }
311                     /** {@inheritDoc} */
312                     @Override
313                     public void valueChanged(final double previousValue, final ParameterDriver driver, final AbsoluteDate date) {
314                         superSetMu(driver.getValue());
315                     }
316                 });
317             } catch (OrekitException oe) {
318                 // this should never happen
319                 throw new OrekitInternalError(oe);
320             }
321 
322             if (hasNewtonianAttraction()) {
323                 // there is already a central attraction model, replace it
324                 forceModels.set(forceModels.size() - 1, model);
325             } else {
326                 // there are no central attraction model yet, add it at the end of the list
327                 forceModels.add(model);
328             }
329         } else {
330             // we want to add a perturbing force model
331             if (hasNewtonianAttraction()) {
332                 // insert the new force model before Newtonian attraction,
333                 // which should always be the last one in the list
334                 forceModels.add(forceModels.size() - 1, model);
335             } else {
336                 // we only have perturbing force models up to now, just append at the end of the list
337                 forceModels.add(model);
338             }
339         }
340 
341     }
342 
343     /** Remove all force models (except central attraction).
344      * <p>Once all perturbing forces have been removed (and as long as no new force
345      * model is added), the integrated orbit will follow a Keplerian evolution
346      * only.</p>
347      * @see #addForceModel(ForceModel)
348      */
349     public void removeForceModels() {
350         final int last = forceModels.size() - 1;
351         if (hasNewtonianAttraction()) {
352             // preserve the Newtonian attraction model at the end
353             final ForceModel newton = forceModels.get(last);
354             forceModels.clear();
355             forceModels.add(newton);
356         } else {
357             forceModels.clear();
358         }
359     }
360 
361     /** Get all the force models, perturbing forces and Newtonian attraction included.
362      * @return list of perturbing force models, with Newtonian attraction being the
363      * last one
364      * @see #addForceModel(ForceModel)
365      * @see #setMu(double)
366      */
367     public List<ForceModel> getAllForceModels() {
368         return Collections.unmodifiableList(forceModels);
369     }
370 
371     /** Set propagation orbit type.
372      * @param orbitType orbit type to use for propagation, null for
373      * propagating using {@link org.orekit.utils.AbsolutePVCoordinates} rather than {@link Orbit}
374      */
375     @Override
376     public void setOrbitType(final OrbitType orbitType) {
377         super.setOrbitType(orbitType);
378     }
379 
380     /** Get propagation parameter type.
381      * @return orbit type used for propagation, null for
382      * propagating using {@link org.orekit.utils.AbsolutePVCoordinates} rather than {@link Orbit}
383      */
384     @Override
385     public OrbitType getOrbitType() {
386         return super.getOrbitType();
387     }
388 
389     /** Set position angle type.
390      * <p>
391      * The position parameter type is meaningful only if {@link
392      * #getOrbitType() propagation orbit type}
393      * support it. As an example, it is not meaningful for propagation
394      * in {@link OrbitType#CARTESIAN Cartesian} parameters.
395      * </p>
396      * @param positionAngleType angle type to use for propagation
397      */
398     @Override
399     public void setPositionAngleType(final PositionAngleType positionAngleType) {
400         super.setPositionAngleType(positionAngleType);
401     }
402 
403     /** Get propagation parameter type.
404      * @return angle type to use for propagation
405      */
406     @Override
407     public PositionAngleType getPositionAngleType() {
408         return super.getPositionAngleType();
409     }
410 
411     /** Set the initial state.
412      * @param initialState initial state
413      */
414     public void setInitialState(final SpacecraftState initialState) {
415         resetInitialState(initialState);
416     }
417 
418     /** {@inheritDoc} */
419     @Override
420     public void resetInitialState(final SpacecraftState state) {
421         super.resetInitialState(state);
422         if (!hasNewtonianAttraction()) {
423             // use the state to define central attraction
424             setMu(state.isOrbitDefined() ? state.getOrbit().getMu() : Double.NaN);
425         }
426         setStartDate(state.getDate());
427     }
428 
429     /** Get the names of the parameters in the matrix returned by {@link MatricesHarvester#getParametersJacobian}.
430      * @return names of the parameters (i.e. columns) of the Jacobian matrix
431      */
432     List<String> getJacobiansColumnsNames() {
433         final List<String> columnsNames = new ArrayList<>();
434         for (final ForceModel forceModel : getAllForceModels()) {
435             for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
436                 if (driver.isSelected() && !columnsNames.contains(driver.getNamesSpanMap().getFirstSpan().getData())) {
437                     // As driver with same name should have same NamesSpanMap we only check if the first span is present,
438                     // if not we add all span names to columnsNames
439                     for (Span<String> span = driver.getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {
440                         columnsNames.add(span.getData());
441                     }
442                 }
443             }
444         }
445         Collections.sort(columnsNames);
446         return columnsNames;
447     }
448 
449     /** {@inheritDoc}
450      * <p>
451      * Unlike other propagators, the numerical one can consider the mass as a state variable in the transition matrix.
452      * To do so, a 7x7 initial matrix is to be passed instead of 6x6.
453      * </p>
454      * */
455     @Override
456     protected AbstractMatricesHarvester createHarvester(final String stmName, final RealMatrix initialStm,
457                                                         final DoubleArrayDictionary initialJacobianColumns) {
458         return new NumericalPropagationHarvester(this, stmName, initialStm, initialJacobianColumns);
459     }
460 
461     /** {@inheritDoc} */
462     @Override
463     public void clearMatricesComputation() {
464         final List<AdditionalDerivativesProvider> copiedDerivativesProviders = new ArrayList<>(getAdditionalDerivativesProviders());
465         copiedDerivativesProviders.stream().filter(AbstractStateTransitionMatrixGenerator.class::isInstance)
466                 .forEach(provider -> removeAdditionalDerivativesProvider(provider.getName()));
467         final List<AdditionalDataProvider<?>> copiedDataProviders = new ArrayList<>(getAdditionalDataProviders());
468         for (final AdditionalDataProvider<?> additionalDataProvider: copiedDataProviders) {
469             if (additionalDataProvider instanceof TriggerDate) {
470                 final TriggerDate triggerDate = (TriggerDate) additionalDataProvider;
471                 if (triggerDate.getMassDepletionDelay() != null) {
472                     removeAdditionalDerivativesProvider(triggerDate.getMassDepletionDelay().getName());
473                 }
474                 removeAdditionalDataProvider(additionalDataProvider.getName());
475             } else if (additionalDataProvider instanceof MedianDate || additionalDataProvider instanceof Duration) {
476                 removeAdditionalDataProvider(additionalDataProvider.getName());
477             }
478         }
479         super.clearMatricesComputation();
480     }
481 
482     /** {@inheritDoc} */
483     @Override
484     protected void setUpStmAndJacobianGenerators() {
485 
486         final AbstractMatricesHarvester harvester = getHarvester();
487         if (harvester != null) {
488 
489             // set up the additional equations and additional state providers
490             final AbstractStateTransitionMatrixGenerator stmGenerator = setUpStmGenerator();
491             final List<String> triggersDates = setUpTriggerDatesJacobiansColumns(stmGenerator);
492             setUpRegularParametersJacobiansColumns(stmGenerator, triggersDates);
493 
494             // as we are now starting the propagation, everything is configured
495             // we can freeze the names in the harvester
496             harvester.freezeColumnsNames();
497 
498         }
499 
500     }
501 
502     /** Set up the State Transition Matrix Generator.
503      * @return State Transition Matrix Generator
504      * @since 11.1
505      */
506     private AbstractStateTransitionMatrixGenerator setUpStmGenerator() {
507 
508         final AbstractMatricesHarvester harvester = getHarvester();
509 
510         // add the STM generator corresponding to the current settings, and setup state accordingly
511         AbstractStateTransitionMatrixGenerator stmGenerator = null;
512         for (final AdditionalDerivativesProvider equations : getAdditionalDerivativesProviders()) {
513             if (equations instanceof AbstractStateTransitionMatrixGenerator &&
514                 equations.getName().equals(harvester.getStmName())) {
515                 // the STM generator has already been set up in a previous propagation
516                 stmGenerator = (AbstractStateTransitionMatrixGenerator) equations;
517                 break;
518             }
519         }
520         if (stmGenerator == null) {
521             // this is the first time we need the STM generate, create it
522             if (harvester.getStateDimension() > 6) {
523                 stmGenerator = new ExtendedStateTransitionMatrixGenerator(harvester.getStmName(), getAllForceModels(),
524                         getAttitudeProvider());
525             } else {
526                 stmGenerator = new StateTransitionMatrixGenerator(harvester.getStmName(), getAllForceModels(),
527                         getAttitudeProvider());
528             }
529             addAdditionalDerivativesProvider(stmGenerator);
530         }
531 
532         if (!getInitialIntegrationState().hasAdditionalData(harvester.getStmName())) {
533             // add the initial State Transition Matrix if it is not already there
534             // (perhaps due to a previous propagation)
535             setInitialState(stmGenerator.setInitialStateTransitionMatrix(getInitialState(),
536                                                                          harvester.getInitialStateTransitionMatrix(),
537                                                                          getOrbitType(),
538                                                                          getPositionAngleType()));
539         }
540 
541         return stmGenerator;
542 
543     }
544 
545     /** Set up the Jacobians columns generator dedicated to trigger dates.
546      * @param stmGenerator State Transition Matrix generator
547      * @return names of the columns corresponding to trigger dates
548      * @since 13.1
549      */
550     private List<String> setUpTriggerDatesJacobiansColumns(final AbstractStateTransitionMatrixGenerator stmGenerator) {
551 
552         final String stmName = stmGenerator.getName();
553         final boolean isMassInStm = stmGenerator instanceof ExtendedStateTransitionMatrixGenerator;
554         final List<String> names = new ArrayList<>();
555         for (final ForceModel forceModel : getAllForceModels()) {
556             if (forceModel instanceof Maneuver && ((Maneuver) forceModel).getManeuverTriggers() instanceof ResettableManeuverTriggers) {
557                 final Maneuver maneuver = (Maneuver) forceModel;
558                 final ResettableManeuverTriggers maneuverTriggers = (ResettableManeuverTriggers) maneuver.getManeuverTriggers();
559 
560                 final Collection<EventDetector> selectedDetectors = maneuverTriggers.getEventDetectors().
561                         filter(ManeuverTriggerDetector.class::isInstance).
562                         map(triggerDetector -> ((ManeuverTriggerDetector<?>) triggerDetector).getDetector())
563                         .collect(Collectors.toList());
564                 for (final EventDetector detector: selectedDetectors) {
565                     if (detector instanceof ParameterDrivenDateIntervalDetector) {
566                         final ParameterDrivenDateIntervalDetector d = (ParameterDrivenDateIntervalDetector) detector;
567                         TriggerDate start;
568                         TriggerDate stop;
569 
570                         if (d.getStartDriver().isSelected() || d.getMedianDriver().isSelected() || d.getDurationDriver().isSelected()) {
571                             // normally datedriver should have only 1 span but just in case the user defines several span, there will
572                             // be no problem here
573                             for (Span<String> span = d.getStartDriver().getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {
574                                 start = manageTriggerDate(stmName, maneuver, maneuverTriggers, span.getData(), true,
575                                         d.getThreshold(), isMassInStm);
576                                 names.add(start.getName());
577                                 start = null;
578                             }
579                         }
580                         if (d.getStopDriver().isSelected() || d.getMedianDriver().isSelected() || d.getDurationDriver().isSelected()) {
581                             // normally datedriver should have only 1 span but just in case the user defines several span, there will
582                             // be no problem here
583                             for (Span<String> span = d.getStopDriver().getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {
584                                 stop = manageTriggerDate(stmName, maneuver, maneuverTriggers, span.getData(), false,
585                                         d.getThreshold(), isMassInStm);
586                                 names.add(stop.getName());
587                                 stop = null;
588                             }
589                         }
590                         if (d.getMedianDriver().isSelected()) {
591                             // for first span
592                             Span<String> currentMedianNameSpan = d.getMedianDriver().getNamesSpanMap().getFirstSpan();
593                             MedianDate median =
594                                     manageMedianDate(d.getStartDriver().getNamesSpanMap().getFirstSpan().getData(),
595                                             d.getStopDriver().getNamesSpanMap().getFirstSpan().getData(), currentMedianNameSpan.getData());
596                             names.add(median.getName());
597                             // for all span
598                             // normally datedriver should have only 1 span but just in case the user defines several span, there will
599                             // be no problem here. /!\ medianDate driver, startDate driver and stopDate driver must have same span number
600                             for (int spanNumber = 1; spanNumber < d.getMedianDriver().getNamesSpanMap().getSpansNumber(); ++spanNumber) {
601                                 currentMedianNameSpan = d.getMedianDriver().getNamesSpanMap().getSpan(currentMedianNameSpan.getEnd());
602                                 median =
603                                         manageMedianDate(d.getStartDriver().getNamesSpanMap().getSpan(currentMedianNameSpan.getStart()).getData(),
604                                                 d.getStopDriver().getNamesSpanMap().getSpan(currentMedianNameSpan.getStart()).getData(),
605                                                 currentMedianNameSpan.getData());
606                                 names.add(median.getName());
607 
608                             }
609 
610                         }
611                         if (d.getDurationDriver().isSelected()) {
612                             // for first span
613                             Span<String> currentDurationNameSpan = d.getDurationDriver().getNamesSpanMap().getFirstSpan();
614                             Duration duration =
615                                     manageManeuverDuration(d.getStartDriver().getNamesSpanMap().getFirstSpan().getData(),
616                                             d.getStopDriver().getNamesSpanMap().getFirstSpan().getData(), currentDurationNameSpan.getData());
617                             names.add(duration.getName());
618                             // for all span
619                             for (int spanNumber = 1; spanNumber < d.getDurationDriver().getNamesSpanMap().getSpansNumber(); ++spanNumber) {
620                                 currentDurationNameSpan = d.getDurationDriver().getNamesSpanMap().getSpan(currentDurationNameSpan.getEnd());
621                                 duration =
622                                         manageManeuverDuration(d.getStartDriver().getNamesSpanMap().getSpan(currentDurationNameSpan.getStart()).getData(),
623                                                 d.getStopDriver().getNamesSpanMap().getSpan(currentDurationNameSpan.getStart()).getData(),
624                                                 currentDurationNameSpan.getData());
625                                 names.add(duration.getName());
626 
627                             }
628                         }
629                     }
630                 }
631             }
632         }
633 
634         return names;
635 
636     }
637 
638     /** Manage a maneuver trigger date.
639      * @param stmName name of the State Transition Matrix state
640      * @param maneuver maneuver force model
641      * @param mt trigger to which the driver is bound
642      * @param driverName name of the date driver
643      * @param start if true, the driver is a maneuver start
644      * @param threshold event detector threshold
645      * @param isMassInStm flag on presence on mass in STM
646      * @return generator for the date driver
647      * @since 13.1
648      */
649     private TriggerDate manageTriggerDate(final String stmName,
650                                           final Maneuver maneuver,
651                                           final ResettableManeuverTriggers mt,
652                                           final String driverName,
653                                           final boolean start,
654                                           final double threshold,
655                                           final boolean isMassInStm) {
656 
657         TriggerDate triggerGenerator = null;
658 
659         // check if we already have set up the provider
660         for (final AdditionalDataProvider<?> provider : getAdditionalDataProviders()) {
661             if (provider instanceof TriggerDate &&
662                 provider.getName().equals(driverName)) {
663                 // the Jacobian column generator has already been set up in a previous propagation
664                 triggerGenerator = (TriggerDate) provider;
665                 break;
666             }
667         }
668 
669         if (triggerGenerator == null) {
670             // this is the first time we need the Jacobian column generator, create it
671             if (isMassInStm) {
672                 triggerGenerator = new TriggerDate(stmName, driverName, start, maneuver, threshold, true);
673             } else {
674                 final Optional<ForceModel> dragForce = getAllForceModels().stream().filter(AbstractDragForceModel.class::isInstance).findFirst();
675                 final Optional<ForceModel> srpForce = getAllForceModels().stream().filter(RadiationForceModel.class::isInstance).findFirst();
676                 final List<ForceModel> nonGravitationalForces = new ArrayList<>();
677                 dragForce.ifPresent(nonGravitationalForces::add);
678                 srpForce.ifPresent(nonGravitationalForces::add);
679                 triggerGenerator = new TriggerDate(stmName, driverName, start, maneuver, threshold, false,
680                         nonGravitationalForces.toArray(new ForceModel[0]));
681             }
682             mt.addResetter(triggerGenerator);
683             final MassDepletionDelay depletionDelay = triggerGenerator.getMassDepletionDelay();
684             if (depletionDelay != null) {
685                 addAdditionalDerivativesProvider(depletionDelay);
686             }
687             addAdditionalDataProvider(triggerGenerator);
688         }
689 
690         if (!getInitialIntegrationState().hasAdditionalData(driverName)) {
691             // add the initial Jacobian column if it is not already there
692             // (perhaps due to a previous propagation)
693             final MassDepletionDelay depletionDelay = triggerGenerator.getMassDepletionDelay();
694             final double[] zeroes = new double[depletionDelay == null ? 7 : 6];
695             if (depletionDelay != null) {
696                 setInitialColumn(depletionDelay.getName(), zeroes);
697             }
698             setInitialColumn(driverName, getHarvester().getInitialJacobianColumn(driverName));
699         }
700 
701         return triggerGenerator;
702 
703     }
704 
705     /** Manage a maneuver median date.
706      * @param startName name of the start driver
707      * @param stopName name of the stop driver
708      * @param medianName name of the median driver
709      * @return generator for the median driver
710      * @since 11.1
711      */
712     private MedianDate manageMedianDate(final String startName, final String stopName, final String medianName) {
713 
714         MedianDate medianGenerator = null;
715 
716         // check if we already have set up the provider
717         for (final AdditionalDataProvider<?> provider : getAdditionalDataProviders()) {
718             if (provider instanceof MedianDate &&
719                 provider.getName().equals(medianName)) {
720                 // the Jacobian column generator has already been set up in a previous propagation
721                 medianGenerator = (MedianDate) provider;
722                 break;
723             }
724         }
725 
726         if (medianGenerator == null) {
727             // this is the first time we need the Jacobian column generator, create it
728             medianGenerator = new MedianDate(startName, stopName, medianName);
729             addAdditionalDataProvider(medianGenerator);
730         }
731 
732         if (!getInitialIntegrationState().hasAdditionalData(medianName)) {
733             // add the initial Jacobian column if it is not already there
734             // (perhaps due to a previous propagation)
735             setInitialColumn(medianName, getHarvester().getInitialJacobianColumn(medianName));
736         }
737 
738         return medianGenerator;
739 
740     }
741 
742     /** Manage a maneuver duration.
743      * @param startName name of the start driver
744      * @param stopName name of the stop driver
745      * @param durationName name of the duration driver
746      * @return generator for the median driver
747      * @since 11.1
748      */
749     private Duration manageManeuverDuration(final String startName, final String stopName, final String durationName) {
750 
751         Duration durationGenerator = null;
752 
753         // check if we already have set up the provider
754         for (final AdditionalDataProvider<?> provider : getAdditionalDataProviders()) {
755             if (provider instanceof Duration &&
756                 provider.getName().equals(durationName)) {
757                 // the Jacobian column generator has already been set up in a previous propagation
758                 durationGenerator = (Duration) provider;
759                 break;
760             }
761         }
762 
763         if (durationGenerator == null) {
764             // this is the first time we need the Jacobian column generator, create it
765             durationGenerator = new Duration(startName, stopName, durationName);
766             addAdditionalDataProvider(durationGenerator);
767         }
768 
769         if (!getInitialIntegrationState().hasAdditionalData(durationName)) {
770             // add the initial Jacobian column if it is not already there
771             // (perhaps due to a previous propagation)
772             setInitialColumn(durationName, getHarvester().getInitialJacobianColumn(durationName));
773         }
774 
775         return durationGenerator;
776 
777     }
778 
779     /** Set up the Jacobians columns generator for regular parameters.
780      * @param stmGenerator generator for the State Transition Matrix
781      * @param triggerDates names of the columns already managed as trigger dates
782      * @since 11.1
783      */
784     private void setUpRegularParametersJacobiansColumns(final AbstractStateTransitionMatrixGenerator stmGenerator,
785                                                         final List<String> triggerDates) {
786 
787         // first pass: gather all parameters (excluding trigger dates), binding similar names together
788         final ParameterDriversList selected = new ParameterDriversList();
789         for (final ForceModel forceModel : getAllForceModels()) {
790             for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
791                 if (!triggerDates.contains(driver.getNamesSpanMap().getFirstSpan().getData())) {
792                     // if the first span is not in triggerDates,
793                     // it means that the driver is not a trigger date and can be selected here
794                     selected.add(driver);
795                 }
796             }
797         }
798 
799         // second pass: now that shared parameter names are bound together,
800         // their selections status have been synchronized, we can filter them
801         selected.filter(true);
802 
803         // third pass: sort parameters lexicographically
804         selected.sort();
805 
806         // add the Jacobians column generators corresponding to parameters, and setup state accordingly
807         // a new column is needed for each value estimated so for each span of the parameterDriver
808         for (final DelegatingDriver driver : selected.getDrivers()) {
809 
810             for (Span<String> currentNameSpan = driver.getNamesSpanMap().getFirstSpan(); currentNameSpan != null; currentNameSpan = currentNameSpan.next()) {
811 
812                 IntegrableJacobianColumnGenerator generator = null;
813                 // check if we already have set up the providers
814                 for (final AdditionalDerivativesProvider provider : getAdditionalDerivativesProviders()) {
815                     if (provider instanceof IntegrableJacobianColumnGenerator &&
816                         provider.getName().equals(currentNameSpan.getData())) {
817                         // the Jacobian column generator has already been set up in a previous propagation
818                         generator = (IntegrableJacobianColumnGenerator) provider;
819                         break;
820                     }
821 
822                 }
823 
824                 if (generator == null) {
825                     // this is the first time we need the Jacobian column generator, create it
826                     final boolean isMassIncluded = stmGenerator.getStateDimension() == 7;
827                     generator = new IntegrableJacobianColumnGenerator(stmGenerator, currentNameSpan.getData(), isMassIncluded);
828                     addAdditionalDerivativesProvider(generator);
829                 }
830 
831                 if (!getInitialIntegrationState().hasAdditionalData(currentNameSpan.getData())) {
832                     // add the initial Jacobian column if it is not already there
833                     // (perhaps due to a previous propagation)
834                     setInitialColumn(currentNameSpan.getData(), getHarvester().getInitialJacobianColumn(currentNameSpan.getData()));
835                 }
836 
837             }
838 
839         }
840 
841     }
842 
843     /** Add the initial value of the column to the initial state.
844      * <p>
845      * The initial state must already contain the Cartesian State Transition Matrix.
846      * </p>
847      * @param columnName name of the column
848      * @param dYdQ column of the Jacobian ∂Y/∂qₘ with respect to propagation type,
849      * if null (which is the most frequent case), assumed to be 0
850      * @since 11.1
851      */
852     private void setInitialColumn(final String columnName, final double[] dYdQ) {
853 
854         final SpacecraftState state = getInitialState();
855 
856         final AbstractStateTransitionMatrixGenerator generator = (AbstractStateTransitionMatrixGenerator)
857                 getAdditionalDerivativesProviders().stream()
858                         .filter(AbstractStateTransitionMatrixGenerator.class::isInstance)
859                         .collect(Collectors.toList()).get(0);
860         final int expectedSize = generator.getStateDimension();
861         if (dYdQ.length != expectedSize) {
862             throw new OrekitException(LocalizedCoreFormats.DIMENSIONS_MISMATCH, dYdQ.length, expectedSize);
863         }
864 
865         // convert to Cartesian Jacobian
866         final RealMatrix dYdC = MatrixUtils.createRealIdentityMatrix(expectedSize);
867         final double[][] jacobian = new double[6][6];
868         getOrbitType().convertType(state.getOrbit()).getJacobianWrtCartesian(getPositionAngleType(), jacobian);
869         dYdC.setSubMatrix(jacobian, 0, 0);
870         final DecompositionSolver solver = getSolver(dYdC);
871         final double[] column = solver.solve(MatrixUtils.createRealVector(dYdQ)).toArray();
872 
873         // set additional state
874         setInitialState(state.addAdditionalData(columnName, column));
875 
876     }
877 
878     /**
879      * Method to get a linear system solver.
880      * @param matrix matrix involved in linear systems
881      * @return solver
882      * @since 13.1
883      */
884     private DecompositionSolver getSolver(final RealMatrix matrix) {
885         return new QRDecomposition(matrix, THRESHOLD).getSolver();
886     }
887 
888     /** {@inheritDoc} */
889     @Override
890     protected AttitudeProvider initializeAttitudeProviderForDerivatives() {
891         return needFullAttitudeForDerivatives ? getAttitudeProvider() : getFrozenAttitudeProvider();
892     }
893 
894     /** {@inheritDoc} */
895     @Override
896     protected StateMapper createMapper(final AbsoluteDate referenceDate, final double mu,
897                                        final OrbitType orbitType, final PositionAngleType positionAngleType,
898                                        final AttitudeProvider attitudeProvider, final Frame frame) {
899         return new OsculatingMapper(referenceDate, mu, orbitType, positionAngleType, attitudeProvider, frame);
900     }
901 
902     /** Internal mapper using directly osculating parameters. */
903     private static class OsculatingMapper extends StateMapper {
904 
905         /** Simple constructor.
906          * <p>
907          * The position parameter type is meaningful only if {@link
908          * #getOrbitType() propagation orbit type}
909          * support it. As an example, it is not meaningful for propagation
910          * in {@link OrbitType#CARTESIAN Cartesian} parameters.
911          * </p>
912          * @param referenceDate reference date
913          * @param mu central attraction coefficient (m³/s²)
914          * @param orbitType orbit type to use for mapping (can be null for {@link AbsolutePVCoordinates})
915          * @param positionAngleType angle type to use for propagation
916          * @param attitudeProvider attitude provider
917          * @param frame inertial frame
918          */
919         OsculatingMapper(final AbsoluteDate referenceDate, final double mu,
920                          final OrbitType orbitType, final PositionAngleType positionAngleType,
921                          final AttitudeProvider attitudeProvider, final Frame frame) {
922             super(referenceDate, mu, orbitType, positionAngleType, attitudeProvider, frame);
923         }
924 
925         /** {@inheritDoc} */
926         @Override
927         public SpacecraftState mapArrayToState(final AbsoluteDate date, final double[] y, final double[] yDot,
928                                                final PropagationType type) {
929             // the parameter type is ignored for the Numerical Propagator
930 
931             final double mass = y[6];
932             if (mass <= 0.0) {
933                 throw new OrekitException(OrekitMessages.NOT_POSITIVE_SPACECRAFT_MASS, mass);
934             }
935 
936             if (getOrbitType() == null) {
937                 // propagation uses absolute position-velocity-acceleration
938                 final Vector3D p = new Vector3D(y[0],    y[1],    y[2]);
939                 final Vector3D v = new Vector3D(y[3],    y[4],    y[5]);
940                 final Vector3D a;
941                 final AbsolutePVCoordinates absPva;
942                 if (yDot == null) {
943                     absPva = new AbsolutePVCoordinates(getFrame(), new TimeStampedPVCoordinates(date, p, v));
944                 } else {
945                     a = new Vector3D(yDot[3], yDot[4], yDot[5]);
946                     absPva = new AbsolutePVCoordinates(getFrame(), new TimeStampedPVCoordinates(date, p, v, a));
947                 }
948 
949                 final Attitude attitude = getAttitudeProvider().getAttitude(absPva, date, getFrame());
950                 return new SpacecraftState(absPva, attitude).withMass(mass);
951             } else {
952                 // propagation uses regular orbits
953                 final Orbit orbit       = getOrbitType().mapArrayToOrbit(y, yDot, getPositionAngleType(), date, getMu(), getFrame());
954                 final Attitude attitude = getAttitudeProvider().getAttitude(orbit, date, getFrame());
955 
956                 return new SpacecraftState(orbit, attitude).withMass(mass);
957             }
958 
959         }
960 
961         /** {@inheritDoc} */
962         public void mapStateToArray(final SpacecraftState state, final double[] y, final double[] yDot) {
963             if (getOrbitType() == null) {
964                 // propagation uses absolute position-velocity-acceleration
965                 final Vector3D p = state.getAbsPVA().getPosition();
966                 final Vector3D v = state.getAbsPVA().getVelocity();
967                 y[0] = p.getX();
968                 y[1] = p.getY();
969                 y[2] = p.getZ();
970                 y[3] = v.getX();
971                 y[4] = v.getY();
972                 y[5] = v.getZ();
973                 y[6] = state.getMass();
974             }
975             else {
976                 getOrbitType().mapOrbitToArray(state.getOrbit(), getPositionAngleType(), y, yDot);
977                 y[6] = state.getMass();
978             }
979         }
980 
981     }
982 
983     /** {@inheritDoc} */
984     protected MainStateEquations getMainStateEquations(final ODEIntegrator integrator) {
985         return new Main(integrator, getOrbitType(), getPositionAngleType(), getAllForceModels());
986     }
987 
988     /** Internal class for osculating parameters integration. */
989     private class Main extends NumericalTimeDerivativesEquations implements MainStateEquations {
990 
991         /** Flag keeping track whether Jacobian matrix needs to be recomputed or not. */
992         private final boolean recomputingJacobian;
993 
994         /** Simple constructor.
995          * @param integrator numerical integrator to use for propagation.
996          * @param orbitType orbit type
997          * @param positionAngleType angle type
998          * @param forceModelList forces
999          */
1000         Main(final ODEIntegrator integrator, final OrbitType orbitType, final PositionAngleType positionAngleType,
1001              final List<ForceModel> forceModelList) {
1002 
1003             super(orbitType, positionAngleType, forceModelList);
1004             final int numberOfForces = forceModelList.size();
1005             if (orbitType != null && orbitType != OrbitType.CARTESIAN && numberOfForces > 0) {
1006                 if (numberOfForces > 1) {
1007                     recomputingJacobian = true;
1008                 } else {
1009                     recomputingJacobian = !(forceModelList.get(0) instanceof NewtonianAttraction);
1010                 }
1011             } else {
1012                 recomputingJacobian = false;
1013             }
1014 
1015             // feed internal event detectors
1016             setUpInternalDetectors(integrator);
1017 
1018         }
1019 
1020         /** Set up all user defined event detectors.
1021          * @param integrator numerical integrator to use for propagation.
1022          */
1023         private void setUpInternalDetectors(final ODEIntegrator integrator) {
1024             final NumericalTimeDerivativesEquations cartesianEquations = new NumericalTimeDerivativesEquations(OrbitType.CARTESIAN,
1025                     null, forceModels);
1026             for (final ForceModel forceModel : getForceModels()) {
1027                 forceModel.getEventDetectors().forEach(detector -> setUpInternalEventDetector(integrator, detector,
1028                         cartesianEquations));
1029             }
1030             getAttitudeProvider().getEventDetectors().forEach(detector -> setUpInternalEventDetector(integrator,
1031                     detector, cartesianEquations));
1032         }
1033 
1034         /** Set up one internal event detector.
1035          * @param integrator numerical integrator to use for propagation.
1036          * @param eventDetector detector
1037          * @param cartesianEquations Cartesian derivatives model
1038          */
1039         private void setUpInternalEventDetector(final ODEIntegrator integrator,
1040                                                 final EventDetector eventDetector,
1041                                                 final NumericalTimeDerivativesEquations cartesianEquations) {
1042             final Optional<ExtendedStateTransitionMatrixGenerator> generator = getAdditionalDerivativesProviders().stream()
1043                     .filter(ExtendedStateTransitionMatrixGenerator.class::isInstance)
1044                     .map(ExtendedStateTransitionMatrixGenerator.class::cast)
1045                     .findFirst();
1046             final EventDetector internalDetector;
1047             if (generator.isPresent() && !eventDetector.dependsOnTimeOnly()) {
1048                 // need to modify STM at each dynamics discontinuities
1049                 final NumericalPropagationHarvester harvester = (NumericalPropagationHarvester) getHarvester();
1050                 final SwitchEventHandler handler = new SwitchEventHandler(eventDetector.getHandler(), harvester,
1051                         cartesianEquations, getAttitudeProvider());
1052                 internalDetector = getLocalDetector(eventDetector, handler);
1053             } else {
1054                 internalDetector = eventDetector;
1055             }
1056             setUpEventDetector(integrator, internalDetector);
1057         }
1058 
1059         /** {@inheritDoc} */
1060         @Override
1061         public void init(final SpacecraftState initialState, final AbsoluteDate target) {
1062             final List<ForceModel> forceModelList = getForceModels();
1063             needFullAttitudeForDerivatives = forceModelList.stream().anyMatch(ForceModel::dependsOnAttitudeRate);
1064 
1065             forceModelList.forEach(fm -> fm.init(initialState, target));
1066 
1067         }
1068 
1069         /** {@inheritDoc} */
1070         @Override
1071         public double[] computeDerivatives(final SpacecraftState state) {
1072             setCurrentState(state);
1073             if (recomputingJacobian) {
1074                 // propagation uses Jacobian matrix of orbital parameters w.r.t. Cartesian ones
1075                 final double[][] jacobian = new double[6][6];
1076                 state.getOrbit().getJacobianWrtCartesian(getPositionAngleType(), jacobian);
1077                 setCoordinatesJacobian(jacobian);
1078             }
1079             return computeTimeDerivatives(state);
1080 
1081         }
1082 
1083     }
1084 
1085     /** Estimate tolerance vectors for integrators when propagating in absolute position-velocity-acceleration.
1086      * @param dP user specified position error
1087      * @param absPva reference absolute position-velocity-acceleration
1088      * @return a two rows array, row 0 being the absolute tolerance error and row 1
1089      * being the relative tolerance error
1090      * @deprecated since 13.0. Use {@link ToleranceProvider} for default and custom tolerances.
1091      */
1092     @Deprecated
1093     public static double[][] tolerances(final double dP, final AbsolutePVCoordinates absPva) {
1094         return ToleranceProvider.of(CartesianToleranceProvider.of(dP)).getTolerances(absPva);
1095     }
1096 
1097     /** Estimate tolerance vectors for integrators when propagating in orbits.
1098      * <p>
1099      * The errors are estimated from partial derivatives properties of orbits,
1100      * starting from a scalar position error specified by the user.
1101      * Considering the energy conservation equation V = sqrt(mu (2/r - 1/a)),
1102      * we get at constant energy (i.e. on a Keplerian trajectory):
1103      * <pre>
1104      * V r² |dV| = mu |dr|
1105      * </pre>
1106      * <p> So we deduce a scalar velocity error consistent with the position error.
1107      * From here, we apply orbits Jacobians matrices to get consistent errors
1108      * on orbital parameters.
1109      *
1110      * <p>
1111      * The tolerances are only <em>orders of magnitude</em>, and integrator tolerances
1112      * are only local estimates, not global ones. So some care must be taken when using
1113      * these tolerances. Setting 1mm as a position error does NOT mean the tolerances
1114      * will guarantee a 1mm error position after several orbits integration.
1115      * </p>
1116      * @param dP user specified position error
1117      * @param orbit reference orbit
1118      * @param type propagation type for the meaning of the tolerance vectors elements
1119      * (it may be different from {@code orbit.getType()})
1120      * @return a two rows array, row 0 being the absolute tolerance error and row 1
1121      * being the relative tolerance error
1122      * @deprecated since 13.0. Use {@link ToleranceProvider} for default and custom tolerances.
1123      */
1124     @Deprecated
1125     public static double[][] tolerances(final double dP, final Orbit orbit, final OrbitType type) {
1126         return ToleranceProvider.getDefaultToleranceProvider(dP).getTolerances(orbit, type, PositionAngleType.TRUE);
1127     }
1128 
1129     /** Estimate tolerance vectors for integrators when propagating in orbits.
1130      * <p>
1131      * The errors are estimated from partial derivatives properties of orbits,
1132      * starting from scalar position and velocity errors specified by the user.
1133      * <p>
1134      * The tolerances are only <em>orders of magnitude</em>, and integrator tolerances
1135      * are only local estimates, not global ones. So some care must be taken when using
1136      * these tolerances. Setting 1mm as a position error does NOT mean the tolerances
1137      * will guarantee a 1mm error position after several orbits integration.
1138      * </p>
1139      * @param dP user specified position error
1140      * @param dV user specified velocity error
1141      * @param orbit reference orbit
1142      * @param type propagation type for the meaning of the tolerance vectors elements
1143      * (it may be different from {@code orbit.getType()})
1144      * @return a two rows array, row 0 being the absolute tolerance error and row 1
1145      * being the relative tolerance error
1146      * @since 10.3
1147      * @deprecated since 13.0. Use {@link ToleranceProvider} for default and custom tolerances.
1148      */
1149     @Deprecated
1150     public static double[][] tolerances(final double dP, final double dV,
1151                                         final Orbit orbit, final OrbitType type) {
1152 
1153         return ToleranceProvider.of(CartesianToleranceProvider.of(dP, dV, CartesianToleranceProvider.DEFAULT_ABSOLUTE_MASS_TOLERANCE))
1154                 .getTolerances(orbit, type, PositionAngleType.TRUE);
1155     }
1156 
1157     /** {@inheritDoc} */
1158     @Override
1159     protected void beforeIntegration(final SpacecraftState initialState, final AbsoluteDate tEnd) {
1160 
1161         if (!getFrame().isPseudoInertial()) {
1162 
1163             // inspect all force models to find InertialForces
1164             for (ForceModel force : forceModels) {
1165                 if (force instanceof InertialForces) {
1166                     return;
1167                 }
1168             }
1169 
1170             // throw exception if no inertial forces found
1171             throw new OrekitException(OrekitMessages.INERTIAL_FORCE_MODEL_MISSING, getFrame().getName());
1172 
1173         }
1174 
1175     }
1176 
1177     /**
1178      * Creates local detector wrapping input one and using specific handler for dynamics discontinuities and STM.
1179      * @param eventDetector detector
1180      * @param switchEventHandler special handler
1181      * @return wrapped detector
1182      */
1183     private static EventDetector getLocalDetector(final EventDetector eventDetector,
1184                                                   final SwitchEventHandler switchEventHandler) {
1185         return new DetectorModifier() {
1186             @Override
1187             public EventDetector getDetector() {
1188                 return eventDetector;
1189             }
1190 
1191             @Override
1192             public EventHandler getHandler() {
1193                 return switchEventHandler;
1194             }
1195         };
1196     }
1197 }