1   /* Copyright 2002-2026 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 java.util.Collections;
20  import java.util.List;
21  
22  import org.hipparchus.geometry.euclidean.threed.Vector3D;
23  import org.hipparchus.linear.RealMatrix;
24  import org.hipparchus.util.FastMath;
25  import org.hipparchus.util.MathUtils;
26  import org.hipparchus.util.Pair;
27  import org.hipparchus.util.SinCos;
28  import org.orekit.annotation.DefaultDataContext;
29  import org.orekit.attitudes.Attitude;
30  import org.orekit.attitudes.AttitudeProvider;
31  import org.orekit.attitudes.FrameAlignedProvider;
32  import org.orekit.data.DataContext;
33  import org.orekit.errors.OrekitException;
34  import org.orekit.errors.OrekitMessages;
35  import org.orekit.frames.Frame;
36  import org.orekit.orbits.CartesianOrbit;
37  import org.orekit.orbits.Orbit;
38  import org.orekit.propagation.AbstractMatricesHarvester;
39  import org.orekit.propagation.MatricesHarvester;
40  import org.orekit.propagation.SpacecraftState;
41  import org.orekit.propagation.analytical.AbstractAnalyticalPropagator;
42  import org.orekit.propagation.analytical.tle.generation.FixedPointTleGenerationAlgorithm;
43  import org.orekit.propagation.analytical.tle.generation.TleGenerationAlgorithm;
44  import org.orekit.time.AbsoluteDate;
45  import org.orekit.time.TimeInterval;
46  import org.orekit.time.TimeScale;
47  import org.orekit.utils.DoubleArrayDictionary;
48  import org.orekit.utils.PVCoordinates;
49  import org.orekit.utils.drivers.ParameterDriver;
50  import org.orekit.utils.drivers.ParameterDriversProvider;
51  import org.orekit.utils.TimeSpanMap;
52  
53  /** This class provides elements to propagate TLE's.
54   * <p>
55   * The models used are SGP4 and SDP4, initially proposed by NORAD as the unique convenient
56   * propagator for TLE's. Inputs and outputs of this propagator are only suited for
57   * NORAD two lines elements sets, since it uses estimations and mean values appropriate
58   * for TLE's only.
59   * </p>
60   * <p>
61   * Deep- or near- space propagator is selected internally according to NORAD recommendations
62   * so that the user has not to worry about the used computation methods. One instance is created
63   * for each TLE (this instance can only be get using {@link #selectExtrapolator(TLE)} method,
64   * and can compute {@link PVCoordinates position and velocity coordinates} at any
65   * time. Maximum accuracy is guaranteed in a 24h range period before and after the provided
66   * TLE epoch (of course this accuracy is not really measurable nor predictable: according to
67   * <a href="https://www.celestrak.com/">CelesTrak</a>, the precision is close to one kilometer
68   * and error won't probably rise above 2 km).
69   * </p>
70   * <p>This implementation is largely inspired from the paper and source code <a
71   * href="https://www.celestrak.com/publications/AIAA/2006-6753/">Revisiting Spacetrack
72   * Report #3</a> and is fully compliant with its results and tests cases.</p>
73   * @author Felix R. Hoots, Ronald L. Roehrich, December 1980 (original fortran)
74   * @author David A. Vallado, Paul Crawford, Richard Hujsak, T.S. Kelso (C++ translation and improvements)
75   * @author Fabien Maussion (java translation)
76   * @see TLE
77   */
78  public abstract class TLEPropagator extends AbstractAnalyticalPropagator implements ParameterDriversProvider {
79  
80      // CHECKSTYLE: stop VisibilityModifier check
81  
82      /** Initial state. */
83      protected TLE tle;
84  
85      /** UTC time scale. */
86      protected final TimeScale utc;
87  
88      /** final RAAN. */
89      protected double xnode;
90  
91      /** final semi major axis. */
92      protected double a;
93  
94      /** final eccentricity. */
95      protected double e;
96  
97      /** final inclination. */
98      protected double i;
99  
100     /** final periapsis argument. */
101     protected double omega;
102 
103     /** L from SPTRCK #3. */
104     protected double xl;
105 
106     /** original recovered semi major axis. */
107     protected double a0dp;
108 
109     /** original recovered mean motion. */
110     protected double xn0dp;
111 
112     /** cosinus original inclination. */
113     protected double cosi0;
114 
115     /** cos io squared. */
116     protected double theta2;
117 
118     /** sinus original inclination. */
119     protected double sini0;
120 
121     /** common parameter for mean anomaly (M) computation. */
122     protected double xmdot;
123 
124     /** common parameter for periapsis argument (omega) computation. */
125     protected double omgdot;
126 
127     /** common parameter for raan (OMEGA) computation. */
128     protected double xnodot;
129 
130     /** original eccentricity squared. */
131     protected double e0sq;
132     /** 1 - e2. */
133     protected double beta02;
134 
135     /** sqrt (1 - e2). */
136     protected double beta0;
137 
138     /** periapsis, expressed in KM and ALTITUDE. */
139     protected double perige;
140 
141     /** eta squared. */
142     protected double etasq;
143 
144     /** original eccentricity * eta. */
145     protected double eeta;
146 
147     /** s* new value for the contant s. */
148     protected double s4;
149 
150     /** tsi from SPTRCK #3. */
151     protected double tsi;
152 
153     /** eta from SPTRCK #3. */
154     protected double eta;
155 
156     /** coef for SGP C3 computation. */
157     protected double coef;
158 
159     /** coef for SGP C5 computation. */
160     protected double coef1;
161 
162     /** C1 from SPTRCK #3. */
163     protected double c1;
164 
165     /** C2 from SPTRCK #3. */
166     protected double c2;
167 
168     /** C4 from SPTRCK #3. */
169     protected double c4;
170 
171     /** common parameter for raan (OMEGA) computation. */
172     protected double xnodcf;
173 
174     /** 3/2 * C1. */
175     protected double t2cof;
176 
177     // CHECKSTYLE: resume VisibilityModifier check
178 
179     /** TLE frame. */
180     private final Frame teme;
181 
182     /** All TLEs and masses. */
183     private TimeSpanMap<Pair<TLE, Double>> tlesAndMasses;
184 
185     /** Driver for the ballistic parameter.
186      * @since 14.0
187      */
188     private final ParameterDriver bStarDriver;
189 
190     /** TLE generation algorithm used when resetting TLE from state. */
191     private TleGenerationAlgorithm generationAlgorithm;
192 
193      /** Protected constructor for derived classes.
194      *
195      * <p>This constructor uses the {@link DataContext#getDefault() default data context}.
196      *
197      * @param initialTLE the unique TLE to propagate
198      * @param attitudeProvider provider for attitude computation
199      * @param mass spacecraft mass (kg)
200      * @see #TLEPropagator(TLE, AttitudeProvider, double, Frame)
201      */
202     @DefaultDataContext
203     protected TLEPropagator(final TLE initialTLE, final AttitudeProvider attitudeProvider, final double mass) {
204         this(initialTLE, attitudeProvider, mass,
205              DataContext.getDefault().getFrames().getTEME());
206     }
207 
208     /** Protected constructor for derived classes.
209      * @param initialTLE the unique TLE to propagate
210      * @param attitudeProvider provider for attitude computation
211      * @param mass spacecraft mass (kg)
212      * @param teme the TEME frame to use for propagation.
213      * @since 10.1
214      */
215     protected TLEPropagator(final TLE initialTLE,
216                             final AttitudeProvider attitudeProvider,
217                             final double mass,
218                             final Frame teme) {
219         super(attitudeProvider);
220         setStartDate(initialTLE.getDate());
221         this.utc           = initialTLE.getUtc();
222         initializeTle(initialTLE);
223         this.teme          = teme;
224         this.tlesAndMasses = new TimeSpanMap<>(new Pair<>(tle, mass));
225         this.bStarDriver   = new ParameterDriver(TleGenerationAlgorithm.B_STAR,
226                                                  initialTLE.getBStar(),
227                                                  TleGenerationAlgorithm.B_STAR_SCALE,
228                                                  Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY,
229                                                  TimeInterval.UNLIMITED);
230         bStarDriver.addObserver((previousValue, driver) -> resetBStar());
231         this.generationAlgorithm = getDefaultTleGenerationAlgorithm(initialTLE, utc, teme);
232 
233         // set the initial state
234         final Orbit orbit = propagateOrbit(initialTLE.getDate());
235         final Attitude attitude = attitudeProvider.getAttitude(orbit, orbit.getDate(), orbit.getFrame());
236         super.resetInitialState(new SpacecraftState(orbit, attitude).withMass(mass));
237     }
238 
239     /** Selects the extrapolator to use with the selected TLE.
240      *
241      * <p>This method uses the {@link DataContext#getDefault() default data context}.
242      *
243      * @param tle the TLE to propagate.
244      * @return the correct propagator.
245      * @see #selectExtrapolator(TLE, Frame)
246      */
247     @DefaultDataContext
248     public static TLEPropagator selectExtrapolator(final TLE tle) {
249         return selectExtrapolator(tle, DataContext.getDefault().getFrames().getTEME());
250     }
251 
252     /** Selects the extrapolator to use with the selected TLE.
253      * @param tle the TLE to propagate.
254      * @param teme TEME frame.
255      * @return the correct propagator.
256      * @since 10.1
257      * @see #selectExtrapolator(TLE, Frame, AttitudeProvider)
258      */
259     public static TLEPropagator selectExtrapolator(final TLE tle, final Frame teme) {
260         return selectExtrapolator(tle, teme, FrameAlignedProvider.of(teme));
261     }
262 
263     /** Selects the extrapolator to use with the selected TLE.
264      * @param tle the TLE to propagate.
265      * @param teme TEME frame.
266      * @param attitudeProvider provider for attitude computation
267      * @return the correct propagator.
268      * @since 12.2
269      */
270     public static TLEPropagator selectExtrapolator(final TLE tle, final Frame teme, final AttitudeProvider attitudeProvider) {
271         return selectExtrapolator(tle, attitudeProvider, DEFAULT_MASS, teme);
272     }
273 
274     /** Selects the extrapolator to use with the selected TLE.
275      *
276      * <p>This method uses the {@link DataContext#getDefault() default data context}.
277      *
278      * @param tle the TLE to propagate.
279      * @param attitudeProvider provider for attitude computation
280      * @param mass spacecraft mass (kg)
281      * @return the correct propagator.
282      * @see #selectExtrapolator(TLE, AttitudeProvider, double, Frame)
283      */
284     @DefaultDataContext
285     public static TLEPropagator selectExtrapolator(final TLE tle, final AttitudeProvider attitudeProvider,
286                                                    final double mass) {
287         return selectExtrapolator(tle, attitudeProvider, mass,
288                                   DataContext.getDefault().getFrames().getTEME());
289     }
290 
291     /** Selects the extrapolator to use with the selected TLE.
292      * @param tle the TLE to propagate.
293      * @param attitudeProvider provider for attitude computation
294      * @param mass spacecraft mass (kg)
295      * @param teme the TEME frame to use for propagation.
296      * @return the correct propagator.
297      * @since 10.1
298      */
299     public static TLEPropagator selectExtrapolator(final TLE tle,
300                                                    final AttitudeProvider attitudeProvider,
301                                                    final double mass,
302                                                    final Frame teme) {
303 
304         final double xkeOverN = TLEConstants.XKE / (tle.getMeanMotion() * 60.0);
305         final double a1 = FastMath.cbrt(xkeOverN * xkeOverN);
306         final double cosi0 = FastMath.cos(tle.getI());
307         final double oMe2  = 1.0 - tle.getE() * tle.getE();
308         final double temp = TLEConstants.CK2 * 1.5 * (3 * cosi0 * cosi0 - 1.0) / (oMe2 * FastMath.sqrt(oMe2));
309         final double delta1 = temp / (a1 * a1);
310         final double a0 = a1 * (1.0 - delta1 * (TLEConstants.ONE_THIRD + delta1 * (delta1 * 134.0 / 81.0 + 1.0)));
311         final double delta0 = temp / (a0 * a0);
312 
313         // recover original mean motion :
314         final double xn0dp = tle.getMeanMotion() * 60.0 / (delta0 + 1.0);
315 
316         // Period >= 225 minutes is deep space
317         if (MathUtils.TWO_PI / (xn0dp * TLEConstants.MINUTES_PER_DAY) >= (1.0 / 6.4)) {
318             return new DeepSDP4(tle, attitudeProvider, mass, teme);
319         } else {
320             return new SGP4(tle, attitudeProvider, mass, teme);
321         }
322     }
323 
324     /** Get the TLE generation algorithm used when resetting TLE from state.
325      * @return TLE generation algorithm
326      * @since 14.0
327      */
328     TleGenerationAlgorithm getTleGenerationAlgorithm() {
329         return generationAlgorithm;
330     }
331 
332     /** Set the TLE generation algorithm used when resetting TLE from state.
333      * @param tleGenerationAlgorithm TLE generation algorithm
334      * @since 14.0
335      */
336     public void setTleGenerationAlgorithm(final TleGenerationAlgorithm tleGenerationAlgorithm) {
337         this.generationAlgorithm = tleGenerationAlgorithm;
338     }
339 
340     /** Get the Earth gravity coefficient used for TLE propagation.
341      * @return the Earth gravity coefficient.
342      */
343     public static double getMU() {
344         return TLEConstants.MU;
345     }
346 
347     /** Get the extrapolated position and velocity from an initial TLE.
348      * @param date the final date
349      * @return the final PVCoordinates
350      */
351     public PVCoordinates getPVCoordinates(final AbsoluteDate date) {
352 
353         sxpPropagate(date.durationFrom(tle.getDate()) / 60.0, tle.getBStar());
354 
355         // Compute PV with previous calculated parameters
356         return computePVCoordinates();
357     }
358 
359     /** Computation of the first commons parameters.
360      * @param bStar value of the ballistic coefficient to use for propagation
361      */
362     private void initializeCommons(final double bStar) {
363 
364         // Sine and cosine of inclination
365         final SinCos scI0 = FastMath.sinCos(tle.getI());
366 
367         final double xkeOverN = TLEConstants.XKE / (tle.getMeanMotion() * 60.0);
368         final double a1 = FastMath.cbrt(xkeOverN * xkeOverN);
369         cosi0 = scI0.cos();
370         theta2 = cosi0 * cosi0;
371         final double x3thm1 = 3.0 * theta2 - 1.0;
372         e0sq = tle.getE() * tle.getE();
373         beta02 = 1.0 - e0sq;
374         beta0 = FastMath.sqrt(beta02);
375         final double tval = TLEConstants.CK2 * 1.5 * x3thm1 / (beta0 * beta02);
376         final double delta1 = tval / (a1 * a1);
377         final double a0 = a1 * (1.0 - delta1 * (TLEConstants.ONE_THIRD + delta1 * (1.0 + 134.0 / 81.0 * delta1)));
378         final double delta0 = tval / (a0 * a0);
379 
380         // recover original mean motion and semi-major axis :
381         xn0dp = tle.getMeanMotion() * 60.0 / (delta0 + 1.0);
382         a0dp = a0 / (1.0 - delta0);
383 
384         // Values of s and qms2t :
385         s4 = TLEConstants.S;  // unmodified value for s
386         double q0ms24 = TLEConstants.QOMS2T; // unmodified value for q0ms2T
387 
388         perige = (a0dp * (1 - tle.getE()) - TLEConstants.NORMALIZED_EQUATORIAL_RADIUS) * TLEConstants.EARTH_RADIUS; // perige
389 
390         //  For periapsis below 156 km, the values of s and qoms2t are changed :
391         if (perige < 156.0) {
392             if (perige <= 98.0) {
393                 s4 = 20.0;
394             } else {
395                 s4 = perige - 78.0;
396             }
397             final double temp_val = (120.0 - s4) * TLEConstants.NORMALIZED_EQUATORIAL_RADIUS / TLEConstants.EARTH_RADIUS;
398             final double temp_val_squared = temp_val * temp_val;
399             q0ms24 = temp_val_squared * temp_val_squared;
400             s4 = s4 / TLEConstants.EARTH_RADIUS + TLEConstants.NORMALIZED_EQUATORIAL_RADIUS; // new value for q0ms2T and s
401         }
402 
403         final double pinv = 1.0 / (a0dp * beta02);
404         final double pinvsq = pinv * pinv;
405         tsi = 1.0 / (a0dp - s4);
406         eta = a0dp * tle.getE() * tsi;
407         etasq = eta * eta;
408         eeta = tle.getE() * eta;
409 
410         final double psisq = FastMath.abs(1.0 - etasq); // abs because pow 3.5 needs positive value
411         final double tsi_squared = tsi * tsi;
412         coef = q0ms24 * tsi_squared * tsi_squared;
413         coef1 = coef / FastMath.pow(psisq, 3.5);
414 
415         // C2 and C1 coefficients computation :
416         c2 = coef1 * xn0dp * (a0dp * (1.0 + 1.5 * etasq + eeta * (4.0 + etasq)) +
417              0.75 * TLEConstants.CK2 * tsi / psisq * x3thm1 * (8.0 + 3.0 * etasq * (8.0 + etasq)));
418         c1 = bStar * c2;
419         sini0 = scI0.sin();
420 
421         final double x1mth2 = 1.0 - theta2;
422 
423         // C4 coefficient computation :
424         c4 = 2.0 * xn0dp * coef1 * a0dp * beta02 * (eta * (2.0 + 0.5 * etasq) +
425                 tle.getE() * (0.5 + 2.0 * etasq) -
426                 2 * TLEConstants.CK2 * tsi / (a0dp * psisq) *
427                         (-3.0 * x3thm1 * (1.0 - 2.0 * eeta + etasq * (1.5 - 0.5 * eeta)) +
428                                 0.75 * x1mth2 * (2.0 * etasq - eeta * (1.0 + etasq)) * FastMath.cos(2.0 * tle.getPeriapsisArgument())));
429 
430         final double theta4 = theta2 * theta2;
431         final double temp1 = 3 * TLEConstants.CK2 * pinvsq * xn0dp;
432         final double temp2 = temp1 * TLEConstants.CK2 * pinvsq;
433         final double temp3 = 1.25 * TLEConstants.CK4 * pinvsq * pinvsq * xn0dp;
434 
435         // atmospheric and gravitation coefs :(Mdf and OMEGAdf)
436         xmdot = xn0dp +
437                 0.5 * temp1 * beta0 * x3thm1 +
438                 0.0625 * temp2 * beta0 * (13.0 - 78.0 * theta2 + 137.0 * theta4);
439 
440         final double x1m5th = 1.0 - 5.0 * theta2;
441 
442         omgdot = -0.5 * temp1 * x1m5th +
443                 0.0625 * temp2 * (7.0 - 114.0 * theta2 + 395.0 * theta4) +
444                 temp3 * (3.0 - 36.0 * theta2 + 49.0 * theta4);
445 
446         final double xhdot1 = -temp1 * cosi0;
447 
448         xnodot = xhdot1 + (0.5 * temp2 * (4.0 - 19.0 * theta2) + 2.0 * temp3 * (3.0 - 7.0 * theta2)) * cosi0;
449         xnodcf = 3.5 * beta02 * xhdot1 * c1;
450         t2cof = 1.5 * c1;
451 
452     }
453 
454     /** Retrieves the position and velocity.
455      * @return the computed PVCoordinates.
456      */
457     private PVCoordinates computePVCoordinates() {
458 
459         // Sine and cosine of final periapsis argument
460         final SinCos scOmega = FastMath.sinCos(omega);
461 
462         // Long period periodics
463         final double axn = e * scOmega.cos();
464         double temp = 1.0 / (a * (1.0 - e * e));
465         final double xlcof = 0.125 * TLEConstants.A3OVK2 * sini0 * (3.0 + 5.0 * cosi0) / (1.0 + cosi0);
466         final double aycof = 0.25 * TLEConstants.A3OVK2 * sini0;
467         final double xll = temp * xlcof * axn;
468         final double aynl = temp * aycof;
469         final double xlt = xl + xll;
470         final double ayn = e * scOmega.sin() + aynl;
471         final double elsq = axn * axn + ayn * ayn;
472         final double capu = MathUtils.normalizeAngle(xlt - xnode, FastMath.PI);
473         double epw = capu;
474         double ecosE = 0;
475         double esinE = 0;
476         double sinEPW = 0;
477         double cosEPW = 0;
478 
479         // Dundee changes:  items dependent on cosio get recomputed:
480         final double cosi0Sq = cosi0 * cosi0;
481         final double x3thm1 = 3.0 * cosi0Sq - 1.0;
482         final double x1mth2 = 1.0 - cosi0Sq;
483         final double x7thm1 = 7.0 * cosi0Sq - 1.0;
484 
485         if (e > (1 - 1e-6)) {
486             throw new OrekitException(OrekitMessages.TOO_LARGE_ECCENTRICITY_FOR_PROPAGATION_MODEL, e);
487         }
488 
489         // Solve Kepler's' Equation.
490         final double newtonRaphsonEpsilon = 1e-12;
491         for (int j = 0; j < 10; j++) {
492 
493             boolean doSecondOrderNewtonRaphson = true;
494 
495             final SinCos scEPW = FastMath.sinCos(epw);
496             sinEPW = scEPW.sin();
497             cosEPW = scEPW.cos();
498             ecosE = axn * cosEPW + ayn * sinEPW;
499             esinE = axn * sinEPW - ayn * cosEPW;
500             final double f = capu - epw + esinE;
501             if (FastMath.abs(f) < newtonRaphsonEpsilon) {
502                 break;
503             }
504             final double fdot = 1.0 - ecosE;
505             double delta_epw = f / fdot;
506             if (j == 0) {
507                 final double maxNewtonRaphson = 1.25 * FastMath.abs(e);
508                 doSecondOrderNewtonRaphson = false;
509                 if (delta_epw > maxNewtonRaphson) {
510                     delta_epw = maxNewtonRaphson;
511                 } else if (delta_epw < -maxNewtonRaphson) {
512                     delta_epw = -maxNewtonRaphson;
513                 } else {
514                     doSecondOrderNewtonRaphson = true;
515                 }
516             }
517             if (doSecondOrderNewtonRaphson) {
518                 delta_epw = f / (fdot + 0.5 * esinE * delta_epw);
519             }
520             epw += delta_epw;
521         }
522 
523         // Short period preliminary quantities
524         temp = 1.0 - elsq;
525         final double pl = a * temp;
526         final double r = a * (1.0 - ecosE);
527         double temp2 = a / r;
528         final double betal = FastMath.sqrt(temp);
529         temp = esinE / (1.0 + betal);
530         final double cosu = temp2 * (cosEPW - axn + ayn * temp);
531         final double sinu = temp2 * (sinEPW - ayn - axn * temp);
532         final double u = FastMath.atan2(sinu, cosu);
533         final double sin2u = 2.0 * sinu * cosu;
534         final double cos2u = 2.0 * cosu * cosu - 1.0;
535         final double temp1 = TLEConstants.CK2 / pl;
536         temp2 = temp1 / pl;
537 
538         // Update for short periodics
539         final double rk = r * (1.0 - 1.5 * temp2 * betal * x3thm1) + 0.5 * temp1 * x1mth2 * cos2u;
540         final double uk = u - 0.25 * temp2 * x7thm1 * sin2u;
541         final double xnodek = xnode + 1.5 * temp2 * cosi0 * sin2u;
542         final double xinck = i + 1.5 * temp2 * cosi0 * sini0 * cos2u;
543 
544         // Orientation vectors
545         final SinCos scuk   = FastMath.sinCos(uk);
546         final SinCos scik   = FastMath.sinCos(xinck);
547         final SinCos scnok  = FastMath.sinCos(xnodek);
548         final double sinuk  = scuk.sin();
549         final double cosuk  = scuk.cos();
550         final double sinik  = scik.sin();
551         final double cosik  = scik.cos();
552         final double sinnok = scnok.sin();
553         final double cosnok = scnok.cos();
554         final double xmx = -sinnok * cosik;
555         final double xmy = cosnok * cosik;
556         final double ux  = xmx * sinuk + cosnok * cosuk;
557         final double uy  = xmy * sinuk + sinnok * cosuk;
558         final double uz  = sinik * sinuk;
559 
560         // Position and velocity
561         final double cr = 1000 * rk * TLEConstants.EARTH_RADIUS;
562         final Vector3D pos = new Vector3D(cr * ux, cr * uy, cr * uz);
563 
564         final double rdot   = TLEConstants.XKE * FastMath.sqrt(a) * esinE / r;
565         final double rfdot  = TLEConstants.XKE * FastMath.sqrt(pl) / r;
566         final double xn     = TLEConstants.XKE / (a * FastMath.sqrt(a));
567         final double rdotk  = rdot - xn * temp1 * x1mth2 * sin2u;
568         final double rfdotk = rfdot + xn * temp1 * (x1mth2 * cos2u + 1.5 * x3thm1);
569         final double vx     = xmx * cosuk - cosnok * sinuk;
570         final double vy     = xmy * cosuk - sinnok * sinuk;
571         final double vz     = sinik * cosuk;
572 
573         final double cv = 1000.0 * TLEConstants.EARTH_RADIUS / 60.0;
574         final Vector3D vel = new Vector3D(cv * (rdotk * ux + rfdotk * vx),
575                                           cv * (rdotk * uy + rfdotk * vy),
576                                           cv * (rdotk * uz + rfdotk * vz));
577 
578         return new PVCoordinates(pos, vel);
579 
580     }
581 
582     /** {@inheritDoc} */
583     @Override
584     public List<ParameterDriver> getParametersDrivers() {
585         return Collections.singletonList(bStarDriver);
586     }
587 
588     /** Initialization proper to each propagator (SGP or SDP).
589      * @param bStar value of the ballistic coefficient to use for propagation
590      */
591     protected abstract void sxpInitialize(double bStar);
592 
593     /** Propagation proper to each propagator (SGP or SDP).
594      * @param t the offset from initial epoch (min)
595      * @param bStar value of the ballistic coefficient to use for propagation
596      */
597     protected abstract void sxpPropagate(double t, double bStar);
598 
599     /** {@inheritDoc}
600      * <p>
601      * For TLE propagator, calling this method is only recommended
602      * for covariance propagation when the new <code>state</code>
603      * differs from the previous one by only adding the additional
604      * state containing the derivatives.
605      * </p>
606      */
607     public void resetInitialState(final SpacecraftState state) {
608         super.resetInitialState(state);
609         resetTle(state);
610         tlesAndMasses = new TimeSpanMap<>(new Pair<>(tle, state.getMass()));
611     }
612 
613     /** {@inheritDoc} */
614     protected void resetIntermediateState(final SpacecraftState state, final boolean forward) {
615         resetTle(state);
616         final Pair<TLE, Double> tleAndMass = new Pair<>(tle, state.getMass());
617         if (forward) {
618             tlesAndMasses.addValidAfter(tleAndMass, state.getDate(), false);
619         } else {
620             tlesAndMasses.addValidBefore(tleAndMass, state.getDate(), false);
621         }
622         stateChanged(state);
623     }
624 
625     /** Reset internal TLE from a SpacecraftState.
626      * @param state spacecraft state on which to base new TLE
627      */
628     private void resetTle(final SpacecraftState state) {
629 
630         final TLE newTle = generationAlgorithm.createFromDrivers();
631         initializeTle(newTle);
632     }
633 
634     /** Reset the B-star value from the parameter driver.
635      * @since 14.0
636      */
637     private void resetBStar() {
638         tle = new TLE(tle.getSatelliteNumber(), tle.getClassification(),
639                       tle.getLaunchYear(), tle.getLaunchNumber(), tle.getLaunchPiece(),
640                       tle.getEphemerisType(), tle.getElementNumber(), tle.getDate(),
641                       tle.getMeanMotion(), tle.getMeanMotionFirstDerivative(),
642                       tle.getMeanMotionSecondDerivative(),
643                       tle.getE(), tle.getI(), tle.getPeriapsisArgument(), tle.getRaan(),
644                       tle.getMeanAnomaly(), tle.getRevolutionNumberAtEpoch(),
645                       bStarDriver.getValue());
646         initializeTle(tle);
647     }
648 
649     /** Initialize internal TLE.
650      * @param newTle tle to replace current one
651      */
652     private void initializeTle(final TLE newTle) {
653         tle = newTle;
654         initializeCommons(newTle.getBStar());
655         sxpInitialize(newTle.getBStar());
656     }
657 
658     /** {@inheritDoc} */
659     protected double getMass(final AbsoluteDate date) {
660         return tlesAndMasses.get(date).getValue();
661     }
662 
663     /** {@inheritDoc} */
664     public Orbit propagateOrbit(final AbsoluteDate date) {
665         final TLE closestTle = tlesAndMasses.get(date).getKey();
666         if (!tle.equals(closestTle)) {
667             initializeTle(closestTle);
668         }
669         return new CartesianOrbit(getPVCoordinates(date), teme, date, TLEConstants.MU);
670     }
671 
672     /** Get the underlying TLE.
673      * If there has been calls to #resetInitialState or #resetIntermediateState,
674      * it will not be the same as given to the constructor.
675      * @return underlying TLE
676      */
677     public TLE getTLE() {
678         return tle;
679     }
680 
681     /** {@inheritDoc} */
682     public Frame getFrame() {
683         return teme;
684     }
685 
686     /** {@inheritDoc} */
687     @Override
688     protected AbstractMatricesHarvester createHarvester(final String stmName, final RealMatrix initialStm,
689                                                         final DoubleArrayDictionary initialJacobianColumns) {
690         // Create the harvester
691         final TLEHarvester harvester = new TLEHarvester(this, stmName, initialStm, initialJacobianColumns);
692         // Update the list of additional state provider
693         addAdditionalDataProvider(harvester);
694         // Return the configured harvester
695         return harvester;
696     }
697 
698     /**
699      * Get the names of the parameters in the matrix returned by {@link MatricesHarvester#getParametersJacobian}.
700      * @return names of the parameters (i.e. columns) of the Jacobian matrix
701      */
702     protected List<String> getJacobiansColumnsNames() {
703         return bStarDriver.isSelected() ?
704                Collections.singletonList(bStarDriver.getName()) :
705                Collections.emptyList();
706     }
707 
708     /**
709      * Get the default TLE generation algorithm.
710      * @param templateTLE template TLE
711      * @param utc UTC time scale
712      * @param teme TEME frame
713      * @return a TLE generation algorithm
714      */
715     public static TleGenerationAlgorithm getDefaultTleGenerationAlgorithm(final TLE templateTLE,
716                                                                           final TimeScale utc, final Frame teme) {
717         return new FixedPointTleGenerationAlgorithm(templateTLE,
718                                                     FixedPointTleGenerationAlgorithm.EPSILON_DEFAULT,
719                                                     FixedPointTleGenerationAlgorithm.MAX_ITERATIONS_DEFAULT,
720                                                     FixedPointTleGenerationAlgorithm.SCALE_DEFAULT, utc, teme);
721     }
722 
723 }