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