NumericalPropagator.java

  1. /* Copyright 2002-2016 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;

  18. import java.io.NotSerializableException;
  19. import java.io.Serializable;
  20. import java.util.ArrayList;
  21. import java.util.Arrays;
  22. import java.util.List;

  23. import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
  24. import org.apache.commons.math3.ode.AbstractIntegrator;
  25. import org.apache.commons.math3.util.FastMath;
  26. import org.orekit.attitudes.Attitude;
  27. import org.orekit.attitudes.AttitudeProvider;
  28. import org.orekit.errors.OrekitIllegalArgumentException;
  29. import org.orekit.errors.OrekitException;
  30. import org.orekit.errors.OrekitMessages;
  31. import org.orekit.errors.PropagationException;
  32. import org.orekit.forces.ForceModel;
  33. import org.orekit.forces.gravity.NewtonianAttraction;
  34. import org.orekit.frames.Frame;
  35. import org.orekit.frames.Transform;
  36. import org.orekit.orbits.Orbit;
  37. import org.orekit.orbits.OrbitType;
  38. import org.orekit.orbits.PositionAngle;
  39. import org.orekit.propagation.SpacecraftState;
  40. import org.orekit.propagation.events.EventDetector;
  41. import org.orekit.propagation.integration.AbstractIntegratedPropagator;
  42. import org.orekit.propagation.integration.StateMapper;
  43. import org.orekit.time.AbsoluteDate;
  44. import org.orekit.utils.PVCoordinates;
  45. import org.orekit.utils.TimeStampedPVCoordinates;

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

  122.  * @see SpacecraftState
  123.  * @see ForceModel
  124.  * @see org.orekit.propagation.sampling.OrekitStepHandler
  125.  * @see org.orekit.propagation.sampling.OrekitFixedStepHandler
  126.  * @see org.orekit.propagation.integration.IntegratedEphemeris
  127.  * @see TimeDerivativesEquations
  128.  *
  129.  * @author Mathieu Rom&eacute;ro
  130.  * @author Luc Maisonobe
  131.  * @author Guylaine Prat
  132.  * @author Fabien Maussion
  133.  * @author V&eacute;ronique Pommier-Maurussane
  134.  */
  135. public class NumericalPropagator extends AbstractIntegratedPropagator {

  136.     /** Central body attraction. */
  137.     private NewtonianAttraction newtonianAttraction;

  138.     /** Force models used during the extrapolation of the Orbit, without jacobians. */
  139.     private final List<ForceModel> forceModels;

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

  161.      /** Set the central attraction coefficient μ.
  162.      * @param mu central attraction coefficient (m³/s²)
  163.      * @see #addForceModel(ForceModel)
  164.      */
  165.     public void setMu(final double mu) {
  166.         super.setMu(mu);
  167.         newtonianAttraction = new NewtonianAttraction(mu);
  168.     }

  169.     /** Add a force model to the global perturbation model.
  170.      * <p>If this method is not called at all, the integrated orbit will follow
  171.      * a keplerian evolution only.</p>
  172.      * @param model perturbing {@link ForceModel} to add
  173.      * @see #removeForceModels()
  174.      * @see #setMu(double)
  175.      */
  176.     public void addForceModel(final ForceModel model) {
  177.         forceModels.add(model);
  178.     }

  179.     /** Remove all perturbing force models from the global perturbation model.
  180.      * <p>Once all perturbing forces have been removed (and as long as no new force
  181.      * model is added), the integrated orbit will follow a keplerian evolution
  182.      * only.</p>
  183.      * @see #addForceModel(ForceModel)
  184.      */
  185.     public void removeForceModels() {
  186.         forceModels.clear();
  187.     }

  188.     /** Get perturbing force models list.
  189.      * @return list of perturbing force models
  190.      * @see #addForceModel(ForceModel)
  191.      * @see #getNewtonianAttractionForceModel()
  192.      */
  193.     public List<ForceModel> getForceModels() {
  194.         return forceModels;
  195.     }

  196.     /** Get the Newtonian attraction from the central body force model.
  197.      * @return Newtonian attraction force model
  198.      * @see #setMu(double)
  199.      * @see #getForceModels()
  200.      */
  201.     public NewtonianAttraction getNewtonianAttractionForceModel() {
  202.         return newtonianAttraction;
  203.     }

  204.     /** Set propagation orbit type.
  205.      * @param orbitType orbit type to use for propagation
  206.      */
  207.     public void setOrbitType(final OrbitType orbitType) {
  208.         super.setOrbitType(orbitType);
  209.     }

  210.     /** Get propagation parameter type.
  211.      * @return orbit type used for propagation
  212.      */
  213.     public OrbitType getOrbitType() {
  214.         return super.getOrbitType();
  215.     }

  216.     /** Set position angle type.
  217.      * <p>
  218.      * The position parameter type is meaningful only if {@link
  219.      * #getOrbitType() propagation orbit type}
  220.      * support it. As an example, it is not meaningful for propagation
  221.      * in {@link OrbitType#CARTESIAN Cartesian} parameters.
  222.      * </p>
  223.      * @param positionAngleType angle type to use for propagation
  224.      */
  225.     public void setPositionAngleType(final PositionAngle positionAngleType) {
  226.         super.setPositionAngleType(positionAngleType);
  227.     }

  228.     /** Get propagation parameter type.
  229.      * @return angle type to use for propagation
  230.      */
  231.     public PositionAngle getPositionAngleType() {
  232.         return super.getPositionAngleType();
  233.     }

  234.     /** Set the initial state.
  235.      * @param initialState initial state
  236.      * @exception PropagationException if initial state cannot be set
  237.      */
  238.     public void setInitialState(final SpacecraftState initialState) throws PropagationException {
  239.         resetInitialState(initialState);
  240.     }

  241.     /** {@inheritDoc} */
  242.     public void resetInitialState(final SpacecraftState state) throws PropagationException {
  243.         super.resetInitialState(state);
  244.         if (newtonianAttraction == null) {
  245.             setMu(state.getMu());
  246.         }
  247.         setStartDate(null);
  248.     }

  249.     /** {@inheritDoc} */
  250.     public TimeStampedPVCoordinates getPVCoordinates(final AbsoluteDate date, final Frame frame)
  251.         throws OrekitException {
  252.         return propagate(date).getPVCoordinates(frame);
  253.     }

  254.     /** {@inheritDoc} */
  255.     protected StateMapper createMapper(final AbsoluteDate referenceDate, final double mu,
  256.                                        final OrbitType orbitType, final PositionAngle positionAngleType,
  257.                                        final AttitudeProvider attitudeProvider, final Frame frame) {
  258.         return new OsculatingMapper(referenceDate, mu, orbitType, positionAngleType, attitudeProvider, frame);
  259.     }

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

  262.         /** Serializable UID. */
  263.         private static final long serialVersionUID = 20130621L;

  264.         /** Simple constructor.
  265.          * <p>
  266.          * The position parameter type is meaningful only if {@link
  267.          * #getOrbitType() propagation orbit type}
  268.          * support it. As an example, it is not meaningful for propagation
  269.          * in {@link OrbitType#CARTESIAN Cartesian} parameters.
  270.          * </p>
  271.          * @param referenceDate reference date
  272.          * @param mu central attraction coefficient (m³/s²)
  273.          * @param orbitType orbit type to use for mapping
  274.          * @param positionAngleType angle type to use for propagation
  275.          * @param attitudeProvider attitude provider
  276.          * @param frame inertial frame
  277.          */
  278.         OsculatingMapper(final AbsoluteDate referenceDate, final double mu,
  279.                          final OrbitType orbitType, final PositionAngle positionAngleType,
  280.                          final AttitudeProvider attitudeProvider, final Frame frame) {
  281.             super(referenceDate, mu, orbitType, positionAngleType, attitudeProvider, frame);
  282.         }

  283.         /** {@inheritDoc} */
  284.         public SpacecraftState mapArrayToState(final AbsoluteDate date, final double[] y, final boolean meanOnly)
  285.             throws OrekitException {
  286.             // the parameter meanOnly is ignored for the Numerical Propagator

  287.             final double mass = y[6];
  288.             if (mass <= 0.0) {
  289.                 throw new PropagationException(OrekitMessages.SPACECRAFT_MASS_BECOMES_NEGATIVE, mass);
  290.             }

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

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

  294.         }

  295.         /** {@inheritDoc} */
  296.         public void mapStateToArray(final SpacecraftState state, final double[] y) {
  297.             getOrbitType().mapOrbitToArray(state.getOrbit(), getPositionAngleType(), y);
  298.             y[6] = state.getMass();
  299.         }

  300.         /** Replace the instance with a data transfer object for serialization.
  301.          * @return data transfer object that will be serialized
  302.          * @exception NotSerializableException if the state mapper cannot be serialized (typically for DSST propagator)
  303.          */
  304.         private Object writeReplace() throws NotSerializableException {
  305.             return new DataTransferObject(getReferenceDate(), getMu(), getOrbitType(),
  306.                                           getPositionAngleType(), getAttitudeProvider(), getFrame());
  307.         }

  308.         /** Internal class used only for serialization. */
  309.         private static class DataTransferObject implements Serializable {

  310.             /** Serializable UID. */
  311.             private static final long serialVersionUID = 20130621L;

  312.             /** Reference date. */
  313.             private final AbsoluteDate referenceDate;

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

  316.             /** Orbit type to use for mapping. */
  317.             private final OrbitType orbitType;

  318.             /** Angle type to use for propagation. */
  319.             private final PositionAngle positionAngleType;

  320.             /** Attitude provider. */
  321.             private final AttitudeProvider attitudeProvider;

  322.             /** Inertial frame. */
  323.             private final Frame frame;

  324.             /** Simple constructor.
  325.              * @param referenceDate reference date
  326.              * @param mu central attraction coefficient (m³/s²)
  327.              * @param orbitType orbit type to use for mapping
  328.              * @param positionAngleType angle type to use for propagation
  329.              * @param attitudeProvider attitude provider
  330.              * @param frame inertial frame
  331.              */
  332.             DataTransferObject(final AbsoluteDate referenceDate, final double mu,
  333.                                       final OrbitType orbitType, final PositionAngle positionAngleType,
  334.                                       final AttitudeProvider attitudeProvider, final Frame frame) {
  335.                 this.referenceDate     = referenceDate;
  336.                 this.mu                = mu;
  337.                 this.orbitType         = orbitType;
  338.                 this.positionAngleType = positionAngleType;
  339.                 this.attitudeProvider  = attitudeProvider;
  340.                 this.frame             = frame;
  341.             }

  342.             /** Replace the deserialized data transfer object with a {@link OsculatingMapper}.
  343.              * @return replacement {@link OsculatingMapper}
  344.              */
  345.             private Object readResolve() {
  346.                 return new OsculatingMapper(referenceDate, mu, orbitType, positionAngleType, attitudeProvider, frame);
  347.             }
  348.         }

  349.     }

  350.     /** {@inheritDoc} */
  351.     protected MainStateEquations getMainStateEquations(final AbstractIntegrator integrator) {
  352.         return new Main(integrator);
  353.     }

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

  356.         /** Derivatives array. */
  357.         private final double[] yDot;

  358.         /** Current orbit. */
  359.         private Orbit orbit;

  360.         /** Jacobian of the orbital parameters with respect to the cartesian parameters. */
  361.         private double[][] jacobian;

  362.         /** Simple constructor.
  363.          * @param integrator numerical integrator to use for propagation.
  364.          */
  365.         Main(final AbstractIntegrator integrator) {

  366.             this.yDot     = new double[7];
  367.             this.jacobian = new double[6][6];

  368.             for (final ForceModel forceModel : forceModels) {
  369.                 final EventDetector[] modelDetectors = forceModel.getEventsDetectors();
  370.                 if (modelDetectors != null) {
  371.                     for (final EventDetector detector : modelDetectors) {
  372.                         setUpEventDetector(integrator, detector);
  373.                     }
  374.                 }
  375.             }

  376.         }

  377.         /** {@inheritDoc} */
  378.         public double[] computeDerivatives(final SpacecraftState state) throws OrekitException {

  379.             orbit = state.getOrbit();
  380.             Arrays.fill(yDot, 0.0);
  381.             orbit.getJacobianWrtCartesian(getPositionAngleType(), jacobian);

  382.             // compute the contributions of all perturbing forces
  383.             for (final ForceModel forceModel : forceModels) {
  384.                 forceModel.addContribution(state, this);
  385.             }

  386.             // finalize derivatives by adding the Kepler contribution
  387.             newtonianAttraction.addContribution(state, this);

  388.             return yDot.clone();

  389.         }

  390.         /** {@inheritDoc} */
  391.         public void addKeplerContribution(final double mu) {
  392.             orbit.addKeplerContribution(getPositionAngleType(), mu, yDot);
  393.         }

  394.         /** {@inheritDoc} */
  395.         public void addXYZAcceleration(final double x, final double y, final double z) {
  396.             for (int i = 0; i < 6; ++i) {
  397.                 final double[] jRow = jacobian[i];
  398.                 yDot[i] += jRow[3] * x + jRow[4] * y + jRow[5] * z;
  399.             }
  400.         }

  401.         /** {@inheritDoc} */
  402.         public void addAcceleration(final Vector3D gamma, final Frame frame)
  403.             throws OrekitException {
  404.             final Transform t = frame.getTransformTo(orbit.getFrame(), orbit.getDate());
  405.             final Vector3D gammInRefFrame = t.transformVector(gamma);
  406.             addXYZAcceleration(gammInRefFrame.getX(), gammInRefFrame.getY(), gammInRefFrame.getZ());
  407.         }

  408.         /** {@inheritDoc} */
  409.         public void addMassDerivative(final double q) {
  410.             if (q > 0) {
  411.                 throw new OrekitIllegalArgumentException(OrekitMessages.POSITIVE_FLOW_RATE, q);
  412.             }
  413.             yDot[6] += q;
  414.         }

  415.     }

  416.     /** Estimate tolerance vectors for integrators.
  417.      * <p>
  418.      * The errors are estimated from partial derivatives properties of orbits,
  419.      * starting from a scalar position error specified by the user.
  420.      * Considering the energy conservation equation V = sqrt(mu (2/r - 1/a)),
  421.      * we get at constant energy (i.e. on a Keplerian trajectory):
  422.      * <pre>
  423.      * V² r |dV| = mu |dr|
  424.      * </pre>
  425.      * So we deduce a scalar velocity error consistent with the position error.
  426.      * From here, we apply orbits Jacobians matrices to get consistent errors
  427.      * on orbital parameters.
  428.      * </p>
  429.      * <p>
  430.      * The tolerances are only <em>orders of magnitude</em>, and integrator tolerances
  431.      * are only local estimates, not global ones. So some care must be taken when using
  432.      * these tolerances. Setting 1mm as a position error does NOT mean the tolerances
  433.      * will guarantee a 1mm error position after several orbits integration.
  434.      * </p>
  435.      * @param dP user specified position error
  436.      * @param orbit reference orbit
  437.      * @param type propagation type for the meaning of the tolerance vectors elements
  438.      * (it may be different from {@code orbit.getType()})
  439.      * @return a two rows array, row 0 being the absolute tolerance error and row 1
  440.      * being the relative tolerance error
  441.      * @exception PropagationException if Jacobian is singular
  442.      */
  443.     public static double[][] tolerances(final double dP, final Orbit orbit, final OrbitType type)
  444.         throws PropagationException {

  445.         // estimate the scalar velocity error
  446.         final PVCoordinates pv = orbit.getPVCoordinates();
  447.         final double r2 = pv.getPosition().getNormSq();
  448.         final double v  = pv.getVelocity().getNorm();
  449.         final double dV = orbit.getMu() * dP / (v * r2);

  450.         final double[] absTol = new double[7];
  451.         final double[] relTol = new double[7];

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

  455.         if (type == OrbitType.CARTESIAN) {
  456.             absTol[0] = dP;
  457.             absTol[1] = dP;
  458.             absTol[2] = dP;
  459.             absTol[3] = dV;
  460.             absTol[4] = dV;
  461.             absTol[5] = dV;
  462.         } else {

  463.             // convert the orbit to the desired type
  464.             final double[][] jacobian = new double[6][6];
  465.             final Orbit converted = type.convertType(orbit);
  466.             converted.getJacobianWrtCartesian(PositionAngle.TRUE, jacobian);

  467.             for (int i = 0; i < 6; ++i) {
  468.                 final double[] row = jacobian[i];
  469.                 absTol[i] = FastMath.abs(row[0]) * dP +
  470.                             FastMath.abs(row[1]) * dP +
  471.                             FastMath.abs(row[2]) * dP +
  472.                             FastMath.abs(row[3]) * dV +
  473.                             FastMath.abs(row[4]) * dV +
  474.                             FastMath.abs(row[5]) * dV;
  475.                 if (Double.isNaN(absTol[i])) {
  476.                     throw new PropagationException(OrekitMessages.SINGULAR_JACOBIAN_FOR_ORBIT_TYPE, type);
  477.                 }
  478.             }

  479.         }

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

  481.         return new double[][] {
  482.             absTol, relTol
  483.         };

  484.     }

  485. }