FieldDSSTPropagator.java

  1. /* Copyright 2002-2022 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 java.util.ArrayList;
  20. import java.util.Arrays;
  21. import java.util.Collection;
  22. import java.util.Collections;
  23. import java.util.HashSet;
  24. import java.util.List;
  25. import java.util.Set;

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

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

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

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

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

  149.     /** Default value for maxIterations. */
  150.     private static final int MAX_ITERATIONS_DEFAULT = 200;

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

  155.     /** Field used by this class.*/
  156.     private final Field<T> field;

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

  161.     /** State mapper holding the force models. */
  162.     private FieldMeanPlusShortPeriodicMapper mapper;

  164.     /** Generator for the interpolation grid. */
  165.     private FieldInterpolationGrid<T> interpolationgrid;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  413.             try {
  414.                 // ensure we are notified of any mu change
  415.                 force.getParametersDrivers().get(0).addObserver(new ParameterObserver() {
  416.                     /** {@inheritDoc} */
  417.                     @Override
  418.                     public void valueChanged(final double previousValue, final ParameterDriver driver) {
  419.                         superSetMu(field.getZero().add(driver.getValue()));
  420.                     }
  421.                 });
  422.             } catch (OrekitException oe) {
  423.                 // this should never happen
  424.                 throw new OrekitInternalError(oe);
  425.             }

  427.             if (hasNewtonianAttraction()) {
  428.                 // there is already a central attraction model, replace it
  429.                 forceModels.set(forceModels.size() - 1, force);
  430.             } else {
  431.                 // there are no central attraction model yet, add it at the end of the list
  432.                 forceModels.add(force);
  433.             }
  434.         } else {
  435.             // we want to add a perturbing force model
  436.             if (hasNewtonianAttraction()) {
  437.                 // insert the new force model before Newtonian attraction,
  438.                 // which should always be the last one in the list
  439.                 forceModels.add(forceModels.size() - 1, force);
  440.             } else {
  441.                 // we only have perturbing force models up to now, just append at the end of the list
  442.                 forceModels.add(force);
  443.             }
  444.         }

  446.         force.registerAttitudeProvider(getAttitudeProvider());

  448.     }

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

  470.     /** Get all the force models, perturbing forces and Newtonian attraction included.
  471.      * @return list of perturbing force models, with Newtonian attraction being the
  472.      * last one
  473.      * @see #addForceModel(DSSTForceModel)
  474.      * @see #setMu(CalculusFieldElement)
  475.      */
  476.     public List<DSSTForceModel> getAllForceModels() {
  477.         return Collections.unmodifiableList(forceModels);
  478.     }

  480.     /** Get propagation parameter type.
  481.      * @return orbit type used for propagation
  482.      */
  483.     public OrbitType getOrbitType() {
  484.         return super.getOrbitType();
  485.     }

  487.     /** Get propagation parameter type.
  488.      * @return angle type to use for propagation
  489.      */
  490.     public PositionAngle getPositionAngleType() {
  491.         return super.getPositionAngleType();
  492.     }

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

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

  519.         // Set the force models
  520.         final List<FieldShortPeriodTerms<T>> shortPeriodTerms = new ArrayList<FieldShortPeriodTerms<T>>();
  521.         for (final DSSTForceModel force : forces) {
  522.             final T[] parameters = force.getParameters(mean.getDate().getField());
  523.             force.registerAttitudeProvider(attitudeProvider);
  524.             shortPeriodTerms.addAll(force.initializeShortPeriodTerms(aux, PropagationType.OSCULATING, parameters));
  525.             force.updateShortPeriodTerms(parameters, mean);
  526.         }

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

  530.         return new FieldSpacecraftState<>(osculatingOrbit, mean.getAttitude(), mean.getMass(),
  531.                                           mean.getAdditionalStatesValues());

  533.     }

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

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

  592.     /** Override the default value of the parameter.
  593.     *  <p>
  594.     *  By default, if the initial orbit is defined as osculating,
  595.     *  it will be averaged over 2 satellite revolutions.
  596.     *  This can be changed by using this method.
  597.     *  </p>
  598.     *  @param satelliteRevolution number of satellite revolutions to use for converting osculating to mean
  599.     *                             elements
  600.     */
  601.     public void setSatelliteRevolution(final int satelliteRevolution) {
  602.         mapper.setSatelliteRevolution(satelliteRevolution);
  603.     }

  605.    /** Get the number of satellite revolutions to use for converting osculating to mean elements.
  606.     *  @return number of satellite revolutions to use for converting osculating to mean elements
  607.     */
  608.     public int getSatelliteRevolution() {
  609.         return mapper.getSatelliteRevolution();
  610.     }

  612.    /** {@inheritDoc} */
  613.     @Override
  614.     public void setAttitudeProvider(final AttitudeProvider attitudeProvider) {
  615.         super.setAttitudeProvider(attitudeProvider);

  617.         //Register the attitude provider for each force model
  618.         for (final DSSTForceModel force : forceModels) {
  619.             force.registerAttitudeProvider(attitudeProvider);
  620.         }
  621.     }

  623.     /** Method called just before integration.
  624.      * <p>
  625.      * The default implementation does nothing, it may be specialized in subclasses.
  626.      * </p>
  627.      * @param initialState initial state
  628.      * @param tEnd target date at which state should be propagated
  629.      */
  630.     @SuppressWarnings("unchecked")
  631.     @Override
  632.     protected void beforeIntegration(final FieldSpacecraftState<T> initialState,
  633.                                      final FieldAbsoluteDate<T> tEnd) {

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

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

  641.         // initialize all perturbing forces
  642.         final List<FieldShortPeriodTerms<T>> shortPeriodTerms = new ArrayList<FieldShortPeriodTerms<T>>();
  643.         for (final DSSTForceModel force : forceModels) {
  644.             shortPeriodTerms.addAll(force.initializeShortPeriodTerms(aux, type, force.getParameters(field)));
  645.         }
  646.         mapper.setShortPeriodTerms(shortPeriodTerms);

  648.         // if required, insert the special short periodics step handler
  649.         if (type == PropagationType.OSCULATING) {
  650.             final FieldShortPeriodicsHandler spHandler = new FieldShortPeriodicsHandler(forceModels);
  651.             // Compute short periodic coefficients for this point
  652.             for (DSSTForceModel forceModel : forceModels) {
  653.                 forceModel.updateShortPeriodTerms(forceModel.getParameters(field), initialState);

  655.             }
  656.             final Collection<FieldODEStepHandler<T>> stepHandlers = new ArrayList<FieldODEStepHandler<T>>();
  657.             stepHandlers.add(spHandler);
  658.             final FieldODEIntegrator<T> integrator = getIntegrator();
  659.             final Collection<FieldODEStepHandler<T>> existing = integrator.getStepHandlers();
  660.             stepHandlers.addAll(existing);

  662.             integrator.clearStepHandlers();

  664.             // add back the existing handlers after the short periodics one
  665.             for (final FieldODEStepHandler<T> sp : stepHandlers) {
  666.                 integrator.addStepHandler(sp);
  667.             }
  668.         }
  669.     }

  671.     /** {@inheritDoc} */
  672.     @Override
  673.     protected void afterIntegration() {
  674.         // remove the special short periodics step handler if added before
  675.         if (isMeanOrbit() ==  PropagationType.OSCULATING) {
  676.             final List<FieldODEStepHandler<T>> preserved = new ArrayList<FieldODEStepHandler<T>>();
  677.             final FieldODEIntegrator<T> integrator = getIntegrator();

  679.             // clear the list
  680.             integrator.clearStepHandlers();

  682.             // add back the step handlers that were important for the user
  683.             for (final FieldODEStepHandler<T> sp : preserved) {
  684.                 integrator.addStepHandler(sp);
  685.             }
  686.         }
  687.     }

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

  709.         // zero
  710.         final T zero = osculating.getDate().getField().getZero();

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

  715.         // threshold for each parameter
  716.         final T epsilonT   = zero.add(epsilon);
  717.         final T thresholdA = epsilonT.multiply(FastMath.abs(meanOrbit.getA()).add(1.));
  718.         final T thresholdE = epsilonT.multiply(meanOrbit.getE().add(1.));
  719.         final T thresholdI = epsilonT.multiply(meanOrbit.getI().add(1.));
  720.         final T thresholdL = epsilonT.multiply(zero.getPi());

  722.         // ensure all Gaussian force models can rely on attitude
  723.         for (final DSSTForceModel force : forceModel) {
  724.             force.registerAttitudeProvider(attitudeProvider);
  725.         }

  727.         int i = 0;
  728.         while (i++ < maxIterations) {

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

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

  735.             // Set the force models
  736.             final List<FieldShortPeriodTerms<T>> shortPeriodTerms = new ArrayList<FieldShortPeriodTerms<T>>();
  737.             for (final DSSTForceModel force : forceModel) {
  738.                 final T[] parameters = force.getParameters(osculating.getDate().getField());
  739.                 shortPeriodTerms.addAll(force.initializeShortPeriodTerms(aux, PropagationType.OSCULATING, parameters));
  740.                 force.updateShortPeriodTerms(parameters, meanState);
  741.             }

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

  746.             // adapted parameters residuals
  747.             final T deltaA  = osculating.getA().subtract(rebuilt.getA());
  748.             final T deltaEx = osculating.getEquinoctialEx().subtract(rebuilt.getEquinoctialEx());
  749.             final T deltaEy = osculating.getEquinoctialEy().subtract(rebuilt.getEquinoctialEy());
  750.             final T deltaHx = osculating.getHx().subtract(rebuilt.getHx());
  751.             final T deltaHy = osculating.getHy().subtract(rebuilt.getHy());
  752.             final T deltaLv = MathUtils.normalizeAngle(osculating.getLv().subtract(rebuilt.getLv()), zero);

  754.             // check convergence
  755.             if (FastMath.abs(deltaA).getReal()  < thresholdA.getReal() &&
  756.                 FastMath.abs(deltaEx).getReal() < thresholdE.getReal() &&
  757.                 FastMath.abs(deltaEy).getReal() < thresholdE.getReal() &&
  758.                 FastMath.abs(deltaHx).getReal() < thresholdI.getReal() &&
  759.                 FastMath.abs(deltaHy).getReal() < thresholdI.getReal() &&
  760.                 FastMath.abs(deltaLv).getReal() < thresholdL.getReal()) {
  761.                 return meanOrbit;
  762.             }

  764.             // update mean parameters
  765.             meanOrbit = new FieldEquinoctialOrbit<>(meanOrbit.getA().add(deltaA),
  766.                                                     meanOrbit.getEquinoctialEx().add(deltaEx),
  767.                                                     meanOrbit.getEquinoctialEy().add(deltaEy),
  768.                                                     meanOrbit.getHx().add(deltaHx),
  769.                                                     meanOrbit.getHy().add(deltaHy),
  770.                                                     meanOrbit.getLv().add(deltaLv),
  771.                                                     PositionAngle.TRUE, meanOrbit.getFrame(),
  772.                                                     meanOrbit.getDate(), meanOrbit.getMu());
  773.         }

  775.         throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_DSST_MEAN_PARAMETERS, i);
  776.     }

  778.     /** Compute osculating state from mean state.
  779.      * <p>
  780.      * Compute and add the short periodic variation to the mean {@link SpacecraftState}.
  781.      * </p>
  782.      * @param meanState initial mean state
  783.      * @param shortPeriodTerms short period terms
  784.      * @param <T> type of the elements
  785.      * @return osculating state
  786.      */
  787.     private static <T extends CalculusFieldElement<T>> FieldEquinoctialOrbit<T> computeOsculatingOrbit(final FieldSpacecraftState<T> meanState,
  788.                                                                                                    final List<FieldShortPeriodTerms<T>> shortPeriodTerms) {

  790.         final T[] mean = MathArrays.buildArray(meanState.getDate().getField(), 6);
  791.         final T[] meanDot = MathArrays.buildArray(meanState.getDate().getField(), 6);
  792.         OrbitType.EQUINOCTIAL.mapOrbitToArray(meanState.getOrbit(), PositionAngle.MEAN, mean, meanDot);
  793.         final T[] y = mean.clone();
  794.         for (final FieldShortPeriodTerms<T> spt : shortPeriodTerms) {
  795.             final T[] shortPeriodic = spt.value(meanState.getOrbit());
  796.             for (int i = 0; i < shortPeriodic.length; i++) {
  797.                 y[i] = y[i].add(shortPeriodic[i]);
  798.             }
  799.         }
  800.         return (FieldEquinoctialOrbit<T>) OrbitType.EQUINOCTIAL.mapArrayToOrbit(y, meanDot,
  801.                                                                                 PositionAngle.MEAN, meanState.getDate(),
  802.                                                                                 meanState.getMu(), meanState.getFrame());
  803.     }

  805.     /** {@inheritDoc} */
  806.     @Override
  807.     protected FieldSpacecraftState<T> getInitialIntegrationState() {
  808.         if (initialIsOsculating) {
  809.             // the initial state is an osculating state,
  810.             // it must be converted to mean state
  811.             return computeMeanState(getInitialState(), getAttitudeProvider(), forceModels);
  812.         } else {
  813.             // the initial state is already a mean state
  814.             return getInitialState();
  815.         }
  816.     }

  818.     /** {@inheritDoc}
  819.      * <p>
  820.      * Note that for DSST, orbit type is hardcoded to {@link OrbitType#EQUINOCTIAL}
  821.      * and position angle type is hardcoded to {@link PositionAngle#MEAN}, so
  822.      * the corresponding parameters are ignored.
  823.      * </p>
  824.      */
  825.     @Override
  826.     protected FieldStateMapper<T> createMapper(final FieldAbsoluteDate<T> referenceDate, final T mu,
  827.                                                final OrbitType ignoredOrbitType, final PositionAngle ignoredPositionAngleType,
  828.                                                final AttitudeProvider attitudeProvider, final Frame frame) {

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

  834.         // copy the specific settings from the existing mapper
  835.         if (mapper != null) {
  836.             newMapper.setSatelliteRevolution(mapper.getSatelliteRevolution());
  837.             newMapper.setSelectedCoefficients(mapper.getSelectedCoefficients());
  838.             newMapper.setShortPeriodTerms(mapper.getShortPeriodTerms());
  839.         }

  841.         mapper = newMapper;
  842.         return mapper;

  844.     }

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

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

  852.         /** Number of satellite revolutions in the averaging interval. */
  853.         private int                        satelliteRevolution;

  855.         /** Short period terms. */
  856.         private List<FieldShortPeriodTerms<T>>     shortPeriodTerms;

  858.         /** Simple constructor.
  859.          * @param referenceDate reference date
  860.          * @param mu central attraction coefficient (m³/s²)
  861.          * @param attitudeProvider attitude provider
  862.          * @param frame inertial frame
  863.          */
  864.         FieldMeanPlusShortPeriodicMapper(final FieldAbsoluteDate<T> referenceDate, final T mu,
  865.                                          final AttitudeProvider attitudeProvider, final Frame frame) {

  867.             super(referenceDate, mu, OrbitType.EQUINOCTIAL, PositionAngle.MEAN, attitudeProvider, frame);

  869.             this.selectedCoefficients = null;

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

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

  876.         }

  878.         /** {@inheritDoc} */
  879.         @Override
  880.         public FieldSpacecraftState<T> mapArrayToState(final FieldAbsoluteDate<T> date,
  881.                                                        final T[] y, final T[] yDot,
  882.                                                        final PropagationType type) {

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

  910.             final T mass = elements[6];
  911.             if (mass.getReal() <= 0.0) {
  912.                 throw new OrekitException(OrekitMessages.SPACECRAFT_MASS_BECOMES_NEGATIVE, mass);
  913.             }

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

  918.             if (coefficients == null) {
  919.                 return new FieldSpacecraftState<>(orbit, attitude, mass);
  920.             } else {
  921.                 return new FieldSpacecraftState<>(orbit, attitude, mass, coefficients);
  922.             }

  924.         }

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

  930.             OrbitType.EQUINOCTIAL.mapOrbitToArray(state.getOrbit(), PositionAngle.MEAN, y, yDot);
  931.             y[6] = state.getMass();

  933.         }

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

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

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

  963.         /** Get the selected short periodic coefficients that must be stored as additional states.
  964.          * @return 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 Set<String> getSelectedCoefficients() {
  968.             return selectedCoefficients;
  969.         }

  971.         /** Set the short period terms.
  972.          * @param shortPeriodTerms short period terms
  973.          * @since 7.1
  974.          */
  975.         public void setShortPeriodTerms(final List<FieldShortPeriodTerms<T>> shortPeriodTerms) {
  976.             this.shortPeriodTerms = shortPeriodTerms;
  977.         }

  979.         /** Get the short period terms.
  980.          * @return shortPeriodTerms short period terms
  981.          * @since 7.1
  982.          */
  983.         public List<FieldShortPeriodTerms<T>> getShortPeriodTerms() {
  984.             return shortPeriodTerms;
  985.         }

  987.     }

  989.     /** {@inheritDoc} */
  990.     @Override
  991.     protected MainStateEquations<T> getMainStateEquations(final FieldODEIntegrator<T> integrator) {
  992.         return new Main(integrator);
  993.     }

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

  998.         /** Derivatives array. */
  999.         private final T[] yDot;

  1001.         /** Simple constructor.
  1002.          * @param integrator numerical integrator to use for propagation.
  1003.          */
  1004.         Main(final FieldODEIntegrator<T> integrator) {
  1005.             yDot = MathArrays.buildArray(field, 7);

  1007.             for (final DSSTForceModel forceModel : forceModels) {
  1008.                 final FieldEventDetector<T>[] modelDetectors = forceModel.getFieldEventsDetectors(field);
  1009.                 if (modelDetectors != null) {
  1010.                     for (final FieldEventDetector<T> detector : modelDetectors) {
  1011.                         setUpEventDetector(integrator, detector);
  1012.                     }
  1013.                 }
  1014.             }

  1016.         }

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

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

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

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

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

  1042.             return yDot.clone();
  1043.         }

  1045.         /** This method allows to compute the mean equinoctial elements rates da<sub>i</sub> / dt
  1046.          *  for a specific force model.
  1047.          *  @param forceModel force to take into account
  1048.          *  @param state current state
  1049.          *  @param auxiliaryElements auxiliary elements related to the current orbit
  1050.          *  @param parameters force model parameters
  1051.          *  @return the mean equinoctial elements rates da<sub>i</sub> / dt
  1052.          */
  1053.         private T[] elementRates(final DSSTForceModel forceModel,
  1054.                                  final FieldSpacecraftState<T> state,
  1055.                                  final FieldAuxiliaryElements<T> auxiliaryElements,
  1056.                                  final T[] parameters) {
  1057.             return forceModel.getMeanElementRate(state, auxiliaryElements, parameters);
  1058.         }

  1060.     }

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

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

  1116.     /** Step handler used to compute the parameters for the short periodic contributions.
  1117.      * @author Lucian Barbulescu
  1118.      */
  1119.     private class FieldShortPeriodicsHandler implements FieldODEStepHandler<T> {

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

  1124.         /** Constructor.
  1125.          * @param forceModels force models
  1126.          */
  1127.         FieldShortPeriodicsHandler(final List<DSSTForceModel> forceModels) {
  1128.             this.forceModels = forceModels;
  1129.         }

  1131.         /** {@inheritDoc} */
  1132.         @SuppressWarnings("unchecked")
  1133.         @Override
  1134.         public void handleStep(final FieldODEStateInterpolator<T> interpolator) {

  1136.             // Get the grid points to compute
  1137.             final T[] interpolationPoints =
  1138.                             interpolationgrid.getGridPoints(interpolator.getPreviousState().getTime(),
  1139.                                                             interpolator.getCurrentState().getTime());

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

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

  1152.             }

  1154.             // Compute short periodic coefficients for this step
  1155.             for (DSSTForceModel forceModel : forceModels) {
  1156.                 forceModel.updateShortPeriodTerms(forceModel.getParameters(field), meanStates);
  1157.             }

  1159.         }
  1160.     }

  1162. }
Created with JaCoCo 0.8.7.202105040129