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.polynomials.PolynomialFunction;
20 import org.hipparchus.util.FastMath;
21
22 /**
23 * Hansen coefficients K(t,n,s) for t=0 and n > 0.
24 * <p>
25 * Implements Collins 4-254 or Danielson 2.7.3-(7) for Hansen Coefficients and
26 * Danielson 3.2-(3) for derivatives. The recursions are transformed into
27 * composition of linear transformations to obtain the associated polynomials
28 * for coefficients and their derivatives - see Petre's paper
29 *
30 * @author Petre Bazavan
31 * @author Lucian Barbulescu
32 */
33 public class HansenThirdBodyLinear {
34
35 /** The number of coefficients that will be computed with a set of roots. */
36 private static final int SLICE = 10;
37
38 /**
39 * The first vector of polynomials associated to Hansen coefficients and
40 * derivatives.
41 */
42 private PolynomialFunction[][] mpvec;
43
44 /** The second vector of polynomials associated only to derivatives. */
45 private PolynomialFunction[][] mpvecDeriv;
46
47 /** The Hansen coefficients used as roots. */
48 private double[][] hansenRoot;
49
50 /** The derivatives of the Hansen coefficients used as roots. */
51 private double[][] hansenDerivRoot;
52
53 /** The number of slices needed to compute the coefficients. */
54 private int numSlices;
55
56 /** The maximum order of n indexes. */
57 private int nMax;
58
59 /** The index of the initial condition, Petre's paper. */
60 private int N0;
61
62 /** The s index. */
63 private int s;
64
65 /** (-1)<sup>s</sup> * (2*s + 1)!! / (s+1)! */
66 private double twosp1dfosp1f;
67
68 /** (-1)<sup>s</sup> * (2*s + 1)!! / (s+2)! */
69 private double twosp1dfosp2f;
70
71 /** (-1)<sup>s</sup> * 2 * (2*s + 1)!! / (s+2)! */
72 private double two2sp1dfosp2f;
73
74 /** (2*s + 3). */
75 private double twosp3;
76
77 /**
78 * Constructor.
79 *
80 * @param nMax the maximum value of n
81 * @param s the value of s
82 */
83 public HansenThirdBodyLinear(final int nMax, final int s) {
84
85 // initialise fields
86 this.nMax = nMax;
87 N0 = s;
88 this.s = s;
89
90 // initialization of structures for stored data
91 mpvec = new PolynomialFunction[this.nMax + 1][];
92 mpvecDeriv = new PolynomialFunction[this.nMax + 1][];
93
94 //Compute the fields that will be used to determine the initial values for the coefficients
95 this.twosp1dfosp1f = (s % 2 == 0) ? 1.0 : -1.0;
96 for (int i = s; i >= 1; i--) {
97 this.twosp1dfosp1f *= (2.0 * i + 1.0) / (i + 1.0);
98 }
99
100 this.twosp1dfosp2f = this.twosp1dfosp1f / (s + 2.);
101 this.twosp3 = 2 * s + 3;
102 this.two2sp1dfosp2f = 2 * this.twosp1dfosp2f;
103
104 // initialization of structures for stored data
105 mpvec = new PolynomialFunction[this.nMax + 1][];
106 mpvecDeriv = new PolynomialFunction[this.nMax + 1][];
107
108 this.numSlices = FastMath.max(1, (nMax - s + SLICE - 2) / SLICE);
109 hansenRoot = new double[numSlices][2];
110 hansenDerivRoot = new double[numSlices][2];
111
112 // Prepare the database of the associated polynomials
113 generatePolynomials();
114
115 }
116
117 /**
118 * Compute polynomial coefficient a.
119 *
120 * <p>
121 * It is used to generate the coefficient for K₀<sup>n-1, s</sup> when computing K₀<sup>n, s</sup>
122 * and the coefficient for dK₀<sup>n-1, s</sup> / dΧ when computing dK₀<sup>n, s</sup> / dΧ
123 * </p>
124 *
125 * <p>
126 * See Danielson 2.7.3-(7c) and Collins 4-254 and 4-257
127 * </p>
128 *
129 * @param n n value
130 * @return the polynomial
131 */
132 private PolynomialFunction a(final int n) {
133 // from recurrence Danielson 2.7.3-(7c), Collins 4-254
134 final double r1 = 2 * n + 1;
135 final double r2 = n + 1;
136
137 return new PolynomialFunction(new double[] {
138 r1 / r2
139 });
140 }
141
142 /**
143 * Compute polynomial coefficient b.
144 *
145 * <p>
146 * It is used to generate the coefficient for K₀<sup>n-2, s</sup> when computing K₀<sup>n, s</sup>
147 * and the coefficient for dK₀<sup>n-2, s</sup> / dΧ when computing dK₀<sup>n, s</sup> / dΧ
148 * </p>
149 *
150 * <p>
151 * See Danielson 2.7.3-(7c) and Collins 4-254 and 4-257
152 * </p>
153 *
154 * @param n n value
155 * @return the polynomial
156 */
157 private PolynomialFunction b(final int n) {
158 // from recurrence Danielson 2.7.3-(7c), Collins 4-254
159 final double r1 = (n + s) * (n - s);
160 final double r2 = n * (n + 1);
161
162 return new PolynomialFunction(new double[] {
163 0.0, 0.0, -r1 / r2
164 });
165 }
166
167 /**
168 * Compute polynomial coefficient d.
169 *
170 * <p>
171 * It is used to generate the coefficient for K₀<sup>n-2, s</sup> when computing dK₀<sup>n, s</sup> / dΧ
172 * </p>
173 *
174 * <p>
175 * See Danielson 2.7.3-(7c) and Collins 4-254 and 4-257
176 * </p>
177 *
178 * @param n n value
179 * @return the polynomial
180 */
181 private PolynomialFunction d(final int n) {
182 // from Danielson 3.2-(3b)
183 final double r1 = 2 * (n + s) * (n - s);
184 final double r2 = n * (n + 1);
185
186 return new PolynomialFunction(new double[] {
187 0.0, 0.0, 0.0, r1 / r2
188 });
189 }
190
191 /**
192 * Generate the polynomials needed in the linear transformation.
193 *
194 * <p>
195 * See Petre's paper
196 * </p>
197 */
198 private void generatePolynomials() {
199
200 int sliceCounter = 0;
201
202 // Initialization of the matrices for linear transformations
203 // The final configuration of these matrices are obtained by composition
204 // of linear transformations
205
206 // the matrix A for the polynomials associated
207 // to Hansen coefficients, Petre's pater
208 PolynomialFunctionMatrix A = HansenUtilities.buildIdentityMatrix2();
209
210 // the matrix D for the polynomials associated
211 // to derivatives, Petre's paper
212 final PolynomialFunctionMatrix B = HansenUtilities.buildZeroMatrix2();
213 PolynomialFunctionMatrix D = HansenUtilities.buildZeroMatrix2();
214 PolynomialFunctionMatrix E = HansenUtilities.buildIdentityMatrix2();
215
216 // The matrix that contains the coefficients at each step
217 final PolynomialFunctionMatrix a = HansenUtilities.buildZeroMatrix2();
218 a.setElem(0, 1, HansenUtilities.ONE);
219
220 // The generation process
221 for (int i = N0 + 2; i <= nMax; i++) {
222 // Collins 4-254 or Danielson 2.7.3-(7)
223 // Petre's paper
224 // The matrix of the current linear transformation is actualized
225 a.setMatrixLine(1, new PolynomialFunction[] {
226 b(i), a(i)
227 });
228
229 // composition of the linear transformations to calculate
230 // the polynomials associated to Hansen coefficients
231 A = A.multiply(a);
232 // store the polynomials associated to Hansen coefficients
233 this.mpvec[i] = A.getMatrixLine(1);
234 // composition of the linear transformations to calculate
235 // the polynomials associated to derivatives
236 // Danielson 3.2-(3b) and Petre's paper
237 D = D.multiply(a);
238 if (sliceCounter % SLICE != 0) {
239 a.setMatrixLine(1, new PolynomialFunction[] {
240 b(i - 1), a(i - 1)
241 });
242 E = E.multiply(a);
243 }
244
245 B.setElem(1, 0, d(i));
246 // F = E.prod(B);
247 D = D.add(E.multiply(B));
248 // store the polynomials associated to the derivatives
249 this.mpvecDeriv[i] = D.getMatrixLine(1);
250
251 if (++sliceCounter % SLICE == 0) {
252 // Re-Initialization of the matrices for linear transformations
253 // The final configuration of these matrices are obtained by composition
254 // of linear transformations
255 A = HansenUtilities.buildIdentityMatrix2();
256 D = HansenUtilities.buildZeroMatrix2();
257 E = HansenUtilities.buildIdentityMatrix2();
258 }
259 }
260 }
261
262 /**
263 * Compute the initial values (see Collins, 4-255, 4-256 and 4.259)
264 * <p>
265 * K₀<sup>s, s</sup> = (-1)<sup>s</sup> * ( (2*s+1)!! / (s+1)! )
266 * </p>
267 * <p>
268 * K₀<sup>s+1, s</sup> = (-1)<sup>s</sup> * ( (2*s+1)!! / (s+2)!
269 * ) * (2*s+3 - χ<sup>-2</sup>)
270 * </p>
271 * <p>
272 * dK₀<sup>s+1, s</sup> / dχ = = (-1)<sup>s</sup> * 2 * (
273 * (2*s+1)!! / (s+2)! ) * χ<sup>-3</sup>
274 * </p>
275 * @param chitm1 sqrt(1 - e²)
276 * @param chitm2 sqrt(1 - e²)²
277 * @param chitm3 sqrt(1 - e²)³
278 */
279 public void computeInitValues(final double chitm1, final double chitm2, final double chitm3) {
280 this.hansenRoot[0][0] = this.twosp1dfosp1f;
281 this.hansenRoot[0][1] = this.twosp1dfosp2f * (this.twosp3 - chitm2);
282 this.hansenDerivRoot[0][0] = 0;
283 this.hansenDerivRoot[0][1] = this.two2sp1dfosp2f * chitm3;
284
285 for (int i = 1; i < numSlices; i++) {
286 for (int j = 0; j < 2; j++) {
287 // Get the required polynomials
288 final PolynomialFunction[] mv = mpvec[s + (i * SLICE) + j];
289 final PolynomialFunction[] sv = mpvecDeriv[s + (i * SLICE) + j];
290
291 //Compute the root derivatives
292 hansenDerivRoot[i][j] = mv[1].value(chitm1) * hansenDerivRoot[i - 1][1] +
293 mv[0].value(chitm1) * hansenDerivRoot[i - 1][0] +
294 sv[1].value(chitm1) * hansenRoot[i - 1][1] +
295 sv[0].value(chitm1) * hansenRoot[i - 1][0];
296
297 //Compute the root Hansen coefficients
298 hansenRoot[i][j] = mv[1].value(chitm1) * hansenRoot[i - 1][1] +
299 mv[0].value(chitm1) * hansenRoot[i - 1][0];
300 }
301 }
302 }
303
304 /**
305 * Compute the value of the Hansen coefficient K₀<sup>n, s</sup>.
306 *
307 * @param n n value
308 * @param chitm1 χ<sup>-1</sup>
309 * @return the coefficient K₀<sup>n, s</sup>
310 */
311 public double getValue(final int n, final double chitm1) {
312 //Compute the potential slice
313 int sliceNo = (n - s) / SLICE;
314 if (sliceNo < numSlices) {
315 //Compute the index within the slice
316 final int indexInSlice = (n - s) % SLICE;
317
318 //Check if a root must be returned
319 if (indexInSlice <= 1) {
320 return hansenRoot[sliceNo][indexInSlice];
321 }
322 } else {
323 //the value was a potential root for a slice, but that slice was not required
324 //Decrease the slice number
325 sliceNo--;
326 }
327
328 // Danielson 2.7.3-(6c)/Collins 4-242 and Petre's paper
329 final PolynomialFunction[] v = mpvec[n];
330 double ret = v[1].value(chitm1) * hansenRoot[sliceNo][1];
331 if (hansenRoot[sliceNo][0] != 0) {
332 ret += v[0].value(chitm1) * hansenRoot[sliceNo][0];
333 }
334
335 return ret;
336
337 }
338
339 /**
340 * Compute the value of the Hansen coefficient dK₀<sup>n, s</sup> / dΧ.
341 *
342 * @param n n value
343 * @param chitm1 χ<sup>-1</sup>
344 * @return the coefficient dK₀<sup>n, s</sup> / dΧ
345 */
346 public double getDerivative(final int n, final double chitm1) {
347 //Compute the potential slice
348 int sliceNo = (n - s) / SLICE;
349 if (sliceNo < numSlices) {
350 //Compute the index within the slice
351 final int indexInSlice = (n - s) % SLICE;
352
353 //Check if a root must be returned
354 if (indexInSlice <= 1) {
355 return hansenDerivRoot[sliceNo][indexInSlice];
356 }
357 } else {
358 //the value was a potential root for a slice, but that slice was not required
359 //Decrease the slice number
360 sliceNo--;
361 }
362
363 final PolynomialFunction[] v = mpvec[n];
364 double ret = v[1].value(chitm1) * hansenDerivRoot[sliceNo][1];
365 if (hansenDerivRoot[sliceNo][0] != 0) {
366 ret += v[0].value(chitm1) * hansenDerivRoot[sliceNo][0];
367 }
368
369 // Danielson 2.7.3-(7c)/Collins 4-254 and Petre's paper
370 final PolynomialFunction[] v1 = mpvecDeriv[n];
371 ret += v1[1].value(chitm1) * hansenRoot[sliceNo][1];
372 if (hansenRoot[sliceNo][0] != 0) {
373 ret += v1[0].value(chitm1) * hansenRoot[sliceNo][0];
374 }
375 return ret;
376
377 }
378
379 }