AbstractIntegratedPropagator.java

  1. /* Copyright 2002-2017 CS Systèmes d'Information
  2.  * Licensed to CS Systèmes d'Information (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.OrekitExceptionWrapper;
  42. import org.orekit.errors.OrekitIllegalStateException;
  43. import org.orekit.errors.OrekitInternalError;
  44. import org.orekit.errors.OrekitMessages;
  45. import org.orekit.frames.Frame;
  46. import org.orekit.orbits.Orbit;
  47. import org.orekit.orbits.OrbitType;
  48. import org.orekit.orbits.PositionAngle;
  49. import org.orekit.propagation.AbstractPropagator;
  50. import org.orekit.propagation.BoundedPropagator;
  51. import org.orekit.propagation.SpacecraftState;
  52. import org.orekit.propagation.events.EventDetector;
  53. import org.orekit.propagation.events.handlers.EventHandler;
  54. import org.orekit.propagation.sampling.OrekitStepHandler;
  55. import org.orekit.propagation.sampling.OrekitStepInterpolator;
  56. import org.orekit.time.AbsoluteDate;


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

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

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

  66.     /** Mode handler. */
  67.     private ModeHandler modeHandler;

  68.     /** Additional equations. */
  69.     private List<AdditionalEquations> additionalEquations;

  70.     /** Counter for differential equations calls. */
  71.     private int calls;

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

  74.     /** Equations mapper. */
  75.     private EquationsMapper equationsMapper;

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

  78.     /** Output only the mean orbit. <br/>
  79.      * <p>
  80.      * This is used only in the case of semianalitical propagators where there is a clear separation between
  81.      * mean and short periodic elements. It is ignored by the Numerical propagator.
  82.      * </p>
  83.      */
  84.     private boolean meanOrbit;

  85.     /** Build a new instance.
  86.      * @param integrator numerical integrator to use for propagation.
  87.      * @param meanOrbit output only the mean orbit.
  88.      */
  89.     protected AbstractIntegratedPropagator(final ODEIntegrator integrator, final boolean meanOrbit) {
  90.         detectors           = new ArrayList<EventDetector>();
  91.         additionalEquations = new ArrayList<AdditionalEquations>();
  92.         this.integrator     = integrator;
  93.         this.meanOrbit      = meanOrbit;
  94.         this.resetAtEnd     = true;
  95.     }

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

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

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

  122.     /** Set propagation orbit type.
  123.      * @param orbitType orbit type to use for propagation
  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
  132.      */
  133.     protected OrbitType getOrbitType() {
  134.         return stateMapper.getOrbitType();
  135.     }

  136.     /** Check if only the mean elements should be used in a semianalitical propagation.
  137.      * @return true if only mean elements have to be used
  138.      */
  139.     protected boolean isMeanOrbit() {
  140.         return meanOrbit;
  141.     }

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

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

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

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

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

  185.     /** {@inheritDoc} */
  186.     @Override
  187.     public boolean isAdditionalStateManaged(final String name) {

  188.         // first look at already integrated states
  189.         if (super.isAdditionalStateManaged(name)) {
  190.             return true;
  191.         }

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

  198.         return false;
  199.     }

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

  211.     /** Add a set of user-specified equations to be integrated along with the orbit propagation.
  212.      * @param additional additional equations
  213.      * @exception OrekitException if a set of equations with the same name is already present
  214.      */
  215.     public void addAdditionalEquations(final AdditionalEquations additional)
  216.         throws OrekitException {

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

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

  225.     }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  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.      * @exception OrekitException if orbit cannot be propagated
  370.      */
  371.     protected SpacecraftState propagate(final AbsoluteDate tEnd, final boolean activateHandlers)
  372.         throws OrekitException {
  373.         try {

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

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


  382.             // set propagation orbit type
  383.             final Orbit initialOrbit = stateMapper.getOrbitType().convertType(getInitialState().getOrbit());
  384.             if (Double.isNaN(getMu())) {
  385.                 setMu(initialOrbit.getMu());
  386.             }

  387.             if (getInitialState().getMass() <= 0.0) {
  388.                 throw new OrekitException(OrekitMessages.SPACECRAFT_MASS_BECOMES_NEGATIVE,
  389.                                           getInitialState().getMass());
  390.             }

  391.             integrator.clearEventHandlers();

  392.             // set up events added by user
  393.             setUpUserEventDetectors();

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

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

  402.             // mathematical integration
  403.             final ODEStateAndDerivative mathFinalState;
  404.             try {
  405.                 beforeIntegration(getInitialState(), tEnd);
  406.                 mathFinalState = integrator.integrate(mathODE, mathInitialState,
  407.                                                       tEnd.durationFrom(getInitialState().getDate()));
  408.                 afterIntegration();
  409.             } catch (OrekitExceptionWrapper oew) {
  410.                 throw oew.getException();
  411.             }

  412.             // get final state
  413.             SpacecraftState finalState =
  414.                             stateMapper.mapArrayToState(stateMapper.mapDoubleToDate(mathFinalState.getTime(),
  415.                                                                                     tEnd),
  416.                                                         mathFinalState.getPrimaryState(),
  417.                                                         mathFinalState.getPrimaryDerivative(),
  418.                                                         meanOrbit);
  419.             finalState = updateAdditionalStates(finalState);
  420.             for (int i = 0; i < additionalEquations.size(); ++i) {
  421.                 final double[] secondary = mathFinalState.getSecondaryState(i + 1);
  422.                 finalState = finalState.addAdditionalState(additionalEquations.get(i).getName(),
  423.                                                            secondary);
  424.             }
  425.             if (resetAtEnd) {
  426.                 resetInitialState(finalState);
  427.                 setStartDate(finalState.getDate());
  428.             }

  429.             return finalState;

  430.         } catch (MathRuntimeException mre) {
  431.             throw OrekitException.unwrap(mre);
  432.         }
  433.     }

  434.     /** Get the initial state for integration.
  435.      * @return initial state for integration
  436.      * @exception OrekitException if initial state cannot be retrieved
  437.      */
  438.     protected SpacecraftState getInitialIntegrationState() throws OrekitException {
  439.         return getInitialState();
  440.     }

  441.     /** Create an initial state.
  442.      * @param initialState initial state in flight dynamics world
  443.      * @return initial state in mathematics world
  444.      * @exception OrekitException if initial state cannot be mapped
  445.      */
  446.     private ODEState createInitialState(final SpacecraftState initialState)
  447.         throws OrekitException {

  448.         // retrieve initial state
  449.         final double[] primary  = new double[getBasicDimension()];
  450.         stateMapper.mapStateToArray(initialState, primary, null);

  451.         // secondary part of the ODE
  452.         final double[][] secondary = new double[additionalEquations.size()][];
  453.         for (int i = 0; i < additionalEquations.size(); ++i) {
  454.             final AdditionalEquations additional = additionalEquations.get(i);
  455.             secondary[i] = initialState.getAdditionalState(additional.getName());
  456.         }

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

  458.     }

  459.     /** Create an ODE with all equations.
  460.      * @param integ numerical integrator to use for propagation.
  461.      * @param mathInitialState initial state
  462.      * @return a new ode
  463.      * @exception OrekitException if initial state cannot be mapped
  464.      */
  465.     private ExpandableODE createODE(final ODEIntegrator integ,
  466.                                     final ODEState mathInitialState)
  467.         throws OrekitException {

  468.         final ExpandableODE ode =
  469.                 new ExpandableODE(new ConvertedMainStateEquations(getMainStateEquations(integ)));

  470.         // secondary part of the ODE
  471.         for (int i = 0; i < additionalEquations.size(); ++i) {
  472.             final AdditionalEquations additional = additionalEquations.get(i);
  473.             final SecondaryODE secondary =
  474.                     new ConvertedSecondaryStateEquations(additional,
  475.                                                          mathInitialState.getSecondaryStateDimension(i + 1));
  476.             ode.addSecondaryEquations(secondary);
  477.         }

  478.         return ode;

  479.     }

  480.     /** Method called just before integration.
  481.      * <p>
  482.      * The default implementation does nothing, it may be specialized in subclasses.
  483.      * </p>
  484.      * @param initialState initial state
  485.      * @param tEnd target date at which state should be propagated
  486.      * @exception OrekitException if hook cannot be run
  487.      */
  488.     protected void beforeIntegration(final SpacecraftState initialState,
  489.                                      final AbsoluteDate tEnd)
  490.         throws OrekitException {
  491.         // do nothing by default
  492.     }

  493.     /** Method called just after integration.
  494.      * <p>
  495.      * The default implementation does nothing, it may be specialized in subclasses.
  496.      * </p>
  497.      * @exception OrekitException if hook cannot be run
  498.      */
  499.     protected void afterIntegration()
  500.         throws OrekitException {
  501.         // do nothing by default
  502.     }

  503.     /** Get state vector dimension without additional parameters.
  504.      * @return state vector dimension without additional parameters.
  505.      */
  506.     public int getBasicDimension() {
  507.         return 7;

  508.     }

  509.     /** Get the integrator used by the propagator.
  510.      * @return the integrator.
  511.      */
  512.     protected ODEIntegrator getIntegrator() {
  513.         return integrator;
  514.     }

  515.     /** Get a complete state with all additional equations.
  516.      * @param t current value of the independent <I>time</I> variable
  517.      * @param y array containing the current value of the state vector
  518.      * @param yDot array containing the current value of the state vector derivative
  519.      * @return complete state
  520.      * @exception OrekitException if state cannot be mapped
  521.      */
  522.     private SpacecraftState getCompleteState(final double t, final double[] y, final double[] yDot)
  523.         throws OrekitException {

  524.         // main state
  525.         SpacecraftState state = stateMapper.mapArrayToState(t, y, yDot, true);  //not sure of the mean orbit, should be true

  526.         // pre-integrated additional states
  527.         state = updateAdditionalStates(state);

  528.         // additional states integrated here
  529.         if (!additionalEquations.isEmpty()) {

  530.             for (int i = 0; i < additionalEquations.size(); ++i) {
  531.                 state = state.addAdditionalState(additionalEquations.get(i).getName(),
  532.                                                  equationsMapper.extractEquationData(i + 1, y));
  533.             }

  534.         }

  535.         return state;

  536.     }

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

  539.         /**
  540.          * Initialize the equations at the start of propagation. This method will be
  541.          * called before any calls to {@link #computeDerivatives(SpacecraftState)}.
  542.          *
  543.          * <p> The default implementation of this method does nothing.
  544.          *
  545.          * @param initialState initial state information at the start of propagation.
  546.          * @param target       date of propagation. Not equal to {@code
  547.          *                     initialState.getDate()}.
  548.          * @throws OrekitException if there is an Orekit related error during
  549.          *                         initialization.
  550.          */
  551.         default void init(SpacecraftState initialState, AbsoluteDate target)
  552.                 throws OrekitException {
  553.         }

  554.         /** Compute differential equations for main state.
  555.          * @param state current state
  556.          * @return derivatives of main state
  557.          * @throws OrekitException if differentials cannot be computed
  558.          */
  559.         double[] computeDerivatives(SpacecraftState state) throws OrekitException;

  560.     }

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

  563.         /** Main state equations. */
  564.         private final MainStateEquations main;

  565.         /** Simple constructor.
  566.          * @param main main state equations
  567.          */
  568.         ConvertedMainStateEquations(final MainStateEquations main) {
  569.             this.main = main;
  570.             calls = 0;
  571.         }

  572.         /** {@inheritDoc} */
  573.         public int getDimension() {
  574.             return getBasicDimension();
  575.         }

  576.         @Override
  577.         public void init(final double t0, final double[] y0, final double finalTime) {
  578.             try {
  579.                 // update space dynamics view
  580.                 // use only ODE elements
  581.                 SpacecraftState initialState = stateMapper.mapArrayToState(t0, y0, null, true);
  582.                 initialState = updateAdditionalStates(initialState);
  583.                 final AbsoluteDate target = stateMapper.mapDoubleToDate(finalTime);
  584.                 main.init(initialState, target);
  585.             } catch (OrekitException oe) {
  586.                 throw new OrekitExceptionWrapper(oe);
  587.             }
  588.         }

  589.         /** {@inheritDoc} */
  590.         public double[] computeDerivatives(final double t, final double[] y)
  591.             throws OrekitExceptionWrapper {
  592.             try {

  593.                 // increment calls counter
  594.                 ++calls;

  595.                 // update space dynamics view
  596.                 // use only ODE elements
  597.                 SpacecraftState currentState = stateMapper.mapArrayToState(t, y, null, true);
  598.                 currentState = updateAdditionalStates(currentState);

  599.                 // compute main state differentials
  600.                 return main.computeDerivatives(currentState);

  601.             } catch (OrekitException oe) {
  602.                 throw new OrekitExceptionWrapper(oe);
  603.             }
  604.         }

  605.     }

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

  608.         /** Additional equations. */
  609.         private final AdditionalEquations equations;

  610.         /** Dimension of the additional state. */
  611.         private final int dimension;

  612.         /** Simple constructor.
  613.          * @param equations additional equations
  614.          * @param dimension dimension of the additional state
  615.          */
  616.         ConvertedSecondaryStateEquations(final AdditionalEquations equations,
  617.                                          final int dimension) {
  618.             this.equations = equations;
  619.             this.dimension = dimension;
  620.         }

  621.         /** {@inheritDoc} */
  622.         public int getDimension() {
  623.             return dimension;
  624.         }

  625.         /** {@inheritDoc} */
  626.         public double[] computeDerivatives(final double t, final double[] primary,
  627.                                            final double[] primaryDot, final double[] secondary)
  628.             throws OrekitExceptionWrapper {

  629.             try {

  630.                 // update space dynamics view
  631.                 // the state contains only the ODE elements
  632.                 SpacecraftState currentState = stateMapper.mapArrayToState(t, primary, primaryDot, true);
  633.                 currentState = updateAdditionalStates(currentState);
  634.                 currentState = currentState.addAdditionalState(equations.getName(), secondary);

  635.                 // compute additional derivatives
  636.                 final double[] secondaryDot = new double[secondary.length];
  637.                 final double[] additionalMainDot =
  638.                         equations.computeDerivatives(currentState, secondaryDot);
  639.                 if (additionalMainDot != null) {
  640.                     // the additional equations have an effect on main equations
  641.                     for (int i = 0; i < additionalMainDot.length; ++i) {
  642.                         primaryDot[i] += additionalMainDot[i];
  643.                     }
  644.                 }
  645.                 return secondaryDot;

  646.             } catch (OrekitException oe) {
  647.                 throw new OrekitExceptionWrapper(oe);
  648.             }

  649.         }

  650.     }

  651.     /** Adapt an {@link org.orekit.propagation.events.EventDetector}
  652.      * to Hipparchus {@link org.hipparchus.ode.events.ODEEventHandler} interface.
  653.      * @param <T> class type for the generic version
  654.      * @author Fabien Maussion
  655.      */
  656.     private class AdaptedEventDetector implements ODEEventHandler {

  657.         /** Underlying event detector. */
  658.         private final EventDetector detector;

  659.         /** Time of the previous call to g. */
  660.         private double lastT;

  661.         /** Value from the previous call to g. */
  662.         private double lastG;

  663.         /** Build a wrapped event detector.
  664.          * @param detector event detector to wrap
  665.         */
  666.         AdaptedEventDetector(final EventDetector detector) {
  667.             this.detector = detector;
  668.             this.lastT    = Double.NaN;
  669.             this.lastG    = Double.NaN;
  670.         }

  671.         /** {@inheritDoc} */
  672.         public void init(final ODEStateAndDerivative s0, final double t) {
  673.             try {

  674.                 detector.init(getCompleteState(s0.getTime(), s0.getCompleteState(), s0.getCompleteDerivative()),
  675.                               stateMapper.mapDoubleToDate(t));
  676.                 this.lastT = Double.NaN;
  677.                 this.lastG = Double.NaN;

  678.             } catch (OrekitException oe) {
  679.                 throw new OrekitExceptionWrapper(oe);
  680.             }
  681.         }

  682.         /** {@inheritDoc} */
  683.         public double g(final ODEStateAndDerivative s) {
  684.             try {
  685.                 if (!Precision.equals(lastT, s.getTime(), 0)) {
  686.                     lastT = s.getTime();
  687.                     lastG = detector.g(getCompleteState(s.getTime(), s.getCompleteState(), s.getCompleteDerivative()));
  688.                 }
  689.                 return lastG;
  690.             } catch (OrekitException oe) {
  691.                 throw new OrekitExceptionWrapper(oe);
  692.             }
  693.         }

  694.         /** {@inheritDoc} */
  695.         public Action eventOccurred(final ODEStateAndDerivative s, final boolean increasing) {
  696.             try {

  697.                 final EventHandler.Action whatNext = detector.eventOccurred(getCompleteState(s.getTime(),
  698.                                                                                              s.getCompleteState(),
  699.                                                                                              s.getCompleteDerivative()),
  700.                                                                             increasing);

  701.                 switch (whatNext) {
  702.                     case STOP :
  703.                         return Action.STOP;
  704.                     case RESET_STATE :
  705.                         return Action.RESET_STATE;
  706.                     case RESET_DERIVATIVES :
  707.                         return Action.RESET_DERIVATIVES;
  708.                     default :
  709.                         return Action.CONTINUE;
  710.                 }
  711.             } catch (OrekitException oe) {
  712.                 throw new OrekitExceptionWrapper(oe);
  713.             }
  714.         }

  715.         /** {@inheritDoc} */
  716.         public ODEState resetState(final ODEStateAndDerivative s) {
  717.             try {

  718.                 final SpacecraftState oldState = getCompleteState(s.getTime(), s.getCompleteState(), s.getCompleteDerivative());
  719.                 final SpacecraftState newState = detector.resetState(oldState);

  720.                 // main part
  721.                 final double[] primary    = new double[s.getPrimaryStateDimension()];
  722.                 stateMapper.mapStateToArray(newState, primary, null);

  723.                 // secondary part
  724.                 final double[][] secondary    = new double[additionalEquations.size()][];
  725.                 for (int i = 0; i < additionalEquations.size(); ++i) {
  726.                     secondary[i] = newState.getAdditionalState(additionalEquations.get(i).getName());
  727.                 }

  728.                 return new ODEState(newState.getDate().durationFrom(getStartDate()),
  729.                                     primary, secondary);

  730.             } catch (OrekitException oe) {
  731.                 throw new OrekitExceptionWrapper(oe);
  732.             }
  733.         }

  734.     }

  735.     /** Adapt an {@link org.orekit.propagation.sampling.OrekitStepHandler}
  736.      * to Hipparchus {@link ODEStepHandler} interface.
  737.      * @author Luc Maisonobe
  738.      */
  739.     private class AdaptedStepHandler implements ODEStepHandler, ModeHandler {

  740.         /** Underlying handler. */
  741.         private final OrekitStepHandler handler;

  742.         /** Flag for handler . */
  743.         private boolean activate;

  744.         /** Build an instance.
  745.          * @param handler underlying handler to wrap
  746.          */
  747.         AdaptedStepHandler(final OrekitStepHandler handler) {
  748.             this.handler = handler;
  749.         }

  750.         /** {@inheritDoc} */
  751.         public void initialize(final boolean activateHandlers,
  752.                                final AbsoluteDate targetDate) {
  753.             this.activate = activateHandlers;
  754.         }

  755.         /** {@inheritDoc} */
  756.         public void init(final ODEStateAndDerivative s0, final double t) {
  757.             try {
  758.                 if (activate) {
  759.                     handler.init(getCompleteState(s0.getTime(), s0.getCompleteState(), s0.getCompleteDerivative()),
  760.                                  stateMapper.mapDoubleToDate(t));
  761.                 }
  762.             } catch (OrekitException oe) {
  763.                 throw new OrekitExceptionWrapper(oe);
  764.             }
  765.         }

  766.         /** {@inheritDoc} */
  767.         public void handleStep(final ODEStateInterpolator interpolator, final boolean isLast) {
  768.             try {
  769.                 if (activate) {
  770.                     handler.handleStep(new AdaptedStepInterpolator(interpolator), isLast);
  771.                 }
  772.             } catch (OrekitException pe) {
  773.                 throw new OrekitExceptionWrapper(pe);
  774.             }
  775.         }

  776.     }

  777.     /** Adapt an Hipparchus {@link ODEStateInterpolator}
  778.      * to an orekit {@link OrekitStepInterpolator} interface.
  779.      * @author Luc Maisonobe
  780.      */
  781.     private class AdaptedStepInterpolator implements OrekitStepInterpolator {

  782.         /** Underlying raw rawInterpolator. */
  783.         private final ODEStateInterpolator mathInterpolator;

  784.         /** Simple constructor.
  785.          * @param mathInterpolator underlying raw interpolator
  786.          */
  787.         AdaptedStepInterpolator(final ODEStateInterpolator mathInterpolator) {
  788.             this.mathInterpolator = mathInterpolator;
  789.         }

  790.         /** {@inheritDoc}} */
  791.         @Override
  792.         public SpacecraftState getPreviousState()
  793.             throws OrekitException {
  794.             return convert(mathInterpolator.getPreviousState());
  795.         }

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

  801.         /** {@inheritDoc}} */
  802.         @Override
  803.         public SpacecraftState getCurrentState()
  804.             throws OrekitException {
  805.             return convert(mathInterpolator.getCurrentState());
  806.         }

  807.         /** {@inheritDoc}} */
  808.         @Override
  809.         public boolean isCurrentStateInterpolated() {
  810.             return mathInterpolator.isCurrentStateInterpolated();
  811.         }

  812.         /** {@inheritDoc}} */
  813.         @Override
  814.         public SpacecraftState getInterpolatedState(final AbsoluteDate date)
  815.             throws OrekitException {
  816.             return convert(mathInterpolator.getInterpolatedState(date.durationFrom(stateMapper.getReferenceDate())));
  817.         }

  818.         /** Convert a state from mathematical world to space flight dynamics world.
  819.          * @param os mathematical state
  820.          * @return space flight dynamics state
  821.          * @exception OrekitException if arrays are inconsistent
  822.          */
  823.         private SpacecraftState convert(final ODEStateAndDerivative os)
  824.             throws OrekitException {

  825.             SpacecraftState s =
  826.                             stateMapper.mapArrayToState(os.getTime(),
  827.                                                         os.getPrimaryState(),
  828.                                                         os.getPrimaryDerivative(),
  829.                                                         meanOrbit);
  830.             s = updateAdditionalStates(s);
  831.             for (int i = 0; i < additionalEquations.size(); ++i) {
  832.                 final double[] secondary = os.getSecondaryState(i + 1);
  833.                 s = s.addAdditionalState(additionalEquations.get(i).getName(), secondary);
  834.             }

  835.             return s;

  836.         }

  837.         /** Convert a state from space flight dynamics world to mathematical world.
  838.          * @param state space flight dynamics state
  839.          * @return mathematical state
  840.          * @exception OrekitException if arrays are inconsistent
  841.          */
  842.         private ODEStateAndDerivative convert(final SpacecraftState state)
  843.             throws OrekitException {

  844.             // retrieve initial state
  845.             final double[] primary    = new double[getBasicDimension()];
  846.             final double[] primaryDot = new double[getBasicDimension()];
  847.             stateMapper.mapStateToArray(state, primary, primaryDot);

  848.             // secondary part of the ODE
  849.             final double[][] secondary    = new double[additionalEquations.size()][];
  850.             for (int i = 0; i < additionalEquations.size(); ++i) {
  851.                 final AdditionalEquations additional = additionalEquations.get(i);
  852.                 secondary[i] = state.getAdditionalState(additional.getName());
  853.             }

  854.             return new ODEStateAndDerivative(stateMapper.mapDateToDouble(state.getDate()),
  855.                                              primary, primaryDot,
  856.                                              secondary, null);

  857.         }

  858.         /** {@inheritDoc}} */
  859.         @Override
  860.         public boolean isForward() {
  861.             return mathInterpolator.isForward();
  862.         }

  863.         /** {@inheritDoc}} */
  864.         @Override
  865.         public AdaptedStepInterpolator restrictStep(final SpacecraftState newPreviousState,
  866.                                                     final SpacecraftState newCurrentState)
  867.             throws OrekitException {
  868.             try {
  869.                 final AbstractODEStateInterpolator aosi = (AbstractODEStateInterpolator) mathInterpolator;
  870.                 return new AdaptedStepInterpolator(aosi.restrictStep(convert(newPreviousState),
  871.                                                                      convert(newCurrentState)));
  872.             } catch (ClassCastException cce) {
  873.                 // this should never happen
  874.                 throw new OrekitInternalError(cce);
  875.             }
  876.         }

  877.     }

  878.     private class EphemerisModeHandler implements ModeHandler, ODEStepHandler {

  879.         /** Underlying raw mathematical model. */
  880.         private DenseOutputModel model;

  881.         /** Generated ephemeris. */
  882.         private BoundedPropagator ephemeris;

  883.         /** Flag for handler . */
  884.         private boolean activate;

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

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

  889.         /** Creates a new instance of EphemerisModeHandler which must be
  890.          *  filled by the propagator.
  891.          */
  892.         EphemerisModeHandler() {
  893.             this.handler = null;
  894.         }

  895.         /** Creates a new instance of EphemerisModeHandler which must be
  896.          *  filled by the propagator.
  897.          *  @param handler the handler to notify of every integrator step.
  898.          */
  899.         EphemerisModeHandler(final OrekitStepHandler handler) {
  900.             this.handler = new AdaptedStepHandler(handler);
  901.         }

  902.         /** {@inheritDoc} */
  903.         public void initialize(final boolean activateHandlers,
  904.                                final AbsoluteDate targetDate) {
  905.             this.activate = activateHandlers;
  906.             this.model    = new DenseOutputModel();
  907.             this.endDate  = targetDate;

  908.             // ephemeris will be generated when last step is processed
  909.             this.ephemeris = null;
  910.             if (this.handler != null) {
  911.                 this.handler.initialize(activateHandlers, targetDate);
  912.             }
  913.         }

  914.         /** Get the generated ephemeris.
  915.          * @return a new instance of the generated ephemeris
  916.          */
  917.         public BoundedPropagator getEphemeris() {
  918.             return ephemeris;
  919.         }

  920.         /** {@inheritDoc} */
  921.         public void handleStep(final ODEStateInterpolator interpolator, final boolean isLast)
  922.             throws OrekitExceptionWrapper {
  923.             try {
  924.                 if (activate) {
  925.                     if (this.handler != null) {
  926.                         this.handler.handleStep(interpolator, isLast);
  927.                     }

  928.                     model.handleStep(interpolator, isLast);
  929.                     if (isLast) {

  930.                         // set up the boundary dates
  931.                         final double tI = model.getInitialTime();
  932.                         final double tF = model.getFinalTime();
  933.                         // tI is almost? always zero
  934.                         final AbsoluteDate startDate =
  935.                                 stateMapper.mapDoubleToDate(tI);
  936.                         final AbsoluteDate finalDate =
  937.                                 stateMapper.mapDoubleToDate(tF, this.endDate);
  938.                         final AbsoluteDate minDate;
  939.                         final AbsoluteDate maxDate;
  940.                         if (tF < tI) {
  941.                             minDate = finalDate;
  942.                             maxDate = startDate;
  943.                         } else {
  944.                             minDate = startDate;
  945.                             maxDate = finalDate;
  946.                         }

  947.                         // get the initial additional states that are not managed
  948.                         final Map<String, double[]> unmanaged = new HashMap<String, double[]>();
  949.                         for (final Map.Entry<String, double[]> initial : getInitialState().getAdditionalStates().entrySet()) {
  950.                             if (!isAdditionalStateManaged(initial.getKey())) {
  951.                                 // this additional state was in the initial state, but is unknown to the propagator
  952.                                 // we simply copy its initial value as is
  953.                                 unmanaged.put(initial.getKey(), initial.getValue());
  954.                             }
  955.                         }

  956.                         // get the names of additional states managed by differential equations
  957.                         final String[] names = new String[additionalEquations.size()];
  958.                         for (int i = 0; i < names.length; ++i) {
  959.                             names[i] = additionalEquations.get(i).getName();
  960.                         }

  961.                         // create the ephemeris
  962.                         ephemeris = new IntegratedEphemeris(startDate, minDate, maxDate,
  963.                                                             stateMapper, meanOrbit, model, unmanaged,
  964.                                                             getAdditionalStateProviders(), names);

  965.                     }
  966.                 }
  967.             } catch (OrekitException oe) {
  968.                 throw new OrekitExceptionWrapper(oe);
  969.             }
  970.         }

  971.         /** {@inheritDoc} */
  972.         public void init(final ODEStateAndDerivative s0, final double t) {
  973.             model.init(s0, t);
  974.             if (this.handler != null) {
  975.                 this.handler.init(s0, t);
  976.             }
  977.         }

  978.     }

  979. }