FieldDSSTPropagator.java

  1. /* Copyright 2002-2024 CS GROUP
  2.  * Licensed to CS GROUP (CS) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * CS licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *   http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.orekit.propagation.semianalytical.dsst;

  19. import org.hipparchus.CalculusFieldElement;
  20. import org.hipparchus.Field;
  21. import org.hipparchus.ode.FieldODEIntegrator;
  22. import org.hipparchus.ode.FieldODEStateAndDerivative;
  23. import org.hipparchus.ode.sampling.FieldODEStateInterpolator;
  24. import org.hipparchus.ode.sampling.FieldODEStepHandler;
  25. import org.hipparchus.util.FastMath;
  26. import org.hipparchus.util.MathArrays;
  27. import org.hipparchus.util.MathUtils;
  28. import org.orekit.annotation.DefaultDataContext;
  29. import org.orekit.attitudes.AttitudeProvider;
  30. import org.orekit.attitudes.FieldAttitude;
  31. import org.orekit.data.DataContext;
  32. import org.orekit.errors.OrekitException;
  33. import org.orekit.errors.OrekitInternalError;
  34. import org.orekit.errors.OrekitMessages;
  35. import org.orekit.frames.Frame;
  36. import org.orekit.orbits.FieldEquinoctialOrbit;
  37. import org.orekit.orbits.FieldOrbit;
  38. import org.orekit.orbits.OrbitType;
  39. import org.orekit.orbits.PositionAngleType;
  40. import org.orekit.propagation.FieldSpacecraftState;
  41. import org.orekit.propagation.PropagationType;
  42. import org.orekit.propagation.Propagator;
  43. import org.orekit.propagation.SpacecraftState;
  44. import org.orekit.propagation.integration.FieldAbstractIntegratedPropagator;
  45. import org.orekit.propagation.integration.FieldStateMapper;
  46. import org.orekit.propagation.numerical.FieldNumericalPropagator;
  47. import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
  48. import org.orekit.propagation.semianalytical.dsst.forces.DSSTNewtonianAttraction;
  49. import org.orekit.propagation.semianalytical.dsst.forces.FieldShortPeriodTerms;
  50. import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
  51. import org.orekit.propagation.semianalytical.dsst.utilities.FieldFixedNumberInterpolationGrid;
  52. import org.orekit.propagation.semianalytical.dsst.utilities.FieldInterpolationGrid;
  53. import org.orekit.propagation.semianalytical.dsst.utilities.FieldMaxGapInterpolationGrid;
  54. import org.orekit.time.AbsoluteDate;
  55. import org.orekit.time.FieldAbsoluteDate;
  56. import org.orekit.utils.FieldArrayDictionary;
  57. import org.orekit.utils.ParameterDriver;
  58. import org.orekit.utils.ParameterObserver;
  59. import org.orekit.utils.TimeSpanMap;

  61. import java.util.ArrayList;
  62. import java.util.Arrays;
  63. import java.util.Collection;
  64. import java.util.Collections;
  65. import java.util.HashSet;
  66. import java.util.List;
  67. import java.util.Set;

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

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

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

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

  151.     /** Default value for maxIterations. */
  152.     private static final int MAX_ITERATIONS_DEFAULT = 200;

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

  157.     /** Field used by this class.*/
  158.     private final Field<T> field;

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

  163.     /** State mapper holding the force models. */
  164.     private FieldMeanPlusShortPeriodicMapper mapper;

  166.     /** Generator for the interpolation grid. */
  167.     private FieldInterpolationGrid<T> interpolationgrid;

  169.     /** Create a new instance of DSSTPropagator.
  170.      *  <p>
  171.      *  After creation, there are no perturbing forces at all.
  172.      *  This means that if {@link #addForceModel addForceModel}
  173.      *  is not called after creation, the integrated orbit will
  174.      *  follow a Keplerian evolution only.
  175.      *  </p>
  176.      *
  177.      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
  178.      *
  179.      *  @param field field used by default
  180.      *  @param integrator numerical integrator to use for propagation.
  181.      *  @param propagationType type of orbit to output (mean or osculating).
  182.      * @see #FieldDSSTPropagator(Field, FieldODEIntegrator, PropagationType,
  183.      * AttitudeProvider)
  184.      */
  185.     @DefaultDataContext
  186.     public FieldDSSTPropagator(final Field<T> field, final FieldODEIntegrator<T> integrator, final PropagationType propagationType) {
  187.         this(field, integrator, propagationType,
  188.                 Propagator.getDefaultLaw(DataContext.getDefault().getFrames()));
  189.     }

  191.     /** Create a new instance of DSSTPropagator.
  192.      *  <p>
  193.      *  After creation, there are no perturbing forces at all.
  194.      *  This means that if {@link #addForceModel addForceModel}
  195.      *  is not called after creation, the integrated orbit will
  196.      *  follow a Keplerian evolution only.
  197.      *  </p>
  198.      * @param field field used by default
  199.      *  @param integrator numerical integrator to use for propagation.
  200.      * @param propagationType type of orbit to output (mean or osculating).
  201.      * @param attitudeProvider attitude law to use.
  202.      * @since 10.1
  203.      */
  204.     public FieldDSSTPropagator(final Field<T> field,
  205.                                final FieldODEIntegrator<T> integrator,
  206.                                final PropagationType propagationType,
  207.                                final AttitudeProvider attitudeProvider) {
  208.         super(field, integrator, propagationType);
  209.         this.field  = field;
  210.         forceModels = new ArrayList<>();
  211.         initMapper(field);
  212.         // DSST uses only equinoctial orbits and mean longitude argument
  213.         setOrbitType(OrbitType.EQUINOCTIAL);
  214.         setPositionAngleType(PositionAngleType.MEAN);
  215.         setAttitudeProvider(attitudeProvider);
  216.         setInterpolationGridToFixedNumberOfPoints(INTERPOLATION_POINTS_PER_STEP);
  217.     }

  219.     /** Create a new instance of DSSTPropagator.
  220.      *  <p>
  221.      *  After creation, there are no perturbing forces at all.
  222.      *  This means that if {@link #addForceModel addForceModel}
  223.      *  is not called after creation, the integrated orbit will
  224.      *  follow a Keplerian evolution only. Only the mean orbits
  225.      *  will be generated.
  226.      *  </p>
  227.      *
  228.      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
  229.      *
  230.      *  @param field fied used by default
  231.      *  @param integrator numerical integrator to use for propagation.
  232.      * @see #FieldDSSTPropagator(Field, FieldODEIntegrator, AttitudeProvider)
  233.      */
  234.     @DefaultDataContext
  235.     public FieldDSSTPropagator(final Field<T> field, final FieldODEIntegrator<T> integrator) {
  236.         this(field, integrator,
  237.                 Propagator.getDefaultLaw(DataContext.getDefault().getFrames()));
  238.     }

  240.     /** Create a new instance of DSSTPropagator.
  241.      *  <p>
  242.      *  After creation, there are no perturbing forces at all.
  243.      *  This means that if {@link #addForceModel addForceModel}
  244.      *  is not called after creation, the integrated orbit will
  245.      *  follow a Keplerian evolution only. Only the mean orbits
  246.      *  will be generated.
  247.      *  </p>
  248.      * @param field fied used by default
  249.      *  @param integrator numerical integrator to use for propagation.
  250.      * @param attitudeProvider attitude law to use.
  251.      * @since 10.1
  252.      */
  253.     public FieldDSSTPropagator(final Field<T> field,
  254.                                final FieldODEIntegrator<T> integrator,
  255.                                final AttitudeProvider attitudeProvider) {
  256.         super(field, integrator, PropagationType.MEAN);
  257.         this.field  = field;
  258.         forceModels = new ArrayList<>();
  259.         initMapper(field);
  260.         // DSST uses only equinoctial orbits and mean longitude argument
  261.         setOrbitType(OrbitType.EQUINOCTIAL);
  262.         setPositionAngleType(PositionAngleType.MEAN);
  263.         setAttitudeProvider(attitudeProvider);
  264.         setInterpolationGridToFixedNumberOfPoints(INTERPOLATION_POINTS_PER_STEP);
  265.     }

  267.     /** Set the central attraction coefficient μ.
  268.      * <p>
  269.      * Setting the central attraction coefficient is
  270.      * equivalent to {@link #addForceModel(DSSTForceModel) add}
  271.      * a {@link DSSTNewtonianAttraction} force model.
  272.      * </p>
  273.     * @param mu central attraction coefficient (m³/s²)
  274.     * @see #addForceModel(DSSTForceModel)
  275.     * @see #getAllForceModels()
  276.     */
  277.     public void setMu(final T mu) {
  278.         addForceModel(new DSSTNewtonianAttraction(mu.getReal()));
  279.     }

  281.     /** Set the central attraction coefficient μ only in upper class.
  282.      * @param mu central attraction coefficient (m³/s²)
  283.      */
  284.     private void superSetMu(final T mu) {
  285.         super.setMu(mu);
  286.     }

  288.     /** Check if Newtonian attraction force model is available.
  289.      * <p>
  290.      * Newtonian attraction is always the last force model in the list.
  291.      * </p>
  292.      * @return true if Newtonian attraction force model is available
  293.      */
  294.     private boolean hasNewtonianAttraction() {
  295.         final int last = forceModels.size() - 1;
  296.         return last >= 0 && forceModels.get(last) instanceof DSSTNewtonianAttraction;
  297.     }

  299.     /** Set the initial state with osculating orbital elements.
  300.      *  @param initialState initial state (defined with osculating elements)
  301.      */
  302.     public void setInitialState(final FieldSpacecraftState<T> initialState) {
  303.         setInitialState(initialState, PropagationType.OSCULATING);
  304.     }

  306.     /** Set the initial state.
  307.      *  @param initialState initial state
  308.      *  @param stateType defined if the orbital state is defined with osculating or mean elements
  309.      */
  310.     public void setInitialState(final FieldSpacecraftState<T> initialState,
  311.                                 final PropagationType stateType) {
  312.         switch (stateType) {
  313.             case MEAN:
  314.                 initialIsOsculating = false;
  315.                 break;
  316.             case OSCULATING:
  317.                 initialIsOsculating = true;
  318.                 break;
  319.             default:
  320.                 throw new OrekitInternalError(null);
  321.         }
  322.         resetInitialState(initialState);
  323.     }

  325.     /** Reset the initial state.
  326.      *
  327.      *  @param state new initial state
  328.      */
  329.     @Override
  330.     public void resetInitialState(final FieldSpacecraftState<T> state) {
  331.         super.resetInitialState(state);
  332.         if (!hasNewtonianAttraction()) {
  333.             // use the state to define central attraction
  334.             setMu(state.getMu());
  335.         }
  336.         super.setStartDate(state.getDate());
  337.     }

  339.     /** Set the selected short periodic coefficients that must be stored as additional states.
  340.      * @param selectedCoefficients short periodic coefficients that must be stored as additional states
  341.      * (null means no coefficients are selected, empty set means all coefficients are selected)
  342.      */
  343.     public void setSelectedCoefficients(final Set<String> selectedCoefficients) {
  344.         mapper.setSelectedCoefficients(selectedCoefficients == null ?
  345.                                        null : new HashSet<String>(selectedCoefficients));
  346.     }

  348.     /** Get the selected short periodic coefficients that must be stored as additional states.
  349.      * @return short periodic coefficients that must be stored as additional states
  350.      * (null means no coefficients are selected, empty set means all coefficients are selected)
  351.      */
  352.     public Set<String> getSelectedCoefficients() {
  353.         final Set<String> set = mapper.getSelectedCoefficients();
  354.         return set == null ? null : Collections.unmodifiableSet(set);
  355.     }

  357.     /** Check if the initial state is provided in osculating elements.
  358.      * @return true if initial state is provided in osculating elements
  359.      */
  360.     public boolean initialIsOsculating() {
  361.         return initialIsOsculating;
  362.     }

  364.     /** Set the interpolation grid generator.
  365.      * <p>
  366.      * The generator will create an interpolation grid with a fixed
  367.      * number of points for each mean element integration step.
  368.      * </p>
  369.      * <p>
  370.      * If neither {@link #setInterpolationGridToFixedNumberOfPoints(int)}
  371.      * nor {@link #setInterpolationGridToMaxTimeGap(CalculusFieldElement)} has been called,
  372.      * by default the propagator is set as to 3 interpolations points per step.
  373.      * </p>
  374.      * @param interpolationPoints number of interpolation points at
  375.      * each integration step
  376.      * @see #setInterpolationGridToMaxTimeGap(CalculusFieldElement)
  377.      * @since 7.1
  378.      */
  379.     public void setInterpolationGridToFixedNumberOfPoints(final int interpolationPoints) {
  380.         interpolationgrid = new FieldFixedNumberInterpolationGrid<>(field, interpolationPoints);
  381.     }

  383.     /** Set the interpolation grid generator.
  384.      * <p>
  385.      * The generator will create an interpolation grid with a maximum
  386.      * time gap between interpolation points.
  387.      * </p>
  388.      * <p>
  389.      * If neither {@link #setInterpolationGridToFixedNumberOfPoints(int)}
  390.      * nor {@link #setInterpolationGridToMaxTimeGap(CalculusFieldElement)} has been called,
  391.      * by default the propagator is set as to 3 interpolations points per step.
  392.      * </p>
  393.      * @param maxGap maximum time gap between interpolation points (seconds)
  394.      * @see #setInterpolationGridToFixedNumberOfPoints(int)
  395.      * @since 7.1
  396.      */
  397.     public void setInterpolationGridToMaxTimeGap(final T maxGap) {
  398.         interpolationgrid = new FieldMaxGapInterpolationGrid<>(field, maxGap);
  399.     }

  401.     /** Add a force model to the global perturbation model.
  402.      *  <p>
  403.      *  If this method is not called at all,
  404.      *  the integrated orbit will follow a Keplerian evolution only.
  405.      *  </p>
  406.      *  @param force perturbing {@link DSSTForceModel force} to add
  407.      *  @see #removeForceModels()
  408.      *  @see #setMu(CalculusFieldElement)
  409.      */
  410.     public void addForceModel(final DSSTForceModel force) {

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

  415.             try {
  416.                 // ensure we are notified of any mu change
  417.                 force.getParametersDrivers().get(0).addObserver(new ParameterObserver() {
  418.                     /** {@inheritDoc} */
  419.                     @Override
  420.                     public void valueChanged(final double previousValue, final ParameterDriver driver, final AbsoluteDate date) {
  421.                         // mu PDriver should have only 1 span
  422.                         superSetMu(field.getZero().newInstance(driver.getValue()));
  423.                     }
  424.                     /** {@inheritDoc} */
  425.                     @Override
  426.                     public void valueSpanMapChanged(final TimeSpanMap<Double> previousValue, final ParameterDriver driver) {
  427.                         // mu PDriver should have only 1 span
  428.                         superSetMu(field.getZero().newInstance(driver.getValue()));
  429.                     }
  430.                 });
  431.             } catch (OrekitException oe) {
  432.                 // this should never happen
  433.                 throw new OrekitInternalError(oe);
  434.             }

  436.             if (hasNewtonianAttraction()) {
  437.                 // there is already a central attraction model, replace it
  438.                 forceModels.set(forceModels.size() - 1, force);
  439.             } else {
  440.                 // there are no central attraction model yet, add it at the end of the list
  441.                 forceModels.add(force);
  442.             }
  443.         } else {
  444.             // we want to add a perturbing force model
  445.             if (hasNewtonianAttraction()) {
  446.                 // insert the new force model before Newtonian attraction,
  447.                 // which should always be the last one in the list
  448.                 forceModels.add(forceModels.size() - 1, force);
  449.             } else {
  450.                 // we only have perturbing force models up to now, just append at the end of the list
  451.                 forceModels.add(force);
  452.             }
  453.         }

  455.         force.registerAttitudeProvider(getAttitudeProvider());

  457.     }

  459.     /** Remove all perturbing force models from the global perturbation model
  460.      *  (except central attraction).
  461.      *  <p>
  462.      *  Once all perturbing forces have been removed (and as long as no new force model is added),
  463.      *  the integrated orbit will follow a Keplerian evolution only.
  464.      *  </p>
  465.      *  @see #addForceModel(DSSTForceModel)
  466.      */
  467.     public void removeForceModels() {
  468.         final int last = forceModels.size() - 1;
  469.         if (hasNewtonianAttraction()) {
  470.             // preserve the Newtonian attraction model at the end
  471.             final DSSTForceModel newton = forceModels.get(last);
  472.             forceModels.clear();
  473.             forceModels.add(newton);
  474.         } else {
  475.             forceModels.clear();
  476.         }
  477.     }

  479.     /** Get all the force models, perturbing forces and Newtonian attraction included.
  480.      * @return list of perturbing force models, with Newtonian attraction being the
  481.      * last one
  482.      * @see #addForceModel(DSSTForceModel)
  483.      * @see #setMu(CalculusFieldElement)
  484.      */
  485.     public List<DSSTForceModel> getAllForceModels() {
  486.         return Collections.unmodifiableList(forceModels);
  487.     }

  489.     /** Get propagation parameter type.
  490.      * @return orbit type used for propagation
  491.      */
  492.     public OrbitType getOrbitType() {
  493.         return super.getOrbitType();
  494.     }

  496.     /** Get propagation parameter type.
  497.      * @return angle type to use for propagation
  498.      */
  499.     public PositionAngleType getPositionAngleType() {
  500.         return super.getPositionAngleType();
  501.     }

  503.     /** Conversion from mean to osculating orbit.
  504.      * <p>
  505.      * Compute osculating state <b>in a DSST sense</b>, corresponding to the
  506.      * mean SpacecraftState in input, and according to the Force models taken
  507.      * into account.
  508.      * </p><p>
  509.      * Since the osculating state is obtained by adding short-periodic variation
  510.      * of each force model, the resulting output will depend on the
  511.      * force models parameterized in input.
  512.      * </p>
  513.      * @param mean Mean state to convert
  514.      * @param forces Forces to take into account
  515.      * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
  516.      * like atmospheric drag, radiation pressure or specific user-defined models)
  517.      * @param <T> type of the elements
  518.      * @return osculating state in a DSST sense
  519.      */
  520.     @SuppressWarnings("unchecked")
  521.     public static <T extends CalculusFieldElement<T>> FieldSpacecraftState<T> computeOsculatingState(final FieldSpacecraftState<T> mean,
  522.                                                                                                  final AttitudeProvider attitudeProvider,
  523.                                                                                                  final Collection<DSSTForceModel> forces) {

  525.         //Create the auxiliary object
  526.         final FieldAuxiliaryElements<T> aux = new FieldAuxiliaryElements<>(mean.getOrbit(), I);

  528.         // Set the force models
  529.         final List<FieldShortPeriodTerms<T>> shortPeriodTerms = new ArrayList<>();
  530.         for (final DSSTForceModel force : forces) {
  531.             force.registerAttitudeProvider(attitudeProvider);
  532.             shortPeriodTerms.addAll(force.initializeShortPeriodTerms(aux, PropagationType.OSCULATING, force.getParameters(mean.getDate().getField(), mean.getDate())));
  533.             force.updateShortPeriodTerms(force.getParametersAllValues(mean.getDate().getField()), mean);
  534.         }

  536.         final FieldEquinoctialOrbit<T> osculatingOrbit = computeOsculatingOrbit(mean, shortPeriodTerms);

  538.         return new FieldSpacecraftState<>(osculatingOrbit, mean.getAttitude(), mean.getMass(),
  539.                                           mean.getAdditionalStatesValues());

  541.     }

  543.     /** Conversion from osculating to mean orbit.
  544.      * <p>
  545.      * Compute mean state <b>in a DSST sense</b>, corresponding to the
  546.      * osculating SpacecraftState in input, and according to the Force models
  547.      * taken into account.
  548.      * </p><p>
  549.      * Since the osculating state is obtained with the computation of
  550.      * short-periodic variation of each force model, the resulting output will
  551.      * depend on the force models parameterized in input.
  552.      * </p><p>
  553.      * The computation is done through a fixed-point iteration process.
  554.      * </p>
  555.      * @param osculating Osculating state to convert
  556.      * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
  557.      * like atmospheric drag, radiation pressure or specific user-defined models)
  558.      * @param forceModel Forces to take into account
  559.      * @param <T> type of the elements
  560.      * @return mean state in a DSST sense
  561.      */
  562.     public static <T extends CalculusFieldElement<T>> FieldSpacecraftState<T> computeMeanState(final FieldSpacecraftState<T> osculating,
  563.                                                                                            final AttitudeProvider attitudeProvider,
  564.                                                                                            final Collection<DSSTForceModel> forceModel) {
  565.         return computeMeanState(osculating, attitudeProvider, forceModel, EPSILON_DEFAULT, MAX_ITERATIONS_DEFAULT);
  566.     }

  568.     /** Conversion from osculating to mean orbit.
  569.      * <p>
  570.      * Compute mean state <b>in a DSST sense</b>, corresponding to the
  571.      * osculating SpacecraftState in input, and according to the Force models
  572.      * taken into account.
  573.      * </p><p>
  574.      * Since the osculating state is obtained with the computation of
  575.      * short-periodic variation of each force model, the resulting output will
  576.      * depend on the force models parameterized in input.
  577.      * </p><p>
  578.      * The computation is done through a fixed-point iteration process.
  579.      * </p>
  580.      * @param osculating Osculating state to convert
  581.      * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
  582.      * like atmospheric drag, radiation pressure or specific user-defined models)
  583.      * @param forceModel Forces to take into account
  584.      * @param epsilon convergence threshold for mean parameters conversion
  585.      * @param maxIterations maximum iterations for mean parameters conversion
  586.      * @return mean state in a DSST sense
  587.      * @param <T> type of the elements
  588.      * @since 10.1
  589.      */
  590.     public static <T extends CalculusFieldElement<T>> FieldSpacecraftState<T> computeMeanState(final FieldSpacecraftState<T> osculating,
  591.                                                                                            final AttitudeProvider attitudeProvider,
  592.                                                                                            final Collection<DSSTForceModel> forceModel,
  593.                                                                                            final double epsilon,
  594.                                                                                            final int maxIterations) {
  595.         final FieldOrbit<T> meanOrbit = computeMeanOrbit(osculating, attitudeProvider, forceModel, epsilon, maxIterations);
  596.         return new FieldSpacecraftState<>(meanOrbit, osculating.getAttitude(), osculating.getMass(),
  597.                                           osculating.getAdditionalStatesValues());
  598.     }

  600.     /** Override the default value of the parameter.
  601.     *  <p>
  602.     *  By default, if the initial orbit is defined as osculating,
  603.     *  it will be averaged over 2 satellite revolutions.
  604.     *  This can be changed by using this method.
  605.     *  </p>
  606.     *  @param satelliteRevolution number of satellite revolutions to use for converting osculating to mean
  607.     *                             elements
  608.     */
  609.     public void setSatelliteRevolution(final int satelliteRevolution) {
  610.         mapper.setSatelliteRevolution(satelliteRevolution);
  611.     }

  613.    /** Get the number of satellite revolutions to use for converting osculating to mean elements.
  614.     *  @return number of satellite revolutions to use for converting osculating to mean elements
  615.     */
  616.     public int getSatelliteRevolution() {
  617.         return mapper.getSatelliteRevolution();
  618.     }

  620.    /** {@inheritDoc} */
  621.     @Override
  622.     public void setAttitudeProvider(final AttitudeProvider attitudeProvider) {
  623.         super.setAttitudeProvider(attitudeProvider);

  625.         //Register the attitude provider for each force model
  626.         for (final DSSTForceModel force : forceModels) {
  627.             force.registerAttitudeProvider(attitudeProvider);
  628.         }
  629.     }

  631.     /** Method called just before integration.
  632.      * <p>
  633.      * The default implementation does nothing, it may be specialized in subclasses.
  634.      * </p>
  635.      * @param initialState initial state
  636.      * @param tEnd target date at which state should be propagated
  637.      */
  638.     @SuppressWarnings("unchecked")
  639.     @Override
  640.     protected void beforeIntegration(final FieldSpacecraftState<T> initialState,
  641.                                      final FieldAbsoluteDate<T> tEnd) {

  643.         // check if only mean elements must be used
  644.         final PropagationType type = isMeanOrbit();

  646.         // compute common auxiliary elements
  647.         final FieldAuxiliaryElements<T> aux = new FieldAuxiliaryElements<>(initialState.getOrbit(), I);

  649.         // initialize all perturbing forces
  650.         final List<FieldShortPeriodTerms<T>> shortPeriodTerms = new ArrayList<>();
  651.         for (final DSSTForceModel force : forceModels) {
  652.             shortPeriodTerms.addAll(force.initializeShortPeriodTerms(aux, type, force.getParameters(field, initialState.getDate())));
  653.         }
  654.         mapper.setShortPeriodTerms(shortPeriodTerms);

  656.         // if required, insert the special short periodics step handler
  657.         if (type == PropagationType.OSCULATING) {
  658.             final FieldShortPeriodicsHandler spHandler = new FieldShortPeriodicsHandler(forceModels);
  659.             // Compute short periodic coefficients for this point
  660.             for (DSSTForceModel forceModel : forceModels) {
  661.                 forceModel.updateShortPeriodTerms(forceModel.getParametersAllValues(field), initialState);

  663.             }
  664.             final Collection<FieldODEStepHandler<T>> stepHandlers = new ArrayList<>();
  665.             stepHandlers.add(spHandler);
  666.             final FieldODEIntegrator<T> integrator = getIntegrator();
  667.             final Collection<FieldODEStepHandler<T>> existing = integrator.getStepHandlers();
  668.             stepHandlers.addAll(existing);

  670.             integrator.clearStepHandlers();

  672.             // add back the existing handlers after the short periodics one
  673.             for (final FieldODEStepHandler<T> sp : stepHandlers) {
  674.                 integrator.addStepHandler(sp);
  675.             }
  676.         }
  677.     }

  679.     /** {@inheritDoc} */
  680.     @Override
  681.     protected void afterIntegration() {
  682.         // remove the special short periodics step handler if added before
  683.         if (isMeanOrbit() ==  PropagationType.OSCULATING) {
  684.             final List<FieldODEStepHandler<T>> preserved = new ArrayList<>();
  685.             final FieldODEIntegrator<T> integrator = getIntegrator();

  687.             // clear the list
  688.             integrator.clearStepHandlers();

  690.             // add back the step handlers that were important for the user
  691.             for (final FieldODEStepHandler<T> sp : preserved) {
  692.                 integrator.addStepHandler(sp);
  693.             }
  694.         }
  695.     }

  697.     /** Compute mean state from osculating state.
  698.      * <p>
  699.      * Compute in a DSST sense the mean state corresponding to the input osculating state.
  700.      * </p><p>
  701.      * The computing is done through a fixed-point iteration process.
  702.      * </p>
  703.      * @param osculating initial osculating state
  704.      * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
  705.      * like atmospheric drag, radiation pressure or specific user-defined models)
  706.      * @param forceModel force models
  707.      * @param epsilon convergence threshold for mean parameters conversion
  708.      * @param maxIterations maximum iterations for mean parameters conversion
  709.      * @param <T> type of the elements
  710.      * @return mean state
  711.      * @since 10.1
  712.      */
  713.     @SuppressWarnings("unchecked")
  714.     private static <T extends CalculusFieldElement<T>> FieldOrbit<T> computeMeanOrbit(final FieldSpacecraftState<T> osculating, final AttitudeProvider attitudeProvider, final Collection<DSSTForceModel> forceModel,
  715.                                                                                   final double epsilon, final int maxIterations) {

  717.         // zero
  718.         final T zero = osculating.getDate().getField().getZero();

  720.         // rough initialization of the mean parameters
  721.         FieldEquinoctialOrbit<T> meanOrbit = (FieldEquinoctialOrbit<T>) OrbitType.EQUINOCTIAL.convertType(osculating.getOrbit());

  723.         // threshold for each parameter
  724.         final T epsilonT   = zero.newInstance(epsilon);
  725.         final T thresholdA = epsilonT.multiply(FastMath.abs(meanOrbit.getA()).add(1.));
  726.         final T thresholdE = epsilonT.multiply(meanOrbit.getE().add(1.));
  727.         final T thresholdI = epsilonT.multiply(meanOrbit.getI().add(1.));
  728.         final T thresholdL = epsilonT.multiply(zero.getPi());

  730.         // ensure all Gaussian force models can rely on attitude
  731.         for (final DSSTForceModel force : forceModel) {
  732.             force.registerAttitudeProvider(attitudeProvider);
  733.         }

  735.         int i = 0;
  736.         while (i++ < maxIterations) {

  738.             final FieldSpacecraftState<T> meanState = new FieldSpacecraftState<>(meanOrbit, osculating.getAttitude(), osculating.getMass());

  740.             //Create the auxiliary object
  741.             final FieldAuxiliaryElements<T> aux = new FieldAuxiliaryElements<>(meanOrbit, I);

  743.             // Set the force models
  744.             final List<FieldShortPeriodTerms<T>> shortPeriodTerms = new ArrayList<>();
  745.             for (final DSSTForceModel force : forceModel) {
  746.                 shortPeriodTerms.addAll(force.initializeShortPeriodTerms(aux, PropagationType.OSCULATING,
  747.                                  force.getParameters(osculating.getDate().getField(), osculating.getDate())));
  748.                 force.updateShortPeriodTerms(force.getParametersAllValues(osculating.getDate().getField()), meanState);
  749.             }

  751.             // recompute the osculating parameters from the current mean parameters
  752.             final FieldEquinoctialOrbit<T> rebuilt = computeOsculatingOrbit(meanState, shortPeriodTerms);

  754.             // adapted parameters residuals
  755.             final T deltaA  = osculating.getA().subtract(rebuilt.getA());
  756.             final T deltaEx = osculating.getEquinoctialEx().subtract(rebuilt.getEquinoctialEx());
  757.             final T deltaEy = osculating.getEquinoctialEy().subtract(rebuilt.getEquinoctialEy());
  758.             final T deltaHx = osculating.getHx().subtract(rebuilt.getHx());
  759.             final T deltaHy = osculating.getHy().subtract(rebuilt.getHy());
  760.             final T deltaLv = MathUtils.normalizeAngle(osculating.getLv().subtract(rebuilt.getLv()), zero);

  762.             // check convergence
  763.             if (FastMath.abs(deltaA).getReal()  < thresholdA.getReal() &&
  764.                 FastMath.abs(deltaEx).getReal() < thresholdE.getReal() &&
  765.                 FastMath.abs(deltaEy).getReal() < thresholdE.getReal() &&
  766.                 FastMath.abs(deltaHx).getReal() < thresholdI.getReal() &&
  767.                 FastMath.abs(deltaHy).getReal() < thresholdI.getReal() &&
  768.                 FastMath.abs(deltaLv).getReal() < thresholdL.getReal()) {
  769.                 return meanOrbit;
  770.             }

  772.             // update mean parameters
  773.             meanOrbit = new FieldEquinoctialOrbit<>(meanOrbit.getA().add(deltaA),
  774.                                                     meanOrbit.getEquinoctialEx().add(deltaEx),
  775.                                                     meanOrbit.getEquinoctialEy().add(deltaEy),
  776.                                                     meanOrbit.getHx().add(deltaHx),
  777.                                                     meanOrbit.getHy().add(deltaHy),
  778.                                                     meanOrbit.getLv().add(deltaLv),
  779.                                                     PositionAngleType.TRUE, meanOrbit.getFrame(),
  780.                                                     meanOrbit.getDate(), meanOrbit.getMu());
  781.         }

  783.         throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_DSST_MEAN_PARAMETERS, i);
  784.     }

  786.     /** Compute osculating state from mean state.
  787.      * <p>
  788.      * Compute and add the short periodic variation to the mean {@link SpacecraftState}.
  789.      * </p>
  790.      * @param meanState initial mean state
  791.      * @param shortPeriodTerms short period terms
  792.      * @param <T> type of the elements
  793.      * @return osculating state
  794.      */
  795.     private static <T extends CalculusFieldElement<T>> FieldEquinoctialOrbit<T> computeOsculatingOrbit(final FieldSpacecraftState<T> meanState,
  796.                                                                                                    final List<FieldShortPeriodTerms<T>> shortPeriodTerms) {

  798.         final T[] mean = MathArrays.buildArray(meanState.getDate().getField(), 6);
  799.         final T[] meanDot = MathArrays.buildArray(meanState.getDate().getField(), 6);
  800.         OrbitType.EQUINOCTIAL.mapOrbitToArray(meanState.getOrbit(), PositionAngleType.MEAN, mean, meanDot);
  801.         final T[] y = mean.clone();
  802.         for (final FieldShortPeriodTerms<T> spt : shortPeriodTerms) {
  803.             final T[] shortPeriodic = spt.value(meanState.getOrbit());
  804.             for (int i = 0; i < shortPeriodic.length; i++) {
  805.                 y[i] = y[i].add(shortPeriodic[i]);
  806.             }
  807.         }
  808.         return (FieldEquinoctialOrbit<T>) OrbitType.EQUINOCTIAL.mapArrayToOrbit(y, meanDot,
  809.                                                                                 PositionAngleType.MEAN, meanState.getDate(),
  810.                                                                                 meanState.getMu(), meanState.getFrame());
  811.     }

  813.     /** {@inheritDoc} */
  814.     @Override
  815.     protected FieldSpacecraftState<T> getInitialIntegrationState() {
  816.         if (initialIsOsculating) {
  817.             // the initial state is an osculating state,
  818.             // it must be converted to mean state
  819.             return computeMeanState(getInitialState(), getAttitudeProvider(), forceModels);
  820.         } else {
  821.             // the initial state is already a mean state
  822.             return getInitialState();
  823.         }
  824.     }

  826.     /** {@inheritDoc}
  827.      * <p>
  828.      * Note that for DSST, orbit type is hardcoded to {@link OrbitType#EQUINOCTIAL}
  829.      * and position angle type is hardcoded to {@link PositionAngleType#MEAN}, so
  830.      * the corresponding parameters are ignored.
  831.      * </p>
  832.      */
  833.     @Override
  834.     protected FieldStateMapper<T> createMapper(final FieldAbsoluteDate<T> referenceDate, final T mu,
  835.                                                final OrbitType ignoredOrbitType, final PositionAngleType ignoredPositionAngleType,
  836.                                                final AttitudeProvider attitudeProvider, final Frame frame) {

  838.         // create a mapper with the common settings provided as arguments
  839.         final FieldMeanPlusShortPeriodicMapper newMapper =
  840.                 new FieldMeanPlusShortPeriodicMapper(referenceDate, mu, attitudeProvider, frame);

  842.         // copy the specific settings from the existing mapper
  843.         if (mapper != null) {
  844.             newMapper.setSatelliteRevolution(mapper.getSatelliteRevolution());
  845.             newMapper.setSelectedCoefficients(mapper.getSelectedCoefficients());
  846.             newMapper.setShortPeriodTerms(mapper.getShortPeriodTerms());
  847.         }

  849.         mapper = newMapper;
  850.         return mapper;

  852.     }

  854.     /** Internal mapper using mean parameters plus short periodic terms. */
  855.     private class FieldMeanPlusShortPeriodicMapper extends FieldStateMapper<T> {

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

  860.         /** Number of satellite revolutions in the averaging interval. */
  861.         private int                        satelliteRevolution;

  863.         /** Short period terms. */
  864.         private List<FieldShortPeriodTerms<T>>     shortPeriodTerms;

  866.         /** Simple constructor.
  867.          * @param referenceDate reference date
  868.          * @param mu central attraction coefficient (m³/s²)
  869.          * @param attitudeProvider attitude provider
  870.          * @param frame inertial frame
  871.          */
  872.         FieldMeanPlusShortPeriodicMapper(final FieldAbsoluteDate<T> referenceDate, final T mu,
  873.                                          final AttitudeProvider attitudeProvider, final Frame frame) {

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

  877.             this.selectedCoefficients = null;

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

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

  884.         }

  886.         /** {@inheritDoc} */
  887.         @Override
  888.         public FieldSpacecraftState<T> mapArrayToState(final FieldAbsoluteDate<T> date,
  889.                                                        final T[] y, final T[] yDot,
  890.                                                        final PropagationType type) {

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

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

  923.             final FieldOrbit<T> orbit       = OrbitType.EQUINOCTIAL.mapArrayToOrbit(elements, yDot, PositionAngleType.MEAN, date, getMu(), getFrame());
  924.             final FieldAttitude<T> attitude = getAttitudeProvider().getAttitude(orbit, date, getFrame());

  926.             if (coefficients == null) {
  927.                 return new FieldSpacecraftState<>(orbit, attitude, mass);
  928.             } else {
  929.                 return new FieldSpacecraftState<>(orbit, attitude, mass, coefficients);
  930.             }

  932.         }

  934.         /** {@inheritDoc} */
  935.         @Override
  936.         public void mapStateToArray(final FieldSpacecraftState<T> state, final T[] y, final T[] yDot) {

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

  941.         }

  943.         /** Set the number of satellite revolutions to use for converting osculating to mean elements.
  944.          *  <p>
  945.          *  By default, if the initial orbit is defined as osculating,
  946.          *  it will be averaged over 2 satellite revolutions.
  947.          *  This can be changed by using this method.
  948.          *  </p>
  949.          *  @param satelliteRevolution number of satellite revolutions to use for converting osculating to mean
  950.          *                             elements
  951.          */
  952.         public void setSatelliteRevolution(final int satelliteRevolution) {
  953.             this.satelliteRevolution = satelliteRevolution;
  954.         }

  956.         /** Get the number of satellite revolutions to use for converting osculating to mean elements.
  957.          *  @return number of satellite revolutions to use for converting osculating to mean elements
  958.          */
  959.         public int getSatelliteRevolution() {
  960.             return satelliteRevolution;
  961.         }

  963.         /** Set the selected short periodic coefficients that must be stored as additional states.
  964.          * @param selectedCoefficients short periodic coefficients that must be stored as additional states
  965.          * (null means no coefficients are selected, empty set means all coefficients are selected)
  966.          */
  967.         public void setSelectedCoefficients(final Set<String> selectedCoefficients) {
  968.             this.selectedCoefficients = selectedCoefficients;
  969.         }

  971.         /** Get the selected short periodic coefficients that must be stored as additional states.
  972.          * @return short periodic coefficients that must be stored as additional states
  973.          * (null means no coefficients are selected, empty set means all coefficients are selected)
  974.          */
  975.         public Set<String> getSelectedCoefficients() {
  976.             return selectedCoefficients;
  977.         }

  979.         /** Set the short period terms.
  980.          * @param shortPeriodTerms short period terms
  981.          * @since 7.1
  982.          */
  983.         public void setShortPeriodTerms(final List<FieldShortPeriodTerms<T>> shortPeriodTerms) {
  984.             this.shortPeriodTerms = shortPeriodTerms;
  985.         }

  987.         /** Get the short period terms.
  988.          * @return shortPeriodTerms short period terms
  989.          * @since 7.1
  990.          */
  991.         public List<FieldShortPeriodTerms<T>> getShortPeriodTerms() {
  992.             return shortPeriodTerms;
  993.         }

  995.     }

  997.     /** {@inheritDoc} */
  998.     @Override
  999.     protected MainStateEquations<T> getMainStateEquations(final FieldODEIntegrator<T> integrator) {
  1000.         return new Main(integrator);
  1001.     }

  1003.     /** Internal class for mean parameters integration. */
  1004.     private class Main implements MainStateEquations<T> {

  1006.         /** Derivatives array. */
  1007.         private final T[] yDot;

  1009.         /** Simple constructor.
  1010.          * @param integrator numerical integrator to use for propagation.
  1011.          */
  1012.         Main(final FieldODEIntegrator<T> integrator) {
  1013.             yDot = MathArrays.buildArray(field, 7);

  1015.             // Setup event detectors for each force model
  1016.             forceModels.forEach(dsstForceModel -> dsstForceModel.getFieldEventDetectors(field).
  1017.                                 forEach(eventDetector -> setUpEventDetector(integrator, eventDetector)));
  1018.         }

  1020.         /** {@inheritDoc} */
  1021.         @Override
  1022.         public void init(final FieldSpacecraftState<T> initialState, final FieldAbsoluteDate<T> target) {
  1023.             forceModels.forEach(fm -> fm.init(initialState, target));
  1024.         }

  1026.         /** {@inheritDoc} */
  1027.         @Override
  1028.         public T[] computeDerivatives(final FieldSpacecraftState<T> state) {

  1030.             final T zero = state.getDate().getField().getZero();
  1031.             Arrays.fill(yDot, zero);

  1033.             // compute common auxiliary elements
  1034.             final FieldAuxiliaryElements<T> auxiliaryElements = new FieldAuxiliaryElements<>(state.getOrbit(), I);

  1036.             // compute the contributions of all perturbing forces
  1037.             for (final DSSTForceModel forceModel : forceModels) {
  1038.                 final T[] daidt = elementRates(forceModel, state, auxiliaryElements, forceModel.getParametersAllValues(field));
  1039.                 for (int i = 0; i < daidt.length; i++) {
  1040.                     yDot[i] = yDot[i].add(daidt[i]);
  1041.                 }
  1042.             }

  1044.             return yDot.clone();
  1045.         }

  1047.         /** This method allows to compute the mean equinoctial elements rates da<sub>i</sub> / dt
  1048.          *  for a specific force model.
  1049.          *  @param forceModel force to take into account
  1050.          *  @param state current state
  1051.          *  @param auxiliaryElements auxiliary elements related to the current orbit
  1052.          *  @param parameters force model parameters (all span values for each parameters)
  1053.          *  the extract parameter method {@link #extractParameters(double[], AbsoluteDate)} is called in
  1054.          *  the method to select the right parameter.
  1055.          *  @return the mean equinoctial elements rates da<sub>i</sub> / dt
  1056.          */
  1057.         private T[] elementRates(final DSSTForceModel forceModel,
  1058.                                  final FieldSpacecraftState<T> state,
  1059.                                  final FieldAuxiliaryElements<T> auxiliaryElements,
  1060.                                  final T[] parameters) {
  1061.             return forceModel.getMeanElementRate(state, auxiliaryElements, parameters);
  1062.         }

  1064.     }

  1066.     /** Estimate tolerance vectors for an AdaptativeStepsizeIntegrator.
  1067.      *  <p>
  1068.      *  The errors are estimated from partial derivatives properties of orbits,
  1069.      *  starting from a scalar position error specified by the user.
  1070.      *  Considering the energy conservation equation V = sqrt(mu (2/r - 1/a)),
  1071.      *  we get at constant energy (i.e. on a Keplerian trajectory):
  1072.      *
  1073.      *  <pre>
  1074.      *  V r² |dV| = mu |dr|
  1075.      *  </pre>
  1076.      *
  1077.      *  <p> So we deduce a scalar velocity error consistent with the position error. From here, we apply
  1078.      *  orbits Jacobians matrices to get consistent errors on orbital parameters.
  1079.      *
  1080.      *  <p>
  1081.      *  The tolerances are only <em>orders of magnitude</em>, and integrator tolerances are only
  1082.      *  local estimates, not global ones. So some care must be taken when using these tolerances.
  1083.      *  Setting 1mm as a position error does NOT mean the tolerances will guarantee a 1mm error
  1084.      *  position after several orbits integration.
  1085.      *  </p>
  1086.      * @param <T> elements type
  1087.      * @param dP user specified position error (m)
  1088.      * @param orbit reference orbit
  1089.      * @return a two rows array, row 0 being the absolute tolerance error
  1090.      *                       and row 1 being the relative tolerance error
  1091.      */
  1092.     public static <T extends CalculusFieldElement<T>> double[][] tolerances(final T dP, final FieldOrbit<T> orbit) {
  1093.         return FieldNumericalPropagator.tolerances(dP, orbit, OrbitType.EQUINOCTIAL);
  1094.     }

  1096.     /** Estimate tolerance vectors for an AdaptativeStepsizeIntegrator.
  1097.      *  <p>
  1098.      *  The errors are estimated from partial derivatives properties of orbits,
  1099.      *  starting from scalar position and velocity errors specified by the user.
  1100.      *  <p>
  1101.      *  The tolerances are only <em>orders of magnitude</em>, and integrator tolerances are only
  1102.      *  local estimates, not global ones. So some care must be taken when using these tolerances.
  1103.      *  Setting 1mm as a position error does NOT mean the tolerances will guarantee a 1mm error
  1104.      *  position after several orbits integration.
  1105.      *  </p>
  1106.      *
  1107.      * @param <T> elements type
  1108.      * @param dP user specified position error (m)
  1109.      * @param dV user specified velocity error (m/s)
  1110.      * @param orbit reference orbit
  1111.      * @return a two rows array, row 0 being the absolute tolerance error
  1112.      *                       and row 1 being the relative tolerance error
  1113.      * @since 10.3
  1114.      */
  1115.     public static <T extends CalculusFieldElement<T>> double[][] tolerances(final T dP, final T dV,
  1116.                                                                         final FieldOrbit<T> orbit) {
  1117.         return FieldNumericalPropagator.tolerances(dP, dV, orbit, OrbitType.EQUINOCTIAL);
  1118.     }

  1120.     /** Step handler used to compute the parameters for the short periodic contributions.
  1121.      * @author Lucian Barbulescu
  1122.      */
  1123.     private class FieldShortPeriodicsHandler implements FieldODEStepHandler<T> {

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

  1128.         /** Constructor.
  1129.          * @param forceModels force models
  1130.          */
  1131.         FieldShortPeriodicsHandler(final List<DSSTForceModel> forceModels) {
  1132.             this.forceModels = forceModels;
  1133.         }

  1135.         /** {@inheritDoc} */
  1136.         @SuppressWarnings("unchecked")
  1137.         @Override
  1138.         public void handleStep(final FieldODEStateInterpolator<T> interpolator) {

  1140.             // Get the grid points to compute
  1141.             final T[] interpolationPoints =
  1142.                             interpolationgrid.getGridPoints(interpolator.getPreviousState().getTime(),
  1143.                                                             interpolator.getCurrentState().getTime());

  1145.             final FieldSpacecraftState<T>[] meanStates = new FieldSpacecraftState[interpolationPoints.length];
  1146.             for (int i = 0; i < interpolationPoints.length; ++i) {

  1148.                 // Build the mean state interpolated at grid point
  1149.                 final T time = interpolationPoints[i];
  1150.                 final FieldODEStateAndDerivative<T> sd = interpolator.getInterpolatedState(time);
  1151.                 meanStates[i] = mapper.mapArrayToState(time,
  1152.                                                        sd.getPrimaryState(),
  1153.                                                        sd.getPrimaryDerivative(),
  1154.                                                        PropagationType.MEAN);

  1156.             }

  1158.             // Compute short periodic coefficients for this step
  1159.             for (DSSTForceModel forceModel : forceModels) {
  1160.                 forceModel.updateShortPeriodTerms(forceModel.getParametersAllValues(field), meanStates);
  1161.             }

  1163.         }
  1164.     }

  1166. }
Created with JaCoCo 0.8.12.202403310830