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 /** Χ = 1 / sqrt(1 - e²) = 1 / B. */
41 private T X;
42 /** Χ². */
43 private T XX;
44 /** Χ³. */
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 // Χ = 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 Χ = 1 / sqrt(1 - e²) = 1 / B.
156 * @return Χ
157 */
158 public T getX() {
159 return X;
160 }
161
162 /** Get Χ².
163 * @return Χ².
164 */
165 public T getXX() {
166 return XX;
167 }
168
169 /** Get Χ³.
170 * @return Χ³
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 }