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 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.PositionAngleType;
  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.integration.FieldAbstractIntegratedPropagator;
  53. import org.orekit.propagation.integration.FieldStateMapper;
  54. import org.orekit.propagation.numerical.FieldNumericalPropagator;
  55. import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
  56. import org.orekit.propagation.semianalytical.dsst.forces.DSSTNewtonianAttraction;
  57. import org.orekit.propagation.semianalytical.dsst.forces.FieldShortPeriodTerms;
  58. import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
  59. import org.orekit.propagation.semianalytical.dsst.utilities.FieldFixedNumberInterpolationGrid;
  60. import org.orekit.propagation.semianalytical.dsst.utilities.FieldInterpolationGrid;
  61. import org.orekit.propagation.semianalytical.dsst.utilities.FieldMaxGapInterpolationGrid;
  62. import org.orekit.time.AbsoluteDate;
  63. import org.orekit.time.FieldAbsoluteDate;
  64. import org.orekit.utils.FieldArrayDictionary;
  65. import org.orekit.utils.ParameterDriver;
  66. import org.orekit.utils.ParameterObserver;
  67. import org.orekit.utils.TimeSpanMap;

  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.         resetInitialState(initialState, stateType);
  313.     }

  315.     /** Reset the initial state.
  316.      *
  317.      *  @param state new initial state
  318.      */
  319.     @Override
  320.     public void resetInitialState(final FieldSpacecraftState<T> state) {
  321.         super.resetInitialState(state);
  322.         if (!hasNewtonianAttraction()) {
  323.             // use the state to define central attraction
  324.             setMu(state.getMu());
  325.         }
  326.         super.setStartDate(state.getDate());
  327.     }

  329.     /** {@inheritDoc}.
  330.      *
  331.      * <p>Change parameter {@link #initialIsOsculating()} accordingly
  332.      * @since 12.1.3
  333.      */
  334.     @Override
  335.     public void resetInitialState(final FieldSpacecraftState<T> state, final PropagationType stateType) {
  336.         // Reset initial state
  337.         resetInitialState(state);

  339.         // Change state of initial osculating, if needed
  340.         initialIsOsculating = stateType == PropagationType.OSCULATING;
  341.     }

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

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

  361.     /** Check if the initial state is provided in osculating elements.
  362.      * @return true if initial state is provided in osculating elements
  363.      */
  364.     public boolean initialIsOsculating() {
  365.         return initialIsOsculating;
  366.     }

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

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

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

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

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

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

  459.         force.registerAttitudeProvider(getAttitudeProvider());

  461.     }

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

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

  493.     /** Get propagation parameter type.
  494.      * @return orbit type used for propagation
  495.      */
  496.     public OrbitType getOrbitType() {
  497.         return super.getOrbitType();
  498.     }

  500.     /** Get propagation parameter type.
  501.      * @return angle type to use for propagation
  502.      */
  503.     public PositionAngleType getPositionAngleType() {
  504.         return super.getPositionAngleType();
  505.     }

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

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

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

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

  542.         return new FieldSpacecraftState<>(osculatingOrbit, mean.getAttitude(), mean.getMass(),
  543.                                           mean.getAdditionalStatesValues());

  545.     }

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

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

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

  617.    /** Get the number of satellite revolutions to use for converting osculating to mean elements.
  618.     *  @return number of satellite revolutions to use for converting osculating to mean elements
  619.     */
  620.     public int getSatelliteRevolution() {
  621.         return mapper.getSatelliteRevolution();
  622.     }

  624.    /** {@inheritDoc} */
  625.     @Override
  626.     public void setAttitudeProvider(final AttitudeProvider attitudeProvider) {
  627.         super.setAttitudeProvider(attitudeProvider);

  629.         //Register the attitude provider for each force model
  630.         for (final DSSTForceModel force : forceModels) {
  631.             force.registerAttitudeProvider(attitudeProvider);
  632.         }
  633.     }

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

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

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

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

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

  667.             }
  668.             final Collection<FieldODEStepHandler<T>> stepHandlers = new ArrayList<>();
  669.             stepHandlers.add(spHandler);
  670.             final FieldODEIntegrator<T> integrator = getIntegrator();
  671.             final Collection<FieldODEStepHandler<T>> existing = integrator.getStepHandlers();
  672.             stepHandlers.addAll(existing);

  674.             integrator.clearStepHandlers();

  676.             // add back the existing handlers after the short periodics one
  677.             for (final FieldODEStepHandler<T> sp : stepHandlers) {
  678.                 integrator.addStepHandler(sp);
  679.             }
  680.         }
  681.     }

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

  691.             // clear the list
  692.             integrator.clearStepHandlers();

  694.             // add back the step handlers that were important for the user
  695.             for (final FieldODEStepHandler<T> sp : preserved) {
  696.                 integrator.addStepHandler(sp);
  697.             }
  698.         }
  699.     }

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

  721.         // zero
  722.         final T zero = osculating.getDate().getField().getZero();

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

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

  734.         // ensure all Gaussian force models can rely on attitude
  735.         for (final DSSTForceModel force : forceModel) {
  736.             force.registerAttitudeProvider(attitudeProvider);
  737.         }

  739.         int i = 0;
  740.         while (i++ < maxIterations) {

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

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

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

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

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

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

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

  787.         throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_DSST_MEAN_PARAMETERS, i);
  788.     }

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

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

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

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

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

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

  853.         mapper = newMapper;
  854.         return mapper;

  856.     }

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

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

  864.         /** Number of satellite revolutions in the averaging interval. */
  865.         private int                        satelliteRevolution;

  867.         /** Short period terms. */
  868.         private List<FieldShortPeriodTerms<T>>     shortPeriodTerms;

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

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

  881.             this.selectedCoefficients = null;

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

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

  888.         }

  890.         /** {@inheritDoc} */
  891.         @Override
  892.         public FieldSpacecraftState<T> mapArrayToState(final FieldAbsoluteDate<T> date,
  893.                                                        final T[] y, final T[] yDot,
  894.                                                        final PropagationType type) {

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

  922.             final T mass = elements[6];
  923.             if (mass.getReal() <= 0.0) {
  924.                 throw new OrekitException(OrekitMessages.NOT_POSITIVE_SPACECRAFT_MASS, mass);
  925.             }

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

  930.             if (coefficients == null) {
  931.                 return new FieldSpacecraftState<>(orbit, attitude, mass);
  932.             } else {
  933.                 return new FieldSpacecraftState<>(orbit, attitude, mass, coefficients);
  934.             }

  936.         }

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

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

  945.         }

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

  960.         /** Get the number of satellite revolutions to use for converting osculating to mean elements.
  961.          *  @return number of satellite revolutions to use for converting osculating to mean elements
  962.          */
  963.         public int getSatelliteRevolution() {
  964.             return satelliteRevolution;
  965.         }

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

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

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

  991.         /** Get the short period terms.
  992.          * @return shortPeriodTerms short period terms
  993.          * @since 7.1
  994.          */
  995.         public List<FieldShortPeriodTerms<T>> getShortPeriodTerms() {
  996.             return shortPeriodTerms;
  997.         }

  999.     }

  1001.     /** {@inheritDoc} */
  1002.     @Override
  1003.     protected MainStateEquations<T> getMainStateEquations(final FieldODEIntegrator<T> integrator) {
  1004.         return new Main(integrator);
  1005.     }

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

  1010.         /** Derivatives array. */
  1011.         private final T[] yDot;

  1013.         /** Simple constructor.
  1014.          * @param integrator numerical integrator to use for propagation.
  1015.          */
  1016.         Main(final FieldODEIntegrator<T> integrator) {
  1017.             yDot = MathArrays.buildArray(field, 7);

  1019.             // Setup event detectors for each force model
  1020.             forceModels.forEach(dsstForceModel -> dsstForceModel.getFieldEventDetectors(field).
  1021.                                 forEach(eventDetector -> setUpEventDetector(integrator, eventDetector)));
  1022.         }

  1024.         /** {@inheritDoc} */
  1025.         @Override
  1026.         public void init(final FieldSpacecraftState<T> initialState, final FieldAbsoluteDate<T> target) {
  1027.             forceModels.forEach(fm -> fm.init(initialState, target));
  1028.         }

  1030.         /** {@inheritDoc} */
  1031.         @Override
  1032.         public T[] computeDerivatives(final FieldSpacecraftState<T> state) {

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

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

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

  1048.             return yDot.clone();
  1049.         }

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

  1068.     }

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

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

  1124.     /** Step handler used to compute the parameters for the short periodic contributions.
  1125.      * @author Lucian Barbulescu
  1126.      */
  1127.     private class FieldShortPeriodicsHandler implements FieldODEStepHandler<T> {

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

  1132.         /** Constructor.
  1133.          * @param forceModels force models
  1134.          */
  1135.         FieldShortPeriodicsHandler(final List<DSSTForceModel> forceModels) {
  1136.             this.forceModels = forceModels;
  1137.         }

  1139.         /** {@inheritDoc} */
  1140.         @SuppressWarnings("unchecked")
  1141.         @Override
  1142.         public void handleStep(final FieldODEStateInterpolator<T> interpolator) {

  1144.             // Get the grid points to compute
  1145.             final T[] interpolationPoints =
  1146.                             interpolationgrid.getGridPoints(interpolator.getPreviousState().getTime(),
  1147.                                                             interpolator.getCurrentState().getTime());

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

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

  1160.             }

  1162.             // Compute short periodic coefficients for this step
  1163.             for (DSSTForceModel forceModel : forceModels) {
  1164.                 forceModel.updateShortPeriodTerms(forceModel.getParametersAllValues(field), meanStates);
  1165.             }

  1167.         }
  1168.     }

  1170. }
Created with JaCoCo 0.8.12.202403310830