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 org.hipparchus.CalculusFieldElement;
20 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
21 import org.hipparchus.util.FastMath;
22 import org.hipparchus.util.MathUtils;
23 import org.orekit.frames.Frame;
24 import org.orekit.orbits.FieldOrbit;
25 import org.orekit.time.FieldAbsoluteDate;
26
27 /** Container class for common parameters used by all DSST forces.
28 * <p>
29 * Most of them are defined in Danielson paper at § 2.1.
30 * </p>
31 */
32 public class FieldAuxiliaryElements<T extends CalculusFieldElement<T>> {
33
34 /** Orbit date. */
35 private final FieldAbsoluteDate<T> date;
36
37 /** Orbit frame. */
38 private final Frame frame;
39
40 /** Eccentricity. */
41 private final T ecc;
42
43 /** Keplerian mean motion. */
44 private final T n;
45
46 /** Keplerian period. */
47 private final T period;
48
49 /** Semi-major axis. */
50 private final T sma;
51
52 /** x component of eccentricity vector. */
53 private final T k;
54
55 /** y component of eccentricity vector. */
56 private final T h;
57
58 /** x component of inclination vector. */
59 private final T q;
60
61 /** y component of inclination vector. */
62 private final T p;
63
64 /** Mean longitude. */
65 private final T lm;
66
67 /** True longitude. */
68 private final T lv;
69
70 /** Eccentric longitude. */
71 private final T le;
72
73 /** Retrograde factor I.
74 * <p>
75 * DSST model needs equinoctial orbit as internal representation.
76 * Classical equinoctial elements have discontinuities when inclination
77 * is close to zero. In this representation, I = +1. <br>
78 * To avoid this discontinuity, another representation exists and equinoctial
79 * elements can be expressed in a different way, called "retrograde" orbit.
80 * This implies I = -1. <br>
81 * As Orekit doesn't implement the retrograde orbit, I is always set to +1.
82 * But for the sake of consistency with the theory, the retrograde factor
83 * has been kept in the formulas.
84 * </p>
85 */
86 private final int I;
87
88 /** B = sqrt(1 - h² - k²). */
89 private final T B;
90
91 /** C = 1 + p² + q². */
92 private final T C;
93
94 /** Equinoctial frame f vector. */
95 private final FieldVector3D<T> f;
96
97 /** Equinoctial frame g vector. */
98 private final FieldVector3D<T> g;
99
100 /** Equinoctial frame w vector. */
101 private final FieldVector3D<T> w;
102
103 /** Direction cosine α. */
104 private final T alpha;
105
106 /** Direction cosine β. */
107 private final T beta;
108
109 /** Direction cosine γ. */
110 private final T gamma;
111
112 /** Simple constructor.
113 * @param orbit related mean orbit for auxiliary elements
114 * @param retrogradeFactor retrograde factor I [Eq. 2.1.2-(2)]
115 */
116 public FieldAuxiliaryElements(final FieldOrbit<T> orbit, final int retrogradeFactor) {
117
118 final T pi = orbit.getDate().getField().getZero().getPi();
119
120 // Date of the orbit
121 date = orbit.getDate();
122
123 // Orbit definition frame
124 frame = orbit.getFrame();
125
126 // Eccentricity
127 ecc = orbit.getE();
128
129 // Keplerian mean motion
130 n = orbit.getKeplerianMeanMotion();
131
132 // Keplerian period
133 period = orbit.getKeplerianPeriod();
134
135 // Equinoctial elements [Eq. 2.1.2-(1)]
136 sma = orbit.getA();
137 k = orbit.getEquinoctialEx();
138 h = orbit.getEquinoctialEy();
139 q = orbit.getHx();
140 p = orbit.getHy();
141 lm = MathUtils.normalizeAngle(orbit.getLM(), pi);
142 lv = MathUtils.normalizeAngle(orbit.getLv(), pi);
143 le = MathUtils.normalizeAngle(orbit.getLE(), pi);
144
145 // Retrograde factor [Eq. 2.1.2-(2)]
146 I = retrogradeFactor;
147
148 final T k2 = k.multiply(k);
149 final T h2 = h.multiply(h);
150 final T q2 = q.multiply(q);
151 final T p2 = p.multiply(p);
152
153 // A, B, C parameters [Eq. 2.1.6-(1)]
154 B = FastMath.sqrt(k2.add(h2).negate().add(1.));
155 C = q2.add(p2).add(1);
156
157 // Equinoctial reference frame [Eq. 2.1.4-(1)]
158 final T ooC = C.reciprocal();
159 final T px2 = p.multiply(2.);
160 final T qx2 = q.multiply(2.);
161 final T pq2 = px2.multiply(q);
162 f = new FieldVector3D<>(ooC, new FieldVector3D<>(p2.negate().add(1.).add(q2), pq2, px2.multiply(I).negate()));
163 g = new FieldVector3D<>(ooC, new FieldVector3D<>(pq2.multiply(I), (p2.add(1.).subtract(q2)).multiply(I), qx2));
164 w = new FieldVector3D<>(ooC, new FieldVector3D<>(px2, qx2.negate(), (p2.add(q2).negate().add(1.)).multiply(I)));
165
166 // Direction cosines for central body [Eq. 2.1.9-(1)]
167 alpha = (T) f.getZ();
168 beta = (T) g.getZ();
169 gamma = (T) w.getZ();
170 }
171
172 /** Get the date of the orbit.
173 * @return the date
174 */
175 public FieldAbsoluteDate<T> getDate() {
176 return date;
177 }
178
179 /** Get the definition frame of the orbit.
180 * @return the definition frame
181 */
182 public Frame getFrame() {
183 return frame;
184 }
185
186 /** Get the eccentricity.
187 * @return ecc
188 */
189 public T getEcc() {
190 return ecc;
191 }
192
193 /** Get the Keplerian mean motion.
194 * @return n
195 */
196 public T getMeanMotion() {
197 return n;
198 }
199
200 /** Get the Keplerian period.
201 * @return period
202 */
203 public T getKeplerianPeriod() {
204 return period;
205 }
206
207 /** Get the semi-major axis.
208 * @return the semi-major axis a
209 */
210 public T getSma() {
211 return sma;
212 }
213
214 /** Get the x component of eccentricity vector.
215 * <p>
216 * This element called k in DSST corresponds to ex
217 * for the {@link org.orekit.orbits.EquinoctialOrbit}
218 * </p>
219 * @return k
220 */
221 public T getK() {
222 return k;
223 }
224
225 /** Get the y component of eccentricity vector.
226 * <p>
227 * This element called h in DSST corresponds to ey
228 * for the {@link org.orekit.orbits.EquinoctialOrbit}
229 * </p>
230 * @return h
231 */
232 public T getH() {
233 return h;
234 }
235
236 /** Get the x component of inclination vector.
237 * <p>
238 * This element called q in DSST corresponds to hx
239 * for the {@link org.orekit.orbits.EquinoctialOrbit}
240 * </p>
241 * @return q
242 */
243 public T getQ() {
244 return q;
245 }
246
247 /** Get the y component of inclination vector.
248 * <p>
249 * This element called p in DSST corresponds to hy
250 * for the {@link org.orekit.orbits.EquinoctialOrbit}
251 * </p>
252 * @return p
253 */
254 public T getP() {
255 return p;
256 }
257
258 /** Get the mean longitude.
259 * @return lm
260 */
261 public T getLM() {
262 return lm;
263 }
264
265 /** Get the true longitude.
266 * @return lv
267 */
268 public T getLv() {
269 return lv;
270 }
271
272 /** Get the eccentric longitude.
273 * @return le
274 */
275 public T getLe() {
276 return le;
277 }
278
279 /** Get the retrograde factor.
280 * @return the retrograde factor I
281 */
282 public int getRetrogradeFactor() {
283 return I;
284 }
285
286 /** Get B = sqrt(1 - e²).
287 * @return B
288 */
289 public T getB() {
290 return B;
291 }
292
293 /** Get C = 1 + p² + q².
294 * @return C
295 */
296 public T getC() {
297 return C;
298 }
299
300 /** Get equinoctial frame vector f.
301 * @return f vector
302 */
303 public FieldVector3D<T> getVectorF() {
304 return f;
305 }
306
307 /** Get equinoctial frame vector g.
308 * @return g vector
309 */
310 public FieldVector3D<T> getVectorG() {
311 return g;
312 }
313
314 /** Get equinoctial frame vector w.
315 * @return w vector
316 */
317 public FieldVector3D<T> getVectorW() {
318 return w;
319 }
320
321 /** Get direction cosine α for central body.
322 * @return α
323 */
324 public T getAlpha() {
325 return alpha;
326 }
327
328 /** Get direction cosine β for central body.
329 * @return β
330 */
331 public T getBeta() {
332 return beta;
333 }
334
335 /** Get direction cosine γ for central body.
336 * @return γ
337 */
338 public T getGamma() {
339 return gamma;
340 }
341
342 }