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.forces;
18  
19  import org.hipparchus.CalculusFieldElement;
20  import org.hipparchus.util.FastMath;
21  import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
22  import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
23  import org.orekit.time.AbsoluteDate;
24  
25  /**
26   * This class is a container for the common "field" parameters used in {@link DSSTZonal}.
27   * <p>
28   * It performs parameters initialization at each integration step for the Zonal contribution
29   * to the central body gravitational perturbation.
30   * </p>
31   * @author Bryan Cazabonne
32   * @since 10.0
33   * @param <T> type of the field elements
34   */
35  public class FieldDSSTZonalContext<T extends CalculusFieldElement<T>> extends FieldForceModelContext<T> {
36  
37      // Common factors for potential computation
38      /** A = sqrt(μ * a). */
39      private final T A;
40      /** &Chi; = 1 / sqrt(1 - e²) = 1 / B. */
41      private T X;
42      /** &Chi;². */
43      private T XX;
44      /** &Chi;³. */
45      private T XXX;
46      /** 1 / (A * B) . */
47      private T ooAB;
48      /** B / A . */
49      private T BoA;
50      /** B / A(1 + B) . */
51      private T BoABpo;
52      /** -C / (2 * A * B) . */
53      private T mCo2AB;
54      /** -2 * a / A . */
55      private T m2aoA;
56      /** μ / a . */
57      private T muoa;
58      /** R / a . */
59      private T roa;
60  
61      /** Keplerian mean motion. */
62      private final T n;
63  
64      // Short period terms
65      /** h * k. */
66      private T hk;
67      /** k² - h². */
68      private T k2mh2;
69      /** (k² - h²) / 2. */
70      private T k2mh2o2;
71      /** 1 / (n² * a²). */
72      private T oon2a2;
73      /** 1 / (n² * a) . */
74      private T oon2a;
75      /** χ³ / (n² * a). */
76      private T x3on2a;
77      /** χ / (n² * a²). */
78      private T xon2a2;
79      /** (C * χ) / ( 2 * n² * a² ). */
80      private T cxo2n2a2;
81      /** (χ²) / (n² * a² * (χ + 1 ) ). */
82      private T x2on2a2xp1;
83      /** B * B. */
84      private T BB;
85  
86      /**
87       * Simple constructor.
88       *
89       * @param auxiliaryElements auxiliary elements related to the current orbit
90       * @param provider          provider for spherical harmonics
91       * @param parameters        values of the force model parameters (only 1 values
92       * for each parameters corresponding to state date) obtained by calling the extract
93       * parameter method {@link #extractParameters(double[], AbsoluteDate)}
94       * to selected the right value for state date or by getting the parameters for a specific date
95       */
96      FieldDSSTZonalContext(final FieldAuxiliaryElements<T> auxiliaryElements,
97                                   final UnnormalizedSphericalHarmonicsProvider provider,
98                                   final T[] parameters) {
99  
100         super(auxiliaryElements);
101 
102         final T mu = parameters[0];
103 
104         // Keplerian mean motion
105         final T absA = FastMath.abs(auxiliaryElements.getSma());
106         n = FastMath.sqrt(mu.divide(absA)).divide(absA);
107 
108         A = FastMath.sqrt(mu.multiply(auxiliaryElements.getSma()));
109 
110         // &Chi; = 1 / B
111         X = auxiliaryElements.getB().reciprocal();
112         XX = X.square();
113         XXX = X.multiply(XX);
114 
115         // 1 / AB
116         ooAB = (A.multiply(auxiliaryElements.getB())).reciprocal();
117         // B / A
118         BoA = auxiliaryElements.getB().divide(A);
119         // -C / 2AB
120         mCo2AB = auxiliaryElements.getC().multiply(ooAB).divide(2.).negate();
121         // B / A(1 + B)
122         BoABpo = BoA.divide(auxiliaryElements.getB().add(1.));
123         // -2 * a / A
124         m2aoA = auxiliaryElements.getSma().divide(A).multiply(2.).negate();
125         // μ / a
126         muoa = mu.divide(auxiliaryElements.getSma());
127         // R / a
128         roa = auxiliaryElements.getSma().divide(provider.getAe()).reciprocal();
129 
130         // Short period terms
131 
132         // h * k.
133         hk = auxiliaryElements.getH().multiply(auxiliaryElements.getK());
134         // k² - h².
135         k2mh2 = auxiliaryElements.getK().multiply(auxiliaryElements.getK()).subtract(auxiliaryElements.getH().multiply(auxiliaryElements.getH()));
136         // (k² - h²) / 2.
137         k2mh2o2 = k2mh2.divide(2.);
138         // 1 / (n² * a²) = 1 / (n * A)
139         oon2a2 = (A.multiply(n)).reciprocal();
140         // 1 / (n² * a) = a / (n * A)
141         oon2a = auxiliaryElements.getSma().multiply(oon2a2);
142         // χ³ / (n² * a)
143         x3on2a = XXX.multiply(oon2a);
144         // χ / (n² * a²)
145         xon2a2 = X.multiply(oon2a2);
146         // (C * χ) / ( 2 * n² * a² )
147         cxo2n2a2 = xon2a2.multiply(auxiliaryElements.getC()).divide(2.);
148         // (χ²) / (n² * a² * (χ + 1 ) )
149         x2on2a2xp1 = xon2a2.multiply(X).divide(X.add(1.));
150         // B * B
151         BB = auxiliaryElements.getB().multiply(auxiliaryElements.getB());
152 
153     }
154 
155     /** Get &Chi; = 1 / sqrt(1 - e²) = 1 / B.
156      * @return &Chi;
157      */
158     public T getX() {
159         return X;
160     }
161 
162     /** Get &Chi;².
163      * @return &Chi;².
164      */
165     public T getXX() {
166         return XX;
167     }
168 
169     /** Get &Chi;³.
170      * @return &Chi;³
171      */
172     public T getXXX() {
173         return XXX;
174     }
175 
176     /** Get m2aoA = -2 * a / A.
177      * @return m2aoA
178      */
179     public T getM2aoA() {
180         return m2aoA;
181     }
182 
183     /** Get B / A.
184      * @return BoA
185      */
186     public T getBoA() {
187         return BoA;
188     }
189 
190     /** Get ooAB = 1 / (A * B).
191      * @return ooAB
192      */
193     public T getOoAB() {
194         return ooAB;
195     }
196 
197     /** Get mCo2AB = -C / 2AB.
198      * @return mCo2AB
199      */
200     public T getMCo2AB() {
201         return mCo2AB;
202     }
203 
204     /** Get BoABpo = B / A(1 + B).
205      * @return BoABpo
206      */
207     public T getBoABpo() {
208         return BoABpo;
209     }
210 
211     /** Get μ / a .
212      * @return muoa
213      */
214     public T getMuoa() {
215         return muoa;
216     }
217 
218     /** Get roa = R / a.
219      * @return roa
220      */
221     public T getRoa() {
222         return roa;
223     }
224 
225     /** Get the Keplerian mean motion.
226      * <p>The Keplerian mean motion is computed directly from semi major axis
227      * and central acceleration constant.</p>
228      * @return Keplerian mean motion in radians per second
229      */
230     public T getMeanMotion() {
231         return n;
232     }
233 
234     /** Get h * k.
235      * @return hk
236      */
237     public T getHK() {
238         return hk;
239     }
240 
241     /** Get k² - h².
242      * @return k2mh2
243      */
244     public T getK2MH2() {
245         return k2mh2;
246     }
247 
248     /** Get (k² - h²) / 2.
249      * @return k2mh2o2
250      */
251     public T getK2MH2O2() {
252         return k2mh2o2;
253     }
254 
255     /** Get 1 / (n² * a²).
256      * @return oon2a2
257      */
258     public T getOON2A2() {
259         return oon2a2;
260     }
261 
262     /** Get χ³ / (n² * a).
263      * @return x3on2a
264      */
265     public T getX3ON2A() {
266         return x3on2a;
267     }
268 
269     /** Get χ / (n² * a²).
270      * @return xon2a2
271      */
272     public T getXON2A2() {
273         return xon2a2;
274     }
275 
276     /** Get (C * χ) / ( 2 * n² * a² ).
277      * @return cxo2n2a2
278      */
279     public T getCXO2N2A2() {
280         return cxo2n2a2;
281     }
282 
283     /** Get (χ²) / (n² * a² * (χ + 1 ) ).
284      * @return x2on2a2xp1
285      */
286     public T getX2ON2A2XP1() {
287         return x2on2a2xp1;
288     }
289 
290     /** Get B * B.
291      * @return BB
292      */
293     public T getBB() {
294         return BB;
295     }
296 
297 }