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