DSSTPropagator.java

  1. /* Copyright 2002-2023 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.semianalytical.dsst;

  18. import java.util.ArrayList;
  19. import java.util.Arrays;
  20. import java.util.Collection;
  21. import java.util.Collections;
  22. import java.util.HashSet;
  23. import java.util.List;
  24. import java.util.Set;

  25. import org.hipparchus.linear.RealMatrix;
  26. import org.hipparchus.ode.ODEIntegrator;
  27. import org.hipparchus.ode.ODEStateAndDerivative;
  28. import org.hipparchus.ode.sampling.ODEStateInterpolator;
  29. import org.hipparchus.ode.sampling.ODEStepHandler;
  30. import org.hipparchus.util.FastMath;
  31. import org.hipparchus.util.MathUtils;
  32. import org.orekit.annotation.DefaultDataContext;
  33. import org.orekit.attitudes.Attitude;
  34. import org.orekit.attitudes.AttitudeProvider;
  35. import org.orekit.data.DataContext;
  36. import org.orekit.errors.OrekitException;
  37. import org.orekit.errors.OrekitInternalError;
  38. import org.orekit.errors.OrekitMessages;
  39. import org.orekit.frames.Frame;
  40. import org.orekit.orbits.EquinoctialOrbit;
  41. import org.orekit.orbits.Orbit;
  42. import org.orekit.orbits.OrbitType;
  43. import org.orekit.orbits.PositionAngleType;
  44. import org.orekit.propagation.AbstractMatricesHarvester;
  45. import org.orekit.propagation.MatricesHarvester;
  46. import org.orekit.propagation.PropagationType;
  47. import org.orekit.propagation.Propagator;
  48. import org.orekit.propagation.SpacecraftState;
  49. import org.orekit.propagation.integration.AbstractIntegratedPropagator;
  50. import org.orekit.propagation.integration.AdditionalDerivativesProvider;
  51. import org.orekit.propagation.integration.StateMapper;
  52. import org.orekit.propagation.numerical.NumericalPropagator;
  53. import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
  54. import org.orekit.propagation.semianalytical.dsst.forces.DSSTNewtonianAttraction;
  55. import org.orekit.propagation.semianalytical.dsst.forces.ShortPeriodTerms;
  56. import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
  57. import org.orekit.propagation.semianalytical.dsst.utilities.FixedNumberInterpolationGrid;
  58. import org.orekit.propagation.semianalytical.dsst.utilities.InterpolationGrid;
  59. import org.orekit.propagation.semianalytical.dsst.utilities.MaxGapInterpolationGrid;
  60. import org.orekit.time.AbsoluteDate;
  61. import org.orekit.utils.DoubleArrayDictionary;
  62. import org.orekit.utils.ParameterDriver;
  63. import org.orekit.utils.ParameterDriversList;
  64. import org.orekit.utils.ParameterDriversList.DelegatingDriver;
  65. import org.orekit.utils.ParameterObserver;
  66. import org.orekit.utils.TimeSpanMap;
  67. import org.orekit.utils.TimeSpanMap.Span;

  68. /**
  69.  * This class propagates {@link org.orekit.orbits.Orbit orbits} using the DSST theory.
  70.  * <p>
  71.  * Whereas analytical propagators are configured only thanks to their various
  72.  * constructors and can be used immediately after construction, such a semianalytical
  73.  * propagator configuration involves setting several parameters between construction
  74.  * time and propagation time, just as numerical propagators.
  75.  * </p>
  76.  * <p>
  77.  * The configuration parameters that can be set are:
  78.  * </p>
  79.  * <ul>
  80.  * <li>the initial spacecraft state ({@link #setInitialState(SpacecraftState)})</li>
  81.  * <li>the various force models ({@link #addForceModel(DSSTForceModel)},
  82.  * {@link #removeForceModels()})</li>
  83.  * <li>the discrete events that should be triggered during propagation (
  84.  * {@link #addEventDetector(org.orekit.propagation.events.EventDetector)},
  85.  * {@link #clearEventsDetectors()})</li>
  86.  * <li>the binding logic with the rest of the application ({@link #getMultiplexer()})</li>
  87.  * </ul>
  88.  * <p>
  89.  * From these configuration parameters, only the initial state is mandatory.
  90.  * The default propagation settings are in {@link OrbitType#EQUINOCTIAL equinoctial}
  91.  * parameters with {@link PositionAngleType#TRUE true} longitude argument.
  92.  * The central attraction coefficient used to define the initial orbit will be used.
  93.  * However, specifying only the initial state would mean the propagator would use
  94.  * only Keplerian forces. In this case, the simpler
  95.  * {@link org.orekit.propagation.analytical.KeplerianPropagator KeplerianPropagator}
  96.  * class would be more effective.
  97.  * </p>
  98.  * <p>
  99.  * The underlying numerical integrator set up in the constructor may also have
  100.  * its own configuration parameters. Typical configuration parameters for adaptive
  101.  * stepsize integrators are the min, max and perhaps start step size as well as
  102.  * the absolute and/or relative errors thresholds.
  103.  * </p>
  104.  * <p>
  105.  * The state that is seen by the integrator is a simple six elements double array.
  106.  * These six elements are:
  107.  * <ul>
  108.  * <li>the {@link org.orekit.orbits.EquinoctialOrbit equinoctial orbit parameters}
  109.  * (a, e<sub>x</sub>, e<sub>y</sub>, h<sub>x</sub>, h<sub>y</sub>, λ<sub>m</sub>)
  110.  * in meters and radians,</li>
  111.  * </ul>
  112.  *
  113.  * <p>By default, at the end of the propagation, the propagator resets the initial state to the final state,
  114.  * thus allowing a new propagation to be started from there without recomputing the part already performed.
  115.  * This behaviour can be chenged by calling {@link #setResetAtEnd(boolean)}.
  116.  * </p>
  117.  * <p>Beware the same instance cannot be used simultaneously by different threads, the class is <em>not</em>
  118.  * thread-safe.</p>
  119.  *
  120.  * @see SpacecraftState
  121.  * @see DSSTForceModel
  122.  * @author Romain Di Costanzo
  123.  * @author Pascal Parraud
  124.  */
  125. public class DSSTPropagator extends AbstractIntegratedPropagator {

  126.     /** Retrograde factor I.
  127.      *  <p>
  128.      *  DSST model needs equinoctial orbit as internal representation.
  129.      *  Classical equinoctial elements have discontinuities when inclination
  130.      *  is close to zero. In this representation, I = +1. <br>
  131.      *  To avoid this discontinuity, another representation exists and equinoctial
  132.      *  elements can be expressed in a different way, called "retrograde" orbit.
  133.      *  This implies I = -1. <br>
  134.      *  As Orekit doesn't implement the retrograde orbit, I is always set to +1.
  135.      *  But for the sake of consistency with the theory, the retrograde factor
  136.      *  has been kept in the formulas.
  137.      *  </p>
  138.      */
  139.     private static final int I = 1;

  140.     /** Default value for epsilon. */
  141.     private static final double EPSILON_DEFAULT = 1.0e-13;

  142.     /** Default value for maxIterations. */
  143.     private static final int MAX_ITERATIONS_DEFAULT = 200;

  144.     /** Number of grid points per integration step to be used in interpolation of short periodics coefficients.*/
  145.     private static final int INTERPOLATION_POINTS_PER_STEP = 3;

  146.     /** Flag specifying whether the initial orbital state is given with osculating elements. */
  147.     private boolean initialIsOsculating;

  148.     /** Force models used to compute short periodic terms. */
  149.     private final transient List<DSSTForceModel> forceModels;

  150.     /** State mapper holding the force models. */
  151.     private MeanPlusShortPeriodicMapper mapper;

  152.     /** Generator for the interpolation grid. */
  153.     private InterpolationGrid interpolationgrid;

  154.     /** Create a new instance of DSSTPropagator.
  155.      *  <p>
  156.      *  After creation, there are no perturbing forces at all.
  157.      *  This means that if {@link #addForceModel addForceModel}
  158.      *  is not called after creation, the integrated orbit will
  159.      *  follow a Keplerian evolution only.
  160.      *  </p>
  161.      *
  162.      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
  163.      *
  164.      *  @param integrator numerical integrator to use for propagation.
  165.      *  @param propagationType type of orbit to output (mean or osculating).
  166.      * @see #DSSTPropagator(ODEIntegrator, PropagationType, AttitudeProvider)
  167.      */
  168.     @DefaultDataContext
  169.     public DSSTPropagator(final ODEIntegrator integrator, final PropagationType propagationType) {
  170.         this(integrator, propagationType,
  171.                 Propagator.getDefaultLaw(DataContext.getDefault().getFrames()));
  172.     }

  173.     /** Create a new instance of DSSTPropagator.
  174.      *  <p>
  175.      *  After creation, there are no perturbing forces at all.
  176.      *  This means that if {@link #addForceModel addForceModel}
  177.      *  is not called after creation, the integrated orbit will
  178.      *  follow a Keplerian evolution only.
  179.      *  </p>
  180.      * @param integrator numerical integrator to use for propagation.
  181.      * @param propagationType type of orbit to output (mean or osculating).
  182.      * @param attitudeProvider the attitude law.
  183.      * @since 10.1
  184.      */
  185.     public DSSTPropagator(final ODEIntegrator integrator,
  186.                           final PropagationType propagationType,
  187.                           final AttitudeProvider attitudeProvider) {
  188.         super(integrator, propagationType);
  189.         forceModels = new ArrayList<DSSTForceModel>();
  190.         initMapper();
  191.         // DSST uses only equinoctial orbits and mean longitude argument
  192.         setOrbitType(OrbitType.EQUINOCTIAL);
  193.         setPositionAngleType(PositionAngleType.MEAN);
  194.         setAttitudeProvider(attitudeProvider);
  195.         setInterpolationGridToFixedNumberOfPoints(INTERPOLATION_POINTS_PER_STEP);
  196.     }


  197.     /** Create a new instance of DSSTPropagator.
  198.      *  <p>
  199.      *  After creation, there are no perturbing forces at all.
  200.      *  This means that if {@link #addForceModel addForceModel}
  201.      *  is not called after creation, the integrated orbit will
  202.      *  follow a Keplerian evolution only. Only the mean orbits
  203.      *  will be generated.
  204.      *  </p>
  205.      *
  206.      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
  207.      *
  208.      *  @param integrator numerical integrator to use for propagation.
  209.      * @see #DSSTPropagator(ODEIntegrator, PropagationType, AttitudeProvider)
  210.      */
  211.     @DefaultDataContext
  212.     public DSSTPropagator(final ODEIntegrator integrator) {
  213.         this(integrator, PropagationType.MEAN);
  214.     }

  215.     /** Set the central attraction coefficient μ.
  216.      * <p>
  217.      * Setting the central attraction coefficient is
  218.      * equivalent to {@link #addForceModel(DSSTForceModel) add}
  219.      * a {@link DSSTNewtonianAttraction} force model.
  220.      * </p>
  221.     * @param mu central attraction coefficient (m³/s²)
  222.     * @see #addForceModel(DSSTForceModel)
  223.     * @see #getAllForceModels()
  224.     */
  225.     public void setMu(final double mu) {
  226.         addForceModel(new DSSTNewtonianAttraction(mu));
  227.     }

  228.     /** Set the central attraction coefficient μ only in upper class.
  229.      * @param mu central attraction coefficient (m³/s²)
  230.      */
  231.     private void superSetMu(final double mu) {
  232.         super.setMu(mu);
  233.     }

  234.     /** Check if Newtonian attraction force model is available.
  235.      * <p>
  236.      * Newtonian attraction is always the last force model in the list.
  237.      * </p>
  238.      * @return true if Newtonian attraction force model is available
  239.      */
  240.     private boolean hasNewtonianAttraction() {
  241.         final int last = forceModels.size() - 1;
  242.         return last >= 0 && forceModels.get(last) instanceof DSSTNewtonianAttraction;
  243.     }

  244.     /** Set the initial state with osculating orbital elements.
  245.      *  @param initialState initial state (defined with osculating elements)
  246.      */
  247.     public void setInitialState(final SpacecraftState initialState) {
  248.         setInitialState(initialState, PropagationType.OSCULATING);
  249.     }

  250.     /** Set the initial state.
  251.      *  @param initialState initial state
  252.      *  @param stateType defined if the orbital state is defined with osculating or mean elements
  253.      */
  254.     public void setInitialState(final SpacecraftState initialState,
  255.                                 final PropagationType stateType) {
  256.         switch (stateType) {
  257.             case MEAN:
  258.                 initialIsOsculating = false;
  259.                 break;
  260.             case OSCULATING:
  261.                 initialIsOsculating = true;
  262.                 break;
  263.             default:
  264.                 throw new OrekitInternalError(null);
  265.         }
  266.         resetInitialState(initialState);
  267.     }

  268.     /** Reset the initial state.
  269.      *
  270.      *  @param state new initial state
  271.      */
  272.     @Override
  273.     public void resetInitialState(final SpacecraftState state) {
  274.         super.resetInitialState(state);
  275.         if (!hasNewtonianAttraction()) {
  276.             // use the state to define central attraction
  277.             setMu(state.getMu());
  278.         }
  279.         super.setStartDate(state.getDate());
  280.     }

  281.     /** Set the selected short periodic coefficients that must be stored as additional states.
  282.      * @param selectedCoefficients short periodic coefficients that must be stored as additional states
  283.      * (null means no coefficients are selected, empty set means all coefficients are selected)
  284.      */
  285.     public void setSelectedCoefficients(final Set<String> selectedCoefficients) {
  286.         mapper.setSelectedCoefficients(selectedCoefficients == null ?
  287.                                        null : new HashSet<String>(selectedCoefficients));
  288.     }

  289.     /** Get the selected short periodic coefficients that must be stored as additional states.
  290.      * @return short periodic coefficients that must be stored as additional states
  291.      * (null means no coefficients are selected, empty set means all coefficients are selected)
  292.      */
  293.     public Set<String> getSelectedCoefficients() {
  294.         final Set<String> set = mapper.getSelectedCoefficients();
  295.         return set == null ? null : Collections.unmodifiableSet(set);
  296.     }

  297.     /** Get the names of the parameters in the matrix returned by {@link MatricesHarvester#getParametersJacobian}.
  298.      * @return names of the parameters (i.e. columns) of the Jacobian matrix
  299.      */
  300.     protected List<String> getJacobiansColumnsNames() {
  301.         final List<String> columnsNames = new ArrayList<>();
  302.         for (final DSSTForceModel forceModel : getAllForceModels()) {
  303.             for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
  304.                 if (driver.isSelected() && !columnsNames.contains(driver.getNamesSpanMap().getFirstSpan().getData())) {
  305.                     // As driver with same name should have same NamesSpanMap we only check if the first span is present,
  306.                     // if not we add all span names to columnsNames
  307.                     for (Span<String> span = driver.getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {
  308.                         columnsNames.add(span.getData());
  309.                     }
  310.                 }
  311.             }
  312.         }
  313.         Collections.sort(columnsNames);
  314.         return columnsNames;
  315.     }

  316.     /** {@inheritDoc} */
  317.     @Override
  318.     protected AbstractMatricesHarvester createHarvester(final String stmName, final RealMatrix initialStm,
  319.                                                         final DoubleArrayDictionary initialJacobianColumns) {
  320.         return new DSSTHarvester(this, stmName, initialStm, initialJacobianColumns);
  321.     }

  322.     /** {@inheritDoc} */
  323.     @Override
  324.     protected void setUpStmAndJacobianGenerators() {

  325.         final AbstractMatricesHarvester harvester = getHarvester();
  326.         if (harvester != null) {

  327.             // set up the additional equations and additional state providers
  328.             final DSSTStateTransitionMatrixGenerator stmGenerator = setUpStmGenerator();
  329.             setUpRegularParametersJacobiansColumns(stmGenerator);

  330.             // as we are now starting the propagation, everything is configured
  331.             // we can freeze the names in the harvester
  332.             harvester.freezeColumnsNames();

  333.         }

  334.     }

  335.     /** Set up the State Transition Matrix Generator.
  336.      * @return State Transition Matrix Generator
  337.      * @since 11.1
  338.      */
  339.     private DSSTStateTransitionMatrixGenerator setUpStmGenerator() {

  340.         final AbstractMatricesHarvester harvester = getHarvester();

  341.         // add the STM generator corresponding to the current settings, and setup state accordingly
  342.         DSSTStateTransitionMatrixGenerator stmGenerator = null;
  343.         for (final AdditionalDerivativesProvider equations : getAdditionalDerivativesProviders()) {
  344.             if (equations instanceof DSSTStateTransitionMatrixGenerator &&
  345.                 equations.getName().equals(harvester.getStmName())) {
  346.                 // the STM generator has already been set up in a previous propagation
  347.                 stmGenerator = (DSSTStateTransitionMatrixGenerator) equations;
  348.                 break;
  349.             }
  350.         }
  351.         if (stmGenerator == null) {
  352.             // this is the first time we need the STM generate, create it
  353.             stmGenerator = new DSSTStateTransitionMatrixGenerator(harvester.getStmName(),
  354.                                                                   getAllForceModels(),
  355.                                                                   getAttitudeProvider());
  356.             addAdditionalDerivativesProvider(stmGenerator);
  357.         }

  358.         if (!getInitialIntegrationState().hasAdditionalState(harvester.getStmName())) {
  359.             // add the initial State Transition Matrix if it is not already there
  360.             // (perhaps due to a previous propagation)
  361.             setInitialState(stmGenerator.setInitialStateTransitionMatrix(getInitialState(),
  362.                                                                          harvester.getInitialStateTransitionMatrix()),
  363.                             getPropagationType());
  364.         }

  365.         return stmGenerator;

  366.     }

  367.     /** Set up the Jacobians columns generator for regular parameters.
  368.      * @param stmGenerator generator for the State Transition Matrix
  369.      * @since 11.1
  370.      */
  371.     private void setUpRegularParametersJacobiansColumns(final DSSTStateTransitionMatrixGenerator stmGenerator) {

  372.         // first pass: gather all parameters (excluding trigger dates), binding similar names together
  373.         final ParameterDriversList selected = new ParameterDriversList();
  374.         for (final DSSTForceModel forceModel : getAllForceModels()) {
  375.             for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
  376.                 selected.add(driver);
  377.             }
  378.         }

  379.         // second pass: now that shared parameter names are bound together,
  380.         // their selections status have been synchronized, we can filter them
  381.         selected.filter(true);

  382.         // third pass: sort parameters lexicographically
  383.         selected.sort();

  384.         // add the Jacobians column generators corresponding to parameters, and setup state accordingly
  385.         for (final DelegatingDriver driver : selected.getDrivers()) {

  386.             for (Span<String> span = driver.getNamesSpanMap().getFirstSpan(); span != null; span = span.next()) {
  387.                 DSSTIntegrableJacobianColumnGenerator generator = null;

  388.                 // check if we already have set up the providers
  389.                 for (final AdditionalDerivativesProvider provider : getAdditionalDerivativesProviders()) {
  390.                     if (provider instanceof DSSTIntegrableJacobianColumnGenerator &&
  391.                         provider.getName().equals(span.getData())) {
  392.                         // the Jacobian column generator has already been set up in a previous propagation
  393.                         generator = (DSSTIntegrableJacobianColumnGenerator) provider;
  394.                         break;
  395.                     }
  396.                 }

  397.                 if (generator == null) {
  398.                     // this is the first time we need the Jacobian column generator, create it
  399.                     generator = new DSSTIntegrableJacobianColumnGenerator(stmGenerator, span.getData());
  400.                     addAdditionalDerivativesProvider(generator);
  401.                 }

  402.                 if (!getInitialIntegrationState().hasAdditionalState(span.getData())) {
  403.                     // add the initial Jacobian column if it is not already there
  404.                     // (perhaps due to a previous propagation)
  405.                     setInitialState(getInitialState().addAdditionalState(span.getData(),
  406.                                                                          getHarvester().getInitialJacobianColumn(span.getData())),
  407.                                     getPropagationType());
  408.                 }
  409.             }

  410.         }

  411.     }

  412.     /** Check if the initial state is provided in osculating elements.
  413.      * @return true if initial state is provided in osculating elements
  414.      */
  415.     public boolean initialIsOsculating() {
  416.         return initialIsOsculating;
  417.     }

  418.     /** Set the interpolation grid generator.
  419.      * <p>
  420.      * The generator will create an interpolation grid with a fixed
  421.      * number of points for each mean element integration step.
  422.      * </p>
  423.      * <p>
  424.      * If neither {@link #setInterpolationGridToFixedNumberOfPoints(int)}
  425.      * nor {@link #setInterpolationGridToMaxTimeGap(double)} has been called,
  426.      * by default the propagator is set as to 3 interpolations points per step.
  427.      * </p>
  428.      * @param interpolationPoints number of interpolation points at
  429.      * each integration step
  430.      * @see #setInterpolationGridToMaxTimeGap(double)
  431.      * @since 7.1
  432.      */
  433.     public void setInterpolationGridToFixedNumberOfPoints(final int interpolationPoints) {
  434.         interpolationgrid = new FixedNumberInterpolationGrid(interpolationPoints);
  435.     }

  436.     /** Set the interpolation grid generator.
  437.      * <p>
  438.      * The generator will create an interpolation grid with a maximum
  439.      * time gap between interpolation points.
  440.      * </p>
  441.      * <p>
  442.      * If neither {@link #setInterpolationGridToFixedNumberOfPoints(int)}
  443.      * nor {@link #setInterpolationGridToMaxTimeGap(double)} has been called,
  444.      * by default the propagator is set as to 3 interpolations points per step.
  445.      * </p>
  446.      * @param maxGap maximum time gap between interpolation points (seconds)
  447.      * @see #setInterpolationGridToFixedNumberOfPoints(int)
  448.      * @since 7.1
  449.      */
  450.     public void setInterpolationGridToMaxTimeGap(final double maxGap) {
  451.         interpolationgrid = new MaxGapInterpolationGrid(maxGap);
  452.     }

  453.     /** Add a force model to the global perturbation model.
  454.      *  <p>
  455.      *  If this method is not called at all,
  456.      *  the integrated orbit will follow a Keplerian evolution only.
  457.      *  </p>
  458.      *  @param force perturbing {@link DSSTForceModel force} to add
  459.      *  @see #removeForceModels()
  460.      *  @see #setMu(double)
  461.      */
  462.     public void addForceModel(final DSSTForceModel force) {

  463.         if (force instanceof DSSTNewtonianAttraction) {
  464.             // we want to add the central attraction force model

  465.             // ensure we are notified of any mu change
  466.             force.getParametersDrivers().get(0).addObserver(new ParameterObserver() {
  467.                 /** {@inheritDoc} */
  468.                 @Override
  469.                 public void valueChanged(final double previousValue, final ParameterDriver driver, final AbsoluteDate date) {
  470.                     // mu PDriver should have only 1 span
  471.                     superSetMu(driver.getValue());
  472.                 }
  473.                 /** {@inheritDoc} */
  474.                 @Override
  475.                 public void valueSpanMapChanged(final TimeSpanMap<Double> previousValue, final ParameterDriver driver) {
  476.                     // mu PDriver should have only 1 span
  477.                     superSetMu(driver.getValue());
  478.                 }
  479.             });

  480.             if (hasNewtonianAttraction()) {
  481.                 // there is already a central attraction model, replace it
  482.                 forceModels.set(forceModels.size() - 1, force);
  483.             } else {
  484.                 // there are no central attraction model yet, add it at the end of the list
  485.                 forceModels.add(force);
  486.             }
  487.         } else {
  488.             // we want to add a perturbing force model
  489.             if (hasNewtonianAttraction()) {
  490.                 // insert the new force model before Newtonian attraction,
  491.                 // which should always be the last one in the list
  492.                 forceModels.add(forceModels.size() - 1, force);
  493.             } else {
  494.                 // we only have perturbing force models up to now, just append at the end of the list
  495.                 forceModels.add(force);
  496.             }
  497.         }

  498.         force.registerAttitudeProvider(getAttitudeProvider());

  499.     }

  500.     /** Remove all perturbing force models from the global perturbation model
  501.      *  (except central attraction).
  502.      *  <p>
  503.      *  Once all perturbing forces have been removed (and as long as no new force model is added),
  504.      *  the integrated orbit will follow a Keplerian evolution only.
  505.      *  </p>
  506.      *  @see #addForceModel(DSSTForceModel)
  507.      */
  508.     public void removeForceModels() {
  509.         final int last = forceModels.size() - 1;
  510.         if (hasNewtonianAttraction()) {
  511.             // preserve the Newtonian attraction model at the end
  512.             final DSSTForceModel newton = forceModels.get(last);
  513.             forceModels.clear();
  514.             forceModels.add(newton);
  515.         } else {
  516.             forceModels.clear();
  517.         }
  518.     }

  519.     /** Get all the force models, perturbing forces and Newtonian attraction included.
  520.      * @return list of perturbing force models, with Newtonian attraction being the
  521.      * last one
  522.      * @see #addForceModel(DSSTForceModel)
  523.      * @see #setMu(double)
  524.      */
  525.     public List<DSSTForceModel> getAllForceModels() {
  526.         return Collections.unmodifiableList(forceModels);
  527.     }

  528.     /** Get propagation parameter type.
  529.      * @return orbit type used for propagation
  530.      */
  531.     public OrbitType getOrbitType() {
  532.         return super.getOrbitType();
  533.     }

  534.     /** Get propagation parameter type.
  535.      * @return angle type to use for propagation
  536.      */
  537.     public PositionAngleType getPositionAngleType() {
  538.         return super.getPositionAngleType();
  539.     }

  540.     /** Conversion from mean to osculating orbit.
  541.      * <p>
  542.      * Compute osculating state <b>in a DSST sense</b>, corresponding to the
  543.      * mean SpacecraftState in input, and according to the Force models taken
  544.      * into account.
  545.      * </p><p>
  546.      * Since the osculating state is obtained by adding short-periodic variation
  547.      * of each force model, the resulting output will depend on the
  548.      * force models parameterized in input.
  549.      * </p>
  550.      * @param mean Mean state to convert
  551.      * @param forces Forces to take into account
  552.      * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
  553.      * like atmospheric drag, radiation pressure or specific user-defined models)
  554.      * @return osculating state in a DSST sense
  555.      */
  556.     public static SpacecraftState computeOsculatingState(final SpacecraftState mean,
  557.                                                          final AttitudeProvider attitudeProvider,
  558.                                                          final Collection<DSSTForceModel> forces) {

  559.         //Create the auxiliary object
  560.         final AuxiliaryElements aux = new AuxiliaryElements(mean.getOrbit(), I);

  561.         // Set the force models
  562.         final List<ShortPeriodTerms> shortPeriodTerms = new ArrayList<ShortPeriodTerms>();
  563.         for (final DSSTForceModel force : forces) {
  564.             force.registerAttitudeProvider(attitudeProvider);
  565.             shortPeriodTerms.addAll(force.initializeShortPeriodTerms(aux, PropagationType.OSCULATING, force.getParameters(mean.getDate())));
  566.             force.updateShortPeriodTerms(force.getParametersAllValues(), mean);
  567.         }

  568.         final EquinoctialOrbit osculatingOrbit = computeOsculatingOrbit(mean, shortPeriodTerms);

  569.         return new SpacecraftState(osculatingOrbit, mean.getAttitude(), mean.getMass(),
  570.                                    mean.getAdditionalStatesValues(), mean.getAdditionalStatesDerivatives());

  571.     }

  572.     /** Conversion from osculating to mean orbit.
  573.      * <p>
  574.      * Compute mean state <b>in a DSST sense</b>, corresponding to the
  575.      * osculating SpacecraftState in input, and according to the Force models
  576.      * taken into account.
  577.      * </p><p>
  578.      * Since the osculating state is obtained with the computation of
  579.      * short-periodic variation of each force model, the resulting output will
  580.      * depend on the force models parameterized in input.
  581.      * </p><p>
  582.      * The computation is done through a fixed-point iteration process.
  583.      * </p>
  584.      * @param osculating Osculating state to convert
  585.      * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
  586.      * like atmospheric drag, radiation pressure or specific user-defined models)
  587.      * @param forceModels Forces to take into account
  588.      * @return mean state in a DSST sense
  589.      */
  590.     public static SpacecraftState computeMeanState(final SpacecraftState osculating,
  591.                                                    final AttitudeProvider attitudeProvider,
  592.                                                    final Collection<DSSTForceModel> forceModels) {
  593.         return computeMeanState(osculating, attitudeProvider, forceModels, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
  594.     }

  595.     /** Conversion from osculating to mean orbit.
  596.      * <p>
  597.      * Compute mean state <b>in a DSST sense</b>, corresponding to the
  598.      * osculating SpacecraftState in input, and according to the Force models
  599.      * taken into account.
  600.      * </p><p>
  601.      * Since the osculating state is obtained with the computation of
  602.      * short-periodic variation of each force model, the resulting output will
  603.      * depend on the force models parameterized in input.
  604.      * </p><p>
  605.      * The computation is done through a fixed-point iteration process.
  606.      * </p>
  607.      * @param osculating Osculating state to convert
  608.      * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
  609.      * like atmospheric drag, radiation pressure or specific user-defined models)
  610.      * @param forceModels Forces to take into account
  611.      * @param epsilon convergence threshold for mean parameters conversion
  612.      * @param maxIterations maximum iterations for mean parameters conversion
  613.      * @return mean state in a DSST sense
  614.      * @since 10.1
  615.      */
  616.     public static SpacecraftState computeMeanState(final SpacecraftState osculating,
  617.                                                    final AttitudeProvider attitudeProvider,
  618.                                                    final Collection<DSSTForceModel> forceModels,
  619.                                                    final double epsilon,
  620.                                                    final int maxIterations) {
  621.         final Orbit meanOrbit = computeMeanOrbit(osculating, attitudeProvider, forceModels, epsilon, maxIterations);
  622.         return new SpacecraftState(meanOrbit, osculating.getAttitude(), osculating.getMass(),
  623.                                    osculating.getAdditionalStatesValues(), osculating.getAdditionalStatesDerivatives());
  624.     }

  625.      /** Override the default value of the parameter.
  626.      *  <p>
  627.      *  By default, if the initial orbit is defined as osculating,
  628.      *  it will be averaged over 2 satellite revolutions.
  629.      *  This can be changed by using this method.
  630.      *  </p>
  631.      *  @param satelliteRevolution number of satellite revolutions to use for converting osculating to mean
  632.      *                             elements
  633.      */
  634.     public void setSatelliteRevolution(final int satelliteRevolution) {
  635.         mapper.setSatelliteRevolution(satelliteRevolution);
  636.     }

  637.     /** Get the number of satellite revolutions to use for converting osculating to mean elements.
  638.      *  @return number of satellite revolutions to use for converting osculating to mean elements
  639.      */
  640.     public int getSatelliteRevolution() {
  641.         return mapper.getSatelliteRevolution();
  642.     }

  643.     /** Override the default value short periodic terms.
  644.     *  <p>
  645.     *  By default, short periodic terms are initialized before
  646.     *  the numerical integration of the mean orbital elements.
  647.     *  </p>
  648.     *  @param shortPeriodTerms short periodic terms
  649.     */
  650.     public void setShortPeriodTerms(final List<ShortPeriodTerms> shortPeriodTerms) {
  651.         mapper.setShortPeriodTerms(shortPeriodTerms);
  652.     }

  653.    /** Get the short periodic terms.
  654.     *  @return the short periodic terms
  655.     */
  656.     public List<ShortPeriodTerms> getShortPeriodTerms() {
  657.         return mapper.getShortPeriodTerms();
  658.     }

  659.     /** {@inheritDoc} */
  660.     @Override
  661.     public void setAttitudeProvider(final AttitudeProvider attitudeProvider) {
  662.         super.setAttitudeProvider(attitudeProvider);

  663.         //Register the attitude provider for each force model
  664.         for (final DSSTForceModel force : forceModels) {
  665.             force.registerAttitudeProvider(attitudeProvider);
  666.         }
  667.     }

  668.     /** Method called just before integration.
  669.      * <p>
  670.      * The default implementation does nothing, it may be specialized in subclasses.
  671.      * </p>
  672.      * @param initialState initial state
  673.      * @param tEnd target date at which state should be propagated
  674.      */
  675.     @Override
  676.     protected void beforeIntegration(final SpacecraftState initialState,
  677.                                      final AbsoluteDate tEnd) {

  678.         // check if only mean elements must be used
  679.         final PropagationType type = getPropagationType();

  680.         // compute common auxiliary elements
  681.         final AuxiliaryElements aux = new AuxiliaryElements(initialState.getOrbit(), I);

  682.         // initialize all perturbing forces
  683.         final List<ShortPeriodTerms> shortPeriodTerms = new ArrayList<ShortPeriodTerms>();
  684.         for (final DSSTForceModel force : forceModels) {
  685.             shortPeriodTerms.addAll(force.initializeShortPeriodTerms(aux, type, force.getParameters(initialState.getDate())));
  686.         }
  687.         mapper.setShortPeriodTerms(shortPeriodTerms);

  688.         // if required, insert the special short periodics step handler
  689.         if (type == PropagationType.OSCULATING) {
  690.             final ShortPeriodicsHandler spHandler = new ShortPeriodicsHandler(forceModels);
  691.             // Compute short periodic coefficients for this point
  692.             for (DSSTForceModel forceModel : forceModels) {
  693.                 forceModel.updateShortPeriodTerms(forceModel.getParametersAllValues(), initialState);
  694.             }
  695.             final Collection<ODEStepHandler> stepHandlers = new ArrayList<ODEStepHandler>();
  696.             stepHandlers.add(spHandler);
  697.             final ODEIntegrator integrator = getIntegrator();
  698.             final Collection<ODEStepHandler> existing = integrator.getStepHandlers();
  699.             stepHandlers.addAll(existing);

  700.             integrator.clearStepHandlers();

  701.             // add back the existing handlers after the short periodics one
  702.             for (final ODEStepHandler sp : stepHandlers) {
  703.                 integrator.addStepHandler(sp);
  704.             }
  705.         }
  706.     }

  707.     /** {@inheritDoc} */
  708.     @Override
  709.     protected void afterIntegration() {
  710.         // remove the special short periodics step handler if added before
  711.         if (getPropagationType() == PropagationType.OSCULATING) {
  712.             final List<ODEStepHandler> preserved = new ArrayList<ODEStepHandler>();
  713.             final ODEIntegrator integrator = getIntegrator();
  714.             for (final ODEStepHandler sp : integrator.getStepHandlers()) {
  715.                 if (!(sp instanceof ShortPeriodicsHandler)) {
  716.                     preserved.add(sp);
  717.                 }
  718.             }

  719.             // clear the list
  720.             integrator.clearStepHandlers();

  721.             // add back the step handlers that were important for the user
  722.             for (final ODEStepHandler sp : preserved) {
  723.                 integrator.addStepHandler(sp);
  724.             }
  725.         }
  726.     }

  727.     /** Compute mean state from osculating state.
  728.      * <p>
  729.      * Compute in a DSST sense the mean state corresponding to the input osculating state.
  730.      * </p><p>
  731.      * The computing is done through a fixed-point iteration process.
  732.      * </p>
  733.      * @param osculating initial osculating state
  734.      * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
  735.      * like atmospheric drag, radiation pressure or specific user-defined models)
  736.      * @param forceModels force models
  737.      * @param epsilon convergence threshold for mean parameters conversion
  738.      * @param maxIterations maximum iterations for mean parameters conversion
  739.      * @return mean state
  740.      */
  741.     private static Orbit computeMeanOrbit(final SpacecraftState osculating,
  742.                                           final AttitudeProvider attitudeProvider,
  743.                                           final Collection<DSSTForceModel> forceModels, final double epsilon, final int maxIterations) {

  744.         // rough initialization of the mean parameters
  745.         EquinoctialOrbit meanOrbit = (EquinoctialOrbit) OrbitType.EQUINOCTIAL.convertType(osculating.getOrbit());

  746.         // threshold for each parameter
  747.         final double thresholdA = epsilon * (1 + FastMath.abs(meanOrbit.getA()));
  748.         final double thresholdE = epsilon * (1 + meanOrbit.getE());
  749.         final double thresholdI = epsilon * (1 + meanOrbit.getI());
  750.         final double thresholdL = epsilon * FastMath.PI;

  751.         // ensure all Gaussian force models can rely on attitude
  752.         for (final DSSTForceModel force : forceModels) {
  753.             force.registerAttitudeProvider(attitudeProvider);
  754.         }

  755.         int i = 0;
  756.         while (i++ < maxIterations) {

  757.             final SpacecraftState meanState = new SpacecraftState(meanOrbit, osculating.getAttitude(), osculating.getMass());

  758.             //Create the auxiliary object
  759.             final AuxiliaryElements aux = new AuxiliaryElements(meanOrbit, I);

  760.             // Set the force models
  761.             final List<ShortPeriodTerms> shortPeriodTerms = new ArrayList<ShortPeriodTerms>();
  762.             for (final DSSTForceModel force : forceModels) {
  763.                 shortPeriodTerms.addAll(force.initializeShortPeriodTerms(aux, PropagationType.OSCULATING, force.getParameters(meanState.getDate())));
  764.                 force.updateShortPeriodTerms(force.getParametersAllValues(), meanState);
  765.             }

  766.             // recompute the osculating parameters from the current mean parameters
  767.             final EquinoctialOrbit rebuilt = computeOsculatingOrbit(meanState, shortPeriodTerms);

  768.             // adapted parameters residuals
  769.             final double deltaA  = osculating.getA() - rebuilt.getA();
  770.             final double deltaEx = osculating.getEquinoctialEx() - rebuilt.getEquinoctialEx();
  771.             final double deltaEy = osculating.getEquinoctialEy() - rebuilt.getEquinoctialEy();
  772.             final double deltaHx = osculating.getHx() - rebuilt.getHx();
  773.             final double deltaHy = osculating.getHy() - rebuilt.getHy();
  774.             final double deltaLv = MathUtils.normalizeAngle(osculating.getLv() - rebuilt.getLv(), 0.0);

  775.             // check convergence
  776.             if (FastMath.abs(deltaA)  < thresholdA &&
  777.                 FastMath.abs(deltaEx) < thresholdE &&
  778.                 FastMath.abs(deltaEy) < thresholdE &&
  779.                 FastMath.abs(deltaHx) < thresholdI &&
  780.                 FastMath.abs(deltaHy) < thresholdI &&
  781.                 FastMath.abs(deltaLv) < thresholdL) {
  782.                 return meanOrbit;
  783.             }

  784.             // update mean parameters
  785.             meanOrbit = new EquinoctialOrbit(meanOrbit.getA() + deltaA,
  786.                                              meanOrbit.getEquinoctialEx() + deltaEx,
  787.                                              meanOrbit.getEquinoctialEy() + deltaEy,
  788.                                              meanOrbit.getHx() + deltaHx,
  789.                                              meanOrbit.getHy() + deltaHy,
  790.                                              meanOrbit.getLv() + deltaLv,
  791.                                              PositionAngleType.TRUE, meanOrbit.getFrame(),
  792.                                              meanOrbit.getDate(), meanOrbit.getMu());
  793.         }

  794.         throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_DSST_MEAN_PARAMETERS, i);

  795.     }

  796.     /** Compute osculating state from mean state.
  797.      * <p>
  798.      * Compute and add the short periodic variation to the mean {@link SpacecraftState}.
  799.      * </p>
  800.      * @param meanState initial mean state
  801.      * @param shortPeriodTerms short period terms
  802.      * @return osculating state
  803.      */
  804.     private static EquinoctialOrbit computeOsculatingOrbit(final SpacecraftState meanState,
  805.                                                            final List<ShortPeriodTerms> shortPeriodTerms) {

  806.         final double[] mean = new double[6];
  807.         final double[] meanDot = new double[6];
  808.         OrbitType.EQUINOCTIAL.mapOrbitToArray(meanState.getOrbit(), PositionAngleType.MEAN, mean, meanDot);
  809.         final double[] y = mean.clone();
  810.         for (final ShortPeriodTerms spt : shortPeriodTerms) {
  811.             final double[] shortPeriodic = spt.value(meanState.getOrbit());
  812.             for (int i = 0; i < shortPeriodic.length; i++) {
  813.                 y[i] += shortPeriodic[i];
  814.             }
  815.         }
  816.         return (EquinoctialOrbit) OrbitType.EQUINOCTIAL.mapArrayToOrbit(y, meanDot,
  817.                                                                         PositionAngleType.MEAN, meanState.getDate(),
  818.                                                                         meanState.getMu(), meanState.getFrame());
  819.     }

  820.     /** {@inheritDoc} */
  821.     @Override
  822.     protected SpacecraftState getInitialIntegrationState() {
  823.         if (initialIsOsculating) {
  824.             // the initial state is an osculating state,
  825.             // it must be converted to mean state
  826.             return computeMeanState(getInitialState(), getAttitudeProvider(), forceModels);
  827.         } else {
  828.             // the initial state is already a mean state
  829.             return getInitialState();
  830.         }
  831.     }

  832.     /** {@inheritDoc}
  833.      * <p>
  834.      * Note that for DSST, orbit type is hardcoded to {@link OrbitType#EQUINOCTIAL}
  835.      * and position angle type is hardcoded to {@link PositionAngleType#MEAN}, so
  836.      * the corresponding parameters are ignored.
  837.      * </p>
  838.      */
  839.     @Override
  840.     protected StateMapper createMapper(final AbsoluteDate referenceDate, final double mu,
  841.                                        final OrbitType ignoredOrbitType, final PositionAngleType ignoredPositionAngleType,
  842.                                        final AttitudeProvider attitudeProvider, final Frame frame) {

  843.         // create a mapper with the common settings provided as arguments
  844.         final MeanPlusShortPeriodicMapper newMapper =
  845.                 new MeanPlusShortPeriodicMapper(referenceDate, mu, attitudeProvider, frame);

  846.         // copy the specific settings from the existing mapper
  847.         if (mapper != null) {
  848.             newMapper.setSatelliteRevolution(mapper.getSatelliteRevolution());
  849.             newMapper.setSelectedCoefficients(mapper.getSelectedCoefficients());
  850.             newMapper.setShortPeriodTerms(mapper.getShortPeriodTerms());
  851.         }

  852.         mapper = newMapper;
  853.         return mapper;

  854.     }


  855.     /** Get the short period terms value.
  856.      * @param meanState the mean state
  857.      * @return shortPeriodTerms short period terms
  858.      * @since 7.1
  859.      */
  860.     public double[] getShortPeriodTermsValue(final SpacecraftState meanState) {
  861.         final double[] sptValue = new double[6];

  862.         for (ShortPeriodTerms spt : mapper.getShortPeriodTerms()) {
  863.             final double[] shortPeriodic = spt.value(meanState.getOrbit());
  864.             for (int i = 0; i < shortPeriodic.length; i++) {
  865.                 sptValue[i] += shortPeriodic[i];
  866.             }
  867.         }
  868.         return sptValue;
  869.     }


  870.     /** Internal mapper using mean parameters plus short periodic terms. */
  871.     private static class MeanPlusShortPeriodicMapper extends StateMapper {

  872.         /** Short periodic coefficients that must be stored as additional states. */
  873.         private Set<String>                selectedCoefficients;

  874.         /** Number of satellite revolutions in the averaging interval. */
  875.         private int                        satelliteRevolution;

  876.         /** Short period terms. */
  877.         private List<ShortPeriodTerms>     shortPeriodTerms;

  878.         /** Simple constructor.
  879.          * @param referenceDate reference date
  880.          * @param mu central attraction coefficient (m³/s²)
  881.          * @param attitudeProvider attitude provider
  882.          * @param frame inertial frame
  883.          */
  884.         MeanPlusShortPeriodicMapper(final AbsoluteDate referenceDate, final double mu,
  885.                                     final AttitudeProvider attitudeProvider, final Frame frame) {

  886.             super(referenceDate, mu, OrbitType.EQUINOCTIAL, PositionAngleType.MEAN, attitudeProvider, frame);

  887.             this.selectedCoefficients = null;

  888.             // Default averaging period for conversion from osculating to mean elements
  889.             this.satelliteRevolution = 2;

  890.             this.shortPeriodTerms    = Collections.emptyList();

  891.         }

  892.         /** {@inheritDoc} */
  893.         @Override
  894.         public SpacecraftState mapArrayToState(final AbsoluteDate date,
  895.                                                final double[] y, final double[] yDot,
  896.                                                final PropagationType type) {

  897.             // add short periodic variations to mean elements to get osculating elements
  898.             // (the loop may not be performed if there are no force models and in the
  899.             //  case we want to remain in mean parameters only)
  900.             final double[] elements = y.clone();
  901.             final DoubleArrayDictionary coefficients;
  902.             if (type == PropagationType.MEAN) {
  903.                 coefficients = null;
  904.             } else {
  905.                 final Orbit meanOrbit = OrbitType.EQUINOCTIAL.mapArrayToOrbit(elements, yDot, PositionAngleType.MEAN, date, getMu(), getFrame());
  906.                 coefficients = selectedCoefficients == null ? null : new DoubleArrayDictionary();
  907.                 for (final ShortPeriodTerms spt : shortPeriodTerms) {
  908.                     final double[] shortPeriodic = spt.value(meanOrbit);
  909.                     for (int i = 0; i < shortPeriodic.length; i++) {
  910.                         elements[i] += shortPeriodic[i];
  911.                     }
  912.                     if (selectedCoefficients != null) {
  913.                         coefficients.putAll(spt.getCoefficients(date, selectedCoefficients));
  914.                     }
  915.                 }
  916.             }

  917.             final double mass = elements[6];
  918.             if (mass <= 0.0) {
  919.                 throw new OrekitException(OrekitMessages.NOT_POSITIVE_SPACECRAFT_MASS, mass);
  920.             }

  921.             final Orbit orbit       = OrbitType.EQUINOCTIAL.mapArrayToOrbit(elements, yDot, PositionAngleType.MEAN, date, getMu(), getFrame());
  922.             final Attitude attitude = getAttitudeProvider().getAttitude(orbit, date, getFrame());

  923.             if (coefficients == null) {
  924.                 return new SpacecraftState(orbit, attitude, mass);
  925.             } else {
  926.                 return new SpacecraftState(orbit, attitude, mass, coefficients);
  927.             }

  928.         }

  929.         /** {@inheritDoc} */
  930.         @Override
  931.         public void mapStateToArray(final SpacecraftState state, final double[] y, final double[] yDot) {

  932.             OrbitType.EQUINOCTIAL.mapOrbitToArray(state.getOrbit(), PositionAngleType.MEAN, y, yDot);
  933.             y[6] = state.getMass();

  934.         }

  935.         /** Set the number of satellite revolutions to use for converting osculating to mean elements.
  936.          *  <p>
  937.          *  By default, if the initial orbit is defined as osculating,
  938.          *  it will be averaged over 2 satellite revolutions.
  939.          *  This can be changed by using this method.
  940.          *  </p>
  941.          *  @param satelliteRevolution number of satellite revolutions to use for converting osculating to mean
  942.          *                             elements
  943.          */
  944.         public void setSatelliteRevolution(final int satelliteRevolution) {
  945.             this.satelliteRevolution = satelliteRevolution;
  946.         }

  947.         /** Get the number of satellite revolutions to use for converting osculating to mean elements.
  948.          *  @return number of satellite revolutions to use for converting osculating to mean elements
  949.          */
  950.         public int getSatelliteRevolution() {
  951.             return satelliteRevolution;
  952.         }

  953.         /** Set the selected short periodic coefficients that must be stored as additional states.
  954.          * @param selectedCoefficients short periodic coefficients that must be stored as additional states
  955.          * (null means no coefficients are selected, empty set means all coefficients are selected)
  956.          */
  957.         public void setSelectedCoefficients(final Set<String> selectedCoefficients) {
  958.             this.selectedCoefficients = selectedCoefficients;
  959.         }

  960.         /** Get the selected short periodic coefficients that must be stored as additional states.
  961.          * @return short periodic coefficients that must be stored as additional states
  962.          * (null means no coefficients are selected, empty set means all coefficients are selected)
  963.          */
  964.         public Set<String> getSelectedCoefficients() {
  965.             return selectedCoefficients;
  966.         }

  967.         /** Set the short period terms.
  968.          * @param shortPeriodTerms short period terms
  969.          * @since 7.1
  970.          */
  971.         public void setShortPeriodTerms(final List<ShortPeriodTerms> shortPeriodTerms) {
  972.             this.shortPeriodTerms = shortPeriodTerms;
  973.         }

  974.         /** Get the short period terms.
  975.          * @return shortPeriodTerms short period terms
  976.          * @since 7.1
  977.          */
  978.         public List<ShortPeriodTerms> getShortPeriodTerms() {
  979.             return shortPeriodTerms;
  980.         }

  981.     }

  982.     /** {@inheritDoc} */
  983.     @Override
  984.     protected MainStateEquations getMainStateEquations(final ODEIntegrator integrator) {
  985.         return new Main(integrator);
  986.     }

  987.     /** Internal class for mean parameters integration. */
  988.     private class Main implements MainStateEquations {

  989.         /** Derivatives array. */
  990.         private final double[] yDot;

  991.         /** Simple constructor.
  992.          * @param integrator numerical integrator to use for propagation.
  993.          */
  994.         Main(final ODEIntegrator integrator) {
  995.             yDot = new double[7];

  996.             // Setup event detectors for each force model
  997.             forceModels.forEach(dsstForceModel -> dsstForceModel.getEventDetectors().
  998.                                 forEach(eventDetector -> setUpEventDetector(integrator, eventDetector)));
  999.         }

  1000.         /** {@inheritDoc} */
  1001.         @Override
  1002.         public void init(final SpacecraftState initialState, final AbsoluteDate target) {
  1003.             forceModels.forEach(fm -> fm.init(initialState, target));
  1004.         }

  1005.         /** {@inheritDoc} */
  1006.         @Override
  1007.         public double[] computeDerivatives(final SpacecraftState state) {

  1008.             Arrays.fill(yDot, 0.0);

  1009.             // compute common auxiliary elements
  1010.             final AuxiliaryElements auxiliaryElements = new AuxiliaryElements(state.getOrbit(), I);

  1011.             // compute the contributions of all perturbing forces
  1012.             for (final DSSTForceModel forceModel : forceModels) {
  1013.                 final double[] daidt = elementRates(forceModel, state, auxiliaryElements, forceModel.getParameters(state.getDate()));
  1014.                 for (int i = 0; i < daidt.length; i++) {
  1015.                     yDot[i] += daidt[i];
  1016.                 }
  1017.             }

  1018.             return yDot.clone();
  1019.         }

  1020.         /** This method allows to compute the mean equinoctial elements rates da<sub>i</sub> / dt
  1021.          *  for a specific force model.
  1022.          *  @param forceModel force to take into account
  1023.          *  @param state current state
  1024.          *  @param auxiliaryElements auxiliary elements related to the current orbit
  1025.          *  @param parameters force model parameters at state date (only 1 value for
  1026.          *  each parameter
  1027.          *  @return the mean equinoctial elements rates da<sub>i</sub> / dt
  1028.          */
  1029.         private double[] elementRates(final DSSTForceModel forceModel,
  1030.                                       final SpacecraftState state,
  1031.                                       final AuxiliaryElements auxiliaryElements,
  1032.                                       final double[] parameters) {
  1033.             return forceModel.getMeanElementRate(state, auxiliaryElements, parameters);
  1034.         }

  1035.     }

  1036.     /** Estimate tolerance vectors for an AdaptativeStepsizeIntegrator.
  1037.      *  <p>
  1038.      *  The errors are estimated from partial derivatives properties of orbits,
  1039.      *  starting from a scalar position error specified by the user.
  1040.      *  Considering the energy conservation equation V = sqrt(mu (2/r - 1/a)),
  1041.      *  we get at constant energy (i.e. on a Keplerian trajectory):
  1042.      *
  1043.      *  <pre>
  1044.      *  V r² |dV| = mu |dr|
  1045.      *  </pre>
  1046.      *
  1047.      *  <p> So we deduce a scalar velocity error consistent with the position error. From here, we apply
  1048.      *  orbits Jacobians matrices to get consistent errors on orbital parameters.
  1049.      *
  1050.      *  <p>
  1051.      *  The tolerances are only <em>orders of magnitude</em>, and integrator tolerances are only
  1052.      *  local estimates, not global ones. So some care must be taken when using these tolerances.
  1053.      *  Setting 1mm as a position error does NOT mean the tolerances will guarantee a 1mm error
  1054.      *  position after several orbits integration.
  1055.      *  </p>
  1056.      *
  1057.      * @param dP user specified position error (m)
  1058.      * @param orbit reference orbit
  1059.      * @return a two rows array, row 0 being the absolute tolerance error
  1060.      *                       and row 1 being the relative tolerance error
  1061.      */
  1062.     public static double[][] tolerances(final double dP, final Orbit orbit) {

  1063.         return NumericalPropagator.tolerances(dP, orbit, OrbitType.EQUINOCTIAL);

  1064.     }

  1065.     /** Estimate tolerance vectors for an AdaptativeStepsizeIntegrator.
  1066.      *  <p>
  1067.      *  The errors are estimated from partial derivatives properties of orbits,
  1068.      *  starting from scalar position and velocity errors specified by the user.
  1069.      *  <p>
  1070.      *  The tolerances are only <em>orders of magnitude</em>, and integrator tolerances are only
  1071.      *  local estimates, not global ones. So some care must be taken when using these tolerances.
  1072.      *  Setting 1mm as a position error does NOT mean the tolerances will guarantee a 1mm error
  1073.      *  position after several orbits integration.
  1074.      *  </p>
  1075.      *
  1076.      * @param dP user specified position error (m)
  1077.      * @param dV user specified velocity error (m/s)
  1078.      * @param orbit reference orbit
  1079.      * @return a two rows array, row 0 being the absolute tolerance error
  1080.      *                       and row 1 being the relative tolerance error
  1081.      * @since 10.3
  1082.      */
  1083.     public static double[][] tolerances(final double dP, final double dV, final Orbit orbit) {

  1084.         return NumericalPropagator.tolerances(dP, dV, orbit, OrbitType.EQUINOCTIAL);

  1085.     }

  1086.     /** Step handler used to compute the parameters for the short periodic contributions.
  1087.      * @author Lucian Barbulescu
  1088.      */
  1089.     private class ShortPeriodicsHandler implements ODEStepHandler {

  1090.         /** Force models used to compute short periodic terms. */
  1091.         private final List<DSSTForceModel> forceModels;

  1092.         /** Constructor.
  1093.          * @param forceModels force models
  1094.          */
  1095.         ShortPeriodicsHandler(final List<DSSTForceModel> forceModels) {
  1096.             this.forceModels = forceModels;
  1097.         }

  1098.         /** {@inheritDoc} */
  1099.         @Override
  1100.         public void handleStep(final ODEStateInterpolator interpolator) {

  1101.             // Get the grid points to compute
  1102.             final double[] interpolationPoints =
  1103.                     interpolationgrid.getGridPoints(interpolator.getPreviousState().getTime(),
  1104.                                                     interpolator.getCurrentState().getTime());

  1105.             final SpacecraftState[] meanStates = new SpacecraftState[interpolationPoints.length];
  1106.             for (int i = 0; i < interpolationPoints.length; ++i) {

  1107.                 // Build the mean state interpolated at grid point
  1108.                 final double time = interpolationPoints[i];
  1109.                 final ODEStateAndDerivative sd = interpolator.getInterpolatedState(time);
  1110.                 meanStates[i] = mapper.mapArrayToState(time,
  1111.                                                        sd.getPrimaryState(),
  1112.                                                        sd.getPrimaryDerivative(),
  1113.                                                        PropagationType.MEAN);
  1114.             }

  1115.             // Computate short periodic coefficients for this step
  1116.             for (DSSTForceModel forceModel : forceModels) {
  1117.                 forceModel.updateShortPeriodTerms(forceModel.getParametersAllValues(), meanStates);
  1118.             }
  1119.         }
  1120.     }
  1121. }