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