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.analytical.tle;
18  
19  import org.hipparchus.CalculusFieldElement;
20  import org.hipparchus.util.FastMath;
21  import org.hipparchus.util.FieldSinCos;
22  import org.hipparchus.util.MathUtils;
23  import org.orekit.attitudes.AttitudeProvider;
24  import org.orekit.frames.Frame;
25  import org.orekit.time.DateTimeComponents;
26  import org.orekit.time.FieldAbsoluteDate;
27  import org.orekit.utils.Constants;
28  
29  /** This class contains methods to compute propagated coordinates with the SDP4 model.
30   * <p>
31   * The user should not bother in this class since it is handled internally by the
32   * {@link TLEPropagator}.
33   * </p>
34   * <p>This implementation is largely inspired from the paper and source code <a
35   * href="https://www.celestrak.com/publications/AIAA/2006-6753/">Revisiting Spacetrack
36   * Report #3</a> and is fully compliant with its results and tests cases.</p>
37   * @author Felix R. Hoots, Ronald L. Roehrich, December 1980 (original fortran)
38   * @author David A. Vallado, Paul Crawford, Richard Hujsak, T.S. Kelso (C++ translation and improvements)
39   * @author Fabien Maussion (java translation)
40   * @author Thomas Paulet (field translation)
41   * @since 11.0
42   * @param <T> type of the field elements
43   */
44  abstract class FieldSDP4<T extends CalculusFieldElement<T>>  extends FieldTLEPropagator<T> {
45  
46      // CHECKSTYLE: stop VisibilityModifier check
47  
48      /** New perigee argument. */
49      protected T omgadf;
50  
51      /** New mean motion. */
52      protected T xn;
53  
54      /** Parameter for xl computation. */
55      protected T xll;
56  
57      /** New eccentricity. */
58      protected T em;
59  
60      /** New inclination. */
61      protected T xinc;
62  
63      // CHECKSTYLE: resume VisibilityModifier check
64  
65      /** Constructor for a unique initial TLE.
66       * @param initialTLE the TLE to propagate.
67       * @param attitudeProvider provider for attitude computation
68       * @param mass spacecraft mass (kg)
69       * @param teme the TEME frame to use for propagation.
70       * @param parameters SGP4 and SDP4 model parameters
71       */
72      protected FieldSDP4(final FieldTLE<T> initialTLE,
73                     final AttitudeProvider attitudeProvider,
74                     final T mass,
75                     final Frame teme,
76                     final T[] parameters) {
77          super(initialTLE, attitudeProvider, mass, teme, parameters);
78      }
79  
80      /** Initialization proper to each propagator (SGP or SDP).
81       * @param parameters model parameters
82       */
83      protected void sxpInitialize(final T[] parameters) {
84          luniSolarTermsComputation();
85      }  // End of initialization
86  
87      /** Propagation proper to each propagator (SGP or SDP).
88       * @param tSince the offset from initial epoch (minutes)
89       * @param parameters model parameters
90       */
91      protected void sxpPropagate(final T tSince, final T[] parameters) {
92  
93          // Update for secular gravity and atmospheric drag
94          final T bStar = parameters[0];
95          omgadf = tle.getPerigeeArgument().add(omgdot.multiply(tSince));
96          final T xnoddf = tle.getRaan().add(xnodot.multiply(tSince));
97          final T tSinceSq = tSince.square();
98          xnode = xnoddf.add(xnodcf.multiply(tSinceSq));
99          xn = xn0dp;
100 
101         // Update for deep-space secular effects
102         xll = tle.getMeanAnomaly().add(xmdot.multiply(tSince));
103 
104         deepSecularEffects(tSince);
105 
106         final T tempa = c1.multiply(tSince).negate().add(1.0);
107         a  = xn.reciprocal().multiply(TLEConstants.XKE).pow(TLEConstants.TWO_THIRD).multiply(tempa).multiply(tempa);
108         em = em.subtract(bStar.multiply(c4).multiply(tSince));
109 
110         // Update for deep-space periodic effects
111         xll = xll.add(xn0dp.multiply(t2cof).multiply(tSinceSq));
112 
113         deepPeriodicEffects(tSince);
114 
115         xl = xll.add(omgadf).add(xnode);
116 
117         // Dundee change:  Reset cosio,  sinio for new xinc:
118         final FieldSinCos<T> scI0 = FastMath.sinCos(xinc);
119         cosi0 = scI0.cos();
120         sini0 = scI0.sin();
121         e = em;
122         i = xinc;
123         omega = omgadf;
124         // end of calculus, go for PV computation
125     }
126 
127     /** Computes SPACETRACK#3 compliant earth rotation angle.
128      * @param date the current date
129      * @return the ERA (rad)
130      */
131     protected double thetaG(final FieldAbsoluteDate<T> date) {
132 
133         // Reference:  The 1992 Astronomical Almanac, page B6.
134         final double omega_E = 1.00273790934;
135         final double jd = date
136                 .getComponents(utc)
137                 .offsetFrom(DateTimeComponents.JULIAN_EPOCH) /
138                 Constants.JULIAN_DAY;
139 
140         // Earth rotations per sidereal day (non-constant)
141         final double UT = (jd + 0.5) % 1;
142         final double seconds_per_day = Constants.JULIAN_DAY;
143         final double jd_2000 = 2451545.0;   /* 1.5 Jan 2000 = JD 2451545. */
144         final double t_cen = (jd - UT - jd_2000) / 36525.;
145         double GMST = 24110.54841 +
146                       t_cen * (8640184.812866 + t_cen * (0.093104 - t_cen * 6.2E-6));
147         GMST = (GMST + seconds_per_day * omega_E * UT) % seconds_per_day;
148         if (GMST < 0.) {
149             GMST += seconds_per_day;
150         }
151 
152         return MathUtils.TWO_PI * GMST / seconds_per_day;
153 
154     }
155 
156     /** Computes luni - solar terms from initial coordinates and epoch.
157      */
158     protected abstract void luniSolarTermsComputation();
159 
160     /** Computes secular terms from current coordinates and epoch.
161      * @param t offset from initial epoch (min)
162      */
163     protected abstract void deepSecularEffects(T t);
164 
165     /** Computes periodic terms from current coordinates and epoch.
166      * @param t offset from initial epoch (min)
167      */
168     protected abstract void deepPeriodicEffects(T t);
169 
170 }