1 /* Copyright 2002-2025 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.Field;
21 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
22 import org.hipparchus.util.FastMath;
23 import org.hipparchus.util.MathUtils;
24 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
25 import org.orekit.frames.FieldStaticTransform;
26 import org.orekit.frames.Frame;
27 import org.orekit.propagation.semianalytical.dsst.utilities.FieldAuxiliaryElements;
28 import org.orekit.time.AbsoluteDate;
29
30 /**
31 * This class is a container for the common "field" parameters used in {@link DSSTTesseral}.
32 * <p>
33 * It performs parameters initialization at each integration step for the Tesseral contribution
34 * to the central body gravitational perturbation.
35 * </p>
36 * @author Bryan Cazabonne
37 * @since 10.0
38 * @param <T> type of the field elements
39 */
40 public class FieldDSSTTesseralContext<T extends CalculusFieldElement<T>> extends FieldDSSTGravityContext<T> {
41
42 /** Retrograde factor I.
43 * <p>
44 * DSST model needs equinoctial orbit as internal representation.
45 * Classical equinoctial elements have discontinuities when inclination
46 * is close to zero. In this representation, I = +1. <br>
47 * To avoid this discontinuity, another representation exists and equinoctial
48 * elements can be expressed in a different way, called "retrograde" orbit.
49 * This implies I = -1. <br>
50 * As Orekit doesn't implement the retrograde orbit, I is always set to +1.
51 * But for the sake of consistency with the theory, the retrograde factor
52 * has been kept in the formulas.
53 * </p>
54 */
55 private static final int I = 1;
56
57 /** Central body rotation angle θ. */
58 private T theta;
59
60 /** ecc². */
61 private T e2;
62
63 /** Keplerian period. */
64 private T period;
65
66 /** Ratio of satellite period to central body rotation period. */
67 private T ratio;
68
69 /**
70 * Simple constructor.
71 *
72 * @param auxiliaryElements auxiliary elements related to the current orbit
73 * @param centralBodyFrame rotating body frame
74 * @param provider provider for spherical harmonics
75 * @param maxFrequencyShortPeriodics maximum value for j
76 * @param bodyPeriod central body rotation period (seconds)
77 * @param parameters values of the force model parameters (only 1 values
78 * for each parameters corresponding to state date) obtained by calling
79 * the extract parameter method {@link #extractParameters(double[], AbsoluteDate)}
80 * to selected the right value for state date or by getting the parameters for a specific date
81 */
82 FieldDSSTTesseralContext(final FieldAuxiliaryElements<T> auxiliaryElements,
83 final Frame centralBodyFrame,
84 final UnnormalizedSphericalHarmonicsProvider provider,
85 final int maxFrequencyShortPeriodics,
86 final double bodyPeriod,
87 final T[] parameters) {
88
89 super(auxiliaryElements, centralBodyFrame, provider, parameters);
90
91 // Get field and zero
92 final Field<T> field = auxiliaryElements.getDate().getField();
93 final T zero = field.getZero();
94
95 // Keplerian period
96 final T a = auxiliaryElements.getSma();
97 period = (a.getReal() < 0) ? zero.newInstance(Double.POSITIVE_INFINITY) : getMeanMotion().reciprocal().multiply(MathUtils.TWO_PI);
98
99 // Eccentricity square
100 e2 = auxiliaryElements.getEcc().multiply(auxiliaryElements.getEcc());
101
102 // Central body rotation angle from equation 2.7.1-(3)(4).
103 final FieldStaticTransform<T> t = getBodyFixedToInertialTransform();
104 final FieldVector3D<T> xB = t.transformVector(FieldVector3D.getPlusI(field));
105 final FieldVector3D<T> yB = t.transformVector(FieldVector3D.getPlusJ(field));
106 theta = FastMath.atan2(auxiliaryElements.getVectorF().dotProduct(yB).negate().add((auxiliaryElements.getVectorG().dotProduct(xB)).multiply(I)),
107 auxiliaryElements.getVectorF().dotProduct(xB).add(auxiliaryElements.getVectorG().dotProduct(yB).multiply(I)));
108
109 // Ratio of satellite to central body periods to define resonant terms
110 ratio = period.divide(bodyPeriod);
111 }
112
113 /** Get ecc².
114 * @return e2
115 */
116 public T getE2() {
117 return e2;
118 }
119
120 /** Get Central body rotation angle θ.
121 * @return theta
122 */
123 public T getTheta() {
124 return theta;
125 }
126
127 /** Get μ / a .
128 * @return moa
129 * @deprecated since 12.2 Use getMuoa() instead
130 */
131 @Deprecated
132 public T getMoa() {
133 return getMuoa();
134 }
135
136 /** Get the Keplerian period.
137 * <p>The Keplerian period is computed directly from semi major axis
138 * and central acceleration constant.</p>
139 * @return Keplerian period in seconds, or positive infinity for hyperbolic orbits
140 */
141 public T getOrbitPeriod() {
142 return period;
143 }
144
145 /** Get the ratio of satellite period to central body rotation period.
146 * @return ratio
147 */
148 public T getRatio() {
149 return ratio;
150 }
151 }