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