FieldNumericalPropagator.java

  1. /* Copyright 2002-2018 CS Systèmes d'Information
  2.  * Licensed to CS Systèmes d'Information (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.numerical;

  19. import java.util.ArrayList;
  20. import java.util.Arrays;
  21. import java.util.Collections;
  22. import java.util.List;

  24. import org.hipparchus.Field;
  25. import org.hipparchus.RealFieldElement;
  26. import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
  27. import org.hipparchus.ode.FieldODEIntegrator;
  28. import org.hipparchus.util.MathArrays;
  29. import org.orekit.attitudes.AttitudeProvider;
  30. import org.orekit.attitudes.FieldAttitude;
  31. import org.orekit.errors.OrekitException;
  32. import org.orekit.errors.OrekitIllegalArgumentException;
  33. import org.orekit.errors.OrekitInternalError;
  34. import org.orekit.errors.OrekitMessages;
  35. import org.orekit.forces.ForceModel;
  36. import org.orekit.forces.gravity.NewtonianAttraction;
  37. import org.orekit.frames.Frame;
  38. import org.orekit.orbits.FieldOrbit;
  39. import org.orekit.orbits.OrbitType;
  40. import org.orekit.orbits.PositionAngle;
  41. import org.orekit.propagation.FieldSpacecraftState;
  42. import org.orekit.propagation.SpacecraftState;
  43. import org.orekit.propagation.events.FieldEventDetector;
  44. import org.orekit.propagation.integration.FieldAbstractIntegratedPropagator;
  45. import org.orekit.propagation.integration.FieldStateMapper;
  46. import org.orekit.time.AbsoluteDate;
  47. import org.orekit.time.FieldAbsoluteDate;
  48. import org.orekit.utils.FieldPVCoordinates;
  49. import org.orekit.utils.ParameterDriver;
  50. import org.orekit.utils.ParameterObserver;
  51. import org.orekit.utils.TimeStampedFieldPVCoordinates;

  53. /** This class propagates {@link org.orekit.orbits.FieldOrbit orbits} using
  54.  * numerical integration.
  55.  * <p>Numerical propagation is much more accurate than analytical propagation
  56.  * like for example {@link org.orekit.propagation.analytical.KeplerianPropagator
  57.  * Keplerian} or {@link org.orekit.propagation.analytical.EcksteinHechlerPropagator
  58.  * Eckstein-Hechler}, but requires a few more steps to set up to be used properly.
  59.  * Whereas analytical propagators are configured only thanks to their various
  60.  * constructors and can be used immediately after construction, numerical propagators
  61.  * configuration involve setting several parameters between construction time
  62.  * and propagation time.</p>
  63.  * <p>The configuration parameters that can be set are:</p>
  64.  * <ul>
  65.  *   <li>the initial spacecraft state ({@link #setInitialState(FieldSpacecraftState)})</li>
  66.  *   <li>the central attraction coefficient ({@link #setMu(double)})</li>
  67.  *   <li>the various force models ({@link #addForceModel(ForceModel)},
  68.  *   {@link #removeForceModels()})</li>
  69.  *   <li>the {@link OrbitType type} of orbital parameters to be used for propagation
  70.  *   ({@link #setOrbitType(OrbitType)}),
  71.  *   <li>the {@link PositionAngle type} of position angle to be used in orbital parameters
  72.  *   to be used for propagation where it is relevant ({@link
  73.  *   #setPositionAngleType(PositionAngle)}),
  74.  *   <li>whether {@link org.orekit.propagation.integration.FieldAdditionalEquations additional equations}
  75.  *   should be propagated along with orbital state
  76.  *   ({@link #addAdditionalEquations(org.orekit.propagation.integration.FieldAdditionalEquations)}),
  77.  *   <li>the discrete events that should be triggered during propagation
  78.  *   ({@link #addEventDetector(FieldEventDetector)},
  79.  *   {@link #clearEventsDetectors()})</li>
  80.  *   <li>the binding logic with the rest of the application ({@link #setSlaveMode()},
  81.  *   {@link #setMasterMode(RealFieldElement, org.orekit.propagation.sampling.FieldOrekitFixedStepHandler)},
  82.  *   {@link #setMasterMode(org.orekit.propagation.sampling.FieldOrekitStepHandler)},
  83.  *   {@link #setEphemerisMode()}, {@link #getGeneratedEphemeris()})</li>
  84.  * </ul>
  85.  * <p>From these configuration parameters, only the initial state is mandatory. The default
  86.  * propagation settings are in {@link OrbitType#EQUINOCTIAL equinoctial} parameters with
  87.  * {@link PositionAngle#TRUE true} longitude argument. If the central attraction coefficient
  88.  * is not explicitly specified, the one used to define the initial orbit will be used.
  89.  * However, specifying only the initial state and perhaps the central attraction coefficient
  90.  * would mean the propagator would use only Keplerian forces. In this case, the simpler {@link
  91.  * org.orekit.propagation.analytical.KeplerianPropagator KeplerianPropagator} class would
  92.  * perhaps be more effective.</p>
  93.  * <p>The underlying numerical integrator set up in the constructor may also have its own
  94.  * configuration parameters. Typical configuration parameters for adaptive stepsize integrators
  95.  * are the min, max and perhaps start step size as well as the absolute and/or relative errors
  96.  * thresholds.</p>
  97.  * <p>The state that is seen by the integrator is a simple seven elements double array.
  98.  * The six first elements are either:
  99.  * <ul>
  100.  *   <li>the {@link org.orekit.orbits.FieldEquinoctialOrbit equinoctial orbit parameters} (a, e<sub>x</sub>,
  101.  *   e<sub>y</sub>, h<sub>x</sub>, h<sub>y</sub>, λ<sub>M</sub> or λ<sub>E</sub>
  102.  *   or λ<sub>v</sub>) in meters and radians,</li>
  103.  *   <li>the {@link org.orekit.orbits.FieldKeplerianOrbit Keplerian orbit parameters} (a, e, i, ω, Ω,
  104.  *   M or E or v) in meters and radians,</li>
  105.  *   <li>the {@link org.orekit.orbits.FieldCircularOrbit circular orbit parameters} (a, e<sub>x</sub>, e<sub>y</sub>, i,
  106.  *   Ω, α<sub>M</sub> or α<sub>E</sub> or α<sub>v</sub>) in meters
  107.  *   and radians,</li>
  108.  *   <li>the {@link org.orekit.orbits.FieldCartesianOrbit Cartesian orbit parameters} (x, y, z, v<sub>x</sub>,
  109.  *   v<sub>y</sub>, v<sub>z</sub>) in meters and meters per seconds.
  110.  * </ul>
  111.  * The last element is the mass in kilograms.
  112.  * </p>
  113.  * <p>The following code snippet shows a typical setting for Low Earth Orbit propagation in
  114.  * equinoctial parameters and true longitude argument:</p>
  115.  * <pre>
  116.  * final T          zero      = field.getZero();
  117.  * final T          dP        = zero.add(0.001);
  118.  * final T          minStep   = zero.add(0.001);
  119.  * final T          maxStep   = zero.add(500);
  120.  * final T          initStep  = zero.add(60);
  121.  * final double[][] tolerance = FieldNumericalPropagator.tolerances(dP, orbit, OrbitType.EQUINOCTIAL);
  122.  * AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, minStep, maxStep, tolerance[0], tolerance[1]);
  123.  * integrator.setInitialStepSize(initStep);
  124.  * propagator = new FieldNumericalPropagator<>(field, integrator);
  125.  * </pre>
  126.  * <p>By default, at the end of the propagation, the propagator resets the initial state to the final state,
  127.  * thus allowing a new propagation to be started from there without recomputing the part already performed.
  128.  * This behaviour can be chenged by calling {@link #setResetAtEnd(boolean)}.
  129.  * </p>
  130.  * <p>Beware the same instance cannot be used simultaneously by different threads, the class is <em>not</em>
  131.  * thread-safe.</p>

  133.  * @see FieldSpacecraftState
  134.  * @see ForceModel
  135.  * @see org.orekit.propagation.sampling.FieldOrekitStepHandler
  136.  * @see org.orekit.propagation.sampling.FieldOrekitFixedStepHandler
  137.  * @see org.orekit.propagation.integration.FieldIntegratedEphemeris
  138.  * @see FieldTimeDerivativesEquations
  139.  *
  140.  * @author Mathieu Rom&eacute;ro
  141.  * @author Luc Maisonobe
  142.  * @author Guylaine Prat
  143.  * @author Fabien Maussion
  144.  * @author V&eacute;ronique Pommier-Maurussane
  145.  */
  146. public class FieldNumericalPropagator<T extends RealFieldElement<T>> extends FieldAbstractIntegratedPropagator<T> {

  148.     /** Force models used during the extrapolation of the FieldOrbit<T>, without Jacobians. */
  149.     private final List<ForceModel> forceModels;

  151.     /** Create a new instance of NumericalPropagator, based on orbit definition mu.
  152.      * After creation, the instance is empty, i.e. the attitude provider is set to an
  153.      * unspecified default law and there are no perturbing forces at all.
  154.      * This means that if {@link #addForceModel addForceModel} is not
  155.      * called after creation, the integrated orbit will follow a Keplerian
  156.      * evolution only. The defaults are {@link OrbitType#EQUINOCTIAL}
  157.      * for {@link #setOrbitType(OrbitType) propagation
  158.      * orbit type} and {@link PositionAngle#TRUE} for {@link
  159.      * #setPositionAngleType(PositionAngle) position angle type}.
  160.      * @param integrator numerical integrator to use for propagation.
  161.      * @param field Field used by default
  162.      */
  163.     public FieldNumericalPropagator(final Field<T> field, final FieldODEIntegrator<T> integrator) {
  164.         super(field, integrator, true);
  165.         forceModels = new ArrayList<ForceModel>();
  166.         initMapper();
  167.         setAttitudeProvider(DEFAULT_LAW);
  168.         setMu(Double.NaN);
  169.         setSlaveMode();
  170.         setOrbitType(OrbitType.EQUINOCTIAL);
  171.         setPositionAngleType(PositionAngle.TRUE);
  172.     }

  174.      /** Set the central attraction coefficient μ.
  175.      * @param mu central attraction coefficient (m³/s²)
  176.      * @see #addForceModel(ForceModel)
  177.      */
  178.     public void setMu(final double mu) {
  179.         addForceModel(new NewtonianAttraction(mu));
  180.     }

  182.     /** Set the central attraction coefficient μ only in upper class.
  183.      * @param mu central attraction coefficient (m³/s²)
  184.      */
  185.     private void superSetMu(final double mu) {
  186.         super.setMu(mu);
  187.     }

  189.     /** Check if Newtonian attraction force model is available.
  190.      * <p>
  191.      * Newtonian attraction is always the last force model in the list.
  192.      * </p>
  193.      * @return true if Newtonian attraction force model is available
  194.      */
  195.     private boolean hasNewtonianAttraction() {
  196.         final int last = forceModels.size() - 1;
  197.         return last >= 0 && forceModels.get(last) instanceof NewtonianAttraction;
  198.     }

  200.     /** Add a force model to the global perturbation model.
  201.      * <p>If this method is not called at all, the integrated orbit will follow
  202.      * a Keplerian evolution only.</p>
  203.      * @param model perturbing {@link ForceModel} to add
  204.      * @see #removeForceModels()
  205.      * @see #setMu(double)
  206.      */
  207.     public void addForceModel(final ForceModel model) {

  209.         if (model instanceof NewtonianAttraction) {
  210.             // we want to add the central attraction force model

  212.             try {
  213.                 // ensure we are notified of any mu change
  214.                 model.getParametersDrivers()[0].addObserver(new ParameterObserver() {
  215.                     /** {@inheritDoc} */
  216.                     @Override
  217.                     public void valueChanged(final double previousValue, final ParameterDriver driver) {
  218.                         superSetMu(driver.getValue());
  219.                     }
  220.                 });
  221.             } catch (OrekitException oe) {
  222.                 // this should never happen
  223.                 throw new OrekitInternalError(oe);
  224.             }

  226.             if (hasNewtonianAttraction()) {
  227.                 // there is already a central attraction model, replace it
  228.                 forceModels.set(forceModels.size() - 1, model);
  229.             } else {
  230.                 // there are no central attraction model yet, add it at the end of the list
  231.                 forceModels.add(model);
  232.             }
  233.         } else {
  234.             // we want to add a perturbing force model
  235.             if (hasNewtonianAttraction()) {
  236.                 // insert the new force model before Newtonian attraction,
  237.                 // which should always be the last one in the list
  238.                 forceModels.add(forceModels.size() - 1, model);
  239.             } else {
  240.                 // we only have perturbing force models up to now, just append at the end of the list
  241.                 forceModels.add(model);
  242.             }
  243.         }

  245.     }

  247.     /** Remove all perturbing force models from the global perturbation model.
  248.      * <p>Once all perturbing forces have been removed (and as long as no new force
  249.      * model is added), the integrated orbit will follow a Keplerian evolution
  250.      * only.</p>
  251.      * @see #addForceModel(ForceModel)
  252.      */
  253.     public void removeForceModels() {
  254.         forceModels.clear();
  255.     }

  257.     /** Get all the force models, perturbing forces and Newtonian attraction included.
  258.      * @return list of perturbing force models, with Newtonian attraction being the
  259.      * last one
  260.      * @see #addForceModel(ForceModel)
  261.      * @see #setMu(double)
  262.      * @since 9.1
  263.      */
  264.     public List<ForceModel> getAllForceModels() {
  265.         return Collections.unmodifiableList(forceModels);
  266.     }

  268.     /** Get perturbing force models list.
  269.      * @return list of perturbing force models
  270.      * @deprecated as of 9.1, this method is deprecated, the perturbing
  271.      * force models are retrieved together with the Newtonian attraction
  272.      * by calling {@link #getAllForceModels()}
  273.      */
  274.     @Deprecated
  275.     public List<ForceModel> getForceModels() {
  276.         return hasNewtonianAttraction() ? forceModels.subList(0, forceModels.size() - 1) : forceModels;
  277.     }

  279.     /** Get the Newtonian attraction from the central body force model.
  280.      * @return Newtonian attraction force model
  281.      * @deprecated as of 9.1, this method is deprecated, the Newtonian
  282.      * attraction force model (if any) is the last in the {@link #getAllForceModels()}
  283.      */
  284.     @Deprecated
  285.     public NewtonianAttraction getNewtonianAttractionForceModel() {
  286.         final int last = forceModels.size() - 1;
  287.         if (last >= 0 && forceModels.get(last) instanceof NewtonianAttraction) {
  288.             return (NewtonianAttraction) forceModels.get(last);
  289.         } else {
  290.             return null;
  291.         }
  292.     }

  294.     /** Set propagation orbit type.
  295.      * @param orbitType orbit type to use for propagation
  296.      */
  297.     public void setOrbitType(final OrbitType orbitType) {
  298.         super.setOrbitType(orbitType);
  299.     }

  301.     /** Get propagation parameter type.
  302.      * @return orbit type used for propagation
  303.      */
  304.     public OrbitType getOrbitType() {
  305.         return super.getOrbitType();
  306.     }

  308.     /** Set position angle type.
  309.      * <p>
  310.      * The position parameter type is meaningful only if {@link
  311.      * #getOrbitType() propagation orbit type}
  312.      * support it. As an example, it is not meaningful for propagation
  313.      * in {@link OrbitType#CARTESIAN Cartesian} parameters.
  314.      * </p>
  315.      * @param positionAngleType angle type to use for propagation
  316.      */
  317.     public void setPositionAngleType(final PositionAngle positionAngleType) {
  318.         super.setPositionAngleType(positionAngleType);
  319.     }

  321.     /** Get propagation parameter type.
  322.      * @return angle type to use for propagation
  323.      */
  324.     public PositionAngle getPositionAngleType() {
  325.         return super.getPositionAngleType();
  326.     }

  328.     /** Set the initial state.
  329.      * @param initialState initial state
  330.      * @exception OrekitException if initial state cannot be set
  331.      */
  332.     public void setInitialState(final FieldSpacecraftState<T> initialState) throws OrekitException {
  333.         resetInitialState(initialState);
  334.     }

  336.     /** {@inheritDoc} */
  337.     public void resetInitialState(final FieldSpacecraftState<T> state) throws OrekitException {
  338.         super.resetInitialState(state);
  339.         if (!hasNewtonianAttraction()) {
  340.             setMu(state.getMu());
  341.         }
  342.         setStartDate(state.getDate());
  343.     }

  345.     /** {@inheritDoc} */
  346.     public TimeStampedFieldPVCoordinates<T> getPVCoordinates(final FieldAbsoluteDate<T> date, final Frame frame)
  347.         throws OrekitException {
  348.         return propagate(date).getPVCoordinates(frame);
  349.     }

  351.     /** {@inheritDoc} */
  352.     protected FieldStateMapper<T> createMapper(final FieldAbsoluteDate<T> referenceDate, final double mu,
  353.                                        final OrbitType orbitType, final PositionAngle positionAngleType,
  354.                                        final AttitudeProvider attitudeProvider, final Frame frame) {
  355.         return new FieldOsculatingMapper(referenceDate, mu, orbitType, positionAngleType, attitudeProvider, frame);
  356.     }

  358.     /** Internal mapper using directly osculating parameters. */
  359.     private class FieldOsculatingMapper extends FieldStateMapper<T> {

  361.         /** Simple constructor.
  362.          * <p>
  363.          * The position parameter type is meaningful only if {@link
  364.          * #getOrbitType() propagation orbit type}
  365.          * support it. As an example, it is not meaningful for propagation
  366.          * in {@link OrbitType#CARTESIAN Cartesian} parameters.
  367.          * </p>
  368.          * @param referenceDate reference date
  369.          * @param mu central attraction coefficient (m³/s²)
  370.          * @param orbitType orbit type to use for mapping
  371.          * @param positionAngleType angle type to use for propagation
  372.          * @param attitudeProvider attitude provider
  373.          * @param frame inertial frame
  374.          */
  375.         FieldOsculatingMapper(final FieldAbsoluteDate<T> referenceDate, final double mu,
  376.                               final OrbitType orbitType, final PositionAngle positionAngleType,
  377.                               final AttitudeProvider attitudeProvider, final Frame frame) {
  378.             super(referenceDate, mu, orbitType, positionAngleType, attitudeProvider, frame);
  379.         }

  381.         /** {@inheritDoc} */
  382.         public FieldSpacecraftState<T> mapArrayToState(final FieldAbsoluteDate<T> date, final T[] y, final T[] yDot,
  383.                                                        final boolean meanOnly)
  384.             throws OrekitException {
  385.             // the parameter meanOnly is ignored for the Numerical Propagator

  387.             final T mass = y[6];
  388.             if (mass.getReal() <= 0.0) {
  389.                 throw new OrekitException(OrekitMessages.SPACECRAFT_MASS_BECOMES_NEGATIVE, mass);
  390.             }
  391.             final FieldOrbit<T> orbit       = super.getOrbitType().mapArrayToOrbit(y, yDot, super.getPositionAngleType(), date, getMu(), getFrame());
  392.             final FieldAttitude<T> attitude = getAttitudeProvider().getAttitude(orbit, date, getFrame());
  393.             return new FieldSpacecraftState<>(orbit, attitude, mass);
  394.         }

  396.         /** {@inheritDoc} */
  397.         public void mapStateToArray(final FieldSpacecraftState<T> state, final T[] y, final T[] yDot) {
  398.             super.getOrbitType().mapOrbitToArray(state.getOrbit(), super.getPositionAngleType(), y, yDot);
  399.             y[6] = state.getMass();
  400.         }

  402.     }

  404.     /** {@inheritDoc} */
  405.     protected MainStateEquations<T> getMainStateEquations(final FieldODEIntegrator<T> integrator) {
  406.         return new Main(integrator);
  407.     }

  409.     /** Internal class for osculating parameters integration. */
  410.     private class Main implements MainStateEquations<T>, FieldTimeDerivativesEquations<T> {

  412.         /** Derivatives array. */
  413.         private final T[] yDot;

  415.         /** Current orbit. */
  416.         private FieldOrbit<T> orbit;

  418.         /** Jacobian of the orbital parameters with respect to the Cartesian parameters. */
  419.         private T[][] jacobian;

  421.         /** Simple constructor.
  422.          * @param integrator numerical integrator to use for propagation.
  423.          */
  424.         Main(final FieldODEIntegrator<T> integrator) {

  426.             this.yDot     = MathArrays.buildArray(getField(),  7);
  427.             this.jacobian = MathArrays.buildArray(getField(),  6, 6);
  428.             for (final ForceModel forceModel : forceModels) {
  429.                 forceModel.getFieldEventsDetectors(getField()).forEach(detector -> setUpEventDetector(integrator, detector));
  430.             }

  432.         }

  434.         /** {@inheritDoc} */
  435.         @Override
  436.         public void init(final FieldSpacecraftState<T> initialState, final FieldAbsoluteDate<T> target)
  437.             throws OrekitException {
  438.             final SpacecraftState stateD  = initialState.toSpacecraftState();
  439.             final AbsoluteDate    targetD = target.toAbsoluteDate();
  440.             for (final ForceModel forceModel : forceModels) {
  441.                 forceModel.init(stateD, targetD);
  442.             }
  443.         }

  445.         /** {@inheritDoc} */
  446.         @Override
  447.         public T[] computeDerivatives(final FieldSpacecraftState<T> state) throws OrekitException {
  448.             final T zero = state.getA().getField().getZero();
  449.             orbit = state.getOrbit();
  450.             Arrays.fill(yDot, zero);
  451.             orbit.getJacobianWrtCartesian(getPositionAngleType(), jacobian);

  453.             // compute the contributions of all perturbing forces,
  454.             // using the Kepler contribution at the end since
  455.             // NewtonianAttraction is always the last instance in the list
  456.             for (final ForceModel forceModel : forceModels) {
  457.                 forceModel.addContribution(state, this);
  458.             }

  460.             return yDot.clone();

  462.         }

  464.         /** {@inheritDoc} */
  465.         @Override
  466.         public void addKeplerContribution(final double mu) {
  467.             orbit.addKeplerContribution(getPositionAngleType(), mu, yDot);
  468.         }

  470.         /** {@inheritDoc} */
  471.         @Override
  472.         public void addNonKeplerianAcceleration(final FieldVector3D<T> gamma)
  473.             throws OrekitException {
  474.             for (int i = 0; i < 6; ++i) {
  475.                 final T[] jRow = jacobian[i];
  476.                 yDot[i] = yDot[i].add(jRow[3].linearCombination(jRow[3], gamma.getX(),
  477.                                                                 jRow[4], gamma.getY(),
  478.                                                                 jRow[5], gamma.getZ()));
  479.             }
  480.         }

  482.         /** {@inheritDoc} */
  483.         @Override
  484.         public void addMassDerivative(final T q) {
  485.             if (q.getReal() > 0) {
  486.                 throw new OrekitIllegalArgumentException(OrekitMessages.POSITIVE_FLOW_RATE, q);
  487.             }
  488.             yDot[6] = yDot[6].add(q);
  489.         }

  491.     }

  493.     /** Estimate tolerance vectors for integrators.
  494.      * <p>
  495.      * The errors are estimated from partial derivatives properties of orbits,
  496.      * starting from a scalar position error specified by the user.
  497.      * Considering the energy conservation equation V = sqrt(mu (2/r - 1/a)),
  498.      * we get at constant energy (i.e. on a Keplerian trajectory):
  499.      * <pre>
  500.      * V² r |dV| = mu |dr|
  501.      * </pre>
  502.      * So we deduce a scalar velocity error consistent with the position error.
  503.      * From here, we apply orbits Jacobians matrices to get consistent errors
  504.      * on orbital parameters.
  505.      * </p>
  506.      * <p>
  507.      * The tolerances are only <em>orders of magnitude</em>, and integrator tolerances
  508.      * are only local estimates, not global ones. So some care must be taken when using
  509.      * these tolerances. Setting 1mm as a position error does NOT mean the tolerances
  510.      * will guarantee a 1mm error position after several orbits integration.
  511.      * </p>
  512.      * @param dP user specified position error
  513.      * @param orbit reference orbit
  514.      * @param type propagation type for the meaning of the tolerance vectors elements
  515.      * (it may be different from {@code orbit.getType()})
  516.      * @return a two rows array, row 0 being the absolute tolerance error and row 1
  517.      * being the relative tolerance error
  518.      * @exception OrekitException if Jacobian is singular
  519.      * @param <T> elements type
  520.      */
  521.     public static <T extends RealFieldElement<T>> double[][] tolerances(final T dP, final FieldOrbit<T> orbit, final OrbitType type)
  522.         throws OrekitException {

  524.         // estimate the scalar velocity error
  525.         final FieldPVCoordinates<T> pv = orbit.getPVCoordinates();
  526.         final T r2 = pv.getPosition().getNormSq();
  527.         final T v  = pv.getVelocity().getNorm();
  528.         final T dV = dP.multiply(orbit.getMu()).divide(v.multiply(r2));

  530.         final double[] absTol = new double[7];
  531.         final double[] relTol = new double[7];

  533.         // we set the mass tolerance arbitrarily to 1.0e-6 kg, as mass evolves linearly
  534.         // with trust, this often has no influence at all on propagation
  535.         absTol[6] = 1.0e-6;

  537.         if (type == OrbitType.CARTESIAN) {
  538.             absTol[0] = dP.getReal();
  539.             absTol[1] = dP.getReal();
  540.             absTol[2] = dP.getReal();
  541.             absTol[3] = dV.getReal();
  542.             absTol[4] = dV.getReal();
  543.             absTol[5] = dV.getReal();
  544.         } else {

  546.             // convert the orbit to the desired type
  547.             final T[][] jacobian = MathArrays.buildArray(dP.getField(), 6, 6);
  548.             final FieldOrbit<T> converted = type.convertType(orbit);
  549.             converted.getJacobianWrtCartesian(PositionAngle.TRUE, jacobian);

  551.             for (int i = 0; i < 6; ++i) {
  552.                 final  T[] row = jacobian[i];
  553.                 absTol[i] =     row[0].abs().multiply(dP).
  554.                             add(row[1].abs().multiply(dP)).
  555.                             add(row[2].abs().multiply(dP)).
  556.                             add(row[3].abs().multiply(dV)).
  557.                             add(row[4].abs().multiply(dV)).
  558.                             add(row[5].abs().multiply(dV)).
  559.                             getReal();
  560.                 if (Double.isNaN(absTol[i])) {
  561.                     throw new OrekitException(OrekitMessages.SINGULAR_JACOBIAN_FOR_ORBIT_TYPE, type);
  562.                 }
  563.             }

  565.         }

  567.         Arrays.fill(relTol, dP.divide(r2.sqrt()).getReal());

  569.         return new double[][]{
  570.             absTol, relTol
  571.         };

  573.     }

  575. }
Created with JaCoCo 0.7.9.201702052155