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