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