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