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 org.hipparchus.analysis.polynomials.PolynomialFunction;
20 import org.hipparchus.util.FastMath;
21
22 import java.util.ArrayList;
23 import java.util.List;
24 import java.util.Map;
25 import java.util.SortedMap;
26 import java.util.TreeMap;
27
28 /**
29 * Implementation of the Modified Newcomb Operators.
30 *
31 * <p> From equations 2.7.3 - (12)(13) of the Danielson paper, those operators
32 * are defined as:
33 *
34 * <p>
35 * 4(ρ + σ)Y<sub>ρ,σ</sub><sup>n,s</sup> = <br>
36 * 2(2s - n)Y<sub>ρ-1,σ</sub><sup>n,s+1</sup> + (s - n)Y<sub>ρ-2,σ</sub><sup>n,s+2</sup> <br>
37 * - 2(2s + n)Y<sub>ρ,σ-1</sub><sup>n,s-1</sup> - (s+n)Y<sub>ρ,σ-2</sub><sup>n,s-2</sup> <br>
38 * + 2(2ρ + 2σ + 2 + 3n)Y<sub>ρ-1,σ-1</sub><sup>n,s</sup>
39 *
40 * <p> Initialization is given by : Y<sub>0,0</sub><sup>n,s</sup> = 1
41 *
42 * <p> Internally, the Modified Newcomb Operators are stored as an array of
43 * {@link PolynomialFunction} :
44 *
45 * <p> Y<sub>ρ,σ</sub><sup>n,s</sup> = P<sub>k0</sub> + P<sub>k1</sub>n + ... +
46 * P<sub>kj</sub>n<sup>j</sup>
47 *
48 * <p> where the P<sub>kj</sub> are given by
49 *
50 * <p> P<sub>kj</sub> = ∑<sub>j=0;ρ</sub> a<sub>j</sub>s<sup>j</sup>
51 *
52 * @author Romain Di Costanzo
53 * @author Pascal Parraud
54 */
55 public class NewcombOperators {
56
57 /** Storage map. */
58 private static final Map<NewKey, Double> MAP = new TreeMap<>((k1, k2) -> {
59 if (k1.n == k2.n) {
60 if (k1.s == k2.s) {
61 if (k1.rho == k2.rho) {
62 return Integer.compare(k1.sigma, k2.sigma);
63 } else if (k1.rho < k2.rho) {
64 return -1;
65 } else {
66 return 1;
67 }
68 } else if (k1.s < k2.s) {
69 return -1;
70 } else {
71 return 1;
72 }
73 } else if (k1.n < k2.n) {
74 return -1;
75 }
76 return 1;
77 });
78
79 /** Private constructor as class is a utility.
80 */
81 private NewcombOperators() {
82 }
83
84 /** Get the Newcomb operator evaluated at n, s, ρ, σ.
85 * <p>
86 * This method is guaranteed to be thread-safe
87 * </p>
88 * @param rho ρ index
89 * @param sigma σ index
90 * @param n n index
91 * @param s s index
92 * @return Y<sub>ρ,σ</sub><sup>n,s</sup>
93 */
94 public static double getValue(final int rho, final int sigma, final int n, final int s) {
95
96 final NewKey key = new NewKey(n, s, rho, sigma);
97 synchronized (MAP) {
98 if (MAP.containsKey(key)) {
99 return MAP.get(key);
100 }
101 }
102
103 // Get the Newcomb polynomials for the given rho and sigma
104 final List<PolynomialFunction> polynomials = PolynomialsGenerator.getPolynomials(rho, sigma);
105
106 // Compute the value from the list of polynomials for the given n and s
107 double nPower = 1.;
108 double value = 0.0;
109 for (final PolynomialFunction polynomial : polynomials) {
110 value += polynomial.value(s) * nPower;
111 nPower = n * nPower;
112 }
113 synchronized (MAP) {
114 MAP.put(key, value);
115 }
116
117 return value;
118
119 }
120
121 /** Generator for Newcomb polynomials. */
122 private static class PolynomialsGenerator {
123
124 /** Polynomials storage. */
125 private static final SortedMap<Couple, List<PolynomialFunction>> POLYNOMIALS =
126 new TreeMap<>((c1, c2) -> {
127 if (c1.rho == c2.rho) {
128 if (c1.sigma < c2.sigma) {
129 return -1;
130 } else if (c1.sigma == c2.sigma) {
131 return 0;
132 }
133 } else if (c1.rho < c2.rho) {
134 return -1;
135 }
136 return 1;
137 });
138
139 /** Private constructor as class is a utility.
140 */
141 private PolynomialsGenerator() {
142 }
143
144 /** Get the list of polynomials representing the Newcomb Operator for the (ρ,σ) couple.
145 * <p>
146 * This method is guaranteed to be thread-safe
147 * </p>
148 * @param rho ρ value
149 * @param sigma σ value
150 * @return Polynomials representing the Newcomb Operator for the (ρ,σ) couple.
151 */
152 private static List<PolynomialFunction> getPolynomials(final int rho, final int sigma) {
153
154 final Couple couple = new Couple(rho, sigma);
155
156 synchronized (POLYNOMIALS) {
157
158 if (POLYNOMIALS.isEmpty()) {
159 // Initialize lists
160 final List<PolynomialFunction> l00 = new ArrayList<>();
161 final List<PolynomialFunction> l01 = new ArrayList<>();
162 final List<PolynomialFunction> l10 = new ArrayList<>();
163 final List<PolynomialFunction> l11 = new ArrayList<>();
164
165 // Y(rho = 0, sigma = 0) = 1
166 l00.add(new PolynomialFunction(new double[] {
167 1.
168 }));
169 // Y(rho = 0, sigma = 1) = -s - n/2
170 l01.add(new PolynomialFunction(new double[] {
171 0, -1.
172 }));
173 l01.add(new PolynomialFunction(new double[] {
174 -0.5
175 }));
176 // Y(rho = 1, sigma = 0) = s - n/2
177 l10.add(new PolynomialFunction(new double[] {
178 0, 1.
179 }));
180 l10.add(new PolynomialFunction(new double[] {
181 -0.5
182 }));
183 // Y(rho = 1, sigma = 1) = 3/2 - s² + 5n/4 + n²/4
184 l11.add(new PolynomialFunction(new double[] {
185 1.5, 0., -1.
186 }));
187 l11.add(new PolynomialFunction(new double[] {
188 1.25
189 }));
190 l11.add(new PolynomialFunction(new double[] {
191 0.25
192 }));
193
194 // Initialize polynomials
195 POLYNOMIALS.put(new Couple(0, 0), l00);
196 POLYNOMIALS.put(new Couple(0, 1), l01);
197 POLYNOMIALS.put(new Couple(1, 0), l10);
198 POLYNOMIALS.put(new Couple(1, 1), l11);
199
200 }
201
202 // If order hasn't been computed yet, update the Newcomb polynomials
203 if (!POLYNOMIALS.containsKey(couple)) {
204 PolynomialsGenerator.computeFor(rho, sigma);
205 }
206
207 return POLYNOMIALS.get(couple);
208
209 }
210 }
211
212 /** Compute the Modified Newcomb Operators up to a given (ρ, σ) couple.
213 * <p>
214 * The recursive computation uses equation 2.7.3-(12) of the Danielson paper.
215 * </p>
216 * @param rho ρ value to reach
217 * @param sigma σ value to reach
218 */
219 private static void computeFor(final int rho, final int sigma) {
220
221 // Initialize result :
222 List<PolynomialFunction> result = new ArrayList<>();
223
224 // Get the coefficient from the recurrence relation
225 final Map<Integer, List<PolynomialFunction>> map = generateRecurrenceCoefficients(rho, sigma);
226
227 // Compute (s - n) * Y[rho - 2, sigma][n, s + 2]
228 if (rho >= 2) {
229 final List<PolynomialFunction> poly = map.get(0);
230 final List<PolynomialFunction> list = getPolynomials(rho - 2, sigma);
231 result = multiplyPolynomialList(poly, shiftList(list, 2));
232 }
233
234 // Compute 2(2rho + 2sigma + 2 + 3n) * Y[rho - 1, sigma - 1][n, s]
235 if (rho >= 1 && sigma >= 1) {
236 final List<PolynomialFunction> poly = map.get(1);
237 final List<PolynomialFunction> list = getPolynomials(rho - 1, sigma - 1);
238 result = sumPolynomialList(result, multiplyPolynomialList(poly, list));
239 }
240
241 // Compute 2(2s - n) * Y[rho - 1, sigma][n, s + 1]
242 if (rho >= 1) {
243 final List<PolynomialFunction> poly = map.get(2);
244 final List<PolynomialFunction> list = getPolynomials(rho - 1, sigma);
245 result = sumPolynomialList(result, multiplyPolynomialList(poly, shiftList(list, 1)));
246 }
247
248 // Compute -(s + n) * Y[rho, sigma - 2][n, s - 2]
249 if (sigma >= 2) {
250 final List<PolynomialFunction> poly = map.get(3);
251 final List<PolynomialFunction> list = getPolynomials(rho, sigma - 2);
252 result = sumPolynomialList(result, multiplyPolynomialList(poly, shiftList(list, -2)));
253 }
254
255 // Compute -2(2s + n) * Y[rho, sigma - 1][n, s - 1]
256 if (sigma >= 1) {
257 final List<PolynomialFunction> poly = map.get(4);
258 final List<PolynomialFunction> list = getPolynomials(rho, sigma - 1);
259 result = sumPolynomialList(result, multiplyPolynomialList(poly, shiftList(list, -1)));
260 }
261
262 // Save polynomials for current (rho, sigma) couple
263 final Couple couple = new Couple(rho, sigma);
264 POLYNOMIALS.put(couple, result);
265 }
266
267 /** Multiply two lists of polynomials defined as the internal representation of the Newcomb Operator.
268 * <p>
269 * Let's call R<sub>s</sub>(n) the result returned by the method :
270 * <pre>
271 * R<sub>s</sub>(n) = (P<sub>s₀</sub> + P<sub>s₁</sub>n + ... + P<sub>s<sub>j</sub></sub>n<sup>j</sup>) *(Q<sub>s₀</sub> + Q<sub>s₁</sub>n + ... + Q<sub>s<sub>k</sub></sub>n<sup>k</sup>)
272 * </pre>
273 *
274 * where the P<sub>s<sub>j</sub></sub> and Q<sub>s<sub>k</sub></sub> are polynomials in s,
275 * s being the index of the Y<sub>ρ,σ</sub><sup>n,s</sup> function
276 *
277 * @param poly1 first list of polynomials
278 * @param poly2 second list of polynomials
279 * @return R<sub>s</sub>(n) as a list of {@link PolynomialFunction}
280 */
281 private static List<PolynomialFunction> multiplyPolynomialList(final List<PolynomialFunction> poly1,
282 final List<PolynomialFunction> poly2) {
283 // Initialize the result list of polynomial function
284 final List<PolynomialFunction> result = new ArrayList<>();
285 initializeListOfPolynomials(poly1.size() + poly2.size() - 1, result);
286
287 int i = 0;
288 // Iterate over first polynomial list
289 for (PolynomialFunction f1 : poly1) {
290 // Iterate over second polynomial list
291 for (int j = i; j < poly2.size() + i; j++) {
292 final PolynomialFunction p2 = poly2.get(j - i);
293 // Get previous polynomial for current 'n' order
294 final PolynomialFunction previousP2 = result.get(j);
295 // Replace the current order by summing the product of both of the polynomials
296 result.set(j, previousP2.add(f1.multiply(p2)));
297 }
298 // shift polynomial order in 'n'
299 i++;
300 }
301 return result;
302 }
303
304 /** Sum two lists of {@link PolynomialFunction}.
305 * @param poly1 first list
306 * @param poly2 second list
307 * @return the summed list
308 */
309 private static List<PolynomialFunction> sumPolynomialList(final List<PolynomialFunction> poly1,
310 final List<PolynomialFunction> poly2) {
311 // identify the lowest degree polynomial
312 final int lowLength = FastMath.min(poly1.size(), poly2.size());
313 final int highLength = FastMath.max(poly1.size(), poly2.size());
314 // Initialize the result list of polynomial function
315 final List<PolynomialFunction> result = new ArrayList<>();
316 initializeListOfPolynomials(highLength, result);
317
318 for (int i = 0; i < lowLength; i++) {
319 // Add polynomials by increasing order of 'n'
320 result.set(i, poly1.get(i).add(poly2.get(i)));
321 }
322 // Complete the list if lists are of different size:
323 for (int i = lowLength; i < highLength; i++) {
324 if (poly1.size() < poly2.size()) {
325 result.set(i, poly2.get(i));
326 } else {
327 result.set(i, poly1.get(i));
328 }
329 }
330 return result;
331 }
332
333 /** Initialize an empty list of polynomials.
334 * @param i order
335 * @param result list into which polynomials should be added
336 */
337 private static void initializeListOfPolynomials(final int i,
338 final List<PolynomialFunction> result) {
339 for (int k = 0; k < i; k++) {
340 result.add(new PolynomialFunction(new double[] {0.}));
341 }
342 }
343
344 /** Shift a list of {@link PolynomialFunction}.
345 *
346 * @param polynomialList list of {@link PolynomialFunction}
347 * @param shift shift value
348 * @return new list of shifted {@link PolynomialFunction}
349 */
350 private static List<PolynomialFunction> shiftList(final List<PolynomialFunction> polynomialList,
351 final int shift) {
352 final List<PolynomialFunction> shiftedList = new ArrayList<>();
353 for (PolynomialFunction function : polynomialList) {
354 shiftedList.add(new PolynomialFunction(shift(function.getCoefficients(), shift)));
355 }
356 return shiftedList;
357 }
358
359 /**
360 * Compute the coefficients of the polynomial \(P_s(x)\)
361 * whose values at point {@code x} will be the same as the those from the
362 * original polynomial \(P(x)\) when computed at {@code x + shift}.
363 * <p>
364 * More precisely, let \(\Delta = \) {@code shift} and let
365 * \(P_s(x) = P(x + \Delta)\). The returned array
366 * consists of the coefficients of \(P_s\). So if \(a_0, ..., a_{n-1}\)
367 * are the coefficients of \(P\), then the returned array
368 * \(b_0, ..., b_{n-1}\) satisfies the identity
369 * \(\sum_{i=0}^{n-1} b_i x^i = \sum_{i=0}^{n-1} a_i (x + \Delta)^i\) for all \(x\).
370 * </p>
371 * <p>
372 * This method is a modified version of the method with the same name
373 * in Hipparchus {@code PolynomialsUtils} class, simply changing
374 * computation of binomial coefficients so degrees higher than 66 can be used.
375 * </p>
376 *
377 * @param coefficients Coefficients of the original polynomial.
378 * @param shift Shift value.
379 * @return the coefficients \(b_i\) of the shifted
380 * polynomial.
381 */
382 public static double[] shift(final double[] coefficients,
383 final double shift) {
384 final int dp1 = coefficients.length;
385 final double[] newCoefficients = new double[dp1];
386
387 // Pascal triangle.
388 final double[][] coeff = new double[dp1][dp1];
389 coeff[0][0] = 1;
390 for (int i = 1; i < dp1; i++) {
391 coeff[i][0] = 1;
392 for (int j = 1; j < i; j++) {
393 coeff[i][j] = coeff[i - 1][j - 1] + coeff[i - 1][j];
394 }
395 coeff[i][i] = 1;
396 }
397
398 // First polynomial coefficient.
399 double shiftI = 1;
400 for (double coefficient : coefficients) {
401 newCoefficients[0] += coefficient * shiftI;
402 shiftI *= shift;
403 }
404
405 // Superior order.
406 final int d = dp1 - 1;
407 for (int i = 0; i < d; i++) {
408 double shiftJmI = 1;
409 for (int j = i; j < d; j++) {
410 newCoefficients[i + 1] += coeff[j + 1][j - i] * coefficients[j + 1] * shiftJmI;
411 shiftJmI *= shift;
412 }
413 }
414
415 return newCoefficients;
416 }
417
418 /** Generate recurrence coefficients.
419 *
420 * @param rho ρ value
421 * @param sigma σ value
422 * @return recurrence coefficients
423 */
424 private static Map<Integer, List<PolynomialFunction>> generateRecurrenceCoefficients(final int rho, final int sigma) {
425 final double den = 1. / (4. * (rho + sigma));
426 final double denx2 = 2. * den;
427 final double denx4 = 4. * den;
428 // Initialization :
429 final Map<Integer, List<PolynomialFunction>> list = new TreeMap<>();
430 final List<PolynomialFunction> poly0 = new ArrayList<>();
431 final List<PolynomialFunction> poly1 = new ArrayList<>();
432 final List<PolynomialFunction> poly2 = new ArrayList<>();
433 final List<PolynomialFunction> poly3 = new ArrayList<>();
434 final List<PolynomialFunction> poly4 = new ArrayList<>();
435 // (s - n)
436 poly0.add(new PolynomialFunction(new double[] {0., den}));
437 poly0.add(new PolynomialFunction(new double[] {-den}));
438 // 2(2 * rho + 2 * sigma + 2 + 3*n)
439 poly1.add(new PolynomialFunction(new double[] {1. + denx4}));
440 poly1.add(new PolynomialFunction(new double[] {denx2 + denx4}));
441 // 2(2s - n)
442 poly2.add(new PolynomialFunction(new double[] {0., denx4}));
443 poly2.add(new PolynomialFunction(new double[] {-denx2}));
444 // - (s + n)
445 poly3.add(new PolynomialFunction(new double[] {0., -den}));
446 poly3.add(new PolynomialFunction(new double[] {-den}));
447 // - 2(2s + n)
448 poly4.add(new PolynomialFunction(new double[] {0., -denx4}));
449 poly4.add(new PolynomialFunction(new double[] {-denx2}));
450
451 // Fill the map :
452 list.put(0, poly0);
453 list.put(1, poly1);
454 list.put(2, poly2);
455 list.put(3, poly3);
456 list.put(4, poly4);
457 return list;
458 }
459
460 }
461
462 /** Private class to define a couple of value. */
463 private static class Couple {
464
465 /** first couple value. */
466 private final int rho;
467
468 /** second couple value. */
469 private final int sigma;
470
471 /** Constructor.
472 * @param rho first couple value
473 * @param sigma second couple value
474 */
475 private Couple(final int rho, final int sigma) {
476 this.rho = rho;
477 this.sigma = sigma;
478 }
479
480 }
481
482 /** Newcomb operator's key. */
483 private static class NewKey {
484
485 /** n value. */
486 private final int n;
487
488 /** s value. */
489 private final int s;
490
491 /** ρ value. */
492 private final int rho;
493
494 /** σ value. */
495 private final int sigma;
496
497 /** Simpleconstructor.
498 * @param n n
499 * @param s s
500 * @param rho ρ
501 * @param sigma σ
502 */
503 NewKey(final int n, final int s, final int rho, final int sigma) {
504 this.n = n;
505 this.s = s;
506 this.rho = rho;
507 this.sigma = sigma;
508 }
509
510 }
511
512 }