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