AbstractIntegratedPropagator.java

  1. /* Copyright 2002-2020 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. import java.util.ArrayList;
  19. import java.util.Collection;
  20. import java.util.Collections;
  21. import java.util.HashMap;
  22. import java.util.List;
  23. import java.util.Map;

  24. import org.hipparchus.exception.MathRuntimeException;
  25. import org.hipparchus.ode.DenseOutputModel;
  26. import org.hipparchus.ode.EquationsMapper;
  27. import org.hipparchus.ode.ExpandableODE;
  28. import org.hipparchus.ode.ODEIntegrator;
  29. import org.hipparchus.ode.ODEState;
  30. import org.hipparchus.ode.ODEStateAndDerivative;
  31. import org.hipparchus.ode.OrdinaryDifferentialEquation;
  32. import org.hipparchus.ode.SecondaryODE;
  33. import org.hipparchus.ode.events.Action;
  34. import org.hipparchus.ode.events.ODEEventHandler;
  35. import org.hipparchus.ode.sampling.AbstractODEStateInterpolator;
  36. import org.hipparchus.ode.sampling.ODEStateInterpolator;
  37. import org.hipparchus.ode.sampling.ODEStepHandler;
  38. import org.hipparchus.util.Precision;
  39. import org.orekit.attitudes.AttitudeProvider;
  40. import org.orekit.errors.OrekitException;
  41. import org.orekit.errors.OrekitIllegalStateException;
  42. import org.orekit.errors.OrekitInternalError;
  43. import org.orekit.errors.OrekitMessages;
  44. import org.orekit.frames.Frame;
  45. import org.orekit.orbits.OrbitType;
  46. import org.orekit.orbits.PositionAngle;
  47. import org.orekit.propagation.AbstractPropagator;
  48. import org.orekit.propagation.BoundedPropagator;
  49. import org.orekit.propagation.PropagationType;
  50. import org.orekit.propagation.SpacecraftState;
  51. import org.orekit.propagation.events.EventDetector;
  52. import org.orekit.propagation.sampling.OrekitStepHandler;
  53. import org.orekit.propagation.sampling.OrekitStepInterpolator;
  54. import org.orekit.time.AbsoluteDate;


  55. /** Common handling of {@link org.orekit.propagation.Propagator Propagator}
  56.  *  methods for both numerical and semi-analytical propagators.
  57.  *  @author Luc Maisonobe
  58.  */
  59. public abstract class AbstractIntegratedPropagator extends AbstractPropagator {

  60.     /** Event detectors not related to force models. */
  61.     private final List<EventDetector> detectors;

  62.     /** Integrator selected by the user for the orbital extrapolation process. */
  63.     private final ODEIntegrator integrator;

  64.     /** Mode handler. */
  65.     private ModeHandler modeHandler;

  66.     /** Additional equations. */
  67.     private List<AdditionalEquations> additionalEquations;

  68.     /** Counter for differential equations calls. */
  69.     private int calls;

  70.     /** Mapper between raw double components and space flight dynamics objects. */
  71.     private StateMapper stateMapper;

  72.     /** Equations mapper. */
  73.     private EquationsMapper equationsMapper;

  74.     /** Flag for resetting the state at end of propagation. */
  75.     private boolean resetAtEnd;

  76.     /** Type of orbit to output (mean or osculating) <br/>
  77.      * <p>
  78.      * This is used only in the case of semianalitical propagators where there is a clear separation between
  79.      * mean and short periodic elements. It is ignored by the Numerical propagator.
  80.      * </p>
  81.      */
  82.     private PropagationType propagationType;

  83.     /** Build a new instance.
  84.      * @param integrator numerical integrator to use for propagation.
  85.      * @param propagationType type of orbit to output (mean or osculating).
  86.      */
  87.     protected AbstractIntegratedPropagator(final ODEIntegrator integrator, final PropagationType propagationType) {
  88.         detectors            = new ArrayList<EventDetector>();
  89.         additionalEquations  = new ArrayList<AdditionalEquations>();
  90.         this.integrator      = integrator;
  91.         this.propagationType = propagationType;
  92.         this.resetAtEnd      = true;
  93.     }

  94.     /** Allow/disallow resetting the initial state at end of propagation.
  95.      * <p>
  96.      * By default, at the end of the propagation, the propagator resets the initial state
  97.      * to the final state, thus allowing a new propagation to be started from there without
  98.      * recomputing the part already performed. Calling this method with {@code resetAtEnd} set
  99.      * to false changes prevents such reset.
  100.      * </p>
  101.      * @param resetAtEnd if true, at end of each propagation, the {@link
  102.      * #getInitialState() initial state} will be reset to the final state of
  103.      * the propagation, otherwise the initial state will be preserved
  104.      * @since 9.0
  105.      */
  106.     public void setResetAtEnd(final boolean resetAtEnd) {
  107.         this.resetAtEnd = resetAtEnd;
  108.     }

  109.     /** Initialize the mapper. */
  110.     protected void initMapper() {
  111.         stateMapper = createMapper(null, Double.NaN, null, null, null, null);
  112.     }

  113.     /**  {@inheritDoc} */
  114.     public void setAttitudeProvider(final AttitudeProvider attitudeProvider) {
  115.         super.setAttitudeProvider(attitudeProvider);
  116.         stateMapper = createMapper(stateMapper.getReferenceDate(), stateMapper.getMu(),
  117.                                    stateMapper.getOrbitType(), stateMapper.getPositionAngleType(),
  118.                                    attitudeProvider, stateMapper.getFrame());
  119.     }

  120.     /** Set propagation orbit type.
  121.      * @param orbitType orbit type to use for propagation, null for
  122.      * propagating using {@link org.orekit.utils.AbsolutePVCoordinates AbsolutePVCoordinates}
  123.      * rather than {@link org.orekit.orbits Orbit}
  124.      */
  125.     protected void setOrbitType(final OrbitType orbitType) {
  126.         stateMapper = createMapper(stateMapper.getReferenceDate(), stateMapper.getMu(),
  127.                                    orbitType, stateMapper.getPositionAngleType(),
  128.                                    stateMapper.getAttitudeProvider(), stateMapper.getFrame());
  129.     }

  130.     /** Get propagation parameter type.
  131.      * @return orbit type used for propagation, null for
  132.      * propagating using {@link org.orekit.utils.AbsolutePVCoordinates AbsolutePVCoordinates}
  133.      * rather than {@link org.orekit.orbits Orbit}
  134.      */
  135.     protected OrbitType getOrbitType() {
  136.         return stateMapper.getOrbitType();
  137.     }

  138.     /** Check if only the mean elements should be used in a semianalitical propagation.
  139.      * @return {@link PropagationType MEAN} if only mean elements have to be used or
  140.      *         {@link PropagationType OSCULATING} if osculating elements have to be also used.
  141.      */
  142.     protected PropagationType isMeanOrbit() {
  143.         return propagationType;
  144.     }

  145.     /** Set position angle type.
  146.      * <p>
  147.      * The position parameter type is meaningful only if {@link
  148.      * #getOrbitType() propagation orbit type}
  149.      * support it. As an example, it is not meaningful for propagation
  150.      * in {@link OrbitType#CARTESIAN Cartesian} parameters.
  151.      * </p>
  152.      * @param positionAngleType angle type to use for propagation
  153.      */
  154.     protected void setPositionAngleType(final PositionAngle positionAngleType) {
  155.         stateMapper = createMapper(stateMapper.getReferenceDate(), stateMapper.getMu(),
  156.                                    stateMapper.getOrbitType(), positionAngleType,
  157.                                    stateMapper.getAttitudeProvider(), stateMapper.getFrame());
  158.     }

  159.     /** Get propagation parameter type.
  160.      * @return angle type to use for propagation
  161.      */
  162.     protected PositionAngle getPositionAngleType() {
  163.         return stateMapper.getPositionAngleType();
  164.     }

  165.     /** Set the central attraction coefficient μ.
  166.      * @param mu central attraction coefficient (m³/s²)
  167.      */
  168.     public void setMu(final double mu) {
  169.         stateMapper = createMapper(stateMapper.getReferenceDate(), mu,
  170.                                    stateMapper.getOrbitType(), stateMapper.getPositionAngleType(),
  171.                                    stateMapper.getAttitudeProvider(), stateMapper.getFrame());
  172.     }

  173.     /** Get the central attraction coefficient μ.
  174.      * @return mu central attraction coefficient (m³/s²)
  175.      * @see #setMu(double)
  176.      */
  177.     public double getMu() {
  178.         return stateMapper.getMu();
  179.     }

  180.     /** Get the number of calls to the differential equations computation method.
  181.      * <p>The number of calls is reset each time the {@link #propagate(AbsoluteDate)}
  182.      * method is called.</p>
  183.      * @return number of calls to the differential equations computation method
  184.      */
  185.     public int getCalls() {
  186.         return calls;
  187.     }

  188.     /** {@inheritDoc} */
  189.     @Override
  190.     public boolean isAdditionalStateManaged(final String name) {

  191.         // first look at already integrated states
  192.         if (super.isAdditionalStateManaged(name)) {
  193.             return true;
  194.         }

  195.         // then look at states we integrate ourselves
  196.         for (final AdditionalEquations equation : additionalEquations) {
  197.             if (equation.getName().equals(name)) {
  198.                 return true;
  199.             }
  200.         }

  201.         return false;
  202.     }

  203.     /** {@inheritDoc} */
  204.     @Override
  205.     public String[] getManagedAdditionalStates() {
  206.         final String[] alreadyIntegrated = super.getManagedAdditionalStates();
  207.         final String[] managed = new String[alreadyIntegrated.length + additionalEquations.size()];
  208.         System.arraycopy(alreadyIntegrated, 0, managed, 0, alreadyIntegrated.length);
  209.         for (int i = 0; i < additionalEquations.size(); ++i) {
  210.             managed[i + alreadyIntegrated.length] = additionalEquations.get(i).getName();
  211.         }
  212.         return managed;
  213.     }

  214.     /** Add a set of user-specified equations to be integrated along with the orbit propagation.
  215.      * @param additional additional equations
  216.      */
  217.     public void addAdditionalEquations(final AdditionalEquations additional) {

  218.         // check if the name is already used
  219.         if (isAdditionalStateManaged(additional.getName())) {
  220.             // this set of equations is already registered, complain
  221.             throw new OrekitException(OrekitMessages.ADDITIONAL_STATE_NAME_ALREADY_IN_USE,
  222.                                       additional.getName());
  223.         }

  224.         // this is really a new set of equations, add it
  225.         additionalEquations.add(additional);

  226.     }

  227.     /** {@inheritDoc} */
  228.     public void addEventDetector(final EventDetector detector) {
  229.         detectors.add(detector);
  230.     }

  231.     /** {@inheritDoc} */
  232.     public Collection<EventDetector> getEventsDetectors() {
  233.         return Collections.unmodifiableCollection(detectors);
  234.     }

  235.     /** {@inheritDoc} */
  236.     public void clearEventsDetectors() {
  237.         detectors.clear();
  238.     }

  239.     /** Set up all user defined event detectors.
  240.      */
  241.     protected void setUpUserEventDetectors() {
  242.         for (final EventDetector detector : detectors) {
  243.             setUpEventDetector(integrator, detector);
  244.         }
  245.     }

  246.     /** Wrap an Orekit event detector and register it to the integrator.
  247.      * @param integ integrator into which event detector should be registered
  248.      * @param detector event detector to wrap
  249.      */
  250.     protected void setUpEventDetector(final ODEIntegrator integ, final EventDetector detector) {
  251.         integ.addEventHandler(new AdaptedEventDetector(detector),
  252.                               detector.getMaxCheckInterval(),
  253.                               detector.getThreshold(),
  254.                               detector.getMaxIterationCount());
  255.     }

  256.     /** {@inheritDoc}
  257.      * <p>Note that this method has the side effect of replacing the step handlers
  258.      * of the underlying integrator set up in the {@link
  259.      * #AbstractIntegratedPropagator(ODEIntegrator, PropagationType) constructor}. So if a specific
  260.      * step handler is needed, it should be added after this method has been callled.</p>
  261.      */
  262.     public void setSlaveMode() {
  263.         super.setSlaveMode();
  264.         if (integrator != null) {
  265.             integrator.clearStepHandlers();
  266.         }
  267.         modeHandler = null;
  268.     }

  269.     /** {@inheritDoc}
  270.      * <p>Note that this method has the side effect of replacing the step handlers
  271.      * of the underlying integrator set up in the {@link
  272.      * #AbstractIntegratedPropagator(ODEIntegrator, PropagationType) constructor}. So if a specific
  273.      * step handler is needed, it should be added after this method has been callled.</p>
  274.      */
  275.     public void setMasterMode(final OrekitStepHandler handler) {
  276.         super.setMasterMode(handler);
  277.         integrator.clearStepHandlers();
  278.         final AdaptedStepHandler wrapped = new AdaptedStepHandler(handler);
  279.         integrator.addStepHandler(wrapped);
  280.         modeHandler = wrapped;
  281.     }

  282.     /** {@inheritDoc}
  283.      * <p>Note that this method has the side effect of replacing the step handlers
  284.      * of the underlying integrator set up in the {@link
  285.      * #AbstractIntegratedPropagator(ODEIntegrator, PropagationType) constructor}. So if a specific
  286.      * step handler is needed, it should be added after this method has been called.</p>
  287.      */
  288.     public void setEphemerisMode() {
  289.         super.setEphemerisMode();
  290.         integrator.clearStepHandlers();
  291.         final EphemerisModeHandler ephemeris = new EphemerisModeHandler();
  292.         modeHandler = ephemeris;
  293.         integrator.addStepHandler(ephemeris);
  294.     }

  295.     /**
  296.      * {@inheritDoc}
  297.      *
  298.      * <p>Note that this method has the side effect of replacing the step handlers of the
  299.      * underlying integrator set up in the {@link #AbstractIntegratedPropagator(ODEIntegrator,
  300.      * PropagationType) constructor}.</p>
  301.      */
  302.     @Override
  303.     public void setEphemerisMode(final OrekitStepHandler handler) {
  304.         super.setEphemerisMode();
  305.         integrator.clearStepHandlers();
  306.         final EphemerisModeHandler ephemeris = new EphemerisModeHandler(handler);
  307.         modeHandler = ephemeris;
  308.         integrator.addStepHandler(ephemeris);
  309.     }

  310.     /** {@inheritDoc} */
  311.     public BoundedPropagator getGeneratedEphemeris()
  312.         throws IllegalStateException {
  313.         if (getMode() != EPHEMERIS_GENERATION_MODE) {
  314.             throw new OrekitIllegalStateException(OrekitMessages.PROPAGATOR_NOT_IN_EPHEMERIS_GENERATION_MODE);
  315.         }
  316.         return ((EphemerisModeHandler) modeHandler).getEphemeris();
  317.     }

  318.     /** Create a mapper between raw double components and spacecraft state.
  319.     /** Simple constructor.
  320.      * <p>
  321.      * The position parameter type is meaningful only if {@link
  322.      * #getOrbitType() propagation orbit type}
  323.      * support it. As an example, it is not meaningful for propagation
  324.      * in {@link OrbitType#CARTESIAN Cartesian} parameters.
  325.      * </p>
  326.      * @param referenceDate reference date
  327.      * @param mu central attraction coefficient (m³/s²)
  328.      * @param orbitType orbit type to use for mapping
  329.      * @param positionAngleType angle type to use for propagation
  330.      * @param attitudeProvider attitude provider
  331.      * @param frame inertial frame
  332.      * @return new mapper
  333.      */
  334.     protected abstract StateMapper createMapper(AbsoluteDate referenceDate, double mu,
  335.                                                 OrbitType orbitType, PositionAngle positionAngleType,
  336.                                                 AttitudeProvider attitudeProvider, Frame frame);

  337.     /** Get the differential equations to integrate (for main state only).
  338.      * @param integ numerical integrator to use for propagation.
  339.      * @return differential equations for main state
  340.      */
  341.     protected abstract MainStateEquations getMainStateEquations(ODEIntegrator integ);

  342.     /** {@inheritDoc} */
  343.     public SpacecraftState propagate(final AbsoluteDate target) {
  344.         if (getStartDate() == null) {
  345.             if (getInitialState() == null) {
  346.                 throw new OrekitException(OrekitMessages.INITIAL_STATE_NOT_SPECIFIED_FOR_ORBIT_PROPAGATION);
  347.             }
  348.             setStartDate(getInitialState().getDate());
  349.         }
  350.         return propagate(getStartDate(), target);
  351.     }

  352.     /** {@inheritDoc} */
  353.     public SpacecraftState propagate(final AbsoluteDate tStart, final AbsoluteDate tEnd) {

  354.         if (getInitialState() == null) {
  355.             throw new OrekitException(OrekitMessages.INITIAL_STATE_NOT_SPECIFIED_FOR_ORBIT_PROPAGATION);
  356.         }

  357.         if (!tStart.equals(getInitialState().getDate())) {
  358.             // if propagation start date is not initial date,
  359.             // propagate from initial to start date without event detection
  360.             propagate(tStart, false);
  361.         }

  362.         // propagate from start date to end date with event detection
  363.         return propagate(tEnd, true);

  364.     }

  365.     /** Propagation with or without event detection.
  366.      * @param tEnd target date to which orbit should be propagated
  367.      * @param activateHandlers if true, step and event handlers should be activated
  368.      * @return state at end of propagation
  369.      */
  370.     protected SpacecraftState propagate(final AbsoluteDate tEnd, final boolean activateHandlers) {
  371.         try {

  372.             initializePropagation();

  373.             if (getInitialState().getDate().equals(tEnd)) {
  374.                 // don't extrapolate
  375.                 return getInitialState();
  376.             }

  377.             // space dynamics view
  378.             stateMapper = createMapper(getInitialState().getDate(), stateMapper.getMu(),
  379.                                        stateMapper.getOrbitType(), stateMapper.getPositionAngleType(),
  380.                                        stateMapper.getAttitudeProvider(), getInitialState().getFrame());


  381.             if (Double.isNaN(getMu())) {
  382.                 setMu(getInitialState().getMu());
  383.             }

  384.             if (getInitialState().getMass() <= 0.0) {
  385.                 throw new OrekitException(OrekitMessages.SPACECRAFT_MASS_BECOMES_NEGATIVE,
  386.                                           getInitialState().getMass());
  387.             }

  388.             integrator.clearEventHandlers();

  389.             // set up events added by user, only if handlers are activated
  390.             if (activateHandlers) {
  391.                 setUpUserEventDetectors();
  392.             }

  393.             // convert space flight dynamics API to math API
  394.             final ODEState mathInitialState = createInitialState(getInitialIntegrationState());
  395.             final ExpandableODE mathODE = createODE(integrator, mathInitialState);
  396.             equationsMapper = mathODE.getMapper();

  397.             // initialize mode handler
  398.             if (modeHandler != null) {
  399.                 modeHandler.initialize(activateHandlers, tEnd);
  400.             }

  401.             // mathematical integration
  402.             final ODEStateAndDerivative mathFinalState;
  403.             beforeIntegration(getInitialState(), tEnd);
  404.             mathFinalState = integrator.integrate(mathODE, mathInitialState,
  405.                                                   tEnd.durationFrom(getInitialState().getDate()));
  406.             afterIntegration();

  407.             // get final state
  408.             SpacecraftState finalState =
  409.                             stateMapper.mapArrayToState(stateMapper.mapDoubleToDate(mathFinalState.getTime(),
  410.                                                                                     tEnd),
  411.                                                         mathFinalState.getPrimaryState(),
  412.                                                         mathFinalState.getPrimaryDerivative(),
  413.                                                         propagationType);

  414.             finalState = updateAdditionalStates(finalState);
  415.             for (int i = 0; i < additionalEquations.size(); ++i) {
  416.                 final double[] secondary = mathFinalState.getSecondaryState(i + 1);
  417.                 finalState = finalState.addAdditionalState(additionalEquations.get(i).getName(),
  418.                                                            secondary);
  419.             }
  420.             if (resetAtEnd) {
  421.                 resetInitialState(finalState);
  422.                 setStartDate(finalState.getDate());
  423.             }

  424.             return finalState;

  425.         } catch (MathRuntimeException mre) {
  426.             throw OrekitException.unwrap(mre);
  427.         }
  428.     }

  429.     /** Get the initial state for integration.
  430.      * @return initial state for integration
  431.      */
  432.     protected SpacecraftState getInitialIntegrationState() {
  433.         return getInitialState();
  434.     }

  435.     /** Create an initial state.
  436.      * @param initialState initial state in flight dynamics world
  437.      * @return initial state in mathematics world
  438.      */
  439.     private ODEState createInitialState(final SpacecraftState initialState) {

  440.         // retrieve initial state
  441.         final double[] primary  = new double[getBasicDimension()];
  442.         stateMapper.mapStateToArray(initialState, primary, null);

  443.         // secondary part of the ODE
  444.         final double[][] secondary = new double[additionalEquations.size()][];
  445.         for (int i = 0; i < additionalEquations.size(); ++i) {
  446.             final AdditionalEquations additional = additionalEquations.get(i);
  447.             secondary[i] = initialState.getAdditionalState(additional.getName());
  448.         }

  449.         return new ODEState(0.0, primary, secondary);

  450.     }

  451.     /** Create an ODE with all equations.
  452.      * @param integ numerical integrator to use for propagation.
  453.      * @param mathInitialState initial state
  454.      * @return a new ode
  455.      */
  456.     private ExpandableODE createODE(final ODEIntegrator integ,
  457.                                     final ODEState mathInitialState) {

  458.         final ExpandableODE ode =
  459.                 new ExpandableODE(new ConvertedMainStateEquations(getMainStateEquations(integ)));

  460.         // secondary part of the ODE
  461.         for (int i = 0; i < additionalEquations.size(); ++i) {
  462.             final AdditionalEquations additional = additionalEquations.get(i);
  463.             final SecondaryODE secondary =
  464.                     new ConvertedSecondaryStateEquations(additional,
  465.                                                          mathInitialState.getSecondaryStateDimension(i + 1));
  466.             ode.addSecondaryEquations(secondary);
  467.         }

  468.         return ode;

  469.     }

  470.     /** Method called just before integration.
  471.      * <p>
  472.      * The default implementation does nothing, it may be specialized in subclasses.
  473.      * </p>
  474.      * @param initialState initial state
  475.      * @param tEnd target date at which state should be propagated
  476.      */
  477.     protected void beforeIntegration(final SpacecraftState initialState,
  478.                                      final AbsoluteDate tEnd) {
  479.         // do nothing by default
  480.     }

  481.     /** Method called just after integration.
  482.      * <p>
  483.      * The default implementation does nothing, it may be specialized in subclasses.
  484.      * </p>
  485.      */
  486.     protected void afterIntegration() {
  487.         // do nothing by default
  488.     }

  489.     /** Get state vector dimension without additional parameters.
  490.      * @return state vector dimension without additional parameters.
  491.      */
  492.     public int getBasicDimension() {
  493.         return 7;

  494.     }

  495.     /** Get the integrator used by the propagator.
  496.      * @return the integrator.
  497.      */
  498.     protected ODEIntegrator getIntegrator() {
  499.         return integrator;
  500.     }

  501.     /** Get a complete state with all additional equations.
  502.      * @param t current value of the independent <I>time</I> variable
  503.      * @param y array containing the current value of the state vector
  504.      * @param yDot array containing the current value of the state vector derivative
  505.      * @return complete state
  506.      */
  507.     private SpacecraftState getCompleteState(final double t, final double[] y, final double[] yDot) {

  508.         // main state
  509.         SpacecraftState state = stateMapper.mapArrayToState(t, y, yDot, propagationType);

  510.         // pre-integrated additional states
  511.         state = updateAdditionalStates(state);

  512.         // additional states integrated here
  513.         if (!additionalEquations.isEmpty()) {

  514.             for (int i = 0; i < additionalEquations.size(); ++i) {
  515.                 state = state.addAdditionalState(additionalEquations.get(i).getName(),
  516.                                                  equationsMapper.extractEquationData(i + 1, y));
  517.             }

  518.         }

  519.         return state;

  520.     }

  521.     /** Differential equations for the main state (orbit, attitude and mass). */
  522.     public interface MainStateEquations {

  523.         /**
  524.          * Initialize the equations at the start of propagation. This method will be
  525.          * called before any calls to {@link #computeDerivatives(SpacecraftState)}.
  526.          *
  527.          * <p> The default implementation of this method does nothing.
  528.          *
  529.          * @param initialState initial state information at the start of propagation.
  530.          * @param target       date of propagation. Not equal to {@code
  531.          *                     initialState.getDate()}.
  532.          */
  533.         default void init(final SpacecraftState initialState, final AbsoluteDate target) {
  534.         }

  535.         /** Compute differential equations for main state.
  536.          * @param state current state
  537.          * @return derivatives of main state
  538.          */
  539.         double[] computeDerivatives(SpacecraftState state);

  540.     }

  541.     /** Differential equations for the main state (orbit, attitude and mass), with converted API. */
  542.     private class ConvertedMainStateEquations implements OrdinaryDifferentialEquation {

  543.         /** Main state equations. */
  544.         private final MainStateEquations main;

  545.         /** Simple constructor.
  546.          * @param main main state equations
  547.          */
  548.         ConvertedMainStateEquations(final MainStateEquations main) {
  549.             this.main = main;
  550.             calls = 0;
  551.         }

  552.         /** {@inheritDoc} */
  553.         public int getDimension() {
  554.             return getBasicDimension();
  555.         }

  556.         @Override
  557.         public void init(final double t0, final double[] y0, final double finalTime) {
  558.             // update space dynamics view
  559.             SpacecraftState initialState = stateMapper.mapArrayToState(t0, y0, null, PropagationType.MEAN);
  560.             initialState = updateAdditionalStates(initialState);
  561.             final AbsoluteDate target = stateMapper.mapDoubleToDate(finalTime);
  562.             main.init(initialState, target);
  563.         }

  564.         /** {@inheritDoc} */
  565.         public double[] computeDerivatives(final double t, final double[] y) {

  566.             // increment calls counter
  567.             ++calls;

  568.             // update space dynamics view
  569.             SpacecraftState currentState = stateMapper.mapArrayToState(t, y, null, PropagationType.MEAN);
  570.             currentState = updateAdditionalStates(currentState);

  571.             // compute main state differentials
  572.             return main.computeDerivatives(currentState);

  573.         }

  574.     }

  575.     /** Differential equations for the secondary state (Jacobians, user variables ...), with converted API. */
  576.     private class ConvertedSecondaryStateEquations implements SecondaryODE {

  577.         /** Additional equations. */
  578.         private final AdditionalEquations equations;

  579.         /** Dimension of the additional state. */
  580.         private final int dimension;

  581.         /** Simple constructor.
  582.          * @param equations additional equations
  583.          * @param dimension dimension of the additional state
  584.          */
  585.         ConvertedSecondaryStateEquations(final AdditionalEquations equations,
  586.                                          final int dimension) {
  587.             this.equations = equations;
  588.             this.dimension = dimension;
  589.         }

  590.         /** {@inheritDoc} */
  591.         @Override
  592.         public int getDimension() {
  593.             return dimension;
  594.         }

  595.         /** {@inheritDoc} */
  596.         @Override
  597.         public void init(final double t0, final double[] primary0,
  598.                          final double[] secondary0, final double finalTime) {
  599.             // update space dynamics view
  600.             SpacecraftState initialState = stateMapper.mapArrayToState(t0, primary0, null, PropagationType.MEAN);
  601.             initialState = updateAdditionalStates(initialState);
  602.             initialState = initialState.addAdditionalState(equations.getName(), secondary0);
  603.             final AbsoluteDate target = stateMapper.mapDoubleToDate(finalTime);
  604.             equations.init(initialState, target);

  605.         }

  606.         /** {@inheritDoc} */
  607.         @Override
  608.         public double[] computeDerivatives(final double t, final double[] primary,
  609.                                            final double[] primaryDot, final double[] secondary) {

  610.             // update space dynamics view
  611.             SpacecraftState currentState = stateMapper.mapArrayToState(t, primary, primaryDot, PropagationType.MEAN);
  612.             currentState = updateAdditionalStates(currentState);
  613.             currentState = currentState.addAdditionalState(equations.getName(), secondary);

  614.             // compute additional derivatives
  615.             final double[] secondaryDot = new double[secondary.length];
  616.             final double[] additionalMainDot =
  617.                             equations.computeDerivatives(currentState, secondaryDot);
  618.             if (additionalMainDot != null) {
  619.                 // the additional equations have an effect on main equations
  620.                 for (int i = 0; i < additionalMainDot.length; ++i) {
  621.                     primaryDot[i] += additionalMainDot[i];
  622.                 }
  623.             }

  624.             return secondaryDot;

  625.         }

  626.     }

  627.     /** Adapt an {@link org.orekit.propagation.events.EventDetector}
  628.      * to Hipparchus {@link org.hipparchus.ode.events.ODEEventHandler} interface.
  629.      * @author Fabien Maussion
  630.      */
  631.     private class AdaptedEventDetector implements ODEEventHandler {

  632.         /** Underlying event detector. */
  633.         private final EventDetector detector;

  634.         /** Time of the previous call to g. */
  635.         private double lastT;

  636.         /** Value from the previous call to g. */
  637.         private double lastG;

  638.         /** Build a wrapped event detector.
  639.          * @param detector event detector to wrap
  640.         */
  641.         AdaptedEventDetector(final EventDetector detector) {
  642.             this.detector = detector;
  643.             this.lastT    = Double.NaN;
  644.             this.lastG    = Double.NaN;
  645.         }

  646.         /** {@inheritDoc} */
  647.         public void init(final ODEStateAndDerivative s0, final double t) {
  648.             detector.init(getCompleteState(s0.getTime(), s0.getCompleteState(), s0.getCompleteDerivative()),
  649.                           stateMapper.mapDoubleToDate(t));
  650.             this.lastT = Double.NaN;
  651.             this.lastG = Double.NaN;
  652.         }

  653.         /** {@inheritDoc} */
  654.         public double g(final ODEStateAndDerivative s) {
  655.             if (!Precision.equals(lastT, s.getTime(), 0)) {
  656.                 lastT = s.getTime();
  657.                 lastG = detector.g(getCompleteState(s.getTime(), s.getCompleteState(), s.getCompleteDerivative()));
  658.             }
  659.             return lastG;
  660.         }

  661.         /** {@inheritDoc} */
  662.         public Action eventOccurred(final ODEStateAndDerivative s, final boolean increasing) {
  663.             return detector.eventOccurred(
  664.                     getCompleteState(
  665.                             s.getTime(),
  666.                             s.getCompleteState(),
  667.                             s.getCompleteDerivative()),
  668.                     increasing);
  669.         }

  670.         /** {@inheritDoc} */
  671.         public ODEState resetState(final ODEStateAndDerivative s) {

  672.             final SpacecraftState oldState = getCompleteState(s.getTime(), s.getCompleteState(), s.getCompleteDerivative());
  673.             final SpacecraftState newState = detector.resetState(oldState);
  674.             stateChanged(newState);

  675.             // main part
  676.             final double[] primary    = new double[s.getPrimaryStateDimension()];
  677.             stateMapper.mapStateToArray(newState, primary, null);

  678.             // secondary part
  679.             final double[][] secondary    = new double[additionalEquations.size()][];
  680.             for (int i = 0; i < additionalEquations.size(); ++i) {
  681.                 secondary[i] = newState.getAdditionalState(additionalEquations.get(i).getName());
  682.             }

  683.             return new ODEState(newState.getDate().durationFrom(getStartDate()),
  684.                                 primary, secondary);

  685.         }

  686.     }

  687.     /** Adapt an {@link org.orekit.propagation.sampling.OrekitStepHandler}
  688.      * to Hipparchus {@link ODEStepHandler} interface.
  689.      * @author Luc Maisonobe
  690.      */
  691.     private class AdaptedStepHandler implements ODEStepHandler, ModeHandler {

  692.         /** Underlying handler. */
  693.         private final OrekitStepHandler handler;

  694.         /** Flag for handler . */
  695.         private boolean activate;

  696.         /** Build an instance.
  697.          * @param handler underlying handler to wrap
  698.          */
  699.         AdaptedStepHandler(final OrekitStepHandler handler) {
  700.             this.handler = handler;
  701.         }

  702.         /** {@inheritDoc} */
  703.         public void initialize(final boolean activateHandlers,
  704.                                final AbsoluteDate targetDate) {
  705.             this.activate = activateHandlers;
  706.         }

  707.         /** {@inheritDoc} */
  708.         public void init(final ODEStateAndDerivative s0, final double t) {
  709.             if (activate) {
  710.                 handler.init(getCompleteState(s0.getTime(), s0.getCompleteState(), s0.getCompleteDerivative()),
  711.                              stateMapper.mapDoubleToDate(t));
  712.             }
  713.         }

  714.         /** {@inheritDoc} */
  715.         public void handleStep(final ODEStateInterpolator interpolator, final boolean isLast) {
  716.             if (activate) {
  717.                 handler.handleStep(new AdaptedStepInterpolator(interpolator), isLast);
  718.             }
  719.         }

  720.     }

  721.     /** Adapt an Hipparchus {@link ODEStateInterpolator}
  722.      * to an orekit {@link OrekitStepInterpolator} interface.
  723.      * @author Luc Maisonobe
  724.      */
  725.     private class AdaptedStepInterpolator implements OrekitStepInterpolator {

  726.         /** Underlying raw rawInterpolator. */
  727.         private final ODEStateInterpolator mathInterpolator;

  728.         /** Simple constructor.
  729.          * @param mathInterpolator underlying raw interpolator
  730.          */
  731.         AdaptedStepInterpolator(final ODEStateInterpolator mathInterpolator) {
  732.             this.mathInterpolator = mathInterpolator;
  733.         }

  734.         /** {@inheritDoc}} */
  735.         @Override
  736.         public SpacecraftState getPreviousState() {
  737.             return convert(mathInterpolator.getPreviousState());
  738.         }

  739.         /** {@inheritDoc}} */
  740.         @Override
  741.         public boolean isPreviousStateInterpolated() {
  742.             return mathInterpolator.isPreviousStateInterpolated();
  743.         }

  744.         /** {@inheritDoc}} */
  745.         @Override
  746.         public SpacecraftState getCurrentState() {
  747.             return convert(mathInterpolator.getCurrentState());
  748.         }

  749.         /** {@inheritDoc}} */
  750.         @Override
  751.         public boolean isCurrentStateInterpolated() {
  752.             return mathInterpolator.isCurrentStateInterpolated();
  753.         }

  754.         /** {@inheritDoc}} */
  755.         @Override
  756.         public SpacecraftState getInterpolatedState(final AbsoluteDate date) {
  757.             return convert(mathInterpolator.getInterpolatedState(date.durationFrom(stateMapper.getReferenceDate())));
  758.         }

  759.         /** Convert a state from mathematical world to space flight dynamics world.
  760.          * @param os mathematical state
  761.          * @return space flight dynamics state
  762.          */
  763.         private SpacecraftState convert(final ODEStateAndDerivative os) {

  764.             SpacecraftState s =
  765.                             stateMapper.mapArrayToState(os.getTime(),
  766.                                                         os.getPrimaryState(),
  767.                                                         os.getPrimaryDerivative(),
  768.                                                         propagationType);
  769.             s = updateAdditionalStates(s);
  770.             for (int i = 0; i < additionalEquations.size(); ++i) {
  771.                 final double[] secondary = os.getSecondaryState(i + 1);
  772.                 s = s.addAdditionalState(additionalEquations.get(i).getName(), secondary);
  773.             }

  774.             return s;

  775.         }

  776.         /** Convert a state from space flight dynamics world to mathematical world.
  777.          * @param state space flight dynamics state
  778.          * @return mathematical state
  779.          */
  780.         private ODEStateAndDerivative convert(final SpacecraftState state) {

  781.             // retrieve initial state
  782.             final double[] primary    = new double[getBasicDimension()];
  783.             final double[] primaryDot = new double[getBasicDimension()];
  784.             stateMapper.mapStateToArray(state, primary, primaryDot);

  785.             // secondary part of the ODE
  786.             final double[][] secondary    = new double[additionalEquations.size()][];
  787.             for (int i = 0; i < additionalEquations.size(); ++i) {
  788.                 final AdditionalEquations additional = additionalEquations.get(i);
  789.                 secondary[i] = state.getAdditionalState(additional.getName());
  790.             }

  791.             return new ODEStateAndDerivative(stateMapper.mapDateToDouble(state.getDate()),
  792.                                              primary, primaryDot,
  793.                                              secondary, null);

  794.         }

  795.         /** {@inheritDoc}} */
  796.         @Override
  797.         public boolean isForward() {
  798.             return mathInterpolator.isForward();
  799.         }

  800.         /** {@inheritDoc}} */
  801.         @Override
  802.         public AdaptedStepInterpolator restrictStep(final SpacecraftState newPreviousState,
  803.                                                     final SpacecraftState newCurrentState) {
  804.             try {
  805.                 final AbstractODEStateInterpolator aosi = (AbstractODEStateInterpolator) mathInterpolator;
  806.                 return new AdaptedStepInterpolator(aosi.restrictStep(convert(newPreviousState),
  807.                                                                      convert(newCurrentState)));
  808.             } catch (ClassCastException cce) {
  809.                 // this should never happen
  810.                 throw new OrekitInternalError(cce);
  811.             }
  812.         }

  813.     }

  814.     private class EphemerisModeHandler implements ModeHandler, ODEStepHandler {

  815.         /** Underlying raw mathematical model. */
  816.         private DenseOutputModel model;

  817.         /** Generated ephemeris. */
  818.         private BoundedPropagator ephemeris;

  819.         /** Flag for handler . */
  820.         private boolean activate;

  821.         /** the user supplied end date. Propagation may not end on this date. */
  822.         private AbsoluteDate endDate;

  823.         /** User's integration step handler. May be null. */
  824.         private final AdaptedStepHandler handler;

  825.         /** Creates a new instance of EphemerisModeHandler which must be
  826.          *  filled by the propagator.
  827.          */
  828.         EphemerisModeHandler() {
  829.             this.handler = null;
  830.         }

  831.         /** Creates a new instance of EphemerisModeHandler which must be
  832.          *  filled by the propagator.
  833.          *  @param handler the handler to notify of every integrator step.
  834.          */
  835.         EphemerisModeHandler(final OrekitStepHandler handler) {
  836.             this.handler = new AdaptedStepHandler(handler);
  837.         }

  838.         /** {@inheritDoc} */
  839.         public void initialize(final boolean activateHandlers,
  840.                                final AbsoluteDate targetDate) {
  841.             this.activate = activateHandlers;
  842.             this.model    = new DenseOutputModel();
  843.             this.endDate  = targetDate;

  844.             // ephemeris will be generated when last step is processed
  845.             this.ephemeris = null;
  846.             if (this.handler != null) {
  847.                 this.handler.initialize(activateHandlers, targetDate);
  848.             }
  849.         }

  850.         /** Get the generated ephemeris.
  851.          * @return a new instance of the generated ephemeris
  852.          */
  853.         public BoundedPropagator getEphemeris() {
  854.             return ephemeris;
  855.         }

  856.         /** {@inheritDoc} */
  857.         public void handleStep(final ODEStateInterpolator interpolator, final boolean isLast) {
  858.             if (activate) {
  859.                 if (this.handler != null) {
  860.                     this.handler.handleStep(interpolator, isLast);
  861.                 }

  862.                 model.handleStep(interpolator, isLast);
  863.                 if (isLast) {

  864.                     // set up the boundary dates
  865.                     final double tI = model.getInitialTime();
  866.                     final double tF = model.getFinalTime();
  867.                     // tI is almost? always zero
  868.                     final AbsoluteDate startDate =
  869.                                     stateMapper.mapDoubleToDate(tI);
  870.                     final AbsoluteDate finalDate =
  871.                                     stateMapper.mapDoubleToDate(tF, this.endDate);
  872.                     final AbsoluteDate minDate;
  873.                     final AbsoluteDate maxDate;
  874.                     if (tF < tI) {
  875.                         minDate = finalDate;
  876.                         maxDate = startDate;
  877.                     } else {
  878.                         minDate = startDate;
  879.                         maxDate = finalDate;
  880.                     }

  881.                     // get the initial additional states that are not managed
  882.                     final Map<String, double[]> unmanaged = new HashMap<String, double[]>();
  883.                     for (final Map.Entry<String, double[]> initial : getInitialState().getAdditionalStates().entrySet()) {
  884.                         if (!isAdditionalStateManaged(initial.getKey())) {
  885.                             // this additional state was in the initial state, but is unknown to the propagator
  886.                             // we simply copy its initial value as is
  887.                             unmanaged.put(initial.getKey(), initial.getValue());
  888.                         }
  889.                     }

  890.                     // get the names of additional states managed by differential equations
  891.                     final String[] names = new String[additionalEquations.size()];
  892.                     for (int i = 0; i < names.length; ++i) {
  893.                         names[i] = additionalEquations.get(i).getName();
  894.                     }

  895.                     // create the ephemeris
  896.                     ephemeris = new IntegratedEphemeris(startDate, minDate, maxDate,
  897.                                                         stateMapper, propagationType, model, unmanaged,
  898.                                                         getAdditionalStateProviders(), names);

  899.                 }
  900.             }

  901.         }

  902.         /** {@inheritDoc} */
  903.         public void init(final ODEStateAndDerivative s0, final double t) {
  904.             model.init(s0, t);
  905.             if (this.handler != null) {
  906.                 this.handler.init(s0, t);
  907.             }
  908.         }

  909.     }

  910. }