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      /** &Chi; = 1 / sqrt(1 - e²) = 1 / B. */
45      private final T chi;
46  
47      /** &Chi;². */
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         // &Chi; = 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 &Chi; = 1 / sqrt(1 - e²) = 1 / B.
160      * @return chi
161      */
162     public T getChi() {
163         return chi;
164     }
165 
166     /** Get &Chi;².
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 }