1   /* Copyright 2002-2018 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.orekit.attitudes.AttitudeProvider;
21  import org.orekit.errors.OrekitException;
22  
23  /** This class contains methods to compute propagated coordinates with the SGP4 model.
24   * <p>
25   * The user should not bother in this class since it is handled internaly by the
26   * {@link TLEPropagator}.
27   * </p>
28   * <p>This implementation is largely inspired from the paper and source code <a
29   * href="http://www.celestrak.com/publications/AIAA/2006-6753/">Revisiting Spacetrack
30   * Report #3</a> and is fully compliant with its results and tests cases.</p>
31   * @author Felix R. Hoots, Ronald L. Roehrich, December 1980 (original fortran)
32   * @author David A. Vallado, Paul Crawford, Richard Hujsak, T.S. Kelso (C++ translation and improvements)
33   * @author Fabien Maussion (java translation)
34   */
35  public class SGP4 extends TLEPropagator {
36  
37      /** If perige is less than 220 km, some calculus are avoided. */
38      private boolean lessThan220;
39  
40      /** (1 + eta * cos(M0))³. */
41      private double delM0;
42  
43      // CHECKSTYLE: stop JavadocVariable check
44      private double d2;
45      private double d3;
46      private double d4;
47      private double t3cof;
48      private double t4cof;
49      private double t5cof;
50      private double sinM0;
51      private double omgcof;
52      private double xmcof;
53      private double c5;
54      // CHECKSTYLE: resume JavadocVariable check
55  
56      /** Constructor for a unique initial TLE.
57       * @param initialTLE the TLE to propagate.
58       * @param attitudeProvider provider for attitude computation
59       * @param mass spacecraft mass (kg)
60       * @exception OrekitException if some specific error occurs
61       */
62      public SGP4(final TLE initialTLE, final AttitudeProvider attitudeProvider,
63                         final double mass) throws OrekitException {
64          super(initialTLE, attitudeProvider, mass);
65      }
66  
67      /** Initialization proper to each propagator (SGP or SDP).
68       */
69      protected void sxpInitialize() {
70  
71          // For perigee less than 220 kilometers, the equations are truncated to
72          // linear variation in sqrt a and quadratic variation in mean anomaly.
73          // Also, the c3 term, the delta omega term, and the delta m term are dropped.
74          lessThan220 = perige < 220;
75          if (!lessThan220) {
76              final double c1sq = c1 * c1;
77              delM0 = 1.0 + eta * FastMath.cos(tle.getMeanAnomaly());
78              delM0 *= delM0 * delM0;
79              d2 = 4 * a0dp * tsi * c1sq;
80              final double temp = d2 * tsi * c1 / 3.0;
81              d3 = (17 * a0dp + s4) * temp;
82              d4 = 0.5 * temp * a0dp * tsi * (221 * a0dp + 31 * s4) * c1;
83              t3cof = d2 + 2 * c1sq;
84              t4cof = 0.25 * (3 * d3 + c1 * (12 * d2 + 10 * c1sq));
85              t5cof = 0.2 * (3 * d4 + 12 * c1 * d3 + 6 * d2 * d2 + 15 * c1sq * (2 * d2 + c1sq));
86              sinM0 = FastMath.sin(tle.getMeanAnomaly());
87              if (tle.getE() < 1e-4) {
88                  omgcof = 0.;
89                  xmcof = 0.;
90              } else  {
91                  final double c3 = coef * tsi * TLEConstants.A3OVK2 * xn0dp *
92                                    TLEConstants.NORMALIZED_EQUATORIAL_RADIUS * sini0 / tle.getE();
93                  xmcof = -TLEConstants.TWO_THIRD * coef * tle.getBStar() *
94                          TLEConstants.NORMALIZED_EQUATORIAL_RADIUS / eeta;
95                  omgcof = tle.getBStar() * c3 * FastMath.cos(tle.getPerigeeArgument());
96              }
97          }
98  
99          c5 = 2 * coef1 * a0dp * beta02 * (1 + 2.75 * (etasq + eeta) + eeta * etasq);
100         // initialized
101     }
102 
103     /** Propagation proper to each propagator (SGP or SDP).
104      * @param tSince the offset from initial epoch (min)
105      */
106     protected void sxpPropagate(final double tSince) {
107 
108         // Update for secular gravity and atmospheric drag.
109         final double xmdf = tle.getMeanAnomaly() + xmdot * tSince;
110         final double omgadf = tle.getPerigeeArgument() + omgdot * tSince;
111         final double xn0ddf = tle.getRaan() + xnodot * tSince;
112         omega = omgadf;
113         double xmp = xmdf;
114         final double tsq = tSince * tSince;
115         xnode = xn0ddf + xnodcf * tsq;
116         double tempa = 1 - c1 * tSince;
117         double tempe = tle.getBStar() * c4 * tSince;
118         double templ = t2cof * tsq;
119 
120         if (!lessThan220) {
121             final double delomg = omgcof * tSince;
122             double delm = 1. + eta * FastMath.cos(xmdf);
123             delm = xmcof * (delm * delm * delm - delM0);
124             final double temp = delomg + delm;
125             xmp = xmdf + temp;
126             omega = omgadf - temp;
127             final double tcube = tsq * tSince;
128             final double tfour = tSince * tcube;
129             tempa = tempa - d2 * tsq - d3 * tcube - d4 * tfour;
130             tempe = tempe + tle.getBStar() * c5 * (FastMath.sin(xmp) - sinM0);
131             templ = templ + t3cof * tcube + tfour * (t4cof + tSince * t5cof);
132         }
133 
134         a = a0dp * tempa * tempa;
135         e = tle.getE() - tempe;
136 
137         // A highly arbitrary lower limit on e,  of 1e-6:
138         if (e < 1e-6) {
139             e = 1e-6;
140         }
141 
142         xl = xmp + omega + xnode + xn0dp * templ;
143 
144         i = tle.getI();
145 
146     }
147 
148 }