1 /* Copyright 2002-2025 CS GROUP
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.semianalytical.dsst.utilities;
18
19 import java.util.Locale;
20
21 import org.hipparchus.analysis.differentiation.Gradient;
22 import org.hipparchus.random.UniformRandomGenerator;
23 import org.hipparchus.random.Well1024a;
24 import org.hipparchus.util.FastMath;
25 import org.hipparchus.util.Precision;
26 import org.junit.jupiter.api.Assertions;
27 import org.junit.jupiter.api.Test;
28
29 /** This class compares and tests different methods for computing value and 1st-order derivative of Jacobi polynomials.
30 *
31 * <p>It was added following issue <a href="https://gitlab.orekit.org/orekit/orekit/-/issues/1098">1098</a> when it was discovered that
32 * the Gradient-based method was considerably slowing down the computation of the tesseral contribution to the osculatin parameters.
33 *
34 * @author Maxime Journot
35 * @since 11.3.3
36 */
37 public class JacobiPolynomialsTest {
38
39 /** Threshold for test on value. */
40 private static final double epsValue = 1.e-20;
41
42 /** Threshold for relative difference on 1st order derivative of Jacobi polynomials. */
43 // Order 8 for tesseral terms
44 //private static final double epsDerivative = 4.44e-13;
45 // Order 32
46 // private static final double epsDerivative = 7.89e-7;
47 // Order 64
48 // private static final double epsDerivative = 78.44;
49
50 // From this threshold values it is obvious that there is still an issue in the computation of the Jacobi
51 // polynomials 1st order derivative at high order of the gravity field.
52 // The problem comes from the fact that some high-order coefficients of the polynomials are big (e.g. 1e18) while the low
53 // order can still be small(e.g. 1e-6)
54 // The numerous additions and multiplications that are done then are not stable. Meaning that, depending on the order of
55 // these operations the result will be different.
56 // In our case, since the "Gradient" method and the "optimized" method don't order the operations the same way, discrepancies arise...
57 // However it is impossible to know which method is "right", both may return results that are wrong.
58 // Some normalization along the biggest absolute value of the coefficients of the polynomial was attempted in
59 // "getValueAndDerivative". Results are different but once again, there is no way to know who's right or wrong here.
60
61 // Order 16 for tesseral terms
62 private static final double epsDerivative = 1.09e-12;
63
64
65 /** Test value and 1st-order derivative of Jacobi polynomials.
66 * <p>This test is designed to reproduce the behavior of a 16th order tesseral gravity field
67 * <p>Values for the polynomials are uniformly drawn between [-1, 1] to reproduce a random cos
68 */
69 @Test
70 public void testValueAndDerivative() {
71
72 // GIVEN
73 // -----
74
75 // Test like a 16x16 gravity field
76 final boolean print = false;
77 final int maxMdailyTesseralSP = 16;
78 final int maxEccPowTesseralSP = 4;
79
80 // Uniform random generator
81 final UniformRandomGenerator gen = new UniformRandomGenerator(new Well1024a(0x366a26b94e520f41l));
82
83 // Normalizing value for random numbers
84 // Using this will guarantee a uniform random number in [-1, 1]
85 // This is consistent with the input "x" of "JacobiPolynomials.getValueAndDerivative" which is:
86 // - γ = cos(2i) in DSSTTesseral (i = inclination in [0, π])
87 // - X = 1/√(1-e²) in DSSTThirdBody (e = eccentricity in [0, 1[)
88 final double sqrt3 = FastMath.sqrt(3.);
89
90 // WHEN
91 // ----
92
93 int nTest = 0;
94 for (int m = 1; m <= maxMdailyTesseralSP; m++) {
95
96 for (int s = 0; s <= maxEccPowTesseralSP; s++) {
97
98 final int nmin = FastMath.max(FastMath.max(2, m), FastMath.abs(s));
99 final int v = FastMath.abs(m - s);
100 final int w = FastMath.abs(m + s);
101
102 for (int n = nmin; n <= maxMdailyTesseralSP; n++) {
103
104 final int l = FastMath.min(n - m, n - FastMath.abs(s));
105 doTestValueAndDerivativeOnePoint(l, v, w, gen.nextNormalizedDouble() / sqrt3, print);
106 nTest++;
107 }
108 }
109 }
110
111 // THEN
112 // ----
113
114 // Check the number of tests performed
115 Assertions.assertEquals(668, nTest);
116 }
117
118 /** Test value and derivative for a given polynomial and input value.
119 * <p>Functions {@link JacobiPolynomials#getValueAndDerivative(int, int, int, double)} is compared to its "Gradient" counterpart.
120 * @param l degree of the polynomial
121 * @param v v value
122 * @param w w value
123 * @param x value to evaluate the polynomial on
124 * @param print print the differences
125 */
126 private void doTestValueAndDerivativeOnePoint(final int l, final int v, final int w, final double x, final boolean print) {
127
128 // WHEN
129 // ----
130
131 final double[] test = JacobiPolynomials.getValueAndDerivative(l, v, w, x);
132 final Gradient ref = JacobiPolynomials.getValue(l, v, w, Gradient.variable(1, 0, x));
133
134 // THEN
135 // ----
136
137 // Value: ref (Gradient) and test
138 final double refValue = ref.getValue();
139 final double testValue = test[0];
140
141 // Derivative: ref (Gradient) and test
142 final double refDer = ref.getPartialDerivative(0);
143 final double testDer = test[1];
144
145 // Threshold and diff: if refDer = 0 → absolute value, else → relative value
146 final double epsDer = FastMath.abs(refDer) > Precision.SAFE_MIN ? epsDerivative : epsValue;
147 final double diffDer = FastMath.abs(refDer) > Precision.SAFE_MIN ? FastMath.abs((testDer - refDer) / refDer) : FastMath.abs(testDer - refDer);
148
149 if (print) {
150 System.out.format(Locale.US,"l = %d%nv = %d%nw = %d%nx = %9.6f%n", l, v, w, x);
151 System.out.format(Locale.US,"\tvalue : %n\t\tref = %15.12e%n\t\ttest = %15.12e%n\t\tdiff = %15.12e%n",
152 refValue, testValue, FastMath.abs(refValue - testValue));
153 System.out.format(Locale.US,"\tderiv : %n\t\tref = %15.12e%n\t\ttest = %15.12e%n\t\tdiff = %15.12e%n",
154 refDer, testDer, diffDer);
155 }
156
157 // Test value directly (absolute difference, should be a clean 0.)
158 Assertions.assertEquals(refValue, testValue, epsValue);
159
160 // Test relative difference for derivative (if possible)
161 Assertions.assertEquals(0., diffDer, epsDer);
162 }
163 }