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éro
170 * @author Luc Maisonobe
171 * @author Guylaine Prat
172 * @author Fabien Maussion
173 * @author Vé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 }