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.integration;
18  
19  import org.hipparchus.analysis.UnivariateFunction;
20  import org.hipparchus.analysis.solvers.BracketedUnivariateSolver;
21  import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
22  import org.hipparchus.exception.MathRuntimeException;
23  import org.hipparchus.ode.DenseOutputModel;
24  import org.hipparchus.ode.ExpandableODE;
25  import org.hipparchus.ode.ODEIntegrator;
26  import org.hipparchus.ode.ODEState;
27  import org.hipparchus.ode.ODEStateAndDerivative;
28  import org.hipparchus.ode.OrdinaryDifferentialEquation;
29  import org.hipparchus.ode.SecondaryODE;
30  import org.hipparchus.ode.events.Action;
31  import org.hipparchus.ode.events.AdaptableInterval;
32  import org.hipparchus.ode.events.ODEEventDetector;
33  import org.hipparchus.ode.events.ODEEventHandler;
34  import org.hipparchus.ode.sampling.AbstractODEStateInterpolator;
35  import org.hipparchus.ode.sampling.ODEStateInterpolator;
36  import org.hipparchus.ode.sampling.ODEStepHandler;
37  import org.hipparchus.util.Precision;
38  import org.orekit.attitudes.AttitudeProvider;
39  import org.orekit.errors.OrekitException;
40  import org.orekit.errors.OrekitInternalError;
41  import org.orekit.errors.OrekitMessages;
42  import org.orekit.frames.Frame;
43  import org.orekit.orbits.OrbitType;
44  import org.orekit.orbits.PositionAngleType;
45  import org.orekit.propagation.AbstractPropagator;
46  import org.orekit.propagation.BoundedPropagator;
47  import org.orekit.propagation.EphemerisGenerator;
48  import org.orekit.propagation.PropagationType;
49  import org.orekit.propagation.SpacecraftState;
50  import org.orekit.propagation.events.EventDetector;
51  import org.orekit.propagation.events.handlers.EventHandler;
52  import org.orekit.propagation.sampling.OrekitStepHandler;
53  import org.orekit.propagation.sampling.OrekitStepInterpolator;
54  import org.orekit.time.AbsoluteDate;
55  import org.orekit.utils.DataDictionary;
56  
57  import java.util.ArrayList;
58  import java.util.Arrays;
59  import java.util.Collection;
60  import java.util.Collections;
61  import java.util.HashMap;
62  import java.util.LinkedList;
63  import java.util.List;
64  import java.util.Map;
65  import java.util.Queue;
66  
67  
68  /** Common handling of {@link org.orekit.propagation.Propagator Propagator}
69   *  methods for both numerical and semi-analytical propagators.
70   *  @author Luc Maisonobe
71   */
72  public abstract class AbstractIntegratedPropagator extends AbstractPropagator {
73  
74      /** Internal name used for complete secondary state dimension.
75       * @since 11.1
76       */
77      private static final String SECONDARY_DIMENSION = "Orekit-secondary-dimension";
78  
79      /** Event detectors not related to force models. */
80      private final List<EventDetector> detectors;
81  
82      /** Step handlers dedicated to ephemeris generation. */
83      private final List<StoringStepHandler> ephemerisGenerators;
84  
85      /** Integrator selected by the user for the orbital extrapolation process. */
86      private final ODEIntegrator integrator;
87  
88      /** Offsets of secondary states managed by {@link AdditionalDerivativesProvider}.
89       * @since 11.1
90       */
91      private final Map<String, Integer> secondaryOffsets;
92  
93      /** Additional derivatives providers.
94       * @since 11.1
95       */
96      private final List<AdditionalDerivativesProvider> additionalDerivativesProviders;
97  
98      /** Map of secondary equation offset in main
99      /** Counter for differential equations calls. */
100     private int calls;
101 
102     /** Mapper between raw double components and space flight dynamics objects. */
103     private StateMapper stateMapper;
104 
105     /**
106      * Attitude provider when evaluating derivatives. Can be a frozen one for performance.
107      * @since 12.1
108      */
109     private AttitudeProvider attitudeProviderForDerivatives;
110 
111     /** Flag for resetting the state at end of propagation. */
112     private boolean resetAtEnd;
113 
114     /** Type of orbit to output (mean or osculating) <br/>
115      * <p>
116      * This is used only in the case of semi-analytical propagators where there is a clear separation between
117      * mean and short periodic elements. It is ignored by the Numerical propagator.
118      * </p>
119      */
120     private final PropagationType propagationType;
121 
122     /** Build a new instance.
123      * @param integrator numerical integrator to use for propagation.
124      * @param propagationType type of orbit to output (mean or osculating).
125      */
126     protected AbstractIntegratedPropagator(final ODEIntegrator integrator, final PropagationType propagationType) {
127         detectors                      = new ArrayList<>();
128         ephemerisGenerators            = new ArrayList<>();
129         additionalDerivativesProviders = new ArrayList<>();
130         this.secondaryOffsets          = new HashMap<>();
131         this.integrator                = integrator;
132         this.propagationType           = propagationType;
133         this.resetAtEnd                = true;
134     }
135 
136     /** Allow/disallow resetting the initial state at end of propagation.
137      * <p>
138      * By default, at the end of the propagation, the propagator resets the initial state
139      * to the final state, thus allowing a new propagation to be started from there without
140      * recomputing the part already performed. Calling this method with {@code resetAtEnd} set
141      * to false changes prevents such reset.
142      * </p>
143      * @param resetAtEnd if true, at end of each propagation, the {@link
144      * #getInitialState() initial state} will be reset to the final state of
145      * the propagation, otherwise the initial state will be preserved
146      * @since 9.0
147      */
148     public void setResetAtEnd(final boolean resetAtEnd) {
149         this.resetAtEnd = resetAtEnd;
150     }
151 
152     /** Getter for the resetting flag regarding initial state.
153      * @return resetting flag
154      * @since 12.0
155      */
156     public boolean getResetAtEnd() {
157         return this.resetAtEnd;
158     }
159 
160     /**
161      * Method called when initializing the attitude provider used when evaluating derivatives.
162      * @return attitude provider for derivatives
163      */
164     protected AttitudeProvider initializeAttitudeProviderForDerivatives() {
165         return getAttitudeProvider();
166     }
167 
168     /** Initialize the mapper. */
169     protected void initMapper() {
170         stateMapper = createMapper(null, Double.NaN, null, null, null, null);
171     }
172 
173     /** Get the integrator's name.
174      * @return name of underlying integrator
175      * @since 12.0
176      */
177     public String getIntegratorName() {
178         return integrator.getName();
179     }
180 
181     /**  {@inheritDoc} */
182     @Override
183     public void setAttitudeProvider(final AttitudeProvider attitudeProvider) {
184         super.setAttitudeProvider(attitudeProvider);
185         stateMapper = createMapper(stateMapper.getReferenceDate(), stateMapper.getMu(),
186                                    stateMapper.getOrbitType(), stateMapper.getPositionAngleType(),
187                                    attitudeProvider, stateMapper.getFrame());
188     }
189 
190     /** Set propagation orbit type.
191      * @param orbitType orbit type to use for propagation, null for
192      * propagating using {@link org.orekit.utils.AbsolutePVCoordinates AbsolutePVCoordinates}
193      * rather than {@link org.orekit.orbits Orbit}
194      */
195     protected void setOrbitType(final OrbitType orbitType) {
196         stateMapper = createMapper(stateMapper.getReferenceDate(), stateMapper.getMu(),
197                                    orbitType, stateMapper.getPositionAngleType(),
198                                    stateMapper.getAttitudeProvider(), stateMapper.getFrame());
199     }
200 
201     /** Get propagation parameter type.
202      * @return orbit type used for propagation, null for
203      * propagating using {@link org.orekit.utils.AbsolutePVCoordinates AbsolutePVCoordinates}
204      * rather than {@link org.orekit.orbits Orbit}
205      */
206     protected OrbitType getOrbitType() {
207         return stateMapper.getOrbitType();
208     }
209 
210     /** Get the propagation type.
211      * @return propagation type.
212      * @since 11.1
213      */
214     public PropagationType getPropagationType() {
215         return propagationType;
216     }
217 
218     /** Set position angle type.
219      * <p>
220      * The position parameter type is meaningful only if {@link
221      * #getOrbitType() propagation orbit type}
222      * support it. As an example, it is not meaningful for propagation
223      * in {@link OrbitType#CARTESIAN Cartesian} parameters.
224      * </p>
225      * @param positionAngleType angle type to use for propagation
226      */
227     protected void setPositionAngleType(final PositionAngleType positionAngleType) {
228         stateMapper = createMapper(stateMapper.getReferenceDate(), stateMapper.getMu(),
229                                    stateMapper.getOrbitType(), positionAngleType,
230                                    stateMapper.getAttitudeProvider(), stateMapper.getFrame());
231     }
232 
233     /** Get propagation parameter type.
234      * @return angle type to use for propagation
235      */
236     protected PositionAngleType getPositionAngleType() {
237         return stateMapper.getPositionAngleType();
238     }
239 
240     /** Set the central attraction coefficient μ.
241      * @param mu central attraction coefficient (m³/s²)
242      */
243     public void setMu(final double mu) {
244         stateMapper = createMapper(stateMapper.getReferenceDate(), mu,
245                                    stateMapper.getOrbitType(), stateMapper.getPositionAngleType(),
246                                    stateMapper.getAttitudeProvider(), stateMapper.getFrame());
247     }
248 
249     /** Get the central attraction coefficient μ.
250      * @return mu central attraction coefficient (m³/s²)
251      * @see #setMu(double)
252      */
253     public double getMu() {
254         return stateMapper.getMu();
255     }
256 
257     /** Get the number of calls to the differential equations computation method.
258      * <p>The number of calls is reset each time the {@link #propagate(AbsoluteDate)}
259      * method is called.</p>
260      * @return number of calls to the differential equations computation method
261      */
262     public int getCalls() {
263         return calls;
264     }
265 
266     /** {@inheritDoc} */
267     @Override
268     public boolean isAdditionalDataManaged(final String name) {
269 
270         // first look at already integrated states
271         if (super.isAdditionalDataManaged(name)) {
272             return true;
273         }
274 
275         // then look at states we integrate ourselves
276         for (final AdditionalDerivativesProvider provider : additionalDerivativesProviders) {
277             if (provider.getName().equals(name)) {
278                 return true;
279             }
280         }
281 
282         return false;
283     }
284 
285     /** {@inheritDoc} */
286     @Override
287     public String[] getManagedAdditionalData() {
288         final String[] alreadyIntegrated = super.getManagedAdditionalData();
289         final String[] managed = new String[alreadyIntegrated.length + additionalDerivativesProviders.size()];
290         System.arraycopy(alreadyIntegrated, 0, managed, 0, alreadyIntegrated.length);
291         for (int i = 0; i < additionalDerivativesProviders.size(); ++i) {
292             managed[i + alreadyIntegrated.length] = additionalDerivativesProviders.get(i).getName();
293         }
294         return managed;
295     }
296 
297     /** Add a provider for user-specified state derivatives to be integrated along with the orbit propagation.
298      * @param provider provider for additional derivatives
299      * @see #addAdditionalDataProvider(org.orekit.propagation.AdditionalDataProvider)
300      * @since 11.1
301      */
302     public void addAdditionalDerivativesProvider(final AdditionalDerivativesProvider provider) {
303 
304         // check if the name is already used
305         if (this.isAdditionalDataManaged(provider.getName())) {
306             // these derivatives are already registered, complain
307             throw new OrekitException(OrekitMessages.ADDITIONAL_STATE_NAME_ALREADY_IN_USE,
308                                       provider.getName());
309         }
310 
311         // this is really a new set of derivatives, add it
312         additionalDerivativesProviders.add(provider);
313 
314         secondaryOffsets.clear();
315 
316     }
317 
318     /** Get an unmodifiable list of providers for additional derivatives.
319      * @return providers for the additional derivatives
320      * @since 11.1
321      */
322     public List<AdditionalDerivativesProvider> getAdditionalDerivativesProviders() {
323         return Collections.unmodifiableList(additionalDerivativesProviders);
324     }
325 
326     /** {@inheritDoc} */
327     public void addEventDetector(final EventDetector detector) {
328         detectors.add(detector);
329     }
330 
331     /** {@inheritDoc} */
332     public Collection<EventDetector> getEventDetectors() {
333         return Collections.unmodifiableCollection(detectors);
334     }
335 
336     /** {@inheritDoc} */
337     public void clearEventsDetectors() {
338         detectors.clear();
339     }
340 
341     /** Set up all user defined event detectors.
342      */
343     protected void setUpUserEventDetectors() {
344         for (final EventDetector detector : detectors) {
345             setUpEventDetector(integrator, detector);
346         }
347     }
348 
349     /** Wrap an Orekit event detector and register it to the integrator.
350      * @param integ integrator into which event detector should be registered
351      * @param detector event detector to wrap
352      */
353     protected void setUpEventDetector(final ODEIntegrator integ, final EventDetector detector) {
354         integ.addEventDetector(new AdaptedEventDetector(detector));
355     }
356 
357     /**
358      * Clear the ephemeris generators.
359      * @since 13.0
360      */
361     public void clearEphemerisGenerators() {
362         ephemerisGenerators.clear();
363     }
364 
365     /** {@inheritDoc} */
366     @Override
367     public EphemerisGenerator getEphemerisGenerator() {
368         final StoringStepHandler storingHandler = new StoringStepHandler();
369         ephemerisGenerators.add(storingHandler);
370         return storingHandler;
371     }
372 
373     /** Create a mapper between raw double components and spacecraft state.
374     /** Simple constructor.
375      * <p>
376      * The position parameter type is meaningful only if {@link
377      * #getOrbitType() propagation orbit type}
378      * support it. As an example, it is not meaningful for propagation
379      * in {@link OrbitType#CARTESIAN Cartesian} parameters.
380      * </p>
381      * @param referenceDate reference date
382      * @param mu central attraction coefficient (m³/s²)
383      * @param orbitType orbit type to use for mapping
384      * @param positionAngleType angle type to use for propagation
385      * @param attitudeProvider attitude provider
386      * @param frame inertial frame
387      * @return new mapper
388      */
389     protected abstract StateMapper createMapper(AbsoluteDate referenceDate, double mu,
390                                                 OrbitType orbitType, PositionAngleType positionAngleType,
391                                                 AttitudeProvider attitudeProvider, Frame frame);
392 
393     /** Get the differential equations to integrate (for main state only).
394      * @param integ numerical integrator to use for propagation.
395      * @return differential equations for main state
396      */
397     protected abstract MainStateEquations getMainStateEquations(ODEIntegrator integ);
398 
399     /** {@inheritDoc} */
400     @Override
401     public SpacecraftState propagate(final AbsoluteDate target) {
402         if (getStartDate() == null) {
403             if (getInitialState() == null) {
404                 throw new OrekitException(OrekitMessages.INITIAL_STATE_NOT_SPECIFIED_FOR_ORBIT_PROPAGATION);
405             }
406             setStartDate(getInitialState().getDate());
407         }
408         return propagate(getStartDate(), target);
409     }
410 
411     /** {@inheritDoc} */
412     public SpacecraftState propagate(final AbsoluteDate tStart, final AbsoluteDate tEnd) {
413 
414         if (getInitialState() == null) {
415             throw new OrekitException(OrekitMessages.INITIAL_STATE_NOT_SPECIFIED_FOR_ORBIT_PROPAGATION);
416         }
417 
418         // make sure the integrator will be reset properly even if we change its events handlers and step handlers
419         try (IntegratorResetter resetter = new IntegratorResetter(integrator)) {
420 
421             // prepare handling of STM and Jacobian matrices
422             setUpStmAndJacobianGenerators();
423 
424             // Initialize additional data
425             initializeAdditionalData(tEnd);
426 
427             if (!tStart.equals(getInitialState().getDate())) {
428                 // if propagation start date is not initial date,
429                 // propagate from initial to start date without event detection
430                 try (IntegratorResetter startResetter = new IntegratorResetter(integrator)) {
431                     integrateDynamics(tStart, true);
432                 }
433             }
434 
435             // set up events added by user
436             setUpUserEventDetectors();
437 
438             // set up step handlers
439             for (final OrekitStepHandler handler : getMultiplexer().getHandlers()) {
440                 integrator.addStepHandler(new AdaptedStepHandler(handler));
441             }
442             for (final StoringStepHandler generator : ephemerisGenerators) {
443                 generator.setEndDate(tEnd);
444                 integrator.addStepHandler(generator);
445             }
446 
447             // propagate from start date to end date with event detection
448             final SpacecraftState finalState = integrateDynamics(tEnd, false);
449 
450             // Finalize event detectors
451             getEventDetectors().forEach(detector -> detector.finish(finalState));
452 
453             return finalState;
454         }
455 
456     }
457 
458     /** Reset initial state with a given propagation type.
459      *
460      * <p> By default this method returns the same as {@link #resetInitialState(SpacecraftState)}
461      * <p> Its purpose is mostly to be derived in DSSTPropagator
462      *
463      * @param state new initial state to consider
464      * @param stateType type of the new state (mean or osculating)
465      * @since 12.1.3
466      */
467     public void resetInitialState(final SpacecraftState state, final PropagationType stateType) {
468         // Default behavior, do not take propagation type into account
469         resetInitialState(state);
470     }
471 
472     /** Set up State Transition Matrix and Jacobian matrix handling.
473      * @since 11.1
474      */
475     protected void setUpStmAndJacobianGenerators() {
476         // nothing to do by default
477     }
478 
479     /** Propagation with or without event detection.
480      * @param tEnd target date to which orbit should be propagated
481      * @param forceResetAtEnd flag to force resetting state and date after integration
482      * @return state at end of propagation
483      */
484     private SpacecraftState integrateDynamics(final AbsoluteDate tEnd, final boolean forceResetAtEnd) {
485         try {
486 
487             initializePropagation();
488 
489             if (getInitialState().getDate().equals(tEnd)) {
490                 // don't extrapolate
491                 return getInitialState();
492             }
493 
494             // space dynamics view
495             stateMapper = createMapper(getInitialState().getDate(), stateMapper.getMu(),
496                                        stateMapper.getOrbitType(), stateMapper.getPositionAngleType(),
497                                        stateMapper.getAttitudeProvider(), getInitialState().getFrame());
498             attitudeProviderForDerivatives = initializeAttitudeProviderForDerivatives();
499 
500             if (Double.isNaN(getMu())) {
501                 setMu(getInitialState().isOrbitDefined() ? getInitialState().getOrbit().getMu() : Double.NaN);
502             }
503 
504             if (getInitialState().getMass() <= 0.0) {
505                 throw new OrekitException(OrekitMessages.NOT_POSITIVE_SPACECRAFT_MASS,
506                                           getInitialState().getMass());
507             }
508 
509             // convert space flight dynamics API to math API
510             final SpacecraftState initialIntegrationState = getInitialIntegrationState();
511             final ODEState mathInitialState = createInitialState(initialIntegrationState);
512             final ExpandableODE mathODE = createODE(integrator);
513 
514             // mathematical integration
515             final ODEStateAndDerivative mathFinalState;
516             beforeIntegration(initialIntegrationState, tEnd);
517             mathFinalState = integrator.integrate(mathODE, mathInitialState,
518                                                   tEnd.durationFrom(getInitialState().getDate()));
519             afterIntegration();
520 
521             // get final state
522             SpacecraftState finalState =
523                             stateMapper.mapArrayToState(stateMapper.mapDoubleToDate(mathFinalState.getTime(),
524                                                                                     tEnd),
525                                                         mathFinalState.getPrimaryState(),
526                                                         mathFinalState.getPrimaryDerivative(),
527                                                         propagationType);
528 
529             finalState = updateAdditionalStatesAndDerivatives(finalState, mathFinalState);
530 
531             if (resetAtEnd || forceResetAtEnd) {
532                 resetInitialState(finalState, propagationType);
533                 setStartDate(finalState.getDate());
534             }
535 
536             return finalState;
537 
538         } catch (MathRuntimeException mre) {
539             throw OrekitException.unwrap(mre);
540         }
541     }
542 
543     /**
544      * Returns an updated version of the inputted state with additional states, including
545      * from derivatives providers.
546      * @param originalState input state
547      * @param os ODE state and derivative
548      * @return new state
549      * @since 12.1
550      */
551     private SpacecraftState updateAdditionalStatesAndDerivatives(final SpacecraftState originalState,
552                                                                  final ODEStateAndDerivative os) {
553         SpacecraftState updatedState = originalState;
554         if (os.getNumberOfSecondaryStates() > 0) {
555             final double[] secondary           = os.getSecondaryState(1);
556             final double[] secondaryDerivative = os.getSecondaryDerivative(1);
557             for (final AdditionalDerivativesProvider provider : additionalDerivativesProviders) {
558                 final String name      = provider.getName();
559                 final int    offset    = secondaryOffsets.get(name);
560                 final int    dimension = provider.getDimension();
561                 updatedState = updatedState.addAdditionalData(name, Arrays.copyOfRange(secondary, offset, offset + dimension));
562                 updatedState = updatedState.addAdditionalStateDerivative(name, Arrays.copyOfRange(secondaryDerivative, offset, offset + dimension));
563             }
564         }
565         return updateAdditionalData(updatedState);
566     }
567 
568     /** Get the initial state for integration.
569      * @return initial state for integration
570      */
571     protected SpacecraftState getInitialIntegrationState() {
572         return getInitialState();
573     }
574 
575     /** Create an initial state.
576      * @param initialState initial state in flight dynamics world
577      * @return initial state in mathematics world
578      */
579     private ODEState createInitialState(final SpacecraftState initialState) {
580 
581         // retrieve initial state
582         final double[] primary = new double[getBasicDimension()];
583         stateMapper.mapStateToArray(initialState, primary, null);
584 
585         if (secondaryOffsets.isEmpty()) {
586             // compute dimension of the secondary state
587             int offset = 0;
588             for (final AdditionalDerivativesProvider provider : additionalDerivativesProviders) {
589                 secondaryOffsets.put(provider.getName(), offset);
590                 offset += provider.getDimension();
591             }
592             secondaryOffsets.put(SECONDARY_DIMENSION, offset);
593         }
594 
595         return new ODEState(0.0, primary, secondary(initialState));
596 
597     }
598 
599     /** Create secondary state.
600      * @param state spacecraft state
601      * @return secondary state
602      * @since 11.1
603      */
604     private double[][] secondary(final SpacecraftState state) {
605 
606         if (secondaryOffsets.isEmpty()) {
607             return null;
608         }
609 
610         final double[][] secondary = new double[1][secondaryOffsets.get(SECONDARY_DIMENSION)];
611         for (final AdditionalDerivativesProvider provider : additionalDerivativesProviders) {
612             final String   name       = provider.getName();
613             final int      offset     = secondaryOffsets.get(name);
614             final double[] additional = state.getAdditionalState(name);
615             System.arraycopy(additional, 0, secondary[0], offset, additional.length);
616         }
617 
618         return secondary;
619 
620     }
621 
622     /** Create secondary state derivative.
623      * @param state spacecraft state
624      * @return secondary state derivative
625      * @since 11.1
626      */
627     private double[][] secondaryDerivative(final SpacecraftState state) {
628 
629         if (secondaryOffsets.isEmpty()) {
630             return null;
631         }
632 
633         final double[][] secondaryDerivative = new double[1][secondaryOffsets.get(SECONDARY_DIMENSION)];
634         for (final AdditionalDerivativesProvider provider : additionalDerivativesProviders) {
635             final String   name       = provider.getName();
636             final int      offset     = secondaryOffsets.get(name);
637             final double[] additionalDerivative = state.getAdditionalStateDerivative(name);
638             System.arraycopy(additionalDerivative, 0, secondaryDerivative[0], offset, additionalDerivative.length);
639         }
640 
641         return secondaryDerivative;
642 
643     }
644 
645     /** Create an ODE with all equations.
646      * @param integ numerical integrator to use for propagation.
647      * @return a new ode
648      */
649     private ExpandableODE createODE(final ODEIntegrator integ) {
650 
651         final ExpandableODE ode =
652                 new ExpandableODE(new ConvertedMainStateEquations(getMainStateEquations(integ)));
653 
654         // secondary part of the ODE
655         if (!additionalDerivativesProviders.isEmpty()) {
656             ode.addSecondaryEquations(new ConvertedSecondaryStateEquations());
657         }
658 
659         return ode;
660 
661     }
662 
663     /** Method called just before integration.
664      * <p>
665      * The default implementation does nothing, it may be specialized in subclasses.
666      * </p>
667      * @param initialState initial state
668      * @param tEnd target date at which state should be propagated
669      */
670     protected void beforeIntegration(final SpacecraftState initialState,
671                                      final AbsoluteDate tEnd) {
672         // do nothing by default
673     }
674 
675     /** Method called just after integration.
676      * <p>
677      * The default implementation does nothing, it may be specialized in subclasses.
678      * </p>
679      */
680     protected void afterIntegration() {
681         // do nothing by default
682     }
683 
684     /** Get state vector dimension without additional parameters.
685      * @return state vector dimension without additional parameters.
686      */
687     public int getBasicDimension() {
688         return 7;
689     }
690 
691     /** Get the integrator used by the propagator.
692      * @return the integrator.
693      */
694     protected ODEIntegrator getIntegrator() {
695         return integrator;
696     }
697 
698     /** Convert a state from mathematical world to space flight dynamics world.
699      * @param os mathematical state
700      * @return space flight dynamics state
701      */
702     private SpacecraftState convert(final ODEStateAndDerivative os) {
703 
704         final SpacecraftState s = stateMapper.mapArrayToState(os.getTime(), os.getPrimaryState(),
705             os.getPrimaryDerivative(), propagationType);
706         return updateAdditionalStatesAndDerivatives(s, os);
707     }
708 
709     /** Convert a state from space flight dynamics world to mathematical world.
710      * @param state space flight dynamics state
711      * @return mathematical state
712      */
713     private ODEStateAndDerivative convert(final SpacecraftState state) {
714 
715         // retrieve initial state
716         final double[] primary    = new double[getBasicDimension()];
717         final double[] primaryDot = new double[getBasicDimension()];
718         stateMapper.mapStateToArray(state, primary, primaryDot);
719 
720         // secondary part of the ODE
721         final double[][] secondary           = secondary(state);
722         final double[][] secondaryDerivative = secondaryDerivative(state);
723 
724         return new ODEStateAndDerivative(stateMapper.mapDateToDouble(state.getDate()),
725                                          primary, primaryDot,
726                                          secondary, secondaryDerivative);
727 
728     }
729 
730     /** Differential equations for the main state (orbit, attitude and mass). */
731     public interface MainStateEquations {
732 
733         /**
734          * Initialize the equations at the start of propagation. This method will be
735          * called before any calls to {@link #computeDerivatives(SpacecraftState)}.
736          *
737          * <p> The default implementation of this method does nothing.
738          *
739          * @param initialState initial state information at the start of propagation.
740          * @param target       date of propagation. Not equal to {@code
741          *                     initialState.getDate()}.
742          */
743         default void init(final SpacecraftState initialState, final AbsoluteDate target) {
744         }
745 
746         /** Compute differential equations for main state.
747          * @param state current state
748          * @return derivatives of main state
749          */
750         double[] computeDerivatives(SpacecraftState state);
751 
752     }
753 
754     /** Differential equations for the main state (orbit, attitude and mass), with converted API. */
755     private class ConvertedMainStateEquations implements OrdinaryDifferentialEquation {
756 
757         /** Main state equations. */
758         private final MainStateEquations main;
759 
760         /** Simple constructor.
761          * @param main main state equations
762          */
763         ConvertedMainStateEquations(final MainStateEquations main) {
764             this.main = main;
765             calls = 0;
766         }
767 
768         /** {@inheritDoc} */
769         public int getDimension() {
770             return getBasicDimension();
771         }
772 
773         @Override
774         public void init(final double t0, final double[] y0, final double finalTime) {
775             // update space dynamics view
776             SpacecraftState initialState = stateMapper.mapArrayToState(t0, y0, null, PropagationType.MEAN);
777             initialState = updateAdditionalData(initialState);
778             initialState = updateStatesFromAdditionalDerivativesIfKnown(initialState);
779             final AbsoluteDate target = stateMapper.mapDoubleToDate(finalTime);
780             main.init(initialState, target);
781             attitudeProviderForDerivatives = initializeAttitudeProviderForDerivatives();
782         }
783 
784         /**
785          * Returns an updated version of the inputted state, with additional states from
786          * derivatives providers as given in the stored initial state.
787          * @param originalState input state
788          * @return new state
789          * @since 12.1
790          */
791         private SpacecraftState updateStatesFromAdditionalDerivativesIfKnown(final SpacecraftState originalState) {
792             SpacecraftState updatedState = originalState;
793             final SpacecraftState storedInitialState = getInitialState();
794             final double originalTime = stateMapper.mapDateToDouble(originalState.getDate());
795             if (storedInitialState != null && stateMapper.mapDateToDouble(storedInitialState.getDate()) == originalTime) {
796                 for (final AdditionalDerivativesProvider provider: additionalDerivativesProviders) {
797                     final String name = provider.getName();
798                     final double[] value = storedInitialState.getAdditionalState(name);
799                     updatedState = updatedState.addAdditionalData(name, value);
800                 }
801             }
802             return updatedState;
803         }
804 
805         /** {@inheritDoc} */
806         public double[] computeDerivatives(final double t, final double[] y) {
807 
808             // increment calls counter
809             ++calls;
810 
811             // update space dynamics view
812             stateMapper.setAttitudeProvider(attitudeProviderForDerivatives);
813             SpacecraftState currentState = stateMapper.mapArrayToState(t, y, null, PropagationType.MEAN);
814             stateMapper.setAttitudeProvider(getAttitudeProvider());
815 
816             currentState = updateAdditionalData(currentState);
817             // compute main state differentials
818             return main.computeDerivatives(currentState);
819 
820         }
821 
822     }
823 
824     /** Differential equations for the secondary state (Jacobians, user variables ...), with converted API. */
825     private class ConvertedSecondaryStateEquations implements SecondaryODE {
826 
827         /** Dimension of the combined additional states. */
828         private final int combinedDimension;
829 
830         /** Simple constructor.
831           */
832         ConvertedSecondaryStateEquations() {
833             this.combinedDimension = secondaryOffsets.get(SECONDARY_DIMENSION);
834         }
835 
836         /** {@inheritDoc} */
837         @Override
838         public int getDimension() {
839             return combinedDimension;
840         }
841 
842         /** {@inheritDoc} */
843         @Override
844         public void init(final double t0, final double[] primary0,
845                          final double[] secondary0, final double finalTime) {
846             // update space dynamics view
847             final SpacecraftState initialState = convert(t0, primary0, null, secondary0);
848 
849             final AbsoluteDate target = stateMapper.mapDoubleToDate(finalTime);
850             for (final AdditionalDerivativesProvider provider : additionalDerivativesProviders) {
851                 provider.init(initialState, target);
852             }
853 
854         }
855 
856         /** {@inheritDoc} */
857         @Override
858         public double[] computeDerivatives(final double t, final double[] primary,
859                                            final double[] primaryDot, final double[] secondary) {
860 
861             // update space dynamics view
862             // the integrable generators generate method will be called here,
863             // according to the generators yield order
864             SpacecraftState updated = convert(t, primary, primaryDot, secondary);
865 
866             // set up queue for equations
867             final Queue<AdditionalDerivativesProvider> pending = new LinkedList<>(additionalDerivativesProviders);
868 
869             // gather the derivatives from all additional equations, taking care of dependencies
870             final double[] secondaryDot = new double[combinedDimension];
871             int yieldCount = 0;
872             while (!pending.isEmpty()) {
873                 final AdditionalDerivativesProvider provider = pending.remove();
874                 if (provider.yields(updated)) {
875                     // this provider has to wait for another one,
876                     // we put it again in the pending queue
877                     pending.add(provider);
878                     if (++yieldCount >= pending.size()) {
879                         // all pending providers yielded!, they probably need data not yet initialized
880                         // we let the propagation proceed, if these data are really needed right now
881                         // an appropriate exception will be triggered when caller tries to access them
882                         break;
883                     }
884                 } else {
885                     // we can use these equations right now
886                     final String              name           = provider.getName();
887                     final int                 offset         = secondaryOffsets.get(name);
888                     final int                 dimension      = provider.getDimension();
889                     final CombinedDerivatives derivatives    = provider.combinedDerivatives(updated);
890                     final double[]            additionalPart = derivatives.getAdditionalDerivatives();
891                     final double[]            mainPart       = derivatives.getMainStateDerivativesIncrements();
892                     System.arraycopy(additionalPart, 0, secondaryDot, offset, dimension);
893                     updated = updated.addAdditionalStateDerivative(name, additionalPart);
894                     if (mainPart != null) {
895                         // this equation does change the main state derivatives
896                         for (int i = 0; i < mainPart.length; ++i) {
897                             primaryDot[i] += mainPart[i];
898                         }
899                     }
900                     yieldCount = 0;
901                 }
902             }
903 
904             return secondaryDot;
905 
906         }
907 
908         /** Convert mathematical view to space view.
909          * @param t current value of the independent <I>time</I> variable
910          * @param primary array containing the current value of the primary state vector
911          * @param primaryDot array containing the derivative of the primary state vector
912          * @param secondary array containing the current value of the secondary state vector
913          * @return space view of the state
914          */
915         private SpacecraftState convert(final double t, final double[] primary,
916                                         final double[] primaryDot, final double[] secondary) {
917 
918             SpacecraftState initialState = stateMapper.mapArrayToState(t, primary, primaryDot, PropagationType.MEAN);
919 
920             for (final AdditionalDerivativesProvider provider : additionalDerivativesProviders) {
921                 final String name      = provider.getName();
922                 final int    offset    = secondaryOffsets.get(name);
923                 final int    dimension = provider.getDimension();
924                 initialState = initialState.addAdditionalData(name, Arrays.copyOfRange(secondary, offset, offset + dimension));
925             }
926 
927             return updateAdditionalData(initialState);
928 
929         }
930 
931     }
932 
933     /** Adapt an {@link org.orekit.propagation.events.EventDetector}
934      * to Hipparchus {@link org.hipparchus.ode.events.ODEEventDetector} interface.
935      * @author Fabien Maussion
936      */
937     private class AdaptedEventDetector implements ODEEventDetector {
938 
939         /** Underlying event detector. */
940         private final EventDetector detector;
941 
942         /** Underlying event handler.
943          * @since 12.0
944          */
945         private final EventHandler handler;
946 
947         /** Time of the previous call to g. */
948         private double lastT;
949 
950         /** Value from the previous call to g. */
951         private double lastG;
952 
953         /** Build a wrapped event detector.
954          * @param detector event detector to wrap
955         */
956         AdaptedEventDetector(final EventDetector detector) {
957             this.detector = detector;
958             this.handler  = detector.getHandler();
959             this.lastT    = Double.NaN;
960             this.lastG    = Double.NaN;
961         }
962 
963         /** {@inheritDoc} */
964         @Override
965         public AdaptableInterval getMaxCheckInterval() {
966             return (state, isForward) -> detector.getMaxCheckInterval().currentInterval(convert(state), isForward);
967         }
968 
969         /** {@inheritDoc} */
970         @Override
971         public int getMaxIterationCount() {
972             return detector.getMaxIterationCount();
973         }
974 
975         /** {@inheritDoc} */
976         @Override
977         public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
978             return new BracketingNthOrderBrentSolver(0, detector.getThreshold(), 0, 5);
979         }
980 
981         /** {@inheritDoc} */
982         @Override
983         public void init(final ODEStateAndDerivative s0, final double t) {
984             detector.init(convert(s0), stateMapper.mapDoubleToDate(t));
985             this.lastT = Double.NaN;
986             this.lastG = Double.NaN;
987         }
988 
989         /** {@inheritDoc} */
990         @Override
991         public void reset(final ODEStateAndDerivative intermediateState, final double finalTime) {
992             detector.reset(convert(intermediateState), stateMapper.mapDoubleToDate(finalTime));
993             this.lastT = Double.NaN;
994             this.lastG = Double.NaN;
995         }
996 
997         /** {@inheritDoc} */
998         public double g(final ODEStateAndDerivative s) {
999             if (!Precision.equals(lastT, s.getTime(), 0)) {
1000                 lastT = s.getTime();
1001                 lastG = detector.g(convert(s));
1002             }
1003             return lastG;
1004         }
1005 
1006         /** {@inheritDoc} */
1007         public ODEEventHandler getHandler() {
1008 
1009             return new ODEEventHandler() {
1010 
1011                 /** {@inheritDoc} */
1012                 public Action eventOccurred(final ODEStateAndDerivative s, final ODEEventDetector d, final boolean increasing) {
1013                     return handler.eventOccurred(convert(s), detector, increasing);
1014                 }
1015 
1016                 /** {@inheritDoc} */
1017                 @Override
1018                 public ODEState resetState(final ODEEventDetector d, final ODEStateAndDerivative s) {
1019 
1020                     final SpacecraftState oldState = convert(s);
1021                     final SpacecraftState newState = handler.resetState(detector, oldState);
1022                     stateChanged(newState);
1023 
1024                     // main part
1025                     final double[] primary    = new double[s.getPrimaryStateDimension()];
1026                     stateMapper.mapStateToArray(newState, primary, null);
1027 
1028                     // secondary part
1029                     final double[][] secondary = new double[1][secondaryOffsets.get(SECONDARY_DIMENSION)];
1030                     for (final AdditionalDerivativesProvider provider : additionalDerivativesProviders) {
1031                         final String name      = provider.getName();
1032                         final int    offset    = secondaryOffsets.get(name);
1033                         final int    dimension = provider.getDimension();
1034                         System.arraycopy(newState.getAdditionalState(name), 0, secondary[0], offset, dimension);
1035                     }
1036 
1037                     return new ODEState(newState.getDate().durationFrom(getStartDate()),
1038                                         primary, secondary);
1039 
1040                 }
1041 
1042             };
1043         }
1044 
1045     }
1046 
1047     /** Adapt an {@link org.orekit.propagation.sampling.OrekitStepHandler}
1048      * to Hipparchus {@link ODEStepHandler} interface.
1049      * @author Luc Maisonobe
1050      */
1051     private class AdaptedStepHandler implements ODEStepHandler {
1052 
1053         /** Underlying handler. */
1054         private final OrekitStepHandler handler;
1055 
1056         /** Build an instance.
1057          * @param handler underlying handler to wrap
1058          */
1059         AdaptedStepHandler(final OrekitStepHandler handler) {
1060             this.handler = handler;
1061         }
1062 
1063         /** {@inheritDoc} */
1064         @Override
1065         public void init(final ODEStateAndDerivative s0, final double t) {
1066             handler.init(convert(s0), stateMapper.mapDoubleToDate(t));
1067         }
1068 
1069         /** {@inheritDoc} */
1070         @Override
1071         public void handleStep(final ODEStateInterpolator interpolator) {
1072             handler.handleStep(new AdaptedStepInterpolator(interpolator));
1073         }
1074 
1075         /** {@inheritDoc} */
1076         @Override
1077         public void finish(final ODEStateAndDerivative finalState) {
1078             handler.finish(convert(finalState));
1079         }
1080 
1081     }
1082 
1083     /** Adapt an Hipparchus {@link ODEStateInterpolator}
1084      * to an orekit {@link OrekitStepInterpolator} interface.
1085      * @author Luc Maisonobe
1086      */
1087     private class AdaptedStepInterpolator implements OrekitStepInterpolator {
1088 
1089         /** Underlying raw rawInterpolator. */
1090         private final ODEStateInterpolator mathInterpolator;
1091 
1092         /** Simple constructor.
1093          * @param mathInterpolator underlying raw interpolator
1094          */
1095         AdaptedStepInterpolator(final ODEStateInterpolator mathInterpolator) {
1096             this.mathInterpolator = mathInterpolator;
1097         }
1098 
1099         /** {@inheritDoc}} */
1100         @Override
1101         public SpacecraftState getPreviousState() {
1102             return convert(mathInterpolator.getPreviousState());
1103         }
1104 
1105         /** {@inheritDoc}} */
1106         @Override
1107         public boolean isPreviousStateInterpolated() {
1108             return mathInterpolator.isPreviousStateInterpolated();
1109         }
1110 
1111         /** {@inheritDoc}} */
1112         @Override
1113         public SpacecraftState getCurrentState() {
1114             return convert(mathInterpolator.getCurrentState());
1115         }
1116 
1117         /** {@inheritDoc}} */
1118         @Override
1119         public boolean isCurrentStateInterpolated() {
1120             return mathInterpolator.isCurrentStateInterpolated();
1121         }
1122 
1123         /** {@inheritDoc}} */
1124         @Override
1125         public SpacecraftState getInterpolatedState(final AbsoluteDate date) {
1126             return convert(mathInterpolator.getInterpolatedState(date.durationFrom(stateMapper.getReferenceDate())));
1127         }
1128 
1129         /** {@inheritDoc}} */
1130         @Override
1131         public boolean isForward() {
1132             return mathInterpolator.isForward();
1133         }
1134 
1135         /** {@inheritDoc}} */
1136         @Override
1137         public AdaptedStepInterpolator restrictStep(final SpacecraftState newPreviousState,
1138                                                     final SpacecraftState newCurrentState) {
1139             try {
1140                 final AbstractODEStateInterpolator aosi = (AbstractODEStateInterpolator) mathInterpolator;
1141                 return new AdaptedStepInterpolator(aosi.restrictStep(convert(newPreviousState),
1142                                                                      convert(newCurrentState)));
1143             } catch (ClassCastException cce) {
1144                 // this should never happen
1145                 throw new OrekitInternalError(cce);
1146             }
1147         }
1148 
1149     }
1150 
1151     /** Specialized step handler storing interpolators for ephemeris generation.
1152      * @since 11.0
1153      */
1154     private class StoringStepHandler implements ODEStepHandler, EphemerisGenerator {
1155 
1156         /** Underlying raw mathematical model. */
1157         private DenseOutputModel model;
1158 
1159         /** the user supplied end date. Propagation may not end on this date. */
1160         private AbsoluteDate endDate;
1161 
1162         /** Generated ephemeris. */
1163         private BoundedPropagator ephemeris;
1164 
1165         /** Last interpolator handled by the object.*/
1166         private  ODEStateInterpolator lastInterpolator;
1167 
1168         /** Set the end date.
1169          * @param endDate end date
1170          */
1171         public void setEndDate(final AbsoluteDate endDate) {
1172             this.endDate = endDate;
1173         }
1174 
1175         /** {@inheritDoc} */
1176         @Override
1177         public void init(final ODEStateAndDerivative s0, final double t) {
1178 
1179             this.model = new DenseOutputModel();
1180             model.init(s0, t);
1181 
1182             // ephemeris will be generated when last step is processed
1183             this.ephemeris = null;
1184 
1185             this.lastInterpolator = null;
1186 
1187         }
1188 
1189         /** {@inheritDoc} */
1190         @Override
1191         public BoundedPropagator getGeneratedEphemeris() {
1192             // Each time we try to get the ephemeris, rebuild it using the last data.
1193             buildEphemeris();
1194             return ephemeris;
1195         }
1196 
1197         /** {@inheritDoc} */
1198         @Override
1199         public void handleStep(final ODEStateInterpolator interpolator) {
1200             model.handleStep(interpolator);
1201             lastInterpolator = interpolator;
1202         }
1203 
1204         /** {@inheritDoc} */
1205         @Override
1206         public void finish(final ODEStateAndDerivative finalState) {
1207             buildEphemeris();
1208         }
1209 
1210         /** Method used to produce ephemeris at a given time.
1211          * Can be used at multiple times, updating the ephemeris to
1212          * its last state.
1213          */
1214         private void buildEphemeris() {
1215             // buildEphemeris was built in order to allow access to what was previously the finish method.
1216             // This now allows to call it through getGeneratedEphemeris, therefore through an external call,
1217             // which was not previously the case.
1218 
1219             // Update the model's finalTime with the last interpolator.
1220             model.finish(lastInterpolator.getCurrentState());
1221 
1222             // set up the boundary dates
1223             final double tI = model.getInitialTime();
1224             final double tF = model.getFinalTime();
1225             // tI is almost? always zero
1226             final AbsoluteDate startDate =
1227                             stateMapper.mapDoubleToDate(tI);
1228             final AbsoluteDate finalDate =
1229                             stateMapper.mapDoubleToDate(tF, this.endDate);
1230             final AbsoluteDate minDate;
1231             final AbsoluteDate maxDate;
1232             if (tF < tI) {
1233                 minDate = finalDate;
1234                 maxDate = startDate;
1235             } else {
1236                 minDate = startDate;
1237                 maxDate = finalDate;
1238             }
1239 
1240             // get the initial additional states that are not managed
1241             final DataDictionary unmanaged = new DataDictionary();
1242             for (final DataDictionary.Entry initial : getInitialState().getAdditionalDataValues().getData()) {
1243                 if (!AbstractIntegratedPropagator.this.isAdditionalDataManaged(initial.getKey())) {
1244                     // this additional state was in the initial state, but is unknown to the propagator
1245                     // we simply copy its initial value as is
1246                     unmanaged.put(initial.getKey(), initial.getValue());
1247                 }
1248             }
1249 
1250             // get the names of additional states managed by differential equations
1251             final String[] names      = new String[additionalDerivativesProviders.size()];
1252             final int[]    dimensions = new int[additionalDerivativesProviders.size()];
1253             for (int i = 0; i < names.length; ++i) {
1254                 names[i] = additionalDerivativesProviders.get(i).getName();
1255                 dimensions[i] = additionalDerivativesProviders.get(i).getDimension();
1256             }
1257 
1258             // create the ephemeris
1259             ephemeris = new IntegratedEphemeris(startDate, minDate, maxDate,
1260                                                 stateMapper, getAttitudeProvider(), propagationType, model,
1261                                                 unmanaged, getAdditionalDataProviders(),
1262                                                 names, dimensions);
1263 
1264         }
1265 
1266     }
1267 
1268     /** Wrapper for resetting an integrator handlers.
1269      * <p>
1270      * This class is intended to be used in a try-with-resource statement.
1271      * If propagator-specific event handlers and step handlers are added to
1272      * the integrator in the try block, they will be removed automatically
1273      * when leaving the block, so the integrator only keeps its own handlers
1274      * between calls to {@link AbstractIntegratedPropagator#propagate(AbsoluteDate, AbsoluteDate).
1275      * </p>
1276      * @since 11.0
1277      */
1278     private static class IntegratorResetter implements AutoCloseable {
1279 
1280         /** Wrapped integrator. */
1281         private final ODEIntegrator integrator;
1282 
1283         /** Initial event detectors list. */
1284         private final List<ODEEventDetector> detectors;
1285 
1286         /** Initial step handlers list. */
1287         private final List<ODEStepHandler> stepHandlers;
1288 
1289         /** Simple constructor.
1290          * @param integrator wrapped integrator
1291          */
1292         IntegratorResetter(final ODEIntegrator integrator) {
1293             this.integrator   = integrator;
1294             this.detectors    = new ArrayList<>(integrator.getEventDetectors());
1295             this.stepHandlers = new ArrayList<>(integrator.getStepHandlers());
1296         }
1297 
1298         /** {@inheritDoc}
1299          * <p>
1300          * Reset event handlers and step handlers back to the initial list
1301          * </p>
1302          */
1303         @Override
1304         public void close() {
1305 
1306             // reset event handlers
1307             integrator.clearEventDetectors();
1308             detectors.forEach(integrator::addEventDetector);
1309 
1310             // reset step handlers
1311             integrator.clearStepHandlers();
1312             stepHandlers.forEach(integrator::addStepHandler);
1313 
1314         }
1315 
1316     }
1317 
1318 }