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