1   /* Copyright 2002-2021 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   */
43  abstract class FieldSDP4<T extends CalculusFieldElement<T>>  extends FieldTLEPropagator<T> {
44  
45      // CHECKSTYLE: stop VisibilityModifier check
46  
47      /** New perigee argument. */
48      protected T omgadf;
49  
50      /** New mean motion. */
51      protected T xn;
52  
53      /** Parameter for xl computation. */
54      protected T xll;
55  
56      /** New eccentricity. */
57      protected T em;
58  
59      /** New inclination. */
60      protected T xinc;
61  
62      // CHECKSTYLE: resume VisibilityModifier check
63  
64      /** Constructor for a unique initial TLE.
65       * @param initialTLE the TLE to propagate.
66       * @param attitudeProvider provider for attitude computation
67       * @param mass spacecraft mass (kg)
68       * @param teme the TEME frame to use for propagation.
69       * @param parameters SGP4 and SDP4 model parameters
70       */
71      protected FieldSDP4(final FieldTLE<T> initialTLE,
72                     final AttitudeProvider attitudeProvider,
73                     final T mass,
74                     final Frame teme,
75                     final T[] parameters) {
76          super(initialTLE, attitudeProvider, mass, teme, parameters);
77      }
78  
79      /** Initialization proper to each propagator (SGP or SDP).
80       * @param parameters model parameters
81       */
82      protected void sxpInitialize(final T[] parameters) {
83          luniSolarTermsComputation();
84      }  // End of initialization
85  
86      /** Propagation proper to each propagator (SGP or SDP).
87       * @param tSince the offset from initial epoch (minutes)
88       * @param parameters model parameters
89       */
90      protected void sxpPropagate(final T tSince, final T[] parameters) {
91  
92          // Update for secular gravity and atmospheric drag
93          final T bStar = parameters[0];
94          omgadf = tle.getPerigeeArgument().add(omgdot.multiply(tSince));
95          final T xnoddf = tle.getRaan().add(xnodot.multiply(tSince));
96          final T tSinceSq = tSince.multiply(tSince);
97          xnode = xnoddf.add(xnodcf.multiply(tSinceSq));
98          xn = xn0dp;
99  
100         // Update for deep-space secular effects
101         xll = tle.getMeanAnomaly().add(xmdot.multiply(tSince));
102 
103         deepSecularEffects(tSince);
104 
105         final T tempa = c1.multiply(tSince).negate().add(1.0);
106         a  = xn.reciprocal().multiply(TLEConstants.XKE).pow(TLEConstants.TWO_THIRD).multiply(tempa).multiply(tempa);
107         em = em.subtract(bStar.multiply(c4).multiply(tSince));
108 
109         // Update for deep-space periodic effects
110         xll = xll.add(xn0dp.multiply(t2cof).multiply(tSinceSq));
111 
112         deepPeriodicEffects(tSince);
113 
114         xl = xll.add(omgadf).add(xnode);
115 
116         // Dundee change:  Reset cosio,  sinio for new xinc:
117         final FieldSinCos<T> scI0 = FastMath.sinCos(xinc);
118         cosi0 = scI0.cos();
119         sini0 = scI0.sin();
120         e = em;
121         i = xinc;
122         omega = omgadf;
123         // end of calculus, go for PV computation
124     }
125 
126     /** Computes SPACETRACK#3 compliant earth rotation angle.
127      * @param date the current date
128      * @return the ERA (rad)
129      */
130     protected double thetaG(final FieldAbsoluteDate<T> date) {
131 
132         // Reference:  The 1992 Astronomical Almanac, page B6.
133         final double omega_E = 1.00273790934;
134         final double jd = date
135                 .getComponents(utc)
136                 .offsetFrom(DateTimeComponents.JULIAN_EPOCH) /
137                 Constants.JULIAN_DAY;
138 
139         // Earth rotations per sidereal day (non-constant)
140         final double UT = (jd + 0.5) % 1;
141         final double seconds_per_day = Constants.JULIAN_DAY;
142         final double jd_2000 = 2451545.0;   /* 1.5 Jan 2000 = JD 2451545. */
143         final double t_cen = (jd - UT - jd_2000) / 36525.;
144         double GMST = 24110.54841 +
145                       t_cen * (8640184.812866 + t_cen * (0.093104 - t_cen * 6.2E-6));
146         GMST = (GMST + seconds_per_day * omega_E * UT) % seconds_per_day;
147         if (GMST < 0.) {
148             GMST += seconds_per_day;
149         }
150 
151         return MathUtils.TWO_PI * GMST / seconds_per_day;
152 
153     }
154 
155     /** Computes luni - solar terms from initial coordinates and epoch.
156      */
157     protected abstract void luniSolarTermsComputation();
158 
159     /** Computes secular terms from current coordinates and epoch.
160      * @param t offset from initial epoch (min)
161      */
162     protected abstract void deepSecularEffects(T t);
163 
164     /** Computes periodic terms from current coordinates and epoch.
165      * @param t offset from initial epoch (min)
166      */
167     protected abstract void deepPeriodicEffects(T t);
168 
169 }