/* Copyright 2002-2022 CS GROUP
    * Licensed to CS GROUP (CS) under one or more
    * contributor license agreements.  See the NOTICE file distributed with
    * this work for additional information regarding copyright ownership.
    * CS licenses this file to You under the Apache License, Version 2.0
    * (the "License"); you may not use this file except in compliance with
    * the License.  You may obtain a copy of the License at
    *
    *   http://www.apache.org/licenses/LICENSE-2.0
   *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   */
  package org.orekit.propagation.semianalytical.dsst;
  
  import java.util.ArrayList;
  import java.util.Arrays;
  import java.util.Collection;
  import java.util.Collections;
  import java.util.HashSet;
  import java.util.List;
  import java.util.Set;
  
  import org.hipparchus.linear.RealMatrix;
  import org.hipparchus.ode.ODEIntegrator;
  import org.hipparchus.ode.ODEStateAndDerivative;
  import org.hipparchus.ode.sampling.ODEStateInterpolator;
  import org.hipparchus.ode.sampling.ODEStepHandler;
  import org.hipparchus.util.FastMath;
  import org.hipparchus.util.MathUtils;
  import org.orekit.annotation.DefaultDataContext;
  import org.orekit.attitudes.Attitude;
  import org.orekit.attitudes.AttitudeProvider;
  import org.orekit.data.DataContext;
  import org.orekit.errors.OrekitException;
  import org.orekit.errors.OrekitInternalError;
  import org.orekit.errors.OrekitMessages;
  import org.orekit.frames.Frame;
  import org.orekit.orbits.EquinoctialOrbit;
  import org.orekit.orbits.Orbit;
  import org.orekit.orbits.OrbitType;
  import org.orekit.orbits.PositionAngle;
  import org.orekit.propagation.AbstractMatricesHarvester;
  import org.orekit.propagation.MatricesHarvester;
  import org.orekit.propagation.PropagationType;
  import org.orekit.propagation.Propagator;
  import org.orekit.propagation.SpacecraftState;
  import org.orekit.propagation.events.EventDetector;
  import org.orekit.propagation.integration.AbstractIntegratedPropagator;
  import org.orekit.propagation.integration.AdditionalDerivativesProvider;
  import org.orekit.propagation.integration.StateMapper;
  import org.orekit.propagation.numerical.NumericalPropagator;
  import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
  import org.orekit.propagation.semianalytical.dsst.forces.DSSTNewtonianAttraction;
  import org.orekit.propagation.semianalytical.dsst.forces.ShortPeriodTerms;
  import org.orekit.propagation.semianalytical.dsst.utilities.AuxiliaryElements;
  import org.orekit.propagation.semianalytical.dsst.utilities.FixedNumberInterpolationGrid;
  import org.orekit.propagation.semianalytical.dsst.utilities.InterpolationGrid;
  import org.orekit.propagation.semianalytical.dsst.utilities.MaxGapInterpolationGrid;
  import org.orekit.time.AbsoluteDate;
  import org.orekit.utils.DoubleArrayDictionary;
  import org.orekit.utils.ParameterDriver;
  import org.orekit.utils.ParameterDriversList;
  import org.orekit.utils.ParameterDriversList.DelegatingDriver;
  import org.orekit.utils.ParameterObserver;
  
  /**
   * This class propagates {@link org.orekit.orbits.Orbit orbits} using the DSST theory.
   * <p>
   * Whereas analytical propagators are configured only thanks to their various
   * constructors and can be used immediately after construction, such a semianalytical
   * propagator configuration involves setting several parameters between construction
   * time and propagation time, just as numerical propagators.
   * </p>
   * <p>
   * The configuration parameters that can be set are:
   * </p>
   * <ul>
   * <li>the initial spacecraft state ({@link #setInitialState(SpacecraftState)})</li>
   * <li>the various force models ({@link #addForceModel(DSSTForceModel)},
   * {@link #removeForceModels()})</li>
   * <li>the discrete events that should be triggered during propagation (
   * {@link #addEventDetector(org.orekit.propagation.events.EventDetector)},
   * {@link #clearEventsDetectors()})</li>
   * <li>the binding logic with the rest of the application ({@link #getMultiplexer()})</li>
   * </ul>
   * <p>
   * From these configuration parameters, only the initial state is mandatory.
   * The default propagation settings are in {@link OrbitType#EQUINOCTIAL equinoctial}
   * parameters with {@link PositionAngle#TRUE true} longitude argument.
   * The central attraction coefficient used to define the initial orbit will be used.
   * However, specifying only the initial state would mean the propagator would use
   * only Keplerian forces. In this case, the simpler
   * {@link org.orekit.propagation.analytical.KeplerianPropagator KeplerianPropagator}
   * class would be more effective.
   * </p>
  * <p>
  * The underlying numerical integrator set up in the constructor may also have
  * its own configuration parameters. Typical configuration parameters for adaptive
  * stepsize integrators are the min, max and perhaps start step size as well as
  * the absolute and/or relative errors thresholds.
  * </p>
  * <p>
  * The state that is seen by the integrator is a simple six elements double array.
  * These six elements are:
  * <ul>
  * <li>the {@link org.orekit.orbits.EquinoctialOrbit equinoctial orbit parameters}
  * (a, e<sub>x</sub>, e<sub>y</sub>, h<sub>x</sub>, h<sub>y</sub>, λ<sub>m</sub>)
  * in meters and radians,</li>
  * </ul>
  *
  * <p>By default, at the end of the propagation, the propagator resets the initial state to the final state,
  * thus allowing a new propagation to be started from there without recomputing the part already performed.
  * This behaviour can be chenged by calling {@link #setResetAtEnd(boolean)}.
  * </p>
  * <p>Beware the same instance cannot be used simultaneously by different threads, the class is <em>not</em>
  * thread-safe.</p>
  *
  * @see SpacecraftState
  * @see DSSTForceModel
  * @author Romain Di Costanzo
  * @author Pascal Parraud
  */
 public class DSSTPropagator extends AbstractIntegratedPropagator {
 
     /** Retrograde factor I.
      *  <p>
      *  DSST model needs equinoctial orbit as internal representation.
      *  Classical equinoctial elements have discontinuities when inclination
      *  is close to zero. In this representation, I = +1. <br>
      *  To avoid this discontinuity, another representation exists and equinoctial
      *  elements can be expressed in a different way, called "retrograde" orbit.
      *  This implies I = -1. <br>
      *  As Orekit doesn't implement the retrograde orbit, I is always set to +1.
      *  But for the sake of consistency with the theory, the retrograde factor
      *  has been kept in the formulas.
      *  </p>
      */
     private static final int I = 1;
 
     /** Default value for epsilon. */
     private static final double EPSILON_DEFAULT = 1.0e-13;
 
     /** Default value for maxIterations. */
     private static final int MAX_ITERATIONS_DEFAULT = 200;
 
     /** Number of grid points per integration step to be used in interpolation of short periodics coefficients.*/
     private static final int INTERPOLATION_POINTS_PER_STEP = 3;
 
     /** Flag specifying whether the initial orbital state is given with osculating elements. */
     private boolean initialIsOsculating;
 
     /** Force models used to compute short periodic terms. */
     private final transient List<DSSTForceModel> forceModels;
 
     /** State mapper holding the force models. */
     private MeanPlusShortPeriodicMapper mapper;
 
     /** Generator for the interpolation grid. */
     private InterpolationGrid interpolationgrid;
 
     /** Create a new instance of DSSTPropagator.
      *  <p>
      *  After creation, there are no perturbing forces at all.
      *  This means that if {@link #addForceModel addForceModel}
      *  is not called after creation, the integrated orbit will
      *  follow a Keplerian evolution only.
      *  </p>
      *
      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
      *
      *  @param integrator numerical integrator to use for propagation.
      *  @param propagationType type of orbit to output (mean or osculating).
      * @see #DSSTPropagator(ODEIntegrator, PropagationType, AttitudeProvider)
      */
     @DefaultDataContext
     public DSSTPropagator(final ODEIntegrator integrator, final PropagationType propagationType) {
         this(integrator, propagationType,
                 Propagator.getDefaultLaw(DataContext.getDefault().getFrames()));
     }
 
     /** Create a new instance of DSSTPropagator.
      *  <p>
      *  After creation, there are no perturbing forces at all.
      *  This means that if {@link #addForceModel addForceModel}
      *  is not called after creation, the integrated orbit will
      *  follow a Keplerian evolution only.
      *  </p>
      * @param integrator numerical integrator to use for propagation.
      * @param propagationType type of orbit to output (mean or osculating).
      * @param attitudeProvider the attitude law.
      * @since 10.1
      */
     public DSSTPropagator(final ODEIntegrator integrator,
                           final PropagationType propagationType,
                           final AttitudeProvider attitudeProvider) {
         super(integrator, propagationType);
         forceModels = new ArrayList<DSSTForceModel>();
         initMapper();
         // DSST uses only equinoctial orbits and mean longitude argument
         setOrbitType(OrbitType.EQUINOCTIAL);
         setPositionAngleType(PositionAngle.MEAN);
         setAttitudeProvider(attitudeProvider);
         setInterpolationGridToFixedNumberOfPoints(INTERPOLATION_POINTS_PER_STEP);
     }
 
 
     /** Create a new instance of DSSTPropagator.
      *  <p>
      *  After creation, there are no perturbing forces at all.
      *  This means that if {@link #addForceModel addForceModel}
      *  is not called after creation, the integrated orbit will
      *  follow a Keplerian evolution only. Only the mean orbits
      *  will be generated.
      *  </p>
      *
      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
      *
      *  @param integrator numerical integrator to use for propagation.
      * @see #DSSTPropagator(ODEIntegrator, PropagationType, AttitudeProvider)
      */
     @DefaultDataContext
     public DSSTPropagator(final ODEIntegrator integrator) {
         this(integrator, PropagationType.MEAN);
     }
 
     /** Set the central attraction coefficient μ.
      * <p>
      * Setting the central attraction coefficient is
      * equivalent to {@link #addForceModel(DSSTForceModel) add}
      * a {@link DSSTNewtonianAttraction} force model.
      * </p>
     * @param mu central attraction coefficient (m³/s²)
     * @see #addForceModel(DSSTForceModel)
     * @see #getAllForceModels()
     */
     public void setMu(final double mu) {
         addForceModel(new DSSTNewtonianAttraction(mu));
     }
 
     /** Set the central attraction coefficient μ only in upper class.
      * @param mu central attraction coefficient (m³/s²)
      */
     private void superSetMu(final double mu) {
         super.setMu(mu);
     }
 
     /** Check if Newtonian attraction force model is available.
      * <p>
      * Newtonian attraction is always the last force model in the list.
      * </p>
      * @return true if Newtonian attraction force model is available
      */
     private boolean hasNewtonianAttraction() {
         final int last = forceModels.size() - 1;
         return last >= 0 && forceModels.get(last) instanceof DSSTNewtonianAttraction;
     }
 
     /** Set the initial state with osculating orbital elements.
      *  @param initialState initial state (defined with osculating elements)
      */
     public void setInitialState(final SpacecraftState initialState) {
         setInitialState(initialState, PropagationType.OSCULATING);
     }
 
     /** Set the initial state.
      *  @param initialState initial state
      *  @param stateType defined if the orbital state is defined with osculating or mean elements
      */
     public void setInitialState(final SpacecraftState initialState,
                                 final PropagationType stateType) {
         switch (stateType) {
             case MEAN:
                 initialIsOsculating = false;
                 break;
             case OSCULATING:
                 initialIsOsculating = true;
                 break;
             default:
                 throw new OrekitInternalError(null);
         }
         resetInitialState(initialState);
     }
 
     /** Reset the initial state.
      *
      *  @param state new initial state
      */
     @Override
     public void resetInitialState(final SpacecraftState state) {
         super.resetInitialState(state);
         if (!hasNewtonianAttraction()) {
             // use the state to define central attraction
             setMu(state.getMu());
         }
         super.setStartDate(state.getDate());
     }
 
     /** Set the selected short periodic coefficients that must be stored as additional states.
      * @param selectedCoefficients short periodic coefficients that must be stored as additional states
      * (null means no coefficients are selected, empty set means all coefficients are selected)
      */
     public void setSelectedCoefficients(final Set<String> selectedCoefficients) {
         mapper.setSelectedCoefficients(selectedCoefficients == null ?
                                        null : new HashSet<String>(selectedCoefficients));
     }
 
     /** Get the selected short periodic coefficients that must be stored as additional states.
      * @return short periodic coefficients that must be stored as additional states
      * (null means no coefficients are selected, empty set means all coefficients are selected)
      */
     public Set<String> getSelectedCoefficients() {
         final Set<String> set = mapper.getSelectedCoefficients();
         return set == null ? null : Collections.unmodifiableSet(set);
     }
 
     /** Get the names of the parameters in the matrix returned by {@link MatricesHarvester#getParametersJacobian}.
      * @return names of the parameters (i.e. columns) of the Jacobian matrix
      */
     protected List<String> getJacobiansColumnsNames() {
         final List<String> columnsNames = new ArrayList<>();
         for (final DSSTForceModel forceModel : getAllForceModels()) {
             for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
                 if (driver.isSelected() && !columnsNames.contains(driver.getName())) {
                     columnsNames.add(driver.getName());
                 }
             }
         }
         Collections.sort(columnsNames);
         return columnsNames;
     }
 
     /** {@inheritDoc} */
     @Override
     protected AbstractMatricesHarvester createHarvester(final String stmName, final RealMatrix initialStm,
                                                         final DoubleArrayDictionary initialJacobianColumns) {
         return new DSSTHarvester(this, stmName, initialStm, initialJacobianColumns);
     }
 
     /** {@inheritDoc} */
     @Override
     protected void setUpStmAndJacobianGenerators() {
 
         final AbstractMatricesHarvester harvester = getHarvester();
         if (harvester != null) {
 
             // set up the additional equations and additional state providers
             final DSSTStateTransitionMatrixGenerator stmGenerator = setUpStmGenerator();
             setUpRegularParametersJacobiansColumns(stmGenerator);
 
             // as we are now starting the propagation, everything is configured
             // we can freeze the names in the harvester
             harvester.freezeColumnsNames();
 
         }
 
     }
 
     /** Set up the State Transition Matrix Generator.
      * @return State Transition Matrix Generator
      * @since 11.1
      */
     private DSSTStateTransitionMatrixGenerator setUpStmGenerator() {
 
         final AbstractMatricesHarvester harvester = getHarvester();
 
         // add the STM generator corresponding to the current settings, and setup state accordingly
         DSSTStateTransitionMatrixGenerator stmGenerator = null;
         for (final AdditionalDerivativesProvider equations : getAdditionalDerivativesProviders()) {
             if (equations instanceof DSSTStateTransitionMatrixGenerator &&
                 equations.getName().equals(harvester.getStmName())) {
                 // the STM generator has already been set up in a previous propagation
                 stmGenerator = (DSSTStateTransitionMatrixGenerator) equations;
                 break;
             }
         }
         if (stmGenerator == null) {
             // this is the first time we need the STM generate, create it
             stmGenerator = new DSSTStateTransitionMatrixGenerator(harvester.getStmName(),
                                                                   getAllForceModels(),
                                                                   getAttitudeProvider());
             addAdditionalDerivativesProvider(stmGenerator);
         }
 
         if (!getInitialIntegrationState().hasAdditionalState(harvester.getStmName())) {
             // add the initial State Transition Matrix if it is not already there
             // (perhaps due to a previous propagation)
             setInitialState(stmGenerator.setInitialStateTransitionMatrix(getInitialState(),
                                                                          harvester.getInitialStateTransitionMatrix()),
                             getPropagationType());
         }
 
         return stmGenerator;
 
     }
 
     /** Set up the Jacobians columns generator for regular parameters.
      * @param stmGenerator generator for the State Transition Matrix
      * @since 11.1
      */
     private void setUpRegularParametersJacobiansColumns(final DSSTStateTransitionMatrixGenerator stmGenerator) {
 
         // first pass: gather all parameters (excluding trigger dates), binding similar names together
         final ParameterDriversList selected = new ParameterDriversList();
         for (final DSSTForceModel forceModel : getAllForceModels()) {
             for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
                 selected.add(driver);
             }
         }
 
         // second pass: now that shared parameter names are bound together,
         // their selections status have been synchronized, we can filter them
         selected.filter(true);
 
         // third pass: sort parameters lexicographically
         selected.sort();
 
         // add the Jacobians column generators corresponding to parameters, and setup state accordingly
         for (final DelegatingDriver driver : selected.getDrivers()) {
 
             DSSTIntegrableJacobianColumnGenerator generator = null;
 
             // check if we already have set up the providers
             for (final AdditionalDerivativesProvider provider : getAdditionalDerivativesProviders()) {
                 if (provider instanceof DSSTIntegrableJacobianColumnGenerator &&
                     provider.getName().equals(driver.getName())) {
                     // the Jacobian column generator has already been set up in a previous propagation
                     generator = (DSSTIntegrableJacobianColumnGenerator) provider;
                     break;
                 }
             }
 
             if (generator == null) {
                 // this is the first time we need the Jacobian column generator, create it
                 generator = new DSSTIntegrableJacobianColumnGenerator(stmGenerator, driver.getName());
                 addAdditionalDerivativesProvider(generator);
             }
 
             if (!getInitialIntegrationState().hasAdditionalState(driver.getName())) {
                 // add the initial Jacobian column if it is not already there
                 // (perhaps due to a previous propagation)
                 setInitialState(getInitialState().addAdditionalState(driver.getName(),
                                                                      getHarvester().getInitialJacobianColumn(driver.getName())),
                                 getPropagationType());
             }
 
         }
 
     }
 
     /** Check if the initial state is provided in osculating elements.
      * @return true if initial state is provided in osculating elements
      */
     public boolean initialIsOsculating() {
         return initialIsOsculating;
     }
 
     /** Set the interpolation grid generator.
      * <p>
      * The generator will create an interpolation grid with a fixed
      * number of points for each mean element integration step.
      * </p>
      * <p>
      * If neither {@link #setInterpolationGridToFixedNumberOfPoints(int)}
      * nor {@link #setInterpolationGridToMaxTimeGap(double)} has been called,
      * by default the propagator is set as to 3 interpolations points per step.
      * </p>
      * @param interpolationPoints number of interpolation points at
      * each integration step
      * @see #setInterpolationGridToMaxTimeGap(double)
      * @since 7.1
      */
     public void setInterpolationGridToFixedNumberOfPoints(final int interpolationPoints) {
         interpolationgrid = new FixedNumberInterpolationGrid(interpolationPoints);
     }
 
     /** Set the interpolation grid generator.
      * <p>
      * The generator will create an interpolation grid with a maximum
      * time gap between interpolation points.
      * </p>
      * <p>
      * If neither {@link #setInterpolationGridToFixedNumberOfPoints(int)}
      * nor {@link #setInterpolationGridToMaxTimeGap(double)} has been called,
      * by default the propagator is set as to 3 interpolations points per step.
      * </p>
      * @param maxGap maximum time gap between interpolation points (seconds)
      * @see #setInterpolationGridToFixedNumberOfPoints(int)
      * @since 7.1
      */
     public void setInterpolationGridToMaxTimeGap(final double maxGap) {
         interpolationgrid = new MaxGapInterpolationGrid(maxGap);
     }
 
     /** Add a force model to the global perturbation model.
      *  <p>
      *  If this method is not called at all,
      *  the integrated orbit will follow a Keplerian evolution only.
      *  </p>
      *  @param force perturbing {@link DSSTForceModel force} to add
      *  @see #removeForceModels()
      *  @see #setMu(double)
      */
     public void addForceModel(final DSSTForceModel force) {
 
         if (force instanceof DSSTNewtonianAttraction) {
             // we want to add the central attraction force model
 
             // ensure we are notified of any mu change
             force.getParametersDrivers().get(0).addObserver(new ParameterObserver() {
                 /** {@inheritDoc} */
                 @Override
                 public void valueChanged(final double previousValue, final ParameterDriver driver) {
                     superSetMu(driver.getValue());
                 }
             });
 
             if (hasNewtonianAttraction()) {
                 // there is already a central attraction model, replace it
                 forceModels.set(forceModels.size() - 1, force);
             } else {
                 // there are no central attraction model yet, add it at the end of the list
                 forceModels.add(force);
             }
         } else {
             // we want to add a perturbing force model
             if (hasNewtonianAttraction()) {
                 // insert the new force model before Newtonian attraction,
                 // which should always be the last one in the list
                 forceModels.add(forceModels.size() - 1, force);
             } else {
                 // we only have perturbing force models up to now, just append at the end of the list
                 forceModels.add(force);
             }
         }
 
         force.registerAttitudeProvider(getAttitudeProvider());
 
     }
 
     /** Remove all perturbing force models from the global perturbation model
      *  (except central attraction).
      *  <p>
      *  Once all perturbing forces have been removed (and as long as no new force model is added),
      *  the integrated orbit will follow a Keplerian evolution only.
      *  </p>
      *  @see #addForceModel(DSSTForceModel)
      */
     public void removeForceModels() {
         final int last = forceModels.size() - 1;
         if (hasNewtonianAttraction()) {
             // preserve the Newtonian attraction model at the end
             final DSSTForceModel newton = forceModels.get(last);
             forceModels.clear();
             forceModels.add(newton);
         } else {
             forceModels.clear();
         }
     }
 
     /** Get all the force models, perturbing forces and Newtonian attraction included.
      * @return list of perturbing force models, with Newtonian attraction being the
      * last one
      * @see #addForceModel(DSSTForceModel)
      * @see #setMu(double)
      */
     public List<DSSTForceModel> getAllForceModels() {
         return Collections.unmodifiableList(forceModels);
     }
 
     /** Get propagation parameter type.
      * @return orbit type used for propagation
      */
     public OrbitType getOrbitType() {
         return super.getOrbitType();
     }
 
     /** Get propagation parameter type.
      * @return angle type to use for propagation
      */
     public PositionAngle getPositionAngleType() {
         return super.getPositionAngleType();
     }
 
     /** Conversion from mean to osculating orbit.
      * <p>
      * Compute osculating state <b>in a DSST sense</b>, corresponding to the
      * mean SpacecraftState in input, and according to the Force models taken
      * into account.
      * </p><p>
      * Since the osculating state is obtained by adding short-periodic variation
      * of each force model, the resulting output will depend on the
      * force models parameterized in input.
      * </p>
      * @param mean Mean state to convert
      * @param forces Forces to take into account
      * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
      * like atmospheric drag, radiation pressure or specific user-defined models)
      * @return osculating state in a DSST sense
      */
     public static SpacecraftState computeOsculatingState(final SpacecraftState mean,
                                                          final AttitudeProvider attitudeProvider,
                                                          final Collection<DSSTForceModel> forces) {
 
         //Create the auxiliary object
         final AuxiliaryElements aux = new AuxiliaryElements(mean.getOrbit(), I);
 
         // Set the force models
         final List<ShortPeriodTerms> shortPeriodTerms = new ArrayList<ShortPeriodTerms>();
         for (final DSSTForceModel force : forces) {
             force.registerAttitudeProvider(attitudeProvider);
             shortPeriodTerms.addAll(force.initializeShortPeriodTerms(aux, PropagationType.OSCULATING, force.getParameters()));
             force.updateShortPeriodTerms(force.getParameters(), mean);
         }
 
         final EquinoctialOrbit osculatingOrbit = computeOsculatingOrbit(mean, shortPeriodTerms);
 
         return new SpacecraftState(osculatingOrbit, mean.getAttitude(), mean.getMass(),
                                    mean.getAdditionalStatesValues(), mean.getAdditionalStatesDerivatives());
 
     }
 
     /** Conversion from osculating to mean orbit.
      * <p>
      * Compute mean state <b>in a DSST sense</b>, corresponding to the
      * osculating SpacecraftState in input, and according to the Force models
      * taken into account.
      * </p><p>
      * Since the osculating state is obtained with the computation of
      * short-periodic variation of each force model, the resulting output will
      * depend on the force models parameterized in input.
      * </p><p>
      * The computation is done through a fixed-point iteration process.
      * </p>
      * @param osculating Osculating state to convert
      * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
      * like atmospheric drag, radiation pressure or specific user-defined models)
      * @param forceModels Forces to take into account
      * @return mean state in a DSST sense
      */
     public static SpacecraftState computeMeanState(final SpacecraftState osculating,
                                                    final AttitudeProvider attitudeProvider,
                                                    final Collection<DSSTForceModel> forceModels) {
         return computeMeanState(osculating, attitudeProvider, forceModels, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
     }
 
     /** Conversion from osculating to mean orbit.
      * <p>
      * Compute mean state <b>in a DSST sense</b>, corresponding to the
      * osculating SpacecraftState in input, and according to the Force models
      * taken into account.
      * </p><p>
      * Since the osculating state is obtained with the computation of
      * short-periodic variation of each force model, the resulting output will
      * depend on the force models parameterized in input.
      * </p><p>
      * The computation is done through a fixed-point iteration process.
      * </p>
      * @param osculating Osculating state to convert
      * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
      * like atmospheric drag, radiation pressure or specific user-defined models)
      * @param forceModels Forces to take into account
      * @param epsilon convergence threshold for mean parameters conversion
      * @param maxIterations maximum iterations for mean parameters conversion
      * @return mean state in a DSST sense
      * @since 10.1
      */
     public static SpacecraftState computeMeanState(final SpacecraftState osculating,
                                                    final AttitudeProvider attitudeProvider,
                                                    final Collection<DSSTForceModel> forceModels,
                                                    final double epsilon,
                                                    final int maxIterations) {
         final Orbit meanOrbit = computeMeanOrbit(osculating, attitudeProvider, forceModels, epsilon, maxIterations);
         return new SpacecraftState(meanOrbit, osculating.getAttitude(), osculating.getMass(),
                                    osculating.getAdditionalStatesValues(), osculating.getAdditionalStatesDerivatives());
     }
 
      /** Override the default value of the parameter.
      *  <p>
      *  By default, if the initial orbit is defined as osculating,
      *  it will be averaged over 2 satellite revolutions.
      *  This can be changed by using this method.
      *  </p>
      *  @param satelliteRevolution number of satellite revolutions to use for converting osculating to mean
      *                             elements
      */
     public void setSatelliteRevolution(final int satelliteRevolution) {
         mapper.setSatelliteRevolution(satelliteRevolution);
     }
 
     /** Get the number of satellite revolutions to use for converting osculating to mean elements.
      *  @return number of satellite revolutions to use for converting osculating to mean elements
      */
     public int getSatelliteRevolution() {
         return mapper.getSatelliteRevolution();
     }
 
     /** Override the default value short periodic terms.
     *  <p>
     *  By default, short periodic terms are initialized before
     *  the numerical integration of the mean orbital elements.
     *  </p>
     *  @param shortPeriodTerms short periodic terms
     */
     public void setShortPeriodTerms(final List<ShortPeriodTerms> shortPeriodTerms) {
         mapper.setShortPeriodTerms(shortPeriodTerms);
     }
 
    /** Get the short periodic terms.
     *  @return the short periodic terms
     */
     public List<ShortPeriodTerms> getShortPeriodTerms() {
         return mapper.getShortPeriodTerms();
     }
 
     /** {@inheritDoc} */
     @Override
     public void setAttitudeProvider(final AttitudeProvider attitudeProvider) {
         super.setAttitudeProvider(attitudeProvider);
 
         //Register the attitude provider for each force model
         for (final DSSTForceModel force : forceModels) {
             force.registerAttitudeProvider(attitudeProvider);
         }
     }
 
     /** Method called just before integration.
      * <p>
      * The default implementation does nothing, it may be specialized in subclasses.
      * </p>
      * @param initialState initial state
      * @param tEnd target date at which state should be propagated
      */
     @Override
     protected void beforeIntegration(final SpacecraftState initialState,
                                      final AbsoluteDate tEnd) {
 
         // check if only mean elements must be used
         final PropagationType type = getPropagationType();
 
         // compute common auxiliary elements
         final AuxiliaryElements aux = new AuxiliaryElements(initialState.getOrbit(), I);
 
         // initialize all perturbing forces
         final List<ShortPeriodTerms> shortPeriodTerms = new ArrayList<ShortPeriodTerms>();
         for (final DSSTForceModel force : forceModels) {
             shortPeriodTerms.addAll(force.initializeShortPeriodTerms(aux, type, force.getParameters()));
         }
         mapper.setShortPeriodTerms(shortPeriodTerms);
 
         // if required, insert the special short periodics step handler
         if (type == PropagationType.OSCULATING) {
             final ShortPeriodicsHandler spHandler = new ShortPeriodicsHandler(forceModels);
             // Compute short periodic coefficients for this point
             for (DSSTForceModel forceModel : forceModels) {
                 forceModel.updateShortPeriodTerms(forceModel.getParameters(), initialState);
             }
             final Collection<ODEStepHandler> stepHandlers = new ArrayList<ODEStepHandler>();
             stepHandlers.add(spHandler);
             final ODEIntegrator integrator = getIntegrator();
             final Collection<ODEStepHandler> existing = integrator.getStepHandlers();
             stepHandlers.addAll(existing);
 
             integrator.clearStepHandlers();
 
             // add back the existing handlers after the short periodics one
             for (final ODEStepHandler sp : stepHandlers) {
                 integrator.addStepHandler(sp);
             }
         }
     }
 
     /** {@inheritDoc} */
     @Override
     protected void afterIntegration() {
         // remove the special short periodics step handler if added before
         if (getPropagationType() == PropagationType.OSCULATING) {
             final List<ODEStepHandler> preserved = new ArrayList<ODEStepHandler>();
             final ODEIntegrator integrator = getIntegrator();
             for (final ODEStepHandler sp : integrator.getStepHandlers()) {
                 if (!(sp instanceof ShortPeriodicsHandler)) {
                     preserved.add(sp);
                 }
             }
 
             // clear the list
             integrator.clearStepHandlers();
 
             // add back the step handlers that were important for the user
             for (final ODEStepHandler sp : preserved) {
                 integrator.addStepHandler(sp);
             }
         }
     }
 
     /** Compute mean state from osculating state.
      * <p>
      * Compute in a DSST sense the mean state corresponding to the input osculating state.
      * </p><p>
      * The computing is done through a fixed-point iteration process.
      * </p>
      * @param osculating initial osculating state
      * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
      * like atmospheric drag, radiation pressure or specific user-defined models)
      * @param forceModels force models
      * @param epsilon convergence threshold for mean parameters conversion
      * @param maxIterations maximum iterations for mean parameters conversion
      * @return mean state
      */
     private static Orbit computeMeanOrbit(final SpacecraftState osculating,
                                           final AttitudeProvider attitudeProvider,
                                           final Collection<DSSTForceModel> forceModels, final double epsilon, final int maxIterations) {
 
         // rough initialization of the mean parameters
         EquinoctialOrbit meanOrbit = (EquinoctialOrbit) OrbitType.EQUINOCTIAL.convertType(osculating.getOrbit());
 
         // threshold for each parameter
         final double thresholdA = epsilon * (1 + FastMath.abs(meanOrbit.getA()));
         final double thresholdE = epsilon * (1 + meanOrbit.getE());
         final double thresholdI = epsilon * (1 + meanOrbit.getI());
         final double thresholdL = epsilon * FastMath.PI;
 
         // ensure all Gaussian force models can rely on attitude
         for (final DSSTForceModel force : forceModels) {
             force.registerAttitudeProvider(attitudeProvider);
         }
 
         int i = 0;
         while (i++ < maxIterations) {
 
             final SpacecraftState meanState = new SpacecraftState(meanOrbit, osculating.getAttitude(), osculating.getMass());
 
             //Create the auxiliary object
             final AuxiliaryElements aux = new AuxiliaryElements(meanOrbit, I);
 
             // Set the force models
             final List<ShortPeriodTerms> shortPeriodTerms = new ArrayList<ShortPeriodTerms>();
             for (final DSSTForceModel force : forceModels) {
                 shortPeriodTerms.addAll(force.initializeShortPeriodTerms(aux, PropagationType.OSCULATING, force.getParameters()));
                 force.updateShortPeriodTerms(force.getParameters(), meanState);
             }
 
             // recompute the osculating parameters from the current mean parameters
             final EquinoctialOrbit rebuilt = computeOsculatingOrbit(meanState, shortPeriodTerms);
 
             // adapted parameters residuals
             final double deltaA  = osculating.getA() - rebuilt.getA();
             final double deltaEx = osculating.getEquinoctialEx() - rebuilt.getEquinoctialEx();
             final double deltaEy = osculating.getEquinoctialEy() - rebuilt.getEquinoctialEy();
             final double deltaHx = osculating.getHx() - rebuilt.getHx();
             final double deltaHy = osculating.getHy() - rebuilt.getHy();
             final double deltaLv = MathUtils.normalizeAngle(osculating.getLv() - rebuilt.getLv(), 0.0);
 
             // check convergence
             if (FastMath.abs(deltaA)  < thresholdA &&
                 FastMath.abs(deltaEx) < thresholdE &&
                 FastMath.abs(deltaEy) < thresholdE &&
                 FastMath.abs(deltaHx) < thresholdI &&
                 FastMath.abs(deltaHy) < thresholdI &&
                 FastMath.abs(deltaLv) < thresholdL) {
                 return meanOrbit;
             }
 
             // update mean parameters
             meanOrbit = new EquinoctialOrbit(meanOrbit.getA() + deltaA,
                                              meanOrbit.getEquinoctialEx() + deltaEx,
                                              meanOrbit.getEquinoctialEy() + deltaEy,
                                              meanOrbit.getHx() + deltaHx,
                                              meanOrbit.getHy() + deltaHy,
                                              meanOrbit.getLv() + deltaLv,
                                              PositionAngle.TRUE, meanOrbit.getFrame(),
                                              meanOrbit.getDate(), meanOrbit.getMu());
         }
 
         throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_DSST_MEAN_PARAMETERS, i);
 
     }
 
     /** Compute osculating state from mean state.
      * <p>
      * Compute and add the short periodic variation to the mean {@link SpacecraftState}.
      * </p>
      * @param meanState initial mean state
      * @param shortPeriodTerms short period terms
      * @return osculating state
      */
     private static EquinoctialOrbit computeOsculatingOrbit(final SpacecraftState meanState,
                                                            final List<ShortPeriodTerms> shortPeriodTerms) {
 
         final double[] mean = new double[6];
         final double[] meanDot = new double[6];
         OrbitType.EQUINOCTIAL.mapOrbitToArray(meanState.getOrbit(), PositionAngle.MEAN, mean, meanDot);
         final double[] y = mean.clone();
         for (final ShortPeriodTerms spt : shortPeriodTerms) {
             final double[] shortPeriodic = spt.value(meanState.getOrbit());
             for (int i = 0; i < shortPeriodic.length; i++) {
                 y[i] += shortPeriodic[i];
             }
         }
         return (EquinoctialOrbit) OrbitType.EQUINOCTIAL.mapArrayToOrbit(y, meanDot,
                                                                         PositionAngle.MEAN, meanState.getDate(),
                                                                         meanState.getMu(), meanState.getFrame());
     }
 
     /** {@inheritDoc} */
     @Override
     protected SpacecraftState getInitialIntegrationState() {
         if (initialIsOsculating) {
             // the initial state is an osculating state,
             // it must be converted to mean state
             return computeMeanState(getInitialState(), getAttitudeProvider(), forceModels);
         } else {
             // the initial state is already a mean state
             return getInitialState();
         }
     }
 
     /** {@inheritDoc}
      * <p>
      * Note that for DSST, orbit type is hardcoded to {@link OrbitType#EQUINOCTIAL}
      * and position angle type is hardcoded to {@link PositionAngle#MEAN}, so
      * the corresponding parameters are ignored.
      * </p>
      */
     @Override
     protected StateMapper createMapper(final AbsoluteDate referenceDate, final double mu,
                                        final OrbitType ignoredOrbitType, final PositionAngle ignoredPositionAngleType,
                                        final AttitudeProvider attitudeProvider, final Frame frame) {
 
         // create a mapper with the common settings provided as arguments
         final MeanPlusShortPeriodicMapper newMapper =
                 new MeanPlusShortPeriodicMapper(referenceDate, mu, attitudeProvider, frame);
 
         // copy the specific settings from the existing mapper
         if (mapper != null) {
             newMapper.setSatelliteRevolution(mapper.getSatelliteRevolution());
             newMapper.setSelectedCoefficients(mapper.getSelectedCoefficients());
             newMapper.setShortPeriodTerms(mapper.getShortPeriodTerms());
         }
 
         mapper = newMapper;
         return mapper;
 
     }
 
 
     /** Get the short period terms value.
      * @param meanState the mean state
      * @return shortPeriodTerms short period terms
      * @since 7.1
      */
     public double[] getShortPeriodTermsValue(final SpacecraftState meanState) {
         final double[] sptValue = new double[6];
 
         for (ShortPeriodTerms spt : mapper.getShortPeriodTerms()) {
             final double[] shortPeriodic = spt.value(meanState.getOrbit());
             for (int i = 0; i < shortPeriodic.length; i++) {
                 sptValue[i] += shortPeriodic[i];
             }
         }
         return sptValue;
     }
 
 
     /** Internal mapper using mean parameters plus short periodic terms. */
     private static class MeanPlusShortPeriodicMapper extends StateMapper {
 
         /** Short periodic coefficients that must be stored as additional states. */
         private Set<String>                selectedCoefficients;
 
         /** Number of satellite revolutions in the averaging interval. */
         private int                        satelliteRevolution;
 
         /** Short period terms. */
         private List<ShortPeriodTerms>     shortPeriodTerms;
 
         /** Simple constructor.
          * @param referenceDate reference date
          * @param mu central attraction coefficient (m³/s²)
          * @param attitudeProvider attitude provider
          * @param frame inertial frame
          */
         MeanPlusShortPeriodicMapper(final AbsoluteDate referenceDate, final double mu,
                                     final AttitudeProvider attitudeProvider, final Frame frame) {
 
             super(referenceDate, mu, OrbitType.EQUINOCTIAL, PositionAngle.MEAN, attitudeProvider, frame);
 
             this.selectedCoefficients = null;
 
             // Default averaging period for conversion from osculating to mean elements
             this.satelliteRevolution = 2;
 
             this.shortPeriodTerms    = Collections.emptyList();
 
         }
 
         /** {@inheritDoc} */
         @Override
         public SpacecraftState mapArrayToState(final AbsoluteDate date,
                                                final double[] y, final double[] yDot,
                                                final PropagationType type) {
 
             // add short periodic variations to mean elements to get osculating elements
             // (the loop may not be performed if there are no force models and in the
             //  case we want to remain in mean parameters only)
             final double[] elements = y.clone();
             final DoubleArrayDictionary coefficients;
             if (type == PropagationType.MEAN) {
                 coefficients = null;
             } else {
                 final Orbit meanOrbit = OrbitType.EQUINOCTIAL.mapArrayToOrbit(elements, yDot, PositionAngle.MEAN, date, getMu(), getFrame());
                 coefficients = selectedCoefficients == null ? null : new DoubleArrayDictionary();
                 for (final ShortPeriodTerms spt : shortPeriodTerms) {
                     final double[] shortPeriodic = spt.value(meanOrbit);
                     for (int i = 0; i < shortPeriodic.length; i++) {
                         elements[i] += shortPeriodic[i];
                     }
                     if (selectedCoefficients != null) {
                         coefficients.putAll(spt.getCoefficients(date, selectedCoefficients));
                     }
                 }
             }
 
             final double mass = elements[6];
             if (mass <= 0.0) {
                 throw new OrekitException(OrekitMessages.SPACECRAFT_MASS_BECOMES_NEGATIVE, mass);
             }
 
             final Orbit orbit       = OrbitType.EQUINOCTIAL.mapArrayToOrbit(elements, yDot, PositionAngle.MEAN, date, getMu(), getFrame());
             final Attitude attitude = getAttitudeProvider().getAttitude(orbit, date, getFrame());
 
             if (coefficients == null) {
                 return new SpacecraftState(orbit, attitude, mass);
             } else {
                 return new SpacecraftState(orbit, attitude, mass, coefficients);
             }
 
         }
 
         /** {@inheritDoc} */
         @Override
         public void mapStateToArray(final SpacecraftState state, final double[] y, final double[] yDot) {
 
             OrbitType.EQUINOCTIAL.mapOrbitToArray(state.getOrbit(), PositionAngle.MEAN, y, yDot);
             y[6] = state.getMass();
 
         }
 
         /** Set the number of satellite revolutions to use for converting osculating to mean elements.
          *  <p>
          *  By default, if the initial orbit is defined as osculating,
          *  it will be averaged over 2 satellite revolutions.
          *  This can be changed by using this method.
          *  </p>
          *  @param satelliteRevolution number of satellite revolutions to use for converting osculating to mean
          *                             elements
          */
         public void setSatelliteRevolution(final int satelliteRevolution) {
             this.satelliteRevolution = satelliteRevolution;
         }
 
         /** Get the number of satellite revolutions to use for converting osculating to mean elements.
          *  @return number of satellite revolutions to use for converting osculating to mean elements
          */
         public int getSatelliteRevolution() {
             return satelliteRevolution;
         }
 
         /** Set the selected short periodic coefficients that must be stored as additional states.
          * @param selectedCoefficients short periodic coefficients that must be stored as additional states
          * (null means no coefficients are selected, empty set means all coefficients are selected)
          */
         public void setSelectedCoefficients(final Set<String> selectedCoefficients) {
             this.selectedCoefficients = selectedCoefficients;
         }
 
         /** Get the selected short periodic coefficients that must be stored as additional states.
          * @return short periodic coefficients that must be stored as additional states
          * (null means no coefficients are selected, empty set means all coefficients are selected)
          */
         public Set<String> getSelectedCoefficients() {
             return selectedCoefficients;
         }
 
         /** Set the short period terms.
          * @param shortPeriodTerms short period terms
          * @since 7.1
          */
         public void setShortPeriodTerms(final List<ShortPeriodTerms> shortPeriodTerms) {
             this.shortPeriodTerms = shortPeriodTerms;
         }
 
         /** Get the short period terms.
          * @return shortPeriodTerms short period terms
          * @since 7.1
          */
         public List<ShortPeriodTerms> getShortPeriodTerms() {
             return shortPeriodTerms;
         }
 
     }
 
     /** {@inheritDoc} */
     @Override
     protected MainStateEquations getMainStateEquations(final ODEIntegrator integrator) {
         return new Main(integrator);
     }
 
     /** Internal class for mean parameters integration. */
     private class Main implements MainStateEquations {
 
         /** Derivatives array. */
         private final double[] yDot;
 
         /** Simple constructor.
          * @param integrator numerical integrator to use for propagation.
          */
         Main(final ODEIntegrator integrator) {
             yDot = new double[7];
 
             for (final DSSTForceModel forceModel : forceModels) {
                 final EventDetector[] modelDetectors = forceModel.getEventsDetectors();
                 if (modelDetectors != null) {
                     for (final EventDetector detector : modelDetectors) {
                         setUpEventDetector(integrator, detector);
                     }
                 }
             }
 
         }
 
         /** {@inheritDoc} */
         @Override
         public void init(final SpacecraftState initialState, final AbsoluteDate target) {
             forceModels.forEach(fm -> fm.init(initialState, target));
         }
 
         /** {@inheritDoc} */
         @Override
         public double[] computeDerivatives(final SpacecraftState state) {
 
             Arrays.fill(yDot, 0.0);
 
             // compute common auxiliary elements
             final AuxiliaryElements auxiliaryElements = new AuxiliaryElements(state.getOrbit(), I);
 
             // compute the contributions of all perturbing forces
             for (final DSSTForceModel forceModel : forceModels) {
                 final double[] daidt = elementRates(forceModel, state, auxiliaryElements, forceModel.getParameters());
                 for (int i = 0; i < daidt.length; i++) {
                     yDot[i] += daidt[i];
                 }
             }
 
             return yDot.clone();
         }
 
         /** This method allows to compute the mean equinoctial elements rates da<sub>i</sub> / dt
          *  for a specific force model.
          *  @param forceModel force to take into account
          *  @param state current state
          *  @param auxiliaryElements auxiliary elements related to the current orbit
          *  @param parameters force model parameters
          *  @return the mean equinoctial elements rates da<sub>i</sub> / dt
          */
         private double[] elementRates(final DSSTForceModel forceModel,
                                       final SpacecraftState state,
                                       final AuxiliaryElements auxiliaryElements,
                                       final double[] parameters) {
             return forceModel.getMeanElementRate(state, auxiliaryElements, parameters);
         }
 
     }
 
     /** Estimate tolerance vectors for an AdaptativeStepsizeIntegrator.
      *  <p>
      *  The errors are estimated from partial derivatives properties of orbits,
      *  starting from a scalar position error specified by the user.
      *  Considering the energy conservation equation V = sqrt(mu (2/r - 1/a)),
      *  we get at constant energy (i.e. on a Keplerian trajectory):
      *
      *  <pre>
      *  V r² |dV| = mu |dr|
      *  </pre>
      *
      *  <p> So we deduce a scalar velocity error consistent with the position error. From here, we apply
      *  orbits Jacobians matrices to get consistent errors on orbital parameters.
      *
      *  <p>
      *  The tolerances are only <em>orders of magnitude</em>, and integrator tolerances are only
      *  local estimates, not global ones. So some care must be taken when using these tolerances.
      *  Setting 1mm as a position error does NOT mean the tolerances will guarantee a 1mm error
      *  position after several orbits integration.
      *  </p>
      *
      * @param dP user specified position error (m)
      * @param orbit reference orbit
      * @return a two rows array, row 0 being the absolute tolerance error
      *                       and row 1 being the relative tolerance error
      */
     public static double[][] tolerances(final double dP, final Orbit orbit) {
 
         return NumericalPropagator.tolerances(dP, orbit, OrbitType.EQUINOCTIAL);
 
     }
 
     /** Estimate tolerance vectors for an AdaptativeStepsizeIntegrator.
      *  <p>
      *  The errors are estimated from partial derivatives properties of orbits,
      *  starting from scalar position and velocity errors specified by the user.
      *  <p>
      *  The tolerances are only <em>orders of magnitude</em>, and integrator tolerances are only
      *  local estimates, not global ones. So some care must be taken when using these tolerances.
      *  Setting 1mm as a position error does NOT mean the tolerances will guarantee a 1mm error
      *  position after several orbits integration.
      *  </p>
      *
      * @param dP user specified position error (m)
      * @param dV user specified velocity error (m/s)
      * @param orbit reference orbit
      * @return a two rows array, row 0 being the absolute tolerance error
      *                       and row 1 being the relative tolerance error
      * @since 10.3
      */
     public static double[][] tolerances(final double dP, final double dV, final Orbit orbit) {
 
         return NumericalPropagator.tolerances(dP, dV, orbit, OrbitType.EQUINOCTIAL);
 
     }
 
     /** Step handler used to compute the parameters for the short periodic contributions.
      * @author Lucian Barbulescu
      */
     private class ShortPeriodicsHandler implements ODEStepHandler {
 
         /** Force models used to compute short periodic terms. */
         private final List<DSSTForceModel> forceModels;
 
         /** Constructor.
          * @param forceModels force models
          */
         ShortPeriodicsHandler(final List<DSSTForceModel> forceModels) {
             this.forceModels = forceModels;
         }
 
         /** {@inheritDoc} */
         @Override
         public void handleStep(final ODEStateInterpolator interpolator) {
 
             // Get the grid points to compute
             final double[] interpolationPoints =
                     interpolationgrid.getGridPoints(interpolator.getPreviousState().getTime(),
                                                     interpolator.getCurrentState().getTime());
 
             final SpacecraftState[] meanStates = new SpacecraftState[interpolationPoints.length];
             for (int i = 0; i < interpolationPoints.length; ++i) {
 
                 // Build the mean state interpolated at grid point
                 final double time = interpolationPoints[i];
                 final ODEStateAndDerivative sd = interpolator.getInterpolatedState(time);
                 meanStates[i] = mapper.mapArrayToState(time,
                                                        sd.getPrimaryState(),
                                                        sd.getPrimaryDerivative(),
                                                        PropagationType.MEAN);
             }
 
             // Computate short periodic coefficients for this step
             for (DSSTForceModel forceModel : forceModels) {
                 forceModel.updateShortPeriodTerms(forceModel.getParameters(), meanStates);
             }
 
         }
     }
 }