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;
18  
19  import org.hipparchus.CalculusFieldElement;
20  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
21  import org.hipparchus.geometry.euclidean.threed.Vector3D;
22  import org.hipparchus.util.FastMath;
23  import org.orekit.orbits.CartesianOrbit;
24  import org.orekit.orbits.FieldCartesianOrbit;
25  import org.orekit.propagation.numerical.FieldNumericalPropagator;
26  import org.orekit.propagation.numerical.NumericalPropagator;
27  import org.orekit.utils.AbsolutePVCoordinates;
28  import org.orekit.utils.FieldAbsolutePVCoordinates;
29  
30  import java.util.Arrays;
31  
32  
33  /**
34   * Interface to define integration tolerances for adaptive schemes (like the embedded Runge-Kutta ones) propagating
35   * the position-velocity vector and the mass, for a total of 7 primary dependent variables (in that order).
36   * 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.
37   *
38   * @see NumericalPropagator
39   * @see FieldNumericalPropagator
40   * @see CartesianToleranceProvider
41   * @since 13.0
42   * @author Romain Serra
43   */
44  public interface CartesianToleranceProvider {
45  
46      /** Default absolute tolerance for mass integration. */
47      double DEFAULT_ABSOLUTE_MASS_TOLERANCE = 1e-6;
48  
49      /**
50       * Retrieve the integration tolerances given reference position and velocity vectors.
51       * @param position reference position vector
52       * @param velocity reference velocity vector
53       * @return absolute and relative tolerances
54       */
55      double[][] getTolerances(Vector3D position, Vector3D velocity);
56  
57      /**
58       * Retrieve the integration tolerances given reference position and velocity Field vectors.
59       * @param position reference position vector
60       * @param velocity reference velocity vector
61       * @param <T> field type
62       * @return absolute and relative tolerances
63       */
64      default <T extends CalculusFieldElement<T>> double[][] getTolerances(FieldVector3D<T> position,
65                                                                           FieldVector3D<T> velocity) {
66          return getTolerances(position.toVector3D(), velocity.toVector3D());
67      }
68  
69      /**
70       * Retrieve the integration tolerances given a reference Cartesian orbit.
71       * @param cartesianOrbit reference orbit
72       * @return absolute and relative tolerances
73       */
74      default double[][] getTolerances(CartesianOrbit cartesianOrbit) {
75          return getTolerances(cartesianOrbit.getPosition(), cartesianOrbit.getPVCoordinates().getVelocity());
76      }
77  
78      /**
79       * Retrieve the integration tolerances given a reference Cartesian orbit.
80       * @param cartesianOrbit reference orbit
81       * @param <T> field type
82       * @return absolute and relative tolerances
83       */
84      default <T extends CalculusFieldElement<T>> double[][] getTolerances(FieldCartesianOrbit<T> cartesianOrbit) {
85          return getTolerances(cartesianOrbit.toOrbit());
86      }
87  
88      /**
89       * Retrieve the integration tolerances given a reference absolute position-velocity vector.
90       * @param absolutePVCoordinates reference position-velocity
91       * @return absolute and relative tolerances
92       */
93      default double[][] getTolerances(AbsolutePVCoordinates absolutePVCoordinates) {
94          return getTolerances(absolutePVCoordinates.getPosition(), absolutePVCoordinates.getPVCoordinates().getVelocity());
95      }
96  
97      /**
98       * Retrieve the integration tolerances given a reference absolute position-velocity vector.
99       * @param absolutePVCoordinates reference position-velocity
100      * @param <T> field type
101      * @return absolute and relative tolerances
102      */
103     default <T extends CalculusFieldElement<T>> double[][] getTolerances(FieldAbsolutePVCoordinates<T> absolutePVCoordinates) {
104         return getTolerances(absolutePVCoordinates.toAbsolutePVCoordinates());
105     }
106 
107     /**
108      * Build a provider based on expected errors for position, velocity and mass respectively.
109      *
110      * <p>
111      * The tolerances are only <em>orders of magnitude</em>, and integrator tolerances
112      * are only local estimates, not global ones. So some care must be taken when using
113      * these tolerances. Setting 1mm as a position error does NOT mean the tolerances
114      * will guarantee a 1mm error position after several orbits integration.
115      * </p>
116      *
117      * @param dP expected position error
118      * @param dV expected velocity error
119      * @param dM expected mass error
120      * @return tolerance provider
121      */
122     static CartesianToleranceProvider of(final double dP, final double dV, final double dM) {
123         return (position, velocity) -> {
124             final double[] absTol = new double[7];
125             final double[] relTol = absTol.clone();
126             Arrays.fill(absTol, 0, 3, dP);
127             Arrays.fill(absTol, 3, 6, dV);
128             absTol[6] = dM;
129             final double relPos = dP / position.getNorm();
130             Arrays.fill(relTol, 0, 3, relPos);
131             final double relVel = dV / velocity.getNorm();
132             Arrays.fill(relTol, 3, 6, relVel);
133             relTol[6] = FastMath.min(relPos, relVel);
134             return new double[][] { absTol, relTol };
135         };
136     }
137 
138     /**
139      * Build a provider based on expected errors for position only.
140      *
141      * <p>
142      * The tolerances are only <em>orders of magnitude</em>, and integrator tolerances
143      * are only local estimates, not global ones. So some care must be taken when using
144      * these tolerances. Setting 1mm as a position error does NOT mean the tolerances
145      * will guarantee a 1mm error position after several orbits integration.
146      * </p>
147      *
148      * @param dP expected position error
149      * @return tolerance provider
150      */
151     static CartesianToleranceProvider of(final double dP) {
152         return (position, velocity) -> {
153             final double[] absTol = new double[7];
154             final double[] relTol = absTol.clone();
155             final double relPos = dP / position.getNorm();
156             final double dV = relPos * velocity.getNorm();
157             Arrays.fill(absTol, 0, 3, dP);
158             Arrays.fill(absTol, 3, 6, dV);
159             absTol[6] = DEFAULT_ABSOLUTE_MASS_TOLERANCE;
160             Arrays.fill(relTol, 0, 7, relPos);
161             return new double[][] { absTol, relTol };
162         };
163     }
164 }