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.hansen;
18
19 import java.lang.reflect.Array;
20
21 import org.hipparchus.CalculusFieldElement;
22 import org.hipparchus.Field;
23 import org.hipparchus.analysis.differentiation.FieldGradient;
24 import org.hipparchus.analysis.polynomials.PolynomialFunction;
25 import org.hipparchus.util.FastMath;
26 import org.hipparchus.util.MathArrays;
27 import org.orekit.propagation.semianalytical.dsst.utilities.NewcombOperators;
28
29 /**
30 * Hansen coefficients K(t,n,s) for t!=0 and n < 0.
31 * <p>
32 * Implements Collins 4-236 or Danielson 2.7.3-(9) for Hansen Coefficients and
33 * Collins 4-240 for derivatives. The recursions are transformed into
34 * composition of linear transformations to obtain the associated polynomials
35 * for coefficients and their derivatives - see Petre's paper
36 *
37 * @author Petre Bazavan
38 * @author Lucian Barbulescu
39 * @author Bryan Cazabonne
40 */
41 public class FieldHansenTesseralLinear <T extends CalculusFieldElement<T>> {
42
43 /** The number of coefficients that will be computed with a set of roots. */
44 private static final int SLICE = 10;
45
46 /**
47 * The first vector of polynomials associated to Hansen coefficients and
48 * derivatives.
49 */
50 private PolynomialFunction[][] mpvec;
51
52 /** The second vector of polynomials associated only to derivatives. */
53 private PolynomialFunction[][] mpvecDeriv;
54
55 /** The Hansen coefficients used as roots. */
56 private T[][] hansenRoot;
57
58 /** The derivatives of the Hansen coefficients used as roots. */
59 private T[][] hansenDerivRoot;
60
61 /** The minimum value for the order. */
62 private int Nmin;
63
64 /** The index of the initial condition, Petre's paper. */
65 private int N0;
66
67 /** The s coefficient. */
68 private int s;
69
70 /** The j coefficient. */
71 private int j;
72
73 /** The number of slices needed to compute the coefficients. */
74 private int numSlices;
75
76 /**
77 * The offset used to identify the polynomial that corresponds to a negative.
78 * value of n in the internal array that starts at 0
79 */
80 private int offset;
81
82 /** The objects used to calculate initial data by means of Newcomb operators. */
83 private FieldHansenCoefficientsBySeries<T>[] hansenInit;
84
85 /**
86 * Constructor.
87 *
88 * @param nMax the maximum (absolute) value of n parameter
89 * @param s s parameter
90 * @param j j parameter
91 * @param n0 the minimum (absolute) value of n
92 * @param maxHansen maximum power of the eccentricity to use in Hansen coefficient Kernel expansion.
93 * @param field field used by default
94 */
95 @SuppressWarnings("unchecked")
96 public FieldHansenTesseralLinear(final int nMax, final int s, final int j, final int n0,
97 final int maxHansen, final Field<T> field) {
98 //Initialize the fields
99 this.offset = nMax + 1;
100 this.Nmin = -nMax - 1;
101 this.N0 = -n0 - 4;
102 this.s = s;
103 this.j = j;
104
105 final int maxRoots = FastMath.min(4, N0 - Nmin + 4);
106 //Ensure that only the needed terms are computed
107 this.hansenInit = (FieldHansenCoefficientsBySeries<T>[]) Array.newInstance(FieldHansenCoefficientsBySeries.class, maxRoots);
108 for (int i = 0; i < maxRoots; i++) {
109 this.hansenInit[i] = new FieldHansenCoefficientsBySeries<>(N0 - i + 3, s, j, maxHansen, field);
110 }
111
112 // The first 4 values are computed with series. No linear combination is needed.
113 final int size = N0 - Nmin;
114 this.numSlices = (int) FastMath.max(FastMath.ceil(((double) size) / SLICE), 1);
115 hansenRoot = MathArrays.buildArray(field, numSlices, 4);
116 hansenDerivRoot = MathArrays.buildArray(field, numSlices, 4);
117 if (size > 0) {
118 mpvec = new PolynomialFunction[size][];
119 mpvecDeriv = new PolynomialFunction[size][];
120
121 // Prepare the database of the associated polynomials
122 generatePolynomials();
123 }
124
125 }
126
127 /**
128 * Compute polynomial coefficient a.
129 *
130 * <p>
131 * It is used to generate the coefficient for K<sub>j</sub><sup>-n, s</sup> when computing K<sub>j</sub><sup>-n-1, s</sup>
132 * and the coefficient for dK<sub>j</sub><sup>-n, s</sup> / de² when computing dK<sub>j</sub><sup>-n-1, s</sup> / de²
133 * </p>
134 *
135 * <p>
136 * See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
137 * </p>
138 *
139 * @param mnm1 -n-1
140 * @return the polynomial
141 */
142 private PolynomialFunction a(final int mnm1) {
143 // Collins 4-236, Danielson 2.7.3-(9)
144 final double r1 = (mnm1 + 2.) * (2. * mnm1 + 5.);
145 final double r2 = (2. + mnm1 + s) * (2. + mnm1 - s);
146 return new PolynomialFunction(new double[] {
147 0.0, 0.0, r1 / r2
148 });
149 }
150
151 /**
152 * Compute polynomial coefficient b.
153 *
154 * <p>
155 * It is used to generate the coefficient for K<sub>j</sub><sup>-n+1, s</sup> when computing K<sub>j</sub><sup>-n-1, s</sup>
156 * and the coefficient for dK<sub>j</sub><sup>-n+1, s</sup> / de² when computing dK<sub>j</sub><sup>-n-1, s</sup> / de²
157 * </p>
158 *
159 * <p>
160 * See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
161 * </p>
162 *
163 * @param mnm1 -n-1
164 * @return the polynomial
165 */
166 private PolynomialFunction b(final int mnm1) {
167 // Collins 4-236, Danielson 2.7.3-(9)
168 final double r2 = (2. + mnm1 + s) * (2. + mnm1 - s);
169 final double d1 = (mnm1 + 3.) * 2. * j * s / (r2 * (mnm1 + 4.));
170 final double d2 = (mnm1 + 3.) * (mnm1 + 2.) / r2;
171 return new PolynomialFunction(new double[] {
172 0.0, -d1, -d2
173 });
174 }
175
176 /**
177 * Compute polynomial coefficient c.
178 *
179 * <p>
180 * It is used to generate the coefficient for K<sub>j</sub><sup>-n+3, s</sup> when computing K<sub>j</sub><sup>-n-1, s</sup>
181 * and the coefficient for dK<sub>j</sub><sup>-n+3, s</sup> / de² when computing dK<sub>j</sub><sup>-n-1, s</sup> / de²
182 * </p>
183 *
184 * <p>
185 * See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
186 * </p>
187 *
188 * @param mnm1 -n-1
189 * @return the polynomial
190 */
191 private PolynomialFunction c(final int mnm1) {
192 // Collins 4-236, Danielson 2.7.3-(9)
193 final double r1 = j * j * (mnm1 + 2.);
194 final double r2 = (mnm1 + 4.) * (2. + mnm1 + s) * (2. + mnm1 - s);
195
196 return new PolynomialFunction(new double[] {
197 0.0, 0.0, r1 / r2
198 });
199 }
200
201 /**
202 * Compute polynomial coefficient d.
203 *
204 * <p>
205 * It is used to generate the coefficient for K<sub>j</sub><sup>-n-1, s</sup> / dχ when computing dK<sub>j</sub><sup>-n-1, s</sup> / de²
206 * </p>
207 *
208 * <p>
209 * See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
210 * </p>
211 *
212 * @param mnm1 -n-1
213 * @return the polynomial
214 */
215 private PolynomialFunction d(final int mnm1) {
216 // Collins 4-236, Danielson 2.7.3-(9)
217 return new PolynomialFunction(new double[] {
218 0.0, 0.0, 1.0
219 });
220 }
221
222 /**
223 * Compute polynomial coefficient f.
224 *
225 * <p>
226 * It is used to generate the coefficient for K<sub>j</sub><sup>-n+1, s</sup> / dχ when computing dK<sub>j</sub><sup>-n-1, s</sup> / de²
227 * </p>
228 *
229 * <p>
230 * See Danielson 2.7.3-(9) and Collins 4-236 and 4-240
231 * </p>
232 *
233 * @param n index
234 * @return the polynomial
235 */
236 private PolynomialFunction f(final int n) {
237 // Collins 4-236, Danielson 2.7.3-(9)
238 final double r1 = (n + 3.0) * j * s;
239 final double r2 = (n + 4.0) * (2.0 + n + s) * (2.0 + n - s);
240 return new PolynomialFunction(new double[] {
241 0.0, 0.0, 0.0, r1 / r2
242 });
243 }
244
245 /**
246 * Generate the polynomials needed in the linear transformation.
247 *
248 * <p>
249 * See Petre's paper
250 * </p>
251 */
252 private void generatePolynomials() {
253
254
255 // Initialization of the matrices for linear transformations
256 // The final configuration of these matrices are obtained by composition
257 // of linear transformations
258
259 // The matrix of polynomials associated to Hansen coefficients, Petre's
260 // paper
261 PolynomialFunctionMatrix A = HansenUtilities.buildIdentityMatrix4();
262
263 // The matrix of polynomials associated to derivatives, Petre's paper
264 final PolynomialFunctionMatrix B = HansenUtilities.buildZeroMatrix4();
265 PolynomialFunctionMatrix D = HansenUtilities.buildZeroMatrix4();
266 final PolynomialFunctionMatrix a = HansenUtilities.buildZeroMatrix4();
267
268 // The matrix of the current linear transformation
269 a.setMatrixLine(0, new PolynomialFunction[] {
270 HansenUtilities.ZERO, HansenUtilities.ONE, HansenUtilities.ZERO, HansenUtilities.ZERO
271 });
272 a.setMatrixLine(1, new PolynomialFunction[] {
273 HansenUtilities.ZERO, HansenUtilities.ZERO, HansenUtilities.ONE, HansenUtilities.ZERO
274 });
275 a.setMatrixLine(2, new PolynomialFunction[] {
276 HansenUtilities.ZERO, HansenUtilities.ZERO, HansenUtilities.ZERO, HansenUtilities.ONE
277 });
278 // The generation process
279 int index;
280 int sliceCounter = 0;
281 for (int i = N0 - 1; i > Nmin - 1; i--) {
282 index = i + this.offset;
283 // The matrix of the current linear transformation is updated
284 // Petre's paper
285 a.setMatrixLine(3, new PolynomialFunction[] {
286 c(i), HansenUtilities.ZERO, b(i), a(i)
287 });
288
289 // composition of the linear transformations to calculate
290 // the polynomials associated to Hansen coefficients
291 // Petre's paper
292 A = A.multiply(a);
293 // store the polynomials for Hansen coefficients
294 mpvec[index] = A.getMatrixLine(3);
295 // composition of the linear transformations to calculate
296 // the polynomials associated to derivatives
297 // Petre's paper
298 D = D.multiply(a);
299
300 //Update the B matrix
301 B.setMatrixLine(3, new PolynomialFunction[] {
302 HansenUtilities.ZERO, f(i),
303 HansenUtilities.ZERO, d(i)
304 });
305 D = D.add(A.multiply(B));
306
307 // store the polynomials for Hansen coefficients from the
308 // expressions of derivatives
309 mpvecDeriv[index] = D.getMatrixLine(3);
310
311 if (++sliceCounter % SLICE == 0) {
312 // Re-Initialisation of matrix for linear transformmations
313 // The final configuration of these matrix are obtained by composition
314 // of linear transformations
315 A = HansenUtilities.buildIdentityMatrix4();
316 D = HansenUtilities.buildZeroMatrix4();
317 }
318 }
319 }
320
321 /**
322 * Compute the values for the first four coefficients and their derivatives by means of series.
323 *
324 * @param e2 e²
325 * @param chi Χ
326 * @param chi2 Χ²
327 */
328 public void computeInitValues(final T e2, final T chi, final T chi2) {
329 // compute the values for n, n+1, n+2 and n+3 by series
330 // See Danielson 2.7.3-(10)
331 //Ensure that only the needed terms are computed
332 final int maxRoots = FastMath.min(4, N0 - Nmin + 4);
333 for (int i = 0; i < maxRoots; i++) {
334 final FieldGradient<T> hansenKernel = hansenInit[i].getValueGradient(e2, chi, chi2);
335 this.hansenRoot[0][i] = hansenKernel.getValue();
336 this.hansenDerivRoot[0][i] = hansenKernel.getPartialDerivative(0);
337 }
338
339 for (int i = 1; i < numSlices; i++) {
340 for (int k = 0; k < 4; k++) {
341 final PolynomialFunction[] mv = mpvec[N0 - (i * SLICE) - k + 3 + offset];
342 final PolynomialFunction[] sv = mpvecDeriv[N0 - (i * SLICE) - k + 3 + offset];
343
344 hansenDerivRoot[i][k] = mv[3].value(chi).multiply(hansenDerivRoot[i - 1][3]).
345 add(mv[2].value(chi).multiply(hansenDerivRoot[i - 1][2])).
346 add(mv[1].value(chi).multiply(hansenDerivRoot[i - 1][1])).
347 add(mv[0].value(chi).multiply(hansenDerivRoot[i - 1][0])).
348 add(sv[3].value(chi).multiply(hansenRoot[i - 1][3])).
349 add(sv[2].value(chi).multiply(hansenRoot[i - 1][2])).
350 add(sv[1].value(chi).multiply(hansenRoot[i - 1][1])).
351 add(sv[0].value(chi).multiply(hansenRoot[i - 1][0]));
352
353 hansenRoot[i][k] = mv[3].value(chi).multiply(hansenRoot[i - 1][3]).
354 add(mv[2].value(chi).multiply(hansenRoot[i - 1][2])).
355 add(mv[1].value(chi).multiply(hansenRoot[i - 1][1])).
356 add(mv[0].value(chi).multiply(hansenRoot[i - 1][0]));
357 }
358 }
359 }
360
361 /**
362 * Compute the value of the Hansen coefficient K<sub>j</sub><sup>-n-1, s</sup>.
363 *
364 * @param mnm1 -n-1
365 * @param chi χ
366 * @return the coefficient K<sub>j</sub><sup>-n-1, s</sup>
367 */
368 public T getValue(final int mnm1, final T chi) {
369 //Compute n
370 final int n = -mnm1 - 1;
371
372 //Compute the potential slice
373 int sliceNo = (n + N0 + 4) / SLICE;
374 if (sliceNo < numSlices) {
375 //Compute the index within the slice
376 final int indexInSlice = (n + N0 + 4) % SLICE;
377
378 //Check if a root must be returned
379 if (indexInSlice <= 3) {
380 return hansenRoot[sliceNo][indexInSlice];
381 }
382 } else {
383 //the value was a potential root for a slice, but that slice was not required
384 //Decrease the slice number
385 sliceNo--;
386 }
387
388 // Computes the coefficient by linear transformation
389 // Danielson 2.7.3-(9) or Collins 4-236 and Petre's paper
390 final PolynomialFunction[] v = mpvec[mnm1 + offset];
391 return v[3].value(chi).multiply(hansenRoot[sliceNo][3]).
392 add(v[2].value(chi).multiply(hansenRoot[sliceNo][2])).
393 add(v[1].value(chi).multiply(hansenRoot[sliceNo][1])).
394 add(v[0].value(chi).multiply(hansenRoot[sliceNo][0]));
395
396 }
397
398 /**
399 * Compute the value of the derivative dK<sub>j</sub><sup>-n-1, s</sup> / de².
400 *
401 * @param mnm1 -n-1
402 * @param chi χ
403 * @return the derivative dK<sub>j</sub><sup>-n-1, s</sup> / de²
404 */
405 public T getDerivative(final int mnm1, final T chi) {
406
407 //Compute n
408 final int n = -mnm1 - 1;
409
410 //Compute the potential slice
411 int sliceNo = (n + N0 + 4) / SLICE;
412 if (sliceNo < numSlices) {
413 //Compute the index within the slice
414 final int indexInSlice = (n + N0 + 4) % SLICE;
415
416 //Check if a root must be returned
417 if (indexInSlice <= 3) {
418 return hansenDerivRoot[sliceNo][indexInSlice];
419 }
420 } else {
421 //the value was a potential root for a slice, but that slice was not required
422 //Decrease the slice number
423 sliceNo--;
424 }
425
426 // Computes the coefficient by linear transformation
427 // Danielson 2.7.3-(9) or Collins 4-236 and Petre's paper
428 final PolynomialFunction[] v = mpvec[mnm1 + this.offset];
429 final PolynomialFunction[] vv = mpvecDeriv[mnm1 + this.offset];
430
431 return v[3].value(chi).multiply(hansenDerivRoot[sliceNo][3]).
432 add(v[2].value(chi).multiply(hansenDerivRoot[sliceNo][2])).
433 add(v[1].value(chi).multiply(hansenDerivRoot[sliceNo][1])).
434 add(v[0].value(chi).multiply(hansenDerivRoot[sliceNo][0])).
435 add(vv[3].value(chi).multiply(hansenRoot[sliceNo][3])).
436 add(vv[2].value(chi).multiply(hansenRoot[sliceNo][2])).
437 add( vv[1].value(chi).multiply(hansenRoot[sliceNo][1])).
438 add(vv[0].value(chi).multiply(hansenRoot[sliceNo][0]));
439
440 }
441
442 /**
443 * Compute a Hansen coefficient with series.
444 * <p>
445 * This class implements the computation of the Hansen kernels
446 * through a power series in e² and that is using
447 * modified Newcomb operators. The reference formulae can be found
448 * in Danielson 2.7.3-10 and 3.3-5
449 * </p>
450 */
451 private static class FieldHansenCoefficientsBySeries <T extends CalculusFieldElement<T>> {
452
453 /** -n-1 coefficient. */
454 private final int mnm1;
455
456 /** s coefficient. */
457 private final int s;
458
459 /** j coefficient. */
460 private final int j;
461
462 /** Max power in e² for the Newcomb's series expansion. */
463 private final int maxNewcomb;
464
465 /** Polynomial representing the serie. */
466 private PolynomialFunction polynomial;
467
468 /**
469 * Class constructor.
470 *
471 * @param mnm1 -n-1 value
472 * @param s s value
473 * @param j j value
474 * @param maxHansen max power of e² in series expansion
475 * @param field field for the function parameters and value
476 */
477 FieldHansenCoefficientsBySeries(final int mnm1, final int s,
478 final int j, final int maxHansen, final Field<T> field) {
479 this.mnm1 = mnm1;
480 this.s = s;
481 this.j = j;
482 this.maxNewcomb = maxHansen;
483 this.polynomial = generatePolynomial();
484 }
485
486 /** Computes the value of Hansen kernel and its derivative at e².
487 * <p>
488 * The formulae applied are described in Danielson 2.7.3-10 and
489 * 3.3-5
490 * </p>
491 * @param e2 e²
492 * @param chi Χ
493 * @param chi2 Χ²
494 * @return the value of the Hansen coefficient and its derivative for e²
495 */
496 private FieldGradient<T> getValueGradient(final T e2, final T chi, final T chi2) {
497
498 final T zero = e2.getField().getZero();
499 //Estimation of the serie expansion at e2
500 final FieldGradient<T> serie = polynomial.value(FieldGradient.variable(1, 0, e2));
501
502 final T value = FastMath.pow(chi2, -mnm1 - 1).multiply(serie.getValue()).divide(chi);
503 final T coef = zero.subtract(mnm1 + 1.5);
504 final T derivative = coef.multiply(chi2).multiply(value).
505 add(FastMath.pow(chi2, -mnm1 - 1).multiply(serie.getPartialDerivative(0)).divide(chi));
506 return new FieldGradient<T>(value, derivative);
507 }
508
509 /** Generate the serie expansion in e².
510 * <p>
511 * Generate the series expansion in e² used in the formulation
512 * of the Hansen kernel (see Danielson 2.7.3-10):
513 * Σ Y<sup>ns</sup><sub>α+a,α+b</sub>
514 * *e<sup>2α</sup>
515 * </p>
516 * @return polynomial representing the power serie expansion
517 */
518 private PolynomialFunction generatePolynomial() {
519 // Initialization
520 final int aHT = FastMath.max(j - s, 0);
521 final int bHT = FastMath.max(s - j, 0);
522
523 final double[] coefficients = new double[maxNewcomb + 1];
524
525 //Loop for getting the Newcomb operators
526 for (int alphaHT = 0; alphaHT <= maxNewcomb; alphaHT++) {
527 coefficients[alphaHT] =
528 NewcombOperators.getValue(alphaHT + aHT, alphaHT + bHT, mnm1, s);
529 }
530
531 //Creation of the polynomial
532 return new PolynomialFunction(coefficients);
533 }
534 }
535
536 }