1 /* Copyright 2002-2026 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.geometry.euclidean.threed.FieldVector3D;
21 import org.hipparchus.geometry.euclidean.threed.Vector3D;
22 import org.hipparchus.util.FastMath;
23 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
24 import org.orekit.frames.FieldStaticTransform;
25 import org.orekit.frames.Frame;
26 import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
27
28 /**
29 * This class is a container for the common parameters used in {@link DSSTTesseral} and {@link DSSTZonal}.
30 * <p>
31 * It performs parameters initialization at each integration step for the Tesseral and Zonal contribution
32 * to the central body gravitational perturbation.
33 * </p>
34 * @param <T> type of the field elements
35 * @author Bryan Cazabonne
36 * @author Maxime Journot
37 * @since 12.2
38 */
39 public class FieldDSSTGravityContext<T extends CalculusFieldElement<T>> extends FieldForceModelContext<T> {
40
41 /** A = sqrt(μ * a). */
42 private final T A;
43
44 /** Χ = 1 / sqrt(1 - e²) = 1 / B. */
45 private final T chi;
46
47 /** Χ². */
48 private final T chi2;
49
50 // Common factors from equinoctial coefficients
51 /** 2 * a / A . */
52 private final T ax2oA;
53
54 /** 1 / (A * B) . */
55 private final T ooAB;
56
57 /** B / A . */
58 private final T BoA;
59
60 /** B / (A * (1 + B)) . */
61 private final T BoABpo;
62
63 /** C / (2 * A * B) . */
64 private final T Co2AB;
65
66 /** μ / a . */
67 private final T muoa;
68
69 /** R / a . */
70 private final T roa;
71
72 /** Keplerian mean motion. */
73 private final T n;
74
75 /** Direction cosine α. */
76 private final T alpha;
77
78 /** Direction cosine β. */
79 private final T beta;
80
81 /** Direction cosine γ. */
82 private final T gamma;
83
84 /** Transform from body-fixed frame to inertial frame. */
85 private final FieldStaticTransform<T> bodyFixedToInertialTransform;
86
87 /**
88 * Constructor.
89 *
90 * @param auxiliaryElements auxiliary elements related to the current orbit
91 * @param centralBodyFixedFrame rotating body frame
92 * @param provider provider for spherical harmonics
93 * @param parameters values of the force model parameters
94 */
95 FieldDSSTGravityContext(final FieldAuxiliaryElements<T> auxiliaryElements,
96 final Frame centralBodyFixedFrame,
97 final UnnormalizedSphericalHarmonicsProvider provider,
98 final T[] parameters) {
99
100 super(auxiliaryElements);
101
102 // µ
103 final T mu = parameters[0];
104
105 // Semi-major axis
106 final T a = auxiliaryElements.getSma();
107
108 // Keplerian Mean Motion
109 final T absA = FastMath.abs(a);
110 this.n = FastMath.sqrt(mu.divide(absA)).divide(absA);
111
112 // A = sqrt(µ * |a|)
113 this.A = FastMath.sqrt(mu.multiply(absA));
114
115 // Χ = 1 / B
116 final T B = auxiliaryElements.getB();
117 this.chi = auxiliaryElements.getB().reciprocal();
118 this.chi2 = chi.multiply(chi);
119
120 // Common factors from equinoctial coefficients
121 // 2 * a / A
122 this.ax2oA = a.divide(A).multiply(2.);
123 // B / A
124 this.BoA = B.divide(A);
125 // 1 / AB
126 this.ooAB = A.multiply(B).reciprocal();
127 // C / 2AB
128 this.Co2AB = auxiliaryElements.getC().multiply(ooAB).divide(2.);
129 // B / (A * (1 + B))
130 this.BoABpo = BoA.divide(B.add(1.));
131 // &mu / a
132 this.muoa = mu.divide(a);
133 // R / a
134 this.roa = a.divide(provider.getAe()).reciprocal();
135
136
137 // If (centralBodyFrame == null), then centralBodyFrame = orbit frame (see DSSTZonal constructors for more on this).
138 final Frame internalCentralBodyFrame = centralBodyFixedFrame == null ? auxiliaryElements.getFrame() : centralBodyFixedFrame;
139
140 // Transform from body-fixed frame (typically ITRF) to inertial frame
141 bodyFixedToInertialTransform = internalCentralBodyFrame.
142 getStaticTransformTo(auxiliaryElements.getFrame(), auxiliaryElements.getDate());
143
144 final FieldVector3D<T> zB = bodyFixedToInertialTransform.transformVector(Vector3D.PLUS_K);
145
146 // Direction cosines for central body [Eq. 2.1.9-(1)]
147 alpha = FieldVector3D.dotProduct(zB, auxiliaryElements.getVectorF());
148 beta = FieldVector3D.dotProduct(zB, auxiliaryElements.getVectorG());
149 gamma = FieldVector3D.dotProduct(zB, auxiliaryElements.getVectorW());
150 }
151
152 /** A = sqrt(μ * a).
153 * @return A
154 */
155 public T getA() {
156 return A;
157 }
158
159 /** Get Χ = 1 / sqrt(1 - e²) = 1 / B.
160 * @return chi
161 */
162 public T getChi() {
163 return chi;
164 }
165
166 /** Get Χ².
167 * @return chi2
168 */
169 public T getChi2() {
170 return chi2;
171 }
172
173 /** Getter for the ax2oA.
174 * @return the ax2oA
175 */
176 public T getAx2oA() {
177 return ax2oA;
178 }
179
180 /** Get ooAB = 1 / (A * B).
181 * @return ooAB
182 */
183 public T getOoAB() {
184 return ooAB;
185 }
186
187 /** Get B / A.
188 * @return BoA
189 */
190 public T getBoA() {
191 return BoA;
192 }
193
194 /** Get BoABpo = B / A(1 + B).
195 * @return BoABpo
196 */
197 public T getBoABpo() {
198 return BoABpo;
199 }
200
201 /** Get Co2AB = C / 2AB.
202 * @return Co2AB
203 */
204 public T getCo2AB() {
205 return Co2AB;
206 }
207
208 /** Get muoa = μ / a.
209 * @return the muoa
210 */
211 public T getMuoa() {
212 return muoa;
213 }
214
215 /** Get roa = R / a.
216 * @return roa
217 */
218 public T getRoa() {
219 return roa;
220 }
221
222 /** Get the Keplerian mean motion.
223 * <p>The Keplerian mean motion is computed directly from semi major axis
224 * and central acceleration constant.</p>
225 * @return Keplerian mean motion in radians per second
226 */
227 public T getMeanMotion() {
228 return n;
229 }
230
231 /** Get direction cosine α for central body.
232 * @return α
233 */
234 public T getAlpha() {
235 return alpha;
236 }
237
238 /** Get direction cosine β for central body.
239 * @return β
240 */
241 public T getBeta() {
242 return beta;
243 }
244
245 /** Get direction cosine γ for central body.
246 * @return the γ
247 */
248 public T getGamma() {
249 return gamma;
250 }
251
252 /** Getter for the bodyFixedToInertialTransform.
253 * @return the bodyFixedToInertialTransform
254 */
255 public FieldStaticTransform<T> getBodyFixedToInertialTransform() {
256 return bodyFixedToInertialTransform;
257 }
258 }