NumericalPropagator.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.io.NotSerializableException;
  20. import java.io.Serializable;
  21. import java.util.ArrayList;
  22. import java.util.Arrays;
  23. import java.util.Collections;
  24. import java.util.List;

  26. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  27. import org.hipparchus.ode.ODEIntegrator;
  28. import org.hipparchus.util.FastMath;
  29. import org.orekit.attitudes.Attitude;
  30. import org.orekit.attitudes.AttitudeProvider;
  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.Orbit;
  39. import org.orekit.orbits.OrbitType;
  40. import org.orekit.orbits.PositionAngle;
  41. import org.orekit.propagation.SpacecraftState;
  42. import org.orekit.propagation.events.EventDetector;
  43. import org.orekit.propagation.integration.AbstractIntegratedPropagator;
  44. import org.orekit.propagation.integration.StateMapper;
  45. import org.orekit.time.AbsoluteDate;
  46. import org.orekit.utils.PVCoordinates;
  47. import org.orekit.utils.ParameterDriver;
  48. import org.orekit.utils.ParameterObserver;
  49. import org.orekit.utils.TimeStampedPVCoordinates;

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

  145.     /** Force models used during the extrapolation of the orbit. */
  146.     private final List<ForceModel> forceModels;

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

  169.      /** Set the central attraction coefficient μ.
  170.       * <p>
  171.       * Setting the central attraction coefficient is
  172.       * equivalent to {@link #addForceModel(ForceModel) add}
  173.       * a {@link NewtonianAttraction} force model.
  174.       * </p>
  175.      * @param mu central attraction coefficient (m³/s²)
  176.      * @see #addForceModel(ForceModel)
  177.      * @see #getAllForceModels()
  178.      */
  179.     public void setMu(final double mu) {
  180.         addForceModel(new NewtonianAttraction(mu));
  181.     }

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

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

  201.     /** Add a force model.
  202.      * <p>If this method is not called at all, the integrated orbit will follow
  203.      * a Keplerian evolution only.</p>
  204.      * @param model {@link ForceModel} to add (it can be either a perturbing force
  205.      * model or an instance of {@link NewtonianAttraction})
  206.      * @see #removeForceModels()
  207.      * @see #setMu(double)
  208.      */
  209.     public void addForceModel(final ForceModel model) {

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

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

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

  247.     }

  249.     /** Remove all force models (except central attraction).
  250.      * <p>Once all perturbing forces have been removed (and as long as no new force
  251.      * model is added), the integrated orbit will follow a Keplerian evolution
  252.      * only.</p>
  253.      * @see #addForceModel(ForceModel)
  254.      */
  255.     public void removeForceModels() {
  256.         final int last = forceModels.size() - 1;
  257.         if (hasNewtonianAttraction()) {
  258.             // preserve the Newtonian attraction model at the end
  259.             final ForceModel newton = forceModels.get(last);
  260.             forceModels.clear();
  261.             forceModels.add(newton);
  262.         } else {
  263.             forceModels.clear();
  264.         }
  265.     }

  267.     /** Get all the force models, perturbing forces and Newtonian attraction included.
  268.      * @return list of perturbing force models, with Newtonian attraction being the
  269.      * last one
  270.      * @see #addForceModel(ForceModel)
  271.      * @see #setMu(double)
  272.      */
  273.     public List<ForceModel> getAllForceModels() {
  274.         return Collections.unmodifiableList(forceModels);
  275.     }

  277.     /** Set propagation orbit type.
  278.      * @param orbitType orbit type to use for propagation
  279.      */
  280.     public void setOrbitType(final OrbitType orbitType) {
  281.         super.setOrbitType(orbitType);
  282.     }

  284.     /** Get propagation parameter type.
  285.      * @return orbit type used for propagation
  286.      */
  287.     public OrbitType getOrbitType() {
  288.         return super.getOrbitType();
  289.     }

  291.     /** Set position angle type.
  292.      * <p>
  293.      * The position parameter type is meaningful only if {@link
  294.      * #getOrbitType() propagation orbit type}
  295.      * support it. As an example, it is not meaningful for propagation
  296.      * in {@link OrbitType#CARTESIAN Cartesian} parameters.
  297.      * </p>
  298.      * @param positionAngleType angle type to use for propagation
  299.      */
  300.     public void setPositionAngleType(final PositionAngle positionAngleType) {
  301.         super.setPositionAngleType(positionAngleType);
  302.     }

  304.     /** Get propagation parameter type.
  305.      * @return angle type to use for propagation
  306.      */
  307.     public PositionAngle getPositionAngleType() {
  308.         return super.getPositionAngleType();
  309.     }

  311.     /** Set the initial state.
  312.      * @param initialState initial state
  313.      * @exception OrekitException if initial state cannot be set
  314.      */
  315.     public void setInitialState(final SpacecraftState initialState) throws OrekitException {
  316.         resetInitialState(initialState);
  317.     }

  319.     /** {@inheritDoc} */
  320.     public void resetInitialState(final SpacecraftState state) throws OrekitException {
  321.         super.resetInitialState(state);
  322.         if (!hasNewtonianAttraction()) {
  323.             // use the state to define central attraction
  324.             setMu(state.getMu());
  325.         }
  326.         setStartDate(state.getDate());
  327.     }

  329.     /** {@inheritDoc} */
  330.     public TimeStampedPVCoordinates getPVCoordinates(final AbsoluteDate date, final Frame frame)
  331.         throws OrekitException {
  332.         return propagate(date).getPVCoordinates(frame);
  333.     }

  335.     /** {@inheritDoc} */
  336.     protected StateMapper createMapper(final AbsoluteDate referenceDate, final double mu,
  337.                                        final OrbitType orbitType, final PositionAngle positionAngleType,
  338.                                        final AttitudeProvider attitudeProvider, final Frame frame) {
  339.         return new OsculatingMapper(referenceDate, mu, orbitType, positionAngleType, attitudeProvider, frame);
  340.     }

  342.     /** Internal mapper using directly osculating parameters. */
  343.     private static class OsculatingMapper extends StateMapper implements Serializable {

  345.         /** Serializable UID. */
  346.         private static final long serialVersionUID = 20130621L;

  348.         /** Simple constructor.
  349.          * <p>
  350.          * The position parameter type is meaningful only if {@link
  351.          * #getOrbitType() propagation orbit type}
  352.          * support it. As an example, it is not meaningful for propagation
  353.          * in {@link OrbitType#CARTESIAN Cartesian} parameters.
  354.          * </p>
  355.          * @param referenceDate reference date
  356.          * @param mu central attraction coefficient (m³/s²)
  357.          * @param orbitType orbit type to use for mapping
  358.          * @param positionAngleType angle type to use for propagation
  359.          * @param attitudeProvider attitude provider
  360.          * @param frame inertial frame
  361.          */
  362.         OsculatingMapper(final AbsoluteDate referenceDate, final double mu,
  363.                          final OrbitType orbitType, final PositionAngle positionAngleType,
  364.                          final AttitudeProvider attitudeProvider, final Frame frame) {
  365.             super(referenceDate, mu, orbitType, positionAngleType, attitudeProvider, frame);
  366.         }

  368.         /** {@inheritDoc} */
  369.         public SpacecraftState mapArrayToState(final AbsoluteDate date, final double[] y, final double[] yDot,
  370.                                                final boolean meanOnly)
  371.             throws OrekitException {
  372.             // the parameter meanOnly is ignored for the Numerical Propagator

  374.             final double mass = y[6];
  375.             if (mass <= 0.0) {
  376.                 throw new OrekitException(OrekitMessages.SPACECRAFT_MASS_BECOMES_NEGATIVE, mass);
  377.             }

  379.             final Orbit orbit       = getOrbitType().mapArrayToOrbit(y, yDot, getPositionAngleType(), date, getMu(), getFrame());
  380.             final Attitude attitude = getAttitudeProvider().getAttitude(orbit, date, getFrame());

  382.             return new SpacecraftState(orbit, attitude, mass);

  384.         }

  386.         /** {@inheritDoc} */
  387.         public void mapStateToArray(final SpacecraftState state, final double[] y, final double[] yDot) {
  388.             getOrbitType().mapOrbitToArray(state.getOrbit(), getPositionAngleType(), y, yDot);
  389.             y[6] = state.getMass();
  390.         }

  392.         /** Replace the instance with a data transfer object for serialization.
  393.          * @return data transfer object that will be serialized
  394.          * @exception NotSerializableException if the state mapper cannot be serialized (typically for DSST propagator)
  395.          */
  396.         private Object writeReplace() throws NotSerializableException {
  397.             return new DataTransferObject(getReferenceDate(), getMu(), getOrbitType(),
  398.                                           getPositionAngleType(), getAttitudeProvider(), getFrame());
  399.         }

  401.         /** Internal class used only for serialization. */
  402.         private static class DataTransferObject implements Serializable {

  404.             /** Serializable UID. */
  405.             private static final long serialVersionUID = 20130621L;

  407.             /** Reference date. */
  408.             private final AbsoluteDate referenceDate;

  410.             /** Central attraction coefficient (m³/s²). */
  411.             private final double mu;

  413.             /** Orbit type to use for mapping. */
  414.             private final OrbitType orbitType;

  416.             /** Angle type to use for propagation. */
  417.             private final PositionAngle positionAngleType;

  419.             /** Attitude provider. */
  420.             private final AttitudeProvider attitudeProvider;

  422.             /** Inertial frame. */
  423.             private final Frame frame;

  425.             /** Simple constructor.
  426.              * @param referenceDate reference date
  427.              * @param mu central attraction coefficient (m³/s²)
  428.              * @param orbitType orbit type to use for mapping
  429.              * @param positionAngleType angle type to use for propagation
  430.              * @param attitudeProvider attitude provider
  431.              * @param frame inertial frame
  432.              */
  433.             DataTransferObject(final AbsoluteDate referenceDate, final double mu,
  434.                                       final OrbitType orbitType, final PositionAngle positionAngleType,
  435.                                       final AttitudeProvider attitudeProvider, final Frame frame) {
  436.                 this.referenceDate     = referenceDate;
  437.                 this.mu                = mu;
  438.                 this.orbitType         = orbitType;
  439.                 this.positionAngleType = positionAngleType;
  440.                 this.attitudeProvider  = attitudeProvider;
  441.                 this.frame             = frame;
  442.             }

  444.             /** Replace the deserialized data transfer object with a {@link OsculatingMapper}.
  445.              * @return replacement {@link OsculatingMapper}
  446.              */
  447.             private Object readResolve() {
  448.                 return new OsculatingMapper(referenceDate, mu, orbitType, positionAngleType, attitudeProvider, frame);
  449.             }
  450.         }

  452.     }

  454.     /** {@inheritDoc} */
  455.     protected MainStateEquations getMainStateEquations(final ODEIntegrator integrator) {
  456.         return new Main(integrator);
  457.     }

  459.     /** Internal class for osculating parameters integration. */
  460.     private class Main implements MainStateEquations, TimeDerivativesEquations {

  462.         /** Derivatives array. */
  463.         private final double[] yDot;

  465.         /** Current orbit. */
  466.         private Orbit orbit;

  468.         /** Jacobian of the orbital parameters with respect to the Cartesian parameters. */
  469.         private double[][] jacobian;

  471.         /** Simple constructor.
  472.          * @param integrator numerical integrator to use for propagation.
  473.          */
  474.         Main(final ODEIntegrator integrator) {

  476.             this.yDot     = new double[7];
  477.             this.jacobian = new double[6][6];

  479.             for (final ForceModel forceModel : forceModels) {
  480.                 forceModel.getEventsDetectors().forEach(detector -> setUpEventDetector(integrator, detector));
  481.             }

  483.         }

  485.         /** {@inheritDoc} */
  486.         @Override
  487.         public void init(final SpacecraftState initialState, final AbsoluteDate target)
  488.                 throws OrekitException {
  489.             for (final ForceModel forceModel : forceModels) {
  490.                 forceModel.init(initialState, target);
  491.             }
  492.         }

  494.         /** {@inheritDoc} */
  495.         @Override
  496.         public double[] computeDerivatives(final SpacecraftState state) throws OrekitException {

  498.             orbit = state.getOrbit();
  499.             Arrays.fill(yDot, 0.0);
  500.             orbit.getJacobianWrtCartesian(getPositionAngleType(), jacobian);

  502.             // compute the contributions of all perturbing forces,
  503.             // using the Kepler contribution at the end since
  504.             // NewtonianAttraction is always the last instance in the list
  505.             for (final ForceModel forceModel : forceModels) {
  506.                 forceModel.addContribution(state, this);
  507.             }

  509.             return yDot.clone();

  511.         }

  513.         /** {@inheritDoc} */
  514.         @Override
  515.         public void addKeplerContribution(final double mu) {
  516.             orbit.addKeplerContribution(getPositionAngleType(), mu, yDot);
  517.         }

  519.         /** {@inheritDoc} */
  520.         public void addNonKeplerianAcceleration(final Vector3D gamma)
  521.             throws OrekitException {
  522.             for (int i = 0; i < 6; ++i) {
  523.                 final double[] jRow = jacobian[i];
  524.                 yDot[i] += jRow[3] * gamma.getX() + jRow[4] * gamma.getY() + jRow[5] * gamma.getZ();
  525.             }
  526.         }

  528.         /** {@inheritDoc} */
  529.         @Override
  530.         public void addMassDerivative(final double q) {
  531.             if (q > 0) {
  532.                 throw new OrekitIllegalArgumentException(OrekitMessages.POSITIVE_FLOW_RATE, q);
  533.             }
  534.             yDot[6] += q;
  535.         }

  537.     }

  539.     /** Estimate tolerance vectors for integrators.
  540.      * <p>
  541.      * The errors are estimated from partial derivatives properties of orbits,
  542.      * starting from a scalar position error specified by the user.
  543.      * Considering the energy conservation equation V = sqrt(mu (2/r - 1/a)),
  544.      * we get at constant energy (i.e. on a Keplerian trajectory):
  545.      * <pre>
  546.      * V² r |dV| = mu |dr|
  547.      * </pre>
  548.      * <p> So we deduce a scalar velocity error consistent with the position error.
  549.      * From here, we apply orbits Jacobians matrices to get consistent errors
  550.      * on orbital parameters.
  551.      *
  552.      * <p>
  553.      * The tolerances are only <em>orders of magnitude</em>, and integrator tolerances
  554.      * are only local estimates, not global ones. So some care must be taken when using
  555.      * these tolerances. Setting 1mm as a position error does NOT mean the tolerances
  556.      * will guarantee a 1mm error position after several orbits integration.
  557.      * </p>
  558.      * @param dP user specified position error
  559.      * @param orbit reference orbit
  560.      * @param type propagation type for the meaning of the tolerance vectors elements
  561.      * (it may be different from {@code orbit.getType()})
  562.      * @return a two rows array, row 0 being the absolute tolerance error and row 1
  563.      * being the relative tolerance error
  564.      * @exception OrekitException if Jacobian is singular
  565.      */
  566.     public static double[][] tolerances(final double dP, final Orbit orbit, final OrbitType type)
  567.         throws OrekitException {

  569.         // estimate the scalar velocity error
  570.         final PVCoordinates pv = orbit.getPVCoordinates();
  571.         final double r2 = pv.getPosition().getNormSq();
  572.         final double v  = pv.getVelocity().getNorm();
  573.         final double dV = orbit.getMu() * dP / (v * r2);

  575.         final double[] absTol = new double[7];
  576.         final double[] relTol = new double[7];

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

  582.         if (type == OrbitType.CARTESIAN) {
  583.             absTol[0] = dP;
  584.             absTol[1] = dP;
  585.             absTol[2] = dP;
  586.             absTol[3] = dV;
  587.             absTol[4] = dV;
  588.             absTol[5] = dV;
  589.         } else {

  591.             // convert the orbit to the desired type
  592.             final double[][] jacobian = new double[6][6];
  593.             final Orbit converted = type.convertType(orbit);
  594.             converted.getJacobianWrtCartesian(PositionAngle.TRUE, jacobian);

  596.             for (int i = 0; i < 6; ++i) {
  597.                 final double[] row = jacobian[i];
  598.                 absTol[i] = FastMath.abs(row[0]) * dP +
  599.                             FastMath.abs(row[1]) * dP +
  600.                             FastMath.abs(row[2]) * dP +
  601.                             FastMath.abs(row[3]) * dV +
  602.                             FastMath.abs(row[4]) * dV +
  603.                             FastMath.abs(row[5]) * dV;
  604.                 if (Double.isNaN(absTol[i])) {
  605.                     throw new OrekitException(OrekitMessages.SINGULAR_JACOBIAN_FOR_ORBIT_TYPE, type);
  606.                 }
  607.             }

  609.         }

  611.         Arrays.fill(relTol, dP / FastMath.sqrt(r2));

  613.         return new double[][] {
  614.             absTol, relTol
  615.         };

  617.     }

  619. }
Created with JaCoCo 0.7.9.201702052155