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