1   /* Copyright 2002-2024 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 java.util.ArrayList;
20  import java.util.Arrays;
21  import java.util.Collection;
22  import java.util.Collections;
23  import java.util.HashMap;
24  import java.util.LinkedList;
25  import java.util.List;
26  import java.util.Map;
27  import java.util.Queue;
28  
29  import org.hipparchus.analysis.UnivariateFunction;
30  import org.hipparchus.analysis.solvers.BracketedUnivariateSolver;
31  import org.hipparchus.analysis.solvers.BracketingNthOrderBrentSolver;
32  import org.hipparchus.exception.MathRuntimeException;
33  import org.hipparchus.ode.DenseOutputModel;
34  import org.hipparchus.ode.ExpandableODE;
35  import org.hipparchus.ode.ODEIntegrator;
36  import org.hipparchus.ode.ODEState;
37  import org.hipparchus.ode.ODEStateAndDerivative;
38  import org.hipparchus.ode.OrdinaryDifferentialEquation;
39  import org.hipparchus.ode.SecondaryODE;
40  import org.hipparchus.ode.events.Action;
41  import org.hipparchus.ode.events.AdaptableInterval;
42  import org.hipparchus.ode.events.ODEEventDetector;
43  import org.hipparchus.ode.events.ODEEventHandler;
44  import org.hipparchus.ode.sampling.AbstractODEStateInterpolator;
45  import org.hipparchus.ode.sampling.ODEStateInterpolator;
46  import org.hipparchus.ode.sampling.ODEStepHandler;
47  import org.hipparchus.util.Precision;
48  import org.orekit.attitudes.AttitudeProvider;
49  import org.orekit.errors.OrekitException;
50  import org.orekit.errors.OrekitInternalError;
51  import org.orekit.errors.OrekitMessages;
52  import org.orekit.frames.Frame;
53  import org.orekit.orbits.OrbitType;
54  import org.orekit.orbits.PositionAngleType;
55  import org.orekit.propagation.AbstractPropagator;
56  import org.orekit.propagation.BoundedPropagator;
57  import org.orekit.propagation.EphemerisGenerator;
58  import org.orekit.propagation.PropagationType;
59  import org.orekit.propagation.SpacecraftState;
60  import org.orekit.propagation.events.EventDetector;
61  import org.orekit.propagation.events.handlers.EventHandler;
62  import org.orekit.propagation.sampling.OrekitStepHandler;
63  import org.orekit.propagation.sampling.OrekitStepInterpolator;
64  import org.orekit.time.AbsoluteDate;
65  import org.orekit.utils.DoubleArrayDictionary;
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 isAdditionalStateManaged(final String name) {
269 
270         // first look at already integrated states
271         if (super.isAdditionalStateManaged(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[] getManagedAdditionalStates() {
288         final String[] alreadyIntegrated = super.getManagedAdditionalStates();
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 #addAdditionalStateProvider(org.orekit.propagation.AdditionalStateProvider)
300      * @since 11.1
301      */
302     public void addAdditionalDerivativesProvider(final AdditionalDerivativesProvider provider) {
303 
304         // check if the name is already used
305         if (isAdditionalStateManaged(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> getEventsDetectors() {
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     /** {@inheritDoc} */
358     @Override
359     public EphemerisGenerator getEphemerisGenerator() {
360         final StoringStepHandler storingHandler = new StoringStepHandler();
361         ephemerisGenerators.add(storingHandler);
362         return storingHandler;
363     }
364 
365     /** Create a mapper between raw double components and spacecraft state.
366     /** Simple constructor.
367      * <p>
368      * The position parameter type is meaningful only if {@link
369      * #getOrbitType() propagation orbit type}
370      * support it. As an example, it is not meaningful for propagation
371      * in {@link OrbitType#CARTESIAN Cartesian} parameters.
372      * </p>
373      * @param referenceDate reference date
374      * @param mu central attraction coefficient (m³/s²)
375      * @param orbitType orbit type to use for mapping
376      * @param positionAngleType angle type to use for propagation
377      * @param attitudeProvider attitude provider
378      * @param frame inertial frame
379      * @return new mapper
380      */
381     protected abstract StateMapper createMapper(AbsoluteDate referenceDate, double mu,
382                                                 OrbitType orbitType, PositionAngleType positionAngleType,
383                                                 AttitudeProvider attitudeProvider, Frame frame);
384 
385     /** Get the differential equations to integrate (for main state only).
386      * @param integ numerical integrator to use for propagation.
387      * @return differential equations for main state
388      */
389     protected abstract MainStateEquations getMainStateEquations(ODEIntegrator integ);
390 
391     /** {@inheritDoc} */
392     @Override
393     public SpacecraftState propagate(final AbsoluteDate target) {
394         if (getStartDate() == null) {
395             if (getInitialState() == null) {
396                 throw new OrekitException(OrekitMessages.INITIAL_STATE_NOT_SPECIFIED_FOR_ORBIT_PROPAGATION);
397             }
398             setStartDate(getInitialState().getDate());
399         }
400         return propagate(getStartDate(), target);
401     }
402 
403     /** {@inheritDoc} */
404     public SpacecraftState propagate(final AbsoluteDate tStart, final AbsoluteDate tEnd) {
405 
406         if (getInitialState() == null) {
407             throw new OrekitException(OrekitMessages.INITIAL_STATE_NOT_SPECIFIED_FOR_ORBIT_PROPAGATION);
408         }
409 
410         // make sure the integrator will be reset properly even if we change its events handlers and step handlers
411         try (IntegratorResetter resetter = new IntegratorResetter(integrator)) {
412 
413             // prepare handling of STM and Jacobian matrices
414             setUpStmAndJacobianGenerators();
415 
416             // Initialize additional states
417             initializeAdditionalStates(tEnd);
418 
419             if (!tStart.equals(getInitialState().getDate())) {
420                 // if propagation start date is not initial date,
421                 // propagate from initial to start date without event detection
422                 try (IntegratorResetter startResetter = new IntegratorResetter(integrator)) {
423                     integrateDynamics(tStart, true);
424                 }
425             }
426 
427             // set up events added by user
428             setUpUserEventDetectors();
429 
430             // set up step handlers
431             for (final OrekitStepHandler handler : getMultiplexer().getHandlers()) {
432                 integrator.addStepHandler(new AdaptedStepHandler(handler));
433             }
434             for (final StoringStepHandler generator : ephemerisGenerators) {
435                 generator.setEndDate(tEnd);
436                 integrator.addStepHandler(generator);
437             }
438 
439             // propagate from start date to end date with event detection
440             final SpacecraftState finalState = integrateDynamics(tEnd, false);
441 
442             // Finalize event detectors
443             getEventsDetectors().forEach(detector -> detector.finish(finalState));
444 
445             return finalState;
446         }
447 
448     }
449 
450     /** Reset initial state with a given propagation type.
451      *
452      * <p> By default this method returns the same as {@link #resetInitialState(SpacecraftState)}
453      * <p> Its purpose is mostly to be derived in DSSTPropagator
454      *
455      * @param state new initial state to consider
456      * @param stateType type of the new state (mean or osculating)
457      * @since 12.1.3
458      */
459     public void resetInitialState(final SpacecraftState state, final PropagationType stateType) {
460         // Default behavior, do not take propagation type into account
461         resetInitialState(state);
462     }
463 
464     /** Set up State Transition Matrix and Jacobian matrix handling.
465      * @since 11.1
466      */
467     protected void setUpStmAndJacobianGenerators() {
468         // nothing to do by default
469     }
470 
471     /** Propagation with or without event detection.
472      * @param tEnd target date to which orbit should be propagated
473      * @param forceResetAtEnd flag to force resetting state and date after integration
474      * @return state at end of propagation
475      */
476     private SpacecraftState integrateDynamics(final AbsoluteDate tEnd, final boolean forceResetAtEnd) {
477         try {
478 
479             initializePropagation();
480 
481             if (getInitialState().getDate().equals(tEnd)) {
482                 // don't extrapolate
483                 return getInitialState();
484             }
485 
486             // space dynamics view
487             stateMapper = createMapper(getInitialState().getDate(), stateMapper.getMu(),
488                                        stateMapper.getOrbitType(), stateMapper.getPositionAngleType(),
489                                        stateMapper.getAttitudeProvider(), getInitialState().getFrame());
490             attitudeProviderForDerivatives = initializeAttitudeProviderForDerivatives();
491 
492             if (Double.isNaN(getMu())) {
493                 setMu(getInitialState().getMu());
494             }
495 
496             if (getInitialState().getMass() <= 0.0) {
497                 throw new OrekitException(OrekitMessages.NOT_POSITIVE_SPACECRAFT_MASS,
498                                           getInitialState().getMass());
499             }
500 
501             // convert space flight dynamics API to math API
502             final SpacecraftState initialIntegrationState = getInitialIntegrationState();
503             final ODEState mathInitialState = createInitialState(initialIntegrationState);
504             final ExpandableODE mathODE = createODE(integrator);
505 
506             // mathematical integration
507             final ODEStateAndDerivative mathFinalState;
508             beforeIntegration(initialIntegrationState, tEnd);
509             mathFinalState = integrator.integrate(mathODE, mathInitialState,
510                                                   tEnd.durationFrom(getInitialState().getDate()));
511             afterIntegration();
512 
513             // get final state
514             SpacecraftState finalState =
515                             stateMapper.mapArrayToState(stateMapper.mapDoubleToDate(mathFinalState.getTime(),
516                                                                                     tEnd),
517                                                         mathFinalState.getPrimaryState(),
518                                                         mathFinalState.getPrimaryDerivative(),
519                                                         propagationType);
520 
521             finalState = updateAdditionalStatesAndDerivatives(finalState, mathFinalState);
522 
523             if (resetAtEnd || forceResetAtEnd) {
524                 resetInitialState(finalState, propagationType);
525                 setStartDate(finalState.getDate());
526             }
527 
528             return finalState;
529 
530         } catch (MathRuntimeException mre) {
531             throw OrekitException.unwrap(mre);
532         }
533     }
534 
535     /**
536      * Returns an updated version of the inputted state with additional states, including
537      * from derivatives providers.
538      * @param originalState input state
539      * @param os ODE state and derivative
540      * @return new state
541      * @since 12.1
542      */
543     private SpacecraftState updateAdditionalStatesAndDerivatives(final SpacecraftState originalState,
544                                                                  final ODEStateAndDerivative os) {
545         SpacecraftState updatedState = originalState;
546         if (os.getNumberOfSecondaryStates() > 0) {
547             final double[] secondary           = os.getSecondaryState(1);
548             final double[] secondaryDerivative = os.getSecondaryDerivative(1);
549             for (final AdditionalDerivativesProvider provider : additionalDerivativesProviders) {
550                 final String name      = provider.getName();
551                 final int    offset    = secondaryOffsets.get(name);
552                 final int    dimension = provider.getDimension();
553                 updatedState = updatedState.addAdditionalState(name, Arrays.copyOfRange(secondary, offset, offset + dimension));
554                 updatedState = updatedState.addAdditionalStateDerivative(name, Arrays.copyOfRange(secondaryDerivative, offset, offset + dimension));
555             }
556         }
557         return updateAdditionalStates(updatedState);
558     }
559 
560     /** Get the initial state for integration.
561      * @return initial state for integration
562      */
563     protected SpacecraftState getInitialIntegrationState() {
564         return getInitialState();
565     }
566 
567     /** Create an initial state.
568      * @param initialState initial state in flight dynamics world
569      * @return initial state in mathematics world
570      */
571     private ODEState createInitialState(final SpacecraftState initialState) {
572 
573         // retrieve initial state
574         final double[] primary = new double[getBasicDimension()];
575         stateMapper.mapStateToArray(initialState, primary, null);
576 
577         if (secondaryOffsets.isEmpty()) {
578             // compute dimension of the secondary state
579             int offset = 0;
580             for (final AdditionalDerivativesProvider provider : additionalDerivativesProviders) {
581                 secondaryOffsets.put(provider.getName(), offset);
582                 offset += provider.getDimension();
583             }
584             secondaryOffsets.put(SECONDARY_DIMENSION, offset);
585         }
586 
587         return new ODEState(0.0, primary, secondary(initialState));
588 
589     }
590 
591     /** Create secondary state.
592      * @param state spacecraft state
593      * @return secondary state
594      * @since 11.1
595      */
596     private double[][] secondary(final SpacecraftState state) {
597 
598         if (secondaryOffsets.isEmpty()) {
599             return null;
600         }
601 
602         final double[][] secondary = new double[1][secondaryOffsets.get(SECONDARY_DIMENSION)];
603         for (final AdditionalDerivativesProvider provider : additionalDerivativesProviders) {
604             final String   name       = provider.getName();
605             final int      offset     = secondaryOffsets.get(name);
606             final double[] additional = state.getAdditionalState(name);
607             System.arraycopy(additional, 0, secondary[0], offset, additional.length);
608         }
609 
610         return secondary;
611 
612     }
613 
614     /** Create secondary state derivative.
615      * @param state spacecraft state
616      * @return secondary state derivative
617      * @since 11.1
618      */
619     private double[][] secondaryDerivative(final SpacecraftState state) {
620 
621         if (secondaryOffsets.isEmpty()) {
622             return null;
623         }
624 
625         final double[][] secondaryDerivative = new double[1][secondaryOffsets.get(SECONDARY_DIMENSION)];
626         for (final AdditionalDerivativesProvider provider : additionalDerivativesProviders) {
627             final String   name       = provider.getName();
628             final int      offset     = secondaryOffsets.get(name);
629             final double[] additionalDerivative = state.getAdditionalStateDerivative(name);
630             System.arraycopy(additionalDerivative, 0, secondaryDerivative[0], offset, additionalDerivative.length);
631         }
632 
633         return secondaryDerivative;
634 
635     }
636 
637     /** Create an ODE with all equations.
638      * @param integ numerical integrator to use for propagation.
639      * @return a new ode
640      */
641     private ExpandableODE createODE(final ODEIntegrator integ) {
642 
643         final ExpandableODE ode =
644                 new ExpandableODE(new ConvertedMainStateEquations(getMainStateEquations(integ)));
645 
646         // secondary part of the ODE
647         if (!additionalDerivativesProviders.isEmpty()) {
648             ode.addSecondaryEquations(new ConvertedSecondaryStateEquations());
649         }
650 
651         return ode;
652 
653     }
654 
655     /** Method called just before integration.
656      * <p>
657      * The default implementation does nothing, it may be specialized in subclasses.
658      * </p>
659      * @param initialState initial state
660      * @param tEnd target date at which state should be propagated
661      */
662     protected void beforeIntegration(final SpacecraftState initialState,
663                                      final AbsoluteDate tEnd) {
664         // do nothing by default
665     }
666 
667     /** Method called just after integration.
668      * <p>
669      * The default implementation does nothing, it may be specialized in subclasses.
670      * </p>
671      */
672     protected void afterIntegration() {
673         // do nothing by default
674     }
675 
676     /** Get state vector dimension without additional parameters.
677      * @return state vector dimension without additional parameters.
678      */
679     public int getBasicDimension() {
680         return 7;
681     }
682 
683     /** Get the integrator used by the propagator.
684      * @return the integrator.
685      */
686     protected ODEIntegrator getIntegrator() {
687         return integrator;
688     }
689 
690     /** Convert a state from mathematical world to space flight dynamics world.
691      * @param os mathematical state
692      * @return space flight dynamics state
693      */
694     private SpacecraftState convert(final ODEStateAndDerivative os) {
695 
696         final SpacecraftState s = stateMapper.mapArrayToState(os.getTime(), os.getPrimaryState(),
697             os.getPrimaryDerivative(), propagationType);
698         return updateAdditionalStatesAndDerivatives(s, os);
699     }
700 
701     /** Convert a state from space flight dynamics world to mathematical world.
702      * @param state space flight dynamics state
703      * @return mathematical state
704      */
705     private ODEStateAndDerivative convert(final SpacecraftState state) {
706 
707         // retrieve initial state
708         final double[] primary    = new double[getBasicDimension()];
709         final double[] primaryDot = new double[getBasicDimension()];
710         stateMapper.mapStateToArray(state, primary, primaryDot);
711 
712         // secondary part of the ODE
713         final double[][] secondary           = secondary(state);
714         final double[][] secondaryDerivative = secondaryDerivative(state);
715 
716         return new ODEStateAndDerivative(stateMapper.mapDateToDouble(state.getDate()),
717                                          primary, primaryDot,
718                                          secondary, secondaryDerivative);
719 
720     }
721 
722     /** Differential equations for the main state (orbit, attitude and mass). */
723     public interface MainStateEquations {
724 
725         /**
726          * Initialize the equations at the start of propagation. This method will be
727          * called before any calls to {@link #computeDerivatives(SpacecraftState)}.
728          *
729          * <p> The default implementation of this method does nothing.
730          *
731          * @param initialState initial state information at the start of propagation.
732          * @param target       date of propagation. Not equal to {@code
733          *                     initialState.getDate()}.
734          */
735         default void init(final SpacecraftState initialState, final AbsoluteDate target) {
736         }
737 
738         /** Compute differential equations for main state.
739          * @param state current state
740          * @return derivatives of main state
741          */
742         double[] computeDerivatives(SpacecraftState state);
743 
744     }
745 
746     /** Differential equations for the main state (orbit, attitude and mass), with converted API. */
747     private class ConvertedMainStateEquations implements OrdinaryDifferentialEquation {
748 
749         /** Main state equations. */
750         private final MainStateEquations main;
751 
752         /** Simple constructor.
753          * @param main main state equations
754          */
755         ConvertedMainStateEquations(final MainStateEquations main) {
756             this.main = main;
757             calls = 0;
758         }
759 
760         /** {@inheritDoc} */
761         public int getDimension() {
762             return getBasicDimension();
763         }
764 
765         @Override
766         public void init(final double t0, final double[] y0, final double finalTime) {
767             // update space dynamics view
768             SpacecraftState initialState = stateMapper.mapArrayToState(t0, y0, null, PropagationType.MEAN);
769             initialState = updateAdditionalStates(initialState);
770             initialState = updateStatesFromAdditionalDerivativesIfKnown(initialState);
771             final AbsoluteDate target = stateMapper.mapDoubleToDate(finalTime);
772             main.init(initialState, target);
773             attitudeProviderForDerivatives = initializeAttitudeProviderForDerivatives();
774         }
775 
776         /**
777          * Returns an updated version of the inputted state, with additional states from
778          * derivatives providers as given in the stored initial state.
779          * @param originalState input state
780          * @return new state
781          * @since 12.1
782          */
783         private SpacecraftState updateStatesFromAdditionalDerivativesIfKnown(final SpacecraftState originalState) {
784             SpacecraftState updatedState = originalState;
785             final SpacecraftState storedInitialState = getInitialState();
786             final double originalTime = stateMapper.mapDateToDouble(originalState.getDate());
787             if (storedInitialState != null && stateMapper.mapDateToDouble(storedInitialState.getDate()) == originalTime) {
788                 for (final AdditionalDerivativesProvider provider: additionalDerivativesProviders) {
789                     final String name = provider.getName();
790                     final double[] value = storedInitialState.getAdditionalState(name);
791                     updatedState = updatedState.addAdditionalState(name, value);
792                 }
793             }
794             return updatedState;
795         }
796 
797         /** {@inheritDoc} */
798         public double[] computeDerivatives(final double t, final double[] y) {
799 
800             // increment calls counter
801             ++calls;
802 
803             // update space dynamics view
804             stateMapper.setAttitudeProvider(attitudeProviderForDerivatives);
805             SpacecraftState currentState = stateMapper.mapArrayToState(t, y, null, PropagationType.MEAN);
806             stateMapper.setAttitudeProvider(getAttitudeProvider());
807 
808             currentState = updateAdditionalStates(currentState);
809             // compute main state differentials
810             return main.computeDerivatives(currentState);
811 
812         }
813 
814     }
815 
816     /** Differential equations for the secondary state (Jacobians, user variables ...), with converted API. */
817     private class ConvertedSecondaryStateEquations implements SecondaryODE {
818 
819         /** Dimension of the combined additional states. */
820         private final int combinedDimension;
821 
822         /** Simple constructor.
823           */
824         ConvertedSecondaryStateEquations() {
825             this.combinedDimension = secondaryOffsets.get(SECONDARY_DIMENSION);
826         }
827 
828         /** {@inheritDoc} */
829         @Override
830         public int getDimension() {
831             return combinedDimension;
832         }
833 
834         /** {@inheritDoc} */
835         @Override
836         public void init(final double t0, final double[] primary0,
837                          final double[] secondary0, final double finalTime) {
838             // update space dynamics view
839             final SpacecraftState initialState = convert(t0, primary0, null, secondary0);
840 
841             final AbsoluteDate target = stateMapper.mapDoubleToDate(finalTime);
842             for (final AdditionalDerivativesProvider provider : additionalDerivativesProviders) {
843                 provider.init(initialState, target);
844             }
845 
846         }
847 
848         /** {@inheritDoc} */
849         @Override
850         public double[] computeDerivatives(final double t, final double[] primary,
851                                            final double[] primaryDot, final double[] secondary) {
852 
853             // update space dynamics view
854             // the integrable generators generate method will be called here,
855             // according to the generators yield order
856             SpacecraftState updated = convert(t, primary, primaryDot, secondary);
857 
858             // set up queue for equations
859             final Queue<AdditionalDerivativesProvider> pending = new LinkedList<>(additionalDerivativesProviders);
860 
861             // gather the derivatives from all additional equations, taking care of dependencies
862             final double[] secondaryDot = new double[combinedDimension];
863             int yieldCount = 0;
864             while (!pending.isEmpty()) {
865                 final AdditionalDerivativesProvider provider = pending.remove();
866                 if (provider.yields(updated)) {
867                     // this provider has to wait for another one,
868                     // we put it again in the pending queue
869                     pending.add(provider);
870                     if (++yieldCount >= pending.size()) {
871                         // all pending providers yielded!, they probably need data not yet initialized
872                         // we let the propagation proceed, if these data are really needed right now
873                         // an appropriate exception will be triggered when caller tries to access them
874                         break;
875                     }
876                 } else {
877                     // we can use these equations right now
878                     final String              name           = provider.getName();
879                     final int                 offset         = secondaryOffsets.get(name);
880                     final int                 dimension      = provider.getDimension();
881                     final CombinedDerivatives derivatives    = provider.combinedDerivatives(updated);
882                     final double[]            additionalPart = derivatives.getAdditionalDerivatives();
883                     final double[]            mainPart       = derivatives.getMainStateDerivativesIncrements();
884                     System.arraycopy(additionalPart, 0, secondaryDot, offset, dimension);
885                     updated = updated.addAdditionalStateDerivative(name, additionalPart);
886                     if (mainPart != null) {
887                         // this equation does change the main state derivatives
888                         for (int i = 0; i < mainPart.length; ++i) {
889                             primaryDot[i] += mainPart[i];
890                         }
891                     }
892                     yieldCount = 0;
893                 }
894             }
895 
896             return secondaryDot;
897 
898         }
899 
900         /** Convert mathematical view to space view.
901          * @param t current value of the independent <I>time</I> variable
902          * @param primary array containing the current value of the primary state vector
903          * @param primaryDot array containing the derivative of the primary state vector
904          * @param secondary array containing the current value of the secondary state vector
905          * @return space view of the state
906          */
907         private SpacecraftState convert(final double t, final double[] primary,
908                                         final double[] primaryDot, final double[] secondary) {
909 
910             SpacecraftState initialState = stateMapper.mapArrayToState(t, primary, primaryDot, PropagationType.MEAN);
911 
912             for (final AdditionalDerivativesProvider provider : additionalDerivativesProviders) {
913                 final String name      = provider.getName();
914                 final int    offset    = secondaryOffsets.get(name);
915                 final int    dimension = provider.getDimension();
916                 initialState = initialState.addAdditionalState(name, Arrays.copyOfRange(secondary, offset, offset + dimension));
917             }
918 
919             return updateAdditionalStates(initialState);
920 
921         }
922 
923     }
924 
925     /** Adapt an {@link org.orekit.propagation.events.EventDetector}
926      * to Hipparchus {@link org.hipparchus.ode.events.ODEEventDetector} interface.
927      * @author Fabien Maussion
928      */
929     private class AdaptedEventDetector implements ODEEventDetector {
930 
931         /** Underlying event detector. */
932         private final EventDetector detector;
933 
934         /** Underlying event handler.
935          * @since 12.0
936          */
937         private final EventHandler handler;
938 
939         /** Time of the previous call to g. */
940         private double lastT;
941 
942         /** Value from the previous call to g. */
943         private double lastG;
944 
945         /** Build a wrapped event detector.
946          * @param detector event detector to wrap
947         */
948         AdaptedEventDetector(final EventDetector detector) {
949             this.detector = detector;
950             this.handler  = detector.getHandler();
951             this.lastT    = Double.NaN;
952             this.lastG    = Double.NaN;
953         }
954 
955         /** {@inheritDoc} */
956         @Override
957         public AdaptableInterval getMaxCheckInterval() {
958             return s -> detector.getMaxCheckInterval().currentInterval(convert(s));
959         }
960 
961         /** {@inheritDoc} */
962         @Override
963         public int getMaxIterationCount() {
964             return detector.getMaxIterationCount();
965         }
966 
967         /** {@inheritDoc} */
968         @Override
969         public BracketedUnivariateSolver<UnivariateFunction> getSolver() {
970             return new BracketingNthOrderBrentSolver(0, detector.getThreshold(), 0, 5);
971         }
972 
973         /** {@inheritDoc} */
974         @Override
975         public void init(final ODEStateAndDerivative s0, final double t) {
976             detector.init(convert(s0), stateMapper.mapDoubleToDate(t));
977             this.lastT = Double.NaN;
978             this.lastG = Double.NaN;
979         }
980 
981         /** {@inheritDoc} */
982         public double g(final ODEStateAndDerivative s) {
983             if (!Precision.equals(lastT, s.getTime(), 0)) {
984                 lastT = s.getTime();
985                 lastG = detector.g(convert(s));
986             }
987             return lastG;
988         }
989 
990         /** {@inheritDoc} */
991         public ODEEventHandler getHandler() {
992 
993             return new ODEEventHandler() {
994 
995                 /** {@inheritDoc} */
996                 public Action eventOccurred(final ODEStateAndDerivative s, final ODEEventDetector d, final boolean increasing) {
997                     return handler.eventOccurred(convert(s), detector, increasing);
998                 }
999 
1000                 /** {@inheritDoc} */
1001                 @Override
1002                 public ODEState resetState(final ODEEventDetector d, final ODEStateAndDerivative s) {
1003 
1004                     final SpacecraftState oldState = convert(s);
1005                     final SpacecraftState newState = handler.resetState(detector, oldState);
1006                     stateChanged(newState);
1007 
1008                     // main part
1009                     final double[] primary    = new double[s.getPrimaryStateDimension()];
1010                     stateMapper.mapStateToArray(newState, primary, null);
1011 
1012                     // secondary part
1013                     final double[][] secondary = new double[1][secondaryOffsets.get(SECONDARY_DIMENSION)];
1014                     for (final AdditionalDerivativesProvider provider : additionalDerivativesProviders) {
1015                         final String name      = provider.getName();
1016                         final int    offset    = secondaryOffsets.get(name);
1017                         final int    dimension = provider.getDimension();
1018                         System.arraycopy(newState.getAdditionalState(name), 0, secondary[0], offset, dimension);
1019                     }
1020 
1021                     return new ODEState(newState.getDate().durationFrom(getStartDate()),
1022                                         primary, secondary);
1023 
1024                 }
1025 
1026             };
1027         }
1028 
1029     }
1030 
1031     /** Adapt an {@link org.orekit.propagation.sampling.OrekitStepHandler}
1032      * to Hipparchus {@link ODEStepHandler} interface.
1033      * @author Luc Maisonobe
1034      */
1035     private class AdaptedStepHandler implements ODEStepHandler {
1036 
1037         /** Underlying handler. */
1038         private final OrekitStepHandler handler;
1039 
1040         /** Build an instance.
1041          * @param handler underlying handler to wrap
1042          */
1043         AdaptedStepHandler(final OrekitStepHandler handler) {
1044             this.handler = handler;
1045         }
1046 
1047         /** {@inheritDoc} */
1048         @Override
1049         public void init(final ODEStateAndDerivative s0, final double t) {
1050             handler.init(convert(s0), stateMapper.mapDoubleToDate(t));
1051         }
1052 
1053         /** {@inheritDoc} */
1054         @Override
1055         public void handleStep(final ODEStateInterpolator interpolator) {
1056             handler.handleStep(new AdaptedStepInterpolator(interpolator));
1057         }
1058 
1059         /** {@inheritDoc} */
1060         @Override
1061         public void finish(final ODEStateAndDerivative finalState) {
1062             handler.finish(convert(finalState));
1063         }
1064 
1065     }
1066 
1067     /** Adapt an Hipparchus {@link ODEStateInterpolator}
1068      * to an orekit {@link OrekitStepInterpolator} interface.
1069      * @author Luc Maisonobe
1070      */
1071     private class AdaptedStepInterpolator implements OrekitStepInterpolator {
1072 
1073         /** Underlying raw rawInterpolator. */
1074         private final ODEStateInterpolator mathInterpolator;
1075 
1076         /** Simple constructor.
1077          * @param mathInterpolator underlying raw interpolator
1078          */
1079         AdaptedStepInterpolator(final ODEStateInterpolator mathInterpolator) {
1080             this.mathInterpolator = mathInterpolator;
1081         }
1082 
1083         /** {@inheritDoc}} */
1084         @Override
1085         public SpacecraftState getPreviousState() {
1086             return convert(mathInterpolator.getPreviousState());
1087         }
1088 
1089         /** {@inheritDoc}} */
1090         @Override
1091         public boolean isPreviousStateInterpolated() {
1092             return mathInterpolator.isPreviousStateInterpolated();
1093         }
1094 
1095         /** {@inheritDoc}} */
1096         @Override
1097         public SpacecraftState getCurrentState() {
1098             return convert(mathInterpolator.getCurrentState());
1099         }
1100 
1101         /** {@inheritDoc}} */
1102         @Override
1103         public boolean isCurrentStateInterpolated() {
1104             return mathInterpolator.isCurrentStateInterpolated();
1105         }
1106 
1107         /** {@inheritDoc}} */
1108         @Override
1109         public SpacecraftState getInterpolatedState(final AbsoluteDate date) {
1110             return convert(mathInterpolator.getInterpolatedState(date.durationFrom(stateMapper.getReferenceDate())));
1111         }
1112 
1113         /** {@inheritDoc}} */
1114         @Override
1115         public boolean isForward() {
1116             return mathInterpolator.isForward();
1117         }
1118 
1119         /** {@inheritDoc}} */
1120         @Override
1121         public AdaptedStepInterpolator restrictStep(final SpacecraftState newPreviousState,
1122                                                     final SpacecraftState newCurrentState) {
1123             try {
1124                 final AbstractODEStateInterpolator aosi = (AbstractODEStateInterpolator) mathInterpolator;
1125                 return new AdaptedStepInterpolator(aosi.restrictStep(convert(newPreviousState),
1126                                                                      convert(newCurrentState)));
1127             } catch (ClassCastException cce) {
1128                 // this should never happen
1129                 throw new OrekitInternalError(cce);
1130             }
1131         }
1132 
1133     }
1134 
1135     /** Specialized step handler storing interpolators for ephemeris generation.
1136      * @since 11.0
1137      */
1138     private class StoringStepHandler implements ODEStepHandler, EphemerisGenerator {
1139 
1140         /** Underlying raw mathematical model. */
1141         private DenseOutputModel model;
1142 
1143         /** the user supplied end date. Propagation may not end on this date. */
1144         private AbsoluteDate endDate;
1145 
1146         /** Generated ephemeris. */
1147         private BoundedPropagator ephemeris;
1148 
1149         /** Last interpolator handled by the object.*/
1150         private  ODEStateInterpolator lastInterpolator;
1151 
1152         /** Set the end date.
1153          * @param endDate end date
1154          */
1155         public void setEndDate(final AbsoluteDate endDate) {
1156             this.endDate = endDate;
1157         }
1158 
1159         /** {@inheritDoc} */
1160         @Override
1161         public void init(final ODEStateAndDerivative s0, final double t) {
1162 
1163             this.model = new DenseOutputModel();
1164             model.init(s0, t);
1165 
1166             // ephemeris will be generated when last step is processed
1167             this.ephemeris = null;
1168 
1169             this.lastInterpolator = null;
1170 
1171         }
1172 
1173         /** {@inheritDoc} */
1174         @Override
1175         public BoundedPropagator getGeneratedEphemeris() {
1176             // Each time we try to get the ephemeris, rebuild it using the last data.
1177             buildEphemeris();
1178             return ephemeris;
1179         }
1180 
1181         /** {@inheritDoc} */
1182         @Override
1183         public void handleStep(final ODEStateInterpolator interpolator) {
1184             model.handleStep(interpolator);
1185             lastInterpolator = interpolator;
1186         }
1187 
1188         /** {@inheritDoc} */
1189         @Override
1190         public void finish(final ODEStateAndDerivative finalState) {
1191             buildEphemeris();
1192         }
1193 
1194         /** Method used to produce ephemeris at a given time.
1195          * Can be used at multiple times, updating the ephemeris to
1196          * its last state.
1197          */
1198         private void buildEphemeris() {
1199             // buildEphemeris was built in order to allow access to what was previously the finish method.
1200             // This now allows to call it through getGeneratedEphemeris, therefore through an external call,
1201             // which was not previously the case.
1202 
1203             // Update the model's finalTime with the last interpolator.
1204             model.finish(lastInterpolator.getCurrentState());
1205 
1206             // set up the boundary dates
1207             final double tI = model.getInitialTime();
1208             final double tF = model.getFinalTime();
1209             // tI is almost? always zero
1210             final AbsoluteDate startDate =
1211                             stateMapper.mapDoubleToDate(tI);
1212             final AbsoluteDate finalDate =
1213                             stateMapper.mapDoubleToDate(tF, this.endDate);
1214             final AbsoluteDate minDate;
1215             final AbsoluteDate maxDate;
1216             if (tF < tI) {
1217                 minDate = finalDate;
1218                 maxDate = startDate;
1219             } else {
1220                 minDate = startDate;
1221                 maxDate = finalDate;
1222             }
1223 
1224             // get the initial additional states that are not managed
1225             final DoubleArrayDictionary unmanaged = new DoubleArrayDictionary();
1226             for (final DoubleArrayDictionary.Entry initial : getInitialState().getAdditionalStatesValues().getData()) {
1227                 if (!isAdditionalStateManaged(initial.getKey())) {
1228                     // this additional state was in the initial state, but is unknown to the propagator
1229                     // we simply copy its initial value as is
1230                     unmanaged.put(initial.getKey(), initial.getValue());
1231                 }
1232             }
1233 
1234             // get the names of additional states managed by differential equations
1235             final String[] names      = new String[additionalDerivativesProviders.size()];
1236             final int[]    dimensions = new int[additionalDerivativesProviders.size()];
1237             for (int i = 0; i < names.length; ++i) {
1238                 names[i] = additionalDerivativesProviders.get(i).getName();
1239                 dimensions[i] = additionalDerivativesProviders.get(i).getDimension();
1240             }
1241 
1242             // create the ephemeris
1243             ephemeris = new IntegratedEphemeris(startDate, minDate, maxDate,
1244                                                 stateMapper, propagationType, model,
1245                                                 unmanaged, getAdditionalStateProviders(),
1246                                                 names, dimensions);
1247 
1248         }
1249 
1250     }
1251 
1252     /** Wrapper for resetting an integrator handlers.
1253      * <p>
1254      * This class is intended to be used in a try-with-resource statement.
1255      * If propagator-specific event handlers and step handlers are added to
1256      * the integrator in the try block, they will be removed automatically
1257      * when leaving the block, so the integrator only keeps its own handlers
1258      * between calls to {@link AbstractIntegratedPropagator#propagate(AbsoluteDate, AbsoluteDate).
1259      * </p>
1260      * @since 11.0
1261      */
1262     private static class IntegratorResetter implements AutoCloseable {
1263 
1264         /** Wrapped integrator. */
1265         private final ODEIntegrator integrator;
1266 
1267         /** Initial event detectors list. */
1268         private final List<ODEEventDetector> detectors;
1269 
1270         /** Initial step handlers list. */
1271         private final List<ODEStepHandler> stepHandlers;
1272 
1273         /** Simple constructor.
1274          * @param integrator wrapped integrator
1275          */
1276         IntegratorResetter(final ODEIntegrator integrator) {
1277             this.integrator   = integrator;
1278             this.detectors    = new ArrayList<>(integrator.getEventDetectors());
1279             this.stepHandlers = new ArrayList<>(integrator.getStepHandlers());
1280         }
1281 
1282         /** {@inheritDoc}
1283          * <p>
1284          * Reset event handlers and step handlers back to the initial list
1285          * </p>
1286          */
1287         @Override
1288         public void close() {
1289 
1290             // reset event handlers
1291             integrator.clearEventDetectors();
1292             detectors.forEach(integrator::addEventDetector);
1293 
1294             // reset step handlers
1295             integrator.clearStepHandlers();
1296             stepHandlers.forEach(integrator::addStepHandler);
1297 
1298         }
1299 
1300     }
1301 
1302 }