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