ToleranceProvider.java

  1. /* Copyright 2022-2025 Romain Serra
  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;

  19. import org.hipparchus.CalculusFieldElement;
  20. import org.hipparchus.geometry.euclidean.threed.Vector3D;
  21. import org.hipparchus.util.FastMath;
  22. import org.orekit.errors.OrekitException;
  23. import org.orekit.errors.OrekitMessages;
  24. import org.orekit.orbits.FieldOrbit;
  25. import org.orekit.orbits.Orbit;
  26. import org.orekit.orbits.OrbitType;
  27. import org.orekit.orbits.PositionAngleType;
  28. import org.orekit.propagation.numerical.FieldNumericalPropagator;
  29. import org.orekit.propagation.numerical.NumericalPropagator;
  30. import org.orekit.utils.PVCoordinates;

  32. import java.util.Arrays;

  34. /**
  35.  * Interface to define integration tolerances for adaptive schemes (like the embedded Runge-Kutta ones) propagating
  36.  * the position-velocity vector (or an equivalent set of coordinates) and the mass, for a total of 7 primary dependent variables (in that order).
  37.  * The tolerances are given as an array of array: each row has 7 elements, whilst the first column is the absolute tolerances and the second the relative ones.
  38.  *
  39.  * @see NumericalPropagator
  40.  * @see FieldNumericalPropagator
  41.  * @see CartesianToleranceProvider
  42.  * @since 13.0
  43.  * @author Romain Serra
  44.  */
  45. public interface ToleranceProvider extends CartesianToleranceProvider {

  47.     /**
  48.      * Retrieve the integration tolerances given a reference orbit.
  49.      * @param referenceOrbit orbit
  50.      * @param propagationOrbitType orbit type for propagation (can be different from the input orbit one)
  51.      * @param positionAngleType reference position angle type
  52.      * @return absolute and relative tolerances
  53.      */
  54.     double[][] getTolerances(Orbit referenceOrbit, OrbitType propagationOrbitType,
  55.                              PositionAngleType positionAngleType);

  57.     /**
  58.      * Retrieve the integration tolerances given a reference orbit.
  59.      * @param referenceOrbit orbit
  60.      * @param propagationOrbitType orbit type for propagation (can be different from the input orbit one)
  61.      * @return absolute and relative tolerances
  62.      */
  63.     default double[][] getTolerances(Orbit referenceOrbit, OrbitType propagationOrbitType) {
  64.         return getTolerances(referenceOrbit, propagationOrbitType, NumericalPropagator.DEFAULT_POSITION_ANGLE_TYPE);
  65.     }

  67.     /**
  68.      * Retrieve the integration tolerances given a reference Field orbit.
  69.      * @param referenceOrbit orbit
  70.      * @param propagationOrbitType orbit type for propagation (can be different from the input orbit one)
  71.      * @param positionAngleType reference position angle type
  72.      * @param <T> field type
  73.      * @return absolute and relative tolerances
  74.      */
  75.     default <T extends CalculusFieldElement<T>> double[][] getTolerances(FieldOrbit<T> referenceOrbit,
  76.                                                                          OrbitType propagationOrbitType,
  77.                                                                          PositionAngleType positionAngleType) {
  78.         return getTolerances(referenceOrbit.toOrbit(), propagationOrbitType, positionAngleType);
  79.     }

  81.     /**
  82.      * Retrieve the integration tolerances given a reference Field orbit.
  83.      * @param referenceOrbit orbit
  84.      * @param propagationOrbitType orbit type for propagation (can be different from the input orbit one)
  85.      * @param <T> field type
  86.      * @return absolute and relative tolerances
  87.      */
  88.     default <T extends CalculusFieldElement<T>> double[][] getTolerances(FieldOrbit<T> referenceOrbit,
  89.                                                                          OrbitType propagationOrbitType) {
  90.         return getTolerances(referenceOrbit, propagationOrbitType, NumericalPropagator.DEFAULT_POSITION_ANGLE_TYPE);
  91.     }

  93.     /**
  94.      * Build a provider using a single value for absolute and respective tolerance respectively.
  95.      * @param absoluteTolerance absolute tolerance value to be used
  96.      * @param relativeTolerance relative tolerance value to be used
  97.      * @return tolerance provider
  98.      */
  99.     static ToleranceProvider of(final double absoluteTolerance, final double relativeTolerance) {

  101.         return new ToleranceProvider() {
  102.             @Override
  103.             public double[][] getTolerances(final Orbit referenceOrbit, final OrbitType propagationOrbitType,
  104.                                             final PositionAngleType positionAngleType) {
  105.                 return getTolerances();
  106.             }

  108.             @Override
  109.             public double[][] getTolerances(final Vector3D position, final Vector3D velocity) {
  110.                 return getTolerances();
  111.             }

  113.             /**
  114.              * Retrieve constant absolute and respective tolerances.
  115.              * @return tolerances
  116.              */
  117.             double[][] getTolerances() {
  118.                 final double[] absTol = new double[7];
  119.                 Arrays.fill(absTol, absoluteTolerance);
  120.                 final double[] relTol = new double[absTol.length];
  121.                 Arrays.fill(relTol, relativeTolerance);
  122.                 return new double[][] { absTol, relTol };
  123.             }
  124.         };
  125.     }

  127.     /**
  128.      * Build a provider based on a tolerance provider for Cartesian coordinates.
  129.      * <p> Orbits Jacobian matrices are used to get consistent errors on orbital parameters.
  130.      * <p>
  131.      *
  132.      * @param cartesianToleranceProvider tolerance provider dedicated to Cartesian propagation
  133.      * @return tolerance provider
  134.      */
  135.     static ToleranceProvider of(final CartesianToleranceProvider cartesianToleranceProvider) {
  136.         return new ToleranceProvider() {

  138.             @Override
  139.             public double[][] getTolerances(final Vector3D position, final Vector3D velocity) {
  140.                 return cartesianToleranceProvider.getTolerances(position, velocity);
  141.             }

  143.             @Override
  144.             public double[][] getTolerances(final Orbit referenceOrbit, final OrbitType propagationOrbitType,
  145.                                             final PositionAngleType positionAngleType) {
  146.                 // compute Cartesian-related tolerances
  147.                 final double[][] cartesianTolerances = getTolerances(referenceOrbit.getPosition(),
  148.                         referenceOrbit.getPVCoordinates().getVelocity());
  149.                 if (propagationOrbitType == OrbitType.CARTESIAN) {
  150.                     return cartesianTolerances;
  151.                 }

  153.                 final double[] cartAbsTol = cartesianTolerances[0];
  154.                 final double[] cartRelTol = cartesianTolerances[1];
  155.                 final double[] absTol = new double[7];
  156.                 final double[] relTol = absTol.clone();

  158.                 // convert the orbit to the desired type
  159.                 final double[][] jacobian = new double[6][6];
  160.                 final Orbit converted = propagationOrbitType.convertType(referenceOrbit);
  161.                 converted.getJacobianWrtCartesian(positionAngleType, jacobian);

  163.                 double minimumRel = cartRelTol[6];
  164.                 for (int i = 0; i < jacobian.length; ++i) {
  165.                     final double[] row = jacobian[i];
  166.                     absTol[i] = FastMath.abs(row[0]) * cartAbsTol[0] +
  167.                             FastMath.abs(row[1]) * cartAbsTol[1] +
  168.                             FastMath.abs(row[2]) * cartAbsTol[2] +
  169.                             FastMath.abs(row[3]) * cartAbsTol[3] +
  170.                             FastMath.abs(row[4]) * cartAbsTol[4] +
  171.                             FastMath.abs(row[5]) * cartAbsTol[5];
  172.                     if (Double.isNaN(absTol[i])) {
  173.                         throw new OrekitException(OrekitMessages.SINGULAR_JACOBIAN_FOR_ORBIT_TYPE, propagationOrbitType);
  174.                     }
  175.                     minimumRel = FastMath.min(minimumRel, cartRelTol[i]);
  176.                 }
  177.                 absTol[6] = cartAbsTol[6];

  179.                 Arrays.fill(relTol, 0, 6, minimumRel);
  180.                 relTol[6] = cartRelTol[6];
  181.                 return new double[][] {absTol, relTol};
  182.             }
  183.         };
  184.     }

  186.     /**
  187.      * Defines a default tolerance provider. It is consistent with values from previous versions of Orekit obtained via other APIs.
  188.      *
  189.      * <p>
  190.      * The tolerances are only <em>orders of magnitude</em>, and integrator tolerances
  191.      * are only local estimates, not global ones. So some care must be taken when using
  192.      * these tolerances. Setting 1mm as a position error does NOT mean the tolerances
  193.      * will guarantee a 1mm error position after several orbits integration.
  194.      * </p>
  195.      *
  196.      * @param dP expected position error
  197.      * @return tolerances
  198.      */
  199.     static ToleranceProvider getDefaultToleranceProvider(final double dP) {
  200.         return new ToleranceProvider() {
  201.             @Override
  202.             public double[][] getTolerances(final Orbit referenceOrbit, final OrbitType propagationOrbitType,
  203.                                             final PositionAngleType positionAngleType) {
  204.                 // Cartesian case
  205.                 final double[] relTol = new double[7];
  206.                 final double[] cartAbsTol = new double[7];
  207.                 final double relPos = dP / referenceOrbit.getPosition().getNorm();
  208.                 Arrays.fill(relTol, 0, relTol.length, relPos);
  209.                 Arrays.fill(cartAbsTol, 0, 3, dP);
  210.                 // estimate the scalar velocity error
  211.                 final PVCoordinates pv = referenceOrbit.getPVCoordinates();
  212.                 final double r2 = pv.getPosition().getNormSq();
  213.                 final double v  = pv.getVelocity().getNorm();
  214.                 final double dV = referenceOrbit.getMu() * dP / (v * r2);
  215.                 Arrays.fill(cartAbsTol, 3, 6, dV);
  216.                 cartAbsTol[6] = DEFAULT_ABSOLUTE_MASS_TOLERANCE;

  218.                 if (propagationOrbitType == OrbitType.CARTESIAN) {
  219.                     return new double[][] {cartAbsTol, relTol};
  220.                 }

  222.                 // convert the orbit to the desired type
  223.                 final double[] absTol = cartAbsTol.clone();
  224.                 final double[][] jacobian = new double[6][6];
  225.                 final Orbit converted = propagationOrbitType.convertType(referenceOrbit);
  226.                 converted.getJacobianWrtCartesian(PositionAngleType.TRUE, jacobian);

  228.                 for (int i = 0; i < jacobian.length; ++i) {
  229.                     final double[] row = jacobian[i];
  230.                     absTol[i] = FastMath.abs(row[0]) * cartAbsTol[0] +
  231.                             FastMath.abs(row[1]) * cartAbsTol[1] +
  232.                             FastMath.abs(row[2]) * cartAbsTol[2] +
  233.                             FastMath.abs(row[3]) * cartAbsTol[3] +
  234.                             FastMath.abs(row[4]) * cartAbsTol[4] +
  235.                             FastMath.abs(row[5]) * cartAbsTol[5];
  236.                     if (Double.isNaN(absTol[i])) {
  237.                         throw new OrekitException(OrekitMessages.SINGULAR_JACOBIAN_FOR_ORBIT_TYPE, propagationOrbitType);
  238.                     }
  239.                 }

  241.                 return new double[][] {absTol, relTol};
  242.             }

  244.             @Override
  245.             public double[][] getTolerances(final Vector3D position, final Vector3D velocity) {
  246.                 final double[] absTol = new double[7];
  247.                 final double[] relTol = absTol.clone();
  248.                 final double relPos = dP / position.getNorm();
  249.                 Arrays.fill(relTol, 0, relTol.length, relPos);
  250.                 Arrays.fill(absTol, 0, 3, dP);
  251.                 final double dV = relPos * velocity.getNorm();
  252.                 Arrays.fill(absTol, 3, 6, dV);
  253.                 absTol[6] = DEFAULT_ABSOLUTE_MASS_TOLERANCE;
  254.                 return new double[][] {absTol, relTol};
  255.             }
  256.         };
  257.     }
  258. }
Created with JaCoCo 0.8.12.202403310830