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