1 /* Copyright 2002-2025 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 final double[][] hansenRoot;
51
52 /** The derivatives of the Hansen coefficients used as roots. */
53 private final double[][] hansenDerivRoot;
54
55 /** The minimum value for the order. */
56 private final int Nmin;
57
58 /** The index of the initial condition, Petre's paper. */
59 private final int N0;
60
61 /** The number of slices needed to compute the coefficients. */
62 private final int numSlices;
63
64 /**
65 * The offset used to identify the polynomial that corresponds to a negative.
66 * value of n in the internal array that starts at 0
67 */
68 private final int offset;
69
70 /** The objects used to calculate initial data by means of Newcomb operators. */
71 private final HansenCoefficientsBySeries[] hansenInit;
72
73 /**
74 * Constructor.
75 *
76 * @param nMax the maximum (absolute) value of n parameter
77 * @param s s parameter
78 * @param j j parameter
79 * @param n0 the minimum (absolute) value of n
80 * @param maxHansen maximum power of e2 in Hansen expansion
81 */
82 public HansenTesseralLinear(final int nMax, final int s, final int j, final int n0, final int maxHansen) {
83 //Initialize the fields
84 this.offset = nMax + 1;
85 this.Nmin = -nMax - 1;
86 this.N0 = -n0 - 4;
87
88 //Ensure that only the needed terms are computed
89 final int maxRoots = FastMath.min(4, N0 - Nmin + 4);
90 this.hansenInit = new HansenCoefficientsBySeries[maxRoots];
91 for (int i = 0; i < maxRoots; i++) {
92 this.hansenInit[i] = new HansenCoefficientsBySeries(N0 - i + 3, s, j, maxHansen);
93 }
94
95 // The first 4 values are computed with series. No linear combination is needed.
96 final int size = N0 - Nmin;
97 this.numSlices = (int) FastMath.max(FastMath.ceil(((double) size) / SLICE), 1);
98 hansenRoot = new double[numSlices][4];
99 hansenDerivRoot = new double[numSlices][4];
100 if (size > 0) {
101 mpvec = new PolynomialFunction[size][];
102 mpvecDeriv = new PolynomialFunction[size][];
103
104 // Prepare the database of the associated polynomials
105 HansenUtilities.generateTesseralPolynomials(N0, Nmin, offset, SLICE, j, s,
106 mpvec, mpvecDeriv);
107 }
108
109 }
110
111 /**
112 * Compute the values for the first four coefficients and their derivatives by means of series.
113 *
114 * @param e2 e²
115 * @param chi Χ
116 * @param chi2 Χ²
117 */
118 public void computeInitValues(final double e2, final double chi, final double chi2) {
119 // compute the values for n, n+1, n+2 and n+3 by series
120 // See Danielson 2.7.3-(10)
121 //Ensure that only the needed terms are computed
122 final int maxRoots = FastMath.min(4, N0 - Nmin + 4);
123 for (int i = 0; i < maxRoots; i++) {
124 final Gradient hansenKernel = hansenInit[i].getValueGradient(e2, chi, chi2);
125 this.hansenRoot[0][i] = hansenKernel.getValue();
126 this.hansenDerivRoot[0][i] = hansenKernel.getPartialDerivative(0);
127 }
128
129 for (int i = 1; i < numSlices; i++) {
130 for (int k = 0; k < 4; k++) {
131 final PolynomialFunction[] mv = mpvec[N0 - (i * SLICE) - k + 3 + offset];
132 final PolynomialFunction[] sv = mpvecDeriv[N0 - (i * SLICE) - k + 3 + offset];
133
134 hansenDerivRoot[i][k] = mv[3].value(chi) * hansenDerivRoot[i - 1][3] +
135 mv[2].value(chi) * hansenDerivRoot[i - 1][2] +
136 mv[1].value(chi) * hansenDerivRoot[i - 1][1] +
137 mv[0].value(chi) * hansenDerivRoot[i - 1][0] +
138 sv[3].value(chi) * hansenRoot[i - 1][3] +
139 sv[2].value(chi) * hansenRoot[i - 1][2] +
140 sv[1].value(chi) * hansenRoot[i - 1][1] +
141 sv[0].value(chi) * hansenRoot[i - 1][0];
142
143 hansenRoot[i][k] = mv[3].value(chi) * hansenRoot[i - 1][3] +
144 mv[2].value(chi) * hansenRoot[i - 1][2] +
145 mv[1].value(chi) * hansenRoot[i - 1][1] +
146 mv[0].value(chi) * hansenRoot[i - 1][0];
147 }
148 }
149 }
150
151 /**
152 * Compute the value of the Hansen coefficient K<sub>j</sub><sup>-n-1, s</sup>.
153 *
154 * @param mnm1 -n-1
155 * @param chi χ
156 * @return the coefficient K<sub>j</sub><sup>-n-1, s</sup>
157 */
158 public double getValue(final int mnm1, final double chi) {
159 //Compute n
160 final int n = -mnm1 - 1;
161
162 //Compute the potential slice
163 int sliceNo = (n + N0 + 4) / SLICE;
164 if (sliceNo < numSlices) {
165 //Compute the index within the slice
166 final int indexInSlice = (n + N0 + 4) % SLICE;
167
168 //Check if a root must be returned
169 if (indexInSlice <= 3) {
170 return hansenRoot[sliceNo][indexInSlice];
171 }
172 } else {
173 //the value was a potential root for a slice, but that slice was not required
174 //Decrease the slice number
175 sliceNo--;
176 }
177
178 // Computes the coefficient by linear transformation
179 // Danielson 2.7.3-(9) or Collins 4-236 and Petre's paper
180 final PolynomialFunction[] v = mpvec[mnm1 + offset];
181 return v[3].value(chi) * hansenRoot[sliceNo][3] +
182 v[2].value(chi) * hansenRoot[sliceNo][2] +
183 v[1].value(chi) * hansenRoot[sliceNo][1] +
184 v[0].value(chi) * hansenRoot[sliceNo][0];
185
186 }
187
188 /**
189 * Compute the value of the derivative dK<sub>j</sub><sup>-n-1, s</sup> / de².
190 *
191 * @param mnm1 -n-1
192 * @param chi χ
193 * @return the derivative dK<sub>j</sub><sup>-n-1, s</sup> / de²
194 */
195 public double getDerivative(final int mnm1, final double chi) {
196
197 //Compute n
198 final int n = -mnm1 - 1;
199
200 //Compute the potential slice
201 int sliceNo = (n + N0 + 4) / SLICE;
202 if (sliceNo < numSlices) {
203 //Compute the index within the slice
204 final int indexInSlice = (n + N0 + 4) % SLICE;
205
206 //Check if a root must be returned
207 if (indexInSlice <= 3) {
208 return hansenDerivRoot[sliceNo][indexInSlice];
209 }
210 } else {
211 //the value was a potential root for a slice, but that slice was not required
212 //Decrease the slice number
213 sliceNo--;
214 }
215
216 // Computes the coefficient by linear transformation
217 // Danielson 2.7.3-(9) or Collins 4-236 and Petre's paper
218 final PolynomialFunction[] v = mpvec[mnm1 + this.offset];
219 final PolynomialFunction[] vv = mpvecDeriv[mnm1 + this.offset];
220
221 return v[3].value(chi) * hansenDerivRoot[sliceNo][3] +
222 v[2].value(chi) * hansenDerivRoot[sliceNo][2] +
223 v[1].value(chi) * hansenDerivRoot[sliceNo][1] +
224 v[0].value(chi) * hansenDerivRoot[sliceNo][0] +
225 vv[3].value(chi) * hansenRoot[sliceNo][3] +
226 vv[2].value(chi) * hansenRoot[sliceNo][2] +
227 vv[1].value(chi) * hansenRoot[sliceNo][1] +
228 vv[0].value(chi) * hansenRoot[sliceNo][0];
229
230 }
231
232 /**
233 * Compute a Hansen coefficient with series.
234 * <p>
235 * This class implements the computation of the Hansen kernels
236 * through a power series in e² and that is using
237 * modified Newcomb operators. The reference formulae can be found
238 * in Danielson 2.7.3-10 and 3.3-5
239 * </p>
240 */
241 private static class HansenCoefficientsBySeries {
242
243 /** -n-1 coefficient. */
244 private final int mnm1;
245
246 /** s coefficient. */
247 private final int s;
248
249 /** j coefficient. */
250 private final int j;
251
252 /** Max power in e² for the Newcomb's series expansion. */
253 private final int maxNewcomb;
254
255 /** Polynomial representing the serie. */
256 private final PolynomialFunction polynomial;
257
258 /**
259 * Class constructor.
260 *
261 * @param mnm1 -n-1 value
262 * @param s s value
263 * @param j j value
264 * @param maxHansen max power of e² in series expansion
265 */
266 HansenCoefficientsBySeries(final int mnm1, final int s,
267 final int j, final int maxHansen) {
268 this.mnm1 = mnm1;
269 this.s = s;
270 this.j = j;
271 this.maxNewcomb = maxHansen;
272 this.polynomial = generatePolynomial();
273 }
274
275 /** Computes the value of Hansen kernel and its derivative at e².
276 * <p>
277 * The formulae applied are described in Danielson 2.7.3-10 and
278 * 3.3-5
279 * </p>
280 * @param e2 e²
281 * @param chi Χ
282 * @param chi2 Χ²
283 * @return the value of the Hansen coefficient and its derivative for e²
284 */
285 public Gradient getValueGradient(final double e2, final double chi, final double chi2) {
286
287 //Estimation of the serie expansion at e2
288 final Gradient serie = polynomial.value(Gradient.variable(1, 0, e2));
289
290 final double value = FastMath.pow(chi2, -mnm1 - 1) * serie.getValue() / chi;
291 final double coef = -(mnm1 + 1.5);
292 final double derivative = coef * chi2 * value +
293 FastMath.pow(chi2, -mnm1 - 1) * serie.getPartialDerivative(0) / chi;
294 return new Gradient(value, derivative);
295 }
296
297 /** Generate the serie expansion in e².
298 * <p>
299 * Generate the series expansion in e² used in the formulation
300 * of the Hansen kernel (see Danielson 2.7.3-10):
301 * Σ Y<sup>ns</sup><sub>α+a,α+b</sub>
302 * *e<sup>2α</sup>
303 * </p>
304 * @return polynomial representing the power serie expansion
305 */
306 private PolynomialFunction generatePolynomial() {
307 // Initialization
308 final int aHT = FastMath.max(j - s, 0);
309 final int bHT = FastMath.max(s - j, 0);
310
311 final double[] coefficients = new double[maxNewcomb + 1];
312
313 //Loop for getting the Newcomb operators
314 for (int alphaHT = 0; alphaHT <= maxNewcomb; alphaHT++) {
315 coefficients[alphaHT] =
316 NewcombOperators.getValue(alphaHT + aHT, alphaHT + bHT, mnm1, s);
317 }
318
319 //Creation of the polynomial
320 return new PolynomialFunction(coefficients);
321 }
322 }
323
324 }