1 /* Copyright 2002-2024 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;
18
19 import org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.Field;
21 import org.hipparchus.util.CombinatoricsUtils;
22 import org.hipparchus.util.FastMath;
23 import org.hipparchus.util.MathArrays;
24 import org.orekit.errors.OrekitException;
25 import org.orekit.errors.OrekitMessages;
26
27 import java.util.SortedMap;
28 import java.util.concurrent.ConcurrentSkipListMap;
29
30 /**
31 * This class is designed to provide coefficient from the DSST theory.
32 *
33 * @author Romain Di Costanzo
34 */
35 public class CoefficientsFactory {
36
37 /** Internal storage of the polynomial values. Reused for further computation. */
38 private static SortedMap<NSKey, Double> VNS = new ConcurrentSkipListMap<>();
39
40 /** Last computed order for V<sub>ns</sub> coefficients. */
41 private static int LAST_VNS_ORDER = 2;
42
43 /** Static initialization for the V<sub>ns</sub> coefficient. */
44 static {
45 // Initialization
46 VNS.put(new NSKey(0, 0), 1.);
47 VNS.put(new NSKey(1, 0), 0.);
48 VNS.put(new NSKey(1, 1), 0.5);
49 }
50
51 /** Private constructor as the class is a utility class.
52 */
53 private CoefficientsFactory() {
54 }
55
56 /** Compute the Q<sub>n,s</sub> coefficients evaluated at γ from the recurrence formula 2.8.3-(2).
57 * <p>
58 * Q<sub>n,s</sub> coefficients are computed for n = 0 to nMax
59 * and s = 0 to sMax + 1 in order to also get the derivative dQ<sub>n,s</sub>/dγ = Q(n, s + 1)
60 * </p>
61 * @param gamma γ angle
62 * @param nMax n max value
63 * @param sMax s max value
64 * @return Q<sub>n,s</sub> coefficients array
65 */
66 public static double[][] computeQns(final double gamma, final int nMax, final int sMax) {
67
68 // Initialization
69 final int sDim = FastMath.min(sMax + 1, nMax) + 1;
70 final int rows = nMax + 1;
71 final double[][] Qns = new double[rows][];
72 for (int i = 0; i <= nMax; i++) {
73 final int snDim = FastMath.min(i + 1, sDim);
74 Qns[i] = new double[snDim];
75 }
76
77 // first element
78 Qns[0][0] = 1;
79
80 for (int n = 1; n <= nMax; n++) {
81 final int snDim = FastMath.min(n + 1, sDim);
82 for (int s = 0; s < snDim; s++) {
83 if (n == s) {
84 Qns[n][s] = (2. * s - 1.) * Qns[s - 1][s - 1];
85 } else if (n == (s + 1)) {
86 Qns[n][s] = (2. * s + 1.) * gamma * Qns[s][s];
87 } else {
88 Qns[n][s] = (2. * n - 1.) * gamma * Qns[n - 1][s] - (n + s - 1.) * Qns[n - 2][s];
89 Qns[n][s] /= n - s;
90 }
91 }
92 }
93
94 return Qns;
95 }
96
97 /** Compute the Q<sub>n,s</sub> coefficients evaluated at γ from the recurrence formula 2.8.3-(2).
98 * <p>
99 * Q<sub>n,s</sub> coefficients are computed for n = 0 to nMax
100 * and s = 0 to sMax + 1 in order to also get the derivative dQ<sub>n,s</sub>/dγ = Q(n, s + 1)
101 * </p>
102 * @param gamma γ angle
103 * @param nMax n max value
104 * @param sMax s max value
105 * @param <T> the type of the field elements
106 * @return Q<sub>n,s</sub> coefficients array
107 */
108 public static <T extends CalculusFieldElement<T>> T[][] computeQns(final T gamma, final int nMax, final int sMax) {
109
110 // Initialization
111 final int sDim = FastMath.min(sMax + 1, nMax) + 1;
112 final int rows = nMax + 1;
113 final T[][] Qns = MathArrays.buildArray(gamma.getField(), rows, FastMath.min(nMax + 1, sDim) - 1);
114 for (int i = 0; i <= nMax; i++) {
115 final int snDim = FastMath.min(i + 1, sDim);
116 Qns[i] = MathArrays.buildArray(gamma.getField(), snDim);
117 }
118
119 // first element
120 Qns[0][0] = gamma.subtract(gamma).add(1.);
121
122 for (int n = 1; n <= nMax; n++) {
123 final int snDim = FastMath.min(n + 1, sDim);
124 for (int s = 0; s < snDim; s++) {
125 if (n == s) {
126 Qns[n][s] = Qns[s - 1][s - 1].multiply(2. * s - 1.);
127 } else if (n == (s + 1)) {
128 Qns[n][s] = Qns[s][s].multiply(gamma).multiply(2. * s + 1.);
129 } else {
130 Qns[n][s] = Qns[n - 1][s].multiply(gamma).multiply(2. * n - 1.).subtract(Qns[n - 2][s].multiply(n + s - 1.));
131 Qns[n][s] = Qns[n][s].divide(n - s);
132 }
133 }
134 }
135
136 return Qns;
137 }
138
139 /** Compute recursively G<sub>s</sub> and H<sub>s</sub> polynomials from equation 3.1-(5).
140 * @param k x-component of the eccentricity vector
141 * @param h y-component of the eccentricity vector
142 * @param alpha 1st direction cosine
143 * @param beta 2nd direction cosine
144 * @param order development order
145 * @return Array of G<sub>s</sub> and H<sub>s</sub> polynomials for s from 0 to order.<br>
146 * The 1st column contains the G<sub>s</sub> values.
147 * The 2nd column contains the H<sub>s</sub> values.
148 */
149 public static double[][] computeGsHs(final double k, final double h,
150 final double alpha, final double beta,
151 final int order) {
152 // Constant terms
153 final double hamkb = h * alpha - k * beta;
154 final double kaphb = k * alpha + h * beta;
155 // Initialization
156 final double[][] GsHs = new double[2][order + 1];
157 GsHs[0][0] = 1.;
158 GsHs[1][0] = 0.;
159
160 for (int s = 1; s <= order; s++) {
161 // Gs coefficient
162 GsHs[0][s] = kaphb * GsHs[0][s - 1] - hamkb * GsHs[1][s - 1];
163 // Hs coefficient
164 GsHs[1][s] = hamkb * GsHs[0][s - 1] + kaphb * GsHs[1][s - 1];
165 }
166
167 return GsHs;
168 }
169
170 /** Compute recursively G<sub>s</sub> and H<sub>s</sub> polynomials from equation 3.1-(5).
171 * @param k x-component of the eccentricity vector
172 * @param h y-component of the eccentricity vector
173 * @param alpha 1st direction cosine
174 * @param beta 2nd direction cosine
175 * @param order development order
176 * @param field field of elements
177 * @param <T> the type of the field elements
178 * @return Array of G<sub>s</sub> and H<sub>s</sub> polynomials for s from 0 to order.<br>
179 * The 1st column contains the G<sub>s</sub> values.
180 * The 2nd column contains the H<sub>s</sub> values.
181 */
182 public static <T extends CalculusFieldElement<T>> T[][] computeGsHs(final T k, final T h,
183 final T alpha, final T beta,
184 final int order, final Field<T> field) {
185 // Zero for initialization
186 final T zero = field.getZero();
187
188 // Constant terms
189 final T hamkb = h.multiply(alpha).subtract(k.multiply(beta));
190 final T kaphb = k.multiply(alpha).add(h.multiply(beta));
191 // Initialization
192 final T[][] GsHs = MathArrays.buildArray(field, 2, order + 1);
193 GsHs[0][0] = zero.newInstance(1.);
194 GsHs[1][0] = zero;
195
196 for (int s = 1; s <= order; s++) {
197 // Gs coefficient
198 GsHs[0][s] = kaphb.multiply(GsHs[0][s - 1]).subtract(hamkb.multiply(GsHs[1][s - 1]));
199 // Hs coefficient
200 GsHs[1][s] = hamkb.multiply(GsHs[0][s - 1]).add(kaphb.multiply(GsHs[1][s - 1]));
201 }
202
203 return GsHs;
204 }
205
206 /** Compute the V<sub>n,s</sub> coefficients from 2.8.2-(1)(2).
207 * @param order Order of the computation. Computation will be done from 0 to order -1
208 * @return Map of the V<sub>n, s</sub> coefficients
209 * @since 11.3.3
210 */
211 public static SortedMap<NSKey, Double> computeVns(final int order) {
212
213 if (order > LAST_VNS_ORDER) {
214 // Compute coefficient
215 // Need previous computation as recurrence relation is done at s + 1 and n + 2
216 final int min = FastMath.max(LAST_VNS_ORDER - 2, 0);
217 for (int n = min; n < order; n++) {
218 for (int s = 0; s < n + 1; s++) {
219 if ((n - s) % 2 != 0) {
220 VNS.put(new NSKey(n, s), 0.);
221 } else {
222 // s = n
223 if (n == s && (s + 1) < order) {
224 VNS.put(new NSKey(s + 1, s + 1), VNS.get(new NSKey(s, s)) / (2 * s + 2.));
225 }
226 // otherwise
227 if ((n + 2) < order) {
228 VNS.put(new NSKey(n + 2, s), VNS.get(new NSKey(n, s)) * (-n + s - 1.) / (n + s + 2.));
229 }
230 }
231 }
232 }
233 LAST_VNS_ORDER = order;
234 }
235 return new ConcurrentSkipListMap<>(VNS);
236 }
237
238 /** Get the V<sub>n,s</sub><sup>m</sup> coefficient from V<sub>n,s</sub>.
239 * <br>See § 2.8.2 in Danielson paper.
240 * @param m m
241 * @param n n
242 * @param s s
243 * @return The V<sub>n, s</sub> <sup>m</sup> coefficient
244 */
245 public static double getVmns(final int m, final int n, final int s) {
246 if (m > n) {
247 throw new OrekitException(OrekitMessages.DSST_VMNS_COEFFICIENT_ERROR_MS, m, n);
248 }
249 final double fns = CombinatoricsUtils.factorialDouble(n + FastMath.abs(s));
250 final double fnm = CombinatoricsUtils.factorialDouble(n - m);
251
252 double result = 0;
253 // If (n - s) is odd, the Vmsn coefficient is null
254 if ((n - s) % 2 == 0) {
255 // Update the Vns coefficient
256 if ((n + 1) > LAST_VNS_ORDER) {
257 computeVns(n + 1);
258 }
259 if (s >= 0) {
260 result = fns * VNS.get(new NSKey(n, s)) / fnm;
261 } else {
262 // If s < 0 : Vmn-s = (-1)^(-s) Vmns
263 final int mops = (s % 2 == 0) ? 1 : -1;
264 result = mops * fns * VNS.get(new NSKey(n, -s)) / fnm;
265 }
266 }
267 return result;
268 }
269
270 /** Key formed by two integer values. */
271 public static class NSKey implements Comparable<NSKey> {
272
273 /** n value. */
274 private final int n;
275
276 /** s value. */
277 private final int s;
278
279 /** Simple constructor.
280 * @param n n
281 * @param s s
282 */
283 public NSKey(final int n, final int s) {
284 this.n = n;
285 this.s = s;
286 }
287
288 /** Get n.
289 * @return n
290 */
291 public int getN() {
292 return n;
293 }
294
295 /** Get s.
296 * @return s
297 */
298 public int getS() {
299 return s;
300 }
301
302 /** {@inheritDoc} */
303 public int compareTo(final NSKey key) {
304 int result = 1;
305 if (n == key.n) {
306 if (s < key.s) {
307 result = -1;
308 } else if (s == key.s) {
309 result = 0;
310 }
311 } else if (n < key.n) {
312 result = -1;
313 }
314 return result;
315 }
316
317 /** {@inheritDoc} */
318 public boolean equals(final Object key) {
319
320 if (key == this) {
321 // first fast check
322 return true;
323 }
324
325 if (key instanceof NSKey) {
326 return n == ((NSKey) key).n && s == ((NSKey) key).s;
327 }
328
329 return false;
330
331 }
332
333 /** {@inheritDoc} */
334 public int hashCode() {
335 return 0x998493a6 ^ (n << 8) ^ s;
336 }
337
338 }
339
340 }