FieldDSSTPropagator.java

  1. /* Copyright 2002-2020 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.HashMap;
  24. import java.util.HashSet;
  25. import java.util.List;
  26. import java.util.Map;
  27. import java.util.Set;

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

  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 #setSlaveMode()},
  88.  * {@link #setMasterMode(RealFieldElement, org.orekit.propagation.sampling.FieldOrekitFixedStepHandler)},
  89.  * {@link #setMasterMode(org.orekit.propagation.sampling.FieldOrekitStepHandler)},
  90.  * {@link #setEphemerisMode()}, {@link #getGeneratedEphemeris()})</li>
  91.  * </ul>
  92.  * <p>
  93.  * From these configuration parameters, only the initial state is mandatory.
  94.  * The default propagation settings are in {@link OrbitType#EQUINOCTIAL equinoctial}
  95.  * parameters with {@link PositionAngle#TRUE true} longitude argument.
  96.  * The central attraction coefficient used to define the initial orbit will be used.
  97.  * However, specifying only the initial state would mean the propagator would use
  98.  * only Keplerian forces. In this case, the simpler
  99.  * {@link org.orekit.propagation.analytical.KeplerianPropagator KeplerianPropagator}
  100.  * class would be more effective.
  101.  * </p>
  102.  * <p>
  103.  * The underlying numerical integrator set up in the constructor may also have
  104.  * its own configuration parameters. Typical configuration parameters for adaptive
  105.  * stepsize integrators are the min, max and perhaps start step size as well as
  106.  * the absolute and/or relative errors thresholds.
  107.  * </p>
  108.  * <p>
  109.  * The state that is seen by the integrator is a simple six elements double array.
  110.  * These six elements are:
  111.  * <ul>
  112.  * <li>the {@link org.orekit.orbits.FieldEquinoctialOrbit equinoctial orbit parameters}
  113.  * (a, e<sub>x</sub>, e<sub>y</sub>, h<sub>x</sub>, h<sub>y</sub>, λ<sub>m</sub>)
  114.  * in meters and radians,</li>
  115.  * </ul>
  116.  *
  117.  * <p>By default, at the end of the propagation, the propagator resets the initial state to the final state,
  118.  * thus allowing a new propagation to be started from there without recomputing the part already performed.
  119.  * This behaviour can be chenged by calling {@link #setResetAtEnd(boolean)}.
  120.  * </p>
  121.  * <p>Beware the same instance cannot be used simultaneously by different threads, the class is <em>not</em>
  122.  * thread-safe.</p>
  123.  *
  124.  * @see FieldSpacecraftState
  125.  * @see DSSTForceModel
  126.  * @author Romain Di Costanzo
  127.  * @author Pascal Parraud
  128.  * @since 10.0
  129.  */
  130. public class FieldDSSTPropagator<T extends RealFieldElement<T>> extends FieldAbstractIntegratedPropagator<T>  {

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

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

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

  153.     /** Default value for maxIterations. */
  154.     private static final int MAX_ITERATIONS_DEFAULT = 200;

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

  159.     /** Field used by this class.*/
  160.     private final Field<T> field;

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

  165.     /** State mapper holding the force models. */
  166.     private FieldMeanPlusShortPeriodicMapper mapper;

  168.     /** Generator for the interpolation grid. */
  169.     private FieldInterpolationGrid<T> interpolationgrid;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  450.         force.registerAttitudeProvider(getAttitudeProvider());

  452.     }

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

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

  484.     /** Get propagation parameter type.
  485.      * @return orbit type used for propagation
  486.      */
  487.     public OrbitType getOrbitType() {
  488.         return super.getOrbitType();
  489.     }

  491.     /** Get propagation parameter type.
  492.      * @return angle type to use for propagation
  493.      */
  494.     public PositionAngle getPositionAngleType() {
  495.         return super.getPositionAngleType();
  496.     }

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

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

  522.         // Set the force models
  523.         final List<FieldShortPeriodTerms<T>> shortPeriodTerms = new ArrayList<FieldShortPeriodTerms<T>>();
  524.         for (final DSSTForceModel force : forces) {
  525.             force.registerAttitudeProvider(attitudeProvider);
  526.             shortPeriodTerms.addAll(force.initialize(aux, PropagationType.OSCULATING, force.getParameters(field)));
  527.             force.updateShortPeriodTerms(force.getParameters(field), mean);
  528.         }

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

  532.         return new FieldSpacecraftState<>(osculatingOrbit, mean.getAttitude(), mean.getMass(),
  533.                                           mean.getAdditionalStates());

  535.     }

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

  561.     /** Conversion from osculating to mean orbit.
  562.      * <p>
  563.      * Compute mean state <b>in a DSST sense</b>, corresponding to the
  564.      * osculating SpacecraftState in input, and according to the Force models
  565.      * taken into account.
  566.      * </p><p>
  567.      * Since the osculating state is obtained with the computation of
  568.      * short-periodic variation of each force model, the resulting output will
  569.      * depend on the force models parameterized in input.
  570.      * </p><p>
  571.      * The computation is done through a fixed-point iteration process.
  572.      * </p>
  573.      * @param osculating Osculating state to convert
  574.      * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
  575.      * like atmospheric drag, radiation pressure or specific user-defined models)
  576.      * @param forceModel Forces to take into account
  577.      * @param epsilon convergence threshold for mean parameters conversion
  578.      * @param maxIterations maximum iterations for mean parameters conversion
  579.      * @return mean state in a DSST sense
  580.      * @since 10.1
  581.      */
  582.     public 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(), osculating.getAdditionalStates());
  589.     }

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

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

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

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

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

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

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

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

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

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

  661.             integrator.clearStepHandlers();

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

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

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

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

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

  707.         // zero
  708.         final T zero = field.getZero();

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

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

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

  725.         int i = 0;
  726.         while (i++ < maxIterations) {

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

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

  733.             // Set the force models
  734.             final List<FieldShortPeriodTerms<T>> shortPeriodTerms = new ArrayList<FieldShortPeriodTerms<T>>();
  735.             for (final DSSTForceModel force : forceModel) {
  736.                 shortPeriodTerms.addAll(force.initialize(aux, PropagationType.OSCULATING, force.getParameters(field)));
  737.                 force.updateShortPeriodTerms(force.getParameters(field), meanState);
  738.             }

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

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

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

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

  772.         throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_DSST_MEAN_PARAMETERS, i);
  773.     }

  775.     /** Compute osculating state from mean state.
  776.      * <p>
  777.      * Compute and add the short periodic variation to the mean {@link SpacecraftState}.
  778.      * </p>
  779.      * @param meanState initial mean state
  780.      * @param shortPeriodTerms short period terms
  781.      * @return osculating state
  782.      */
  783.     private FieldEquinoctialOrbit<T> computeOsculatingOrbit(final FieldSpacecraftState<T> meanState,
  784.                                                             final List<FieldShortPeriodTerms<T>> shortPeriodTerms) {

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

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

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

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

  830.         // copy the specific settings from the existing mapper
  831.         if (mapper != null) {
  832.             newMapper.setSatelliteRevolution(mapper.getSatelliteRevolution());
  833.             newMapper.setSelectedCoefficients(mapper.getSelectedCoefficients());
  834.             newMapper.setShortPeriodTerms(mapper.getShortPeriodTerms());
  835.         }

  837.         mapper = newMapper;
  838.         return mapper;

  840.     }

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

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

  848.         /** Number of satellite revolutions in the averaging interval. */
  849.         private int                        satelliteRevolution;

  851.         /** Short period terms. */
  852.         private List<FieldShortPeriodTerms<T>>     shortPeriodTerms;

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

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

  865.             this.selectedCoefficients = null;

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

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

  872.         }

  874.         /** {@inheritDoc} */
  875.         @Override
  876.         public FieldSpacecraftState<T> mapArrayToState(final FieldAbsoluteDate<T> date,
  877.                                                        final T[] y, final T[] yDot,
  878.                                                        final PropagationType type) {

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

  906.             final T mass = elements[6];
  907.             if (mass.getReal() <= 0.0) {
  908.                 throw new OrekitException(OrekitMessages.SPACECRAFT_MASS_BECOMES_NEGATIVE, mass);
  909.             }

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

  914.             if (coefficients == null) {
  915.                 return new FieldSpacecraftState<>(orbit, attitude, mass);
  916.             } else {
  917.                 return new FieldSpacecraftState<>(orbit, attitude, mass, coefficients);
  918.             }

  920.         }

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

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

  929.         }

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

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

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

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

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

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

  983.     }

  985.     /** {@inheritDoc} */
  986.     @Override
  987.     protected MainStateEquations<T> getMainStateEquations(final FieldODEIntegrator<T> integrator) {
  988.         return new Main(integrator);
  989.     }

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

  994.         /** Derivatives array. */
  995.         private final T[] yDot;

  997.         /** Simple constructor.
  998.          * @param integrator numerical integrator to use for propagation.
  999.          */
  1000.         Main(final FieldODEIntegrator<T> integrator) {
  1001.             yDot = MathArrays.buildArray(field, 7);

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

  1012.         }

  1014.         /** {@inheritDoc} */
  1015.         @Override
  1016.         public void init(final FieldSpacecraftState<T> initialState, final FieldAbsoluteDate<T> target) {
  1017.         }

  1019.         /** {@inheritDoc} */
  1020.         @Override
  1021.         public T[] computeDerivatives(final FieldSpacecraftState<T> state) {

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

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

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

  1037.             return yDot.clone();
  1038.         }

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

  1055.     }

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

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

  1111.     /** Step handler used to compute the parameters for the short periodic contributions.
  1112.      * @author Lucian Barbulescu
  1113.      */
  1114.     private class FieldShortPeriodicsHandler implements FieldODEStepHandler<T> {

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

  1119.         /** Constructor.
  1120.          * @param forceModels force models
  1121.          */
  1122.         FieldShortPeriodicsHandler(final List<DSSTForceModel> forceModels) {
  1123.             this.forceModels = forceModels;
  1124.         }

  1126.         /** {@inheritDoc} */
  1127.         @SuppressWarnings("unchecked")
  1128.         @Override
  1129.         public void handleStep(final FieldODEStateInterpolator<T> interpolator, final boolean isLast) {

  1131.             // Get the grid points to compute
  1132.             final T[] interpolationPoints =
  1133.                             interpolationgrid.getGridPoints(interpolator.getPreviousState().getTime(),
  1134.                                                             interpolator.getCurrentState().getTime());

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

  1139.                 // Build the mean state interpolated at grid point
  1140.                 final T time = interpolationPoints[i];
  1141.                 final FieldODEStateAndDerivative<T> sd = interpolator.getInterpolatedState(time);
  1142.                 meanStates[i] = mapper.mapArrayToState(time,
  1143.                                                        sd.getPrimaryState(),
  1144.                                                        sd.getPrimaryDerivative(),
  1145.                                                        PropagationType.MEAN);

  1147.             }

  1149.             // Compute short periodic coefficients for this step
  1150.             for (DSSTForceModel forceModel : forceModels) {
  1151.                 forceModel.updateShortPeriodTerms(forceModel.getParameters(field), meanStates);
  1152.             }

  1154.         }
  1155.     }

  1157. }
Created with JaCoCo 0.8.6.202009150832