1   /* Copyright 2002-2026 Bryan Cazabonne
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    * Bryan Cazabonne 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.generation;
18  
19  import java.util.List;
20  
21  import org.hipparchus.CalculusFieldElement;
22  import org.hipparchus.analysis.differentiation.Gradient;
23  import org.hipparchus.analysis.differentiation.GradientField;
24  import org.hipparchus.linear.MatrixUtils;
25  import org.hipparchus.linear.RealMatrix;
26  import org.hipparchus.util.FastMath;
27  import org.hipparchus.util.MathUtils;
28  import org.orekit.frames.Frame;
29  import org.orekit.orbits.AbstractOrbitalStateFactory;
30  import org.orekit.orbits.FieldKeplerianOrbit;
31  import org.orekit.orbits.KeplerianOrbit;
32  import org.orekit.orbits.Orbit;
33  import org.orekit.orbits.OrbitParamsType;
34  import org.orekit.orbits.PositionAngleType;
35  import org.orekit.propagation.FieldSpacecraftState;
36  import org.orekit.propagation.SpacecraftState;
37  import org.orekit.propagation.analytical.tle.FieldTLE;
38  import org.orekit.propagation.analytical.tle.FieldTLEPropagator;
39  import org.orekit.propagation.analytical.tle.TLE;
40  import org.orekit.propagation.analytical.tle.TLEConstants;
41  import org.orekit.propagation.conversion.osc2mean.OsculatingToMeanConverter;
42  import org.orekit.time.FieldAbsoluteDate;
43  import org.orekit.time.TimeInterval;
44  import org.orekit.utils.drivers.ParameterDriver;
45  import org.orekit.utils.drivers.ParameterDriversList;
46  import org.orekit.utils.drivers.ParameterDriversList.DelegatingDriver;
47  import org.orekit.utils.TimeStampedFieldPVCoordinates;
48  
49  /**
50   * Base class for generating a TLE.
51   * @author Bryan Cazabonne
52   * @since 12.0
53   */
54  public abstract class TleGenerationAlgorithm extends AbstractOrbitalStateFactory<TLE> {
55  
56      /** Name for mean motion. */
57      public static final String MEAN_MOTION = "TleMeanMotion";
58  
59      /** Name for eccentricity. */
60      public static final String ECCENTRICITY   = "TleEccentricity";
61  
62      /** Name for inclination. */
63      public static final String INCLINATION   = "TleInclination";
64  
65      /** Name for periapsis argument. */
66      public static final String PERIAPSIS_ARGUMENT = "TlePeriapsisArgument";
67  
68      /** Name for right ascension of ascending node. */
69      public static final String RAAN    = "TleRighAscensionAscendingNode";
70  
71      /** Name for mean anomaly. */
72      public static final String MEAN_ANOM = "TleMeanAnomaly";
73  
74      /** Parameter name for B* coefficient. */
75      public static final String B_STAR = "BSTAR";
76  
77      /** B* scaling factor.
78       * <p>
79       * We use a power of 2 to avoid numeric noise introduction
80       * in the multiplications/divisions sequences.
81       * </p>
82       */
83      public static final double B_STAR_SCALE = FastMath.scalb(1.0, -20);
84  
85      /** Number of orbital parameters, i.e. of both rows and columns of the Jacobians. */
86      private static final int DEFAULT_STATE_DIMENSION = 6;
87  
88      /** Template TLE. */
89      private final TLE templateTLE;
90  
91      /** Non-Keplerian drivers (containing only for ballistic coefficient parameter). */
92      private ParameterDriversList nonKeplerianDrivers;
93  
94      /** Osculating to mean orbit converter. */
95      private final OsculatingToMeanConverter converter;
96  
97      /** Default constructor.
98       * @param templateTLE template TLE
99       * @param teme teme frame
100      * @param converter osculating to mean orbit converter
101      */
102     protected TleGenerationAlgorithm(final TLE templateTLE,  final Frame teme,
103                                      final OsculatingToMeanConverter converter) {
104         super(null, createOrbitalParametersDrivers(templateTLE), teme, PositionAngleType.MEAN,
105               templateTLE.getDate(), TLEConstants.MU);
106         this.templateTLE = templateTLE;
107 
108         // create model parameter drivers
109         nonKeplerianDrivers = new ParameterDriversList();
110         nonKeplerianDrivers.add(new ParameterDriver(B_STAR, templateTLE.getBStar(), B_STAR_SCALE,
111                                                     Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY,
112                                                     TimeInterval.UNLIMITED));
113 
114         // conversion algorithm
115         this.converter = converter;
116 
117     }
118 
119     /** Get the template TLE.
120      * @return template TLE
121      */
122     public TLE getTemplateTLE() {
123         return templateTLE;
124     }
125 
126     /** Get the osculating to mean orbit converter.
127      * @return osculating to mean orbit converter
128      * @since 14.0
129      */
130     public OsculatingToMeanConverter getConverter() {
131         return converter;
132     }
133 
134     /** {@inheritDoc} */
135     @Override
136     public ParameterDriversList getNonKeplerianParametersDrivers() {
137         return nonKeplerianDrivers;
138     }
139 
140     /** Create orbital parameter drivers.
141      * @param tle reference TLE
142      * @return drivers
143      */
144     private static ParameterDriversList createOrbitalParametersDrivers(final TLE tle) {
145         final ParameterDriversList drivers = new ParameterDriversList();
146         drivers.add(new ParameterDriver(MEAN_MOTION, tle.getMeanMotion(),
147                                         FastMath.scalb(1.0, -32),
148                                         0, Double.POSITIVE_INFINITY, TimeInterval.UNLIMITED));
149         drivers.add(new ParameterDriver(ECCENTRICITY, tle.getE(),
150                                         FastMath.scalb(1.0, -22),
151                                         0.0, 1.0, TimeInterval.UNLIMITED));
152         drivers.add(new ParameterDriver(INCLINATION, tle.getI(),
153                                         FastMath.scalb(1.0, -22),
154                                         0, FastMath.PI, TimeInterval.UNLIMITED));
155         drivers.add(new ParameterDriver(PERIAPSIS_ARGUMENT, tle.getPeriapsisArgument(),
156                                         FastMath.scalb(1.0, -22),
157                                         Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, TimeInterval.UNLIMITED));
158         drivers.add(new ParameterDriver(RAAN, tle.getRaan(),
159                                         FastMath.scalb(1.0, -22),
160                                         Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, TimeInterval.UNLIMITED));
161         drivers.add(new ParameterDriver(MEAN_ANOM, tle.getMeanAnomaly(),
162                                         FastMath.scalb(1.0, -22),
163                                         Double.NEGATIVE_INFINITY, Double.POSITIVE_INFINITY, TimeInterval.UNLIMITED));
164         return drivers;
165     }
166 
167     /** {@inheritDoc} */
168     @Override
169     protected double[] toArray(final Orbit orbit) {
170 
171         // fix both frame and type
172         final Orbit mean               = converter.convertToMean(orbit);
173         final Orbit partiallyConverted = orbit.getFrame() == getFrame() ? mean : mean.inFrame(getFrame());
174         final Orbit fullyConverted     = OrbitParamsType.KEPLERIAN.convertType(partiallyConverted);
175 
176         // retrieve orbital parameters
177         final double[] stateVector = new double[6];
178         OrbitParamsType.KEPLERIAN.mapOrbitToArray(fullyConverted, PositionAngleType.MEAN, stateVector, null);
179 
180         // TLE uses mean motion as first parameter, not semi major axis as Keplerian orbit
181         stateVector[0] = fullyConverted.getKeplerianMeanMotion();
182 
183         return stateVector;
184 
185     }
186 
187     /** {@inheritDoc} */
188     @Override
189     public TLE createFromDrivers() {
190 
191         final List<DelegatingDriver> drivers = getOrbitalParametersDrivers().getDrivers();
192 
193         // adjust revolution number
194         // as neither SDP4 nor SGP4 use mean motion derivatives, we intentionally ignore them here
195         final double latArg0 =
196             MathUtils.normalizeAngle(drivers.get(3).getValue() + drivers.get(5).getValue(),
197                                      FastMath.PI);
198         final double deltaT   = getDate().durationFrom(templateTLE.getDate());
199         final double latArg1  = latArg0 + deltaT * drivers.get(0).getValue();
200         final int    deltaRev = (int) FastMath.floor(latArg1 / MathUtils.TWO_PI);
201 
202         return new TLE(templateTLE.getSatelliteNumber(), templateTLE.getClassification(),
203                        templateTLE.getLaunchYear(), templateTLE.getLaunchNumber(), templateTLE.getLaunchPiece(),
204                        templateTLE.getEphemerisType(),
205                        templateTLE.getElementNumber() + 1,
206                        getDate(),
207                        drivers.get(0).getValue(),
208                        templateTLE.getMeanMotionFirstDerivative(), templateTLE.getMeanMotionSecondDerivative(),
209                        drivers.get(1).getValue(),
210                        drivers.get(2).getValue(),
211                        drivers.get(3).getValue(),
212                        drivers.get(4).getValue(),
213                        drivers.get(5).getValue(),
214                        templateTLE.getRevolutionNumberAtEpoch() + deltaRev,
215                        getBStar(),
216                        templateTLE.getUtc());
217     }
218 
219     /** {@inheritDoc}
220      * <p>
221      * The TLE orbital elements are related to the Cartesian coordinates by the SGP4/SDP4
222      * model itself, which has no closed-form derivatives, so the Jacobian is obtained by
223      * automatic differentiation: the six elements of the TLE built from the current drivers
224      * are turned into {@link Gradient} variables, and the resulting TLE is evaluated at its
225      * own epoch.
226      * </p>
227      */
228     @Override
229     public RealMatrix getJacobianWrtParameters() {
230         return getJacobianWrtParameters(createFromDrivers());
231     }
232 
233     /** Get the Jacobian of the Cartesian coordinates with respect to the orbital elements of a TLE.
234      * <p>
235      * The TLE orbital elements are related to the Cartesian coordinates by the SGP4/SDP4 model
236      * itself, which has no closed-form derivatives, so the Jacobian is obtained by automatic
237      * differentiation: the six orbital elements are turned into {@link Gradient} variables and
238      * the TLE is evaluated at its own epoch. All the remaining data (identification, mean motion
239      * derivatives, B*) are held constant.
240      * </p>
241      * @param tle TLE holding the orbital elements the Jacobian is computed with respect to
242      * @return jacobian matrix dC/dB, at the TLE epoch and in the TEME frame
243      */
244     public RealMatrix getJacobianWrtParameters(final TLE tle) {
245 
246         // evaluate the Cartesian coordinates from a TLE whose orbital elements are variables
247         final TimeStampedFieldPVCoordinates<Gradient> pv =
248             FieldTLEPropagator.selectExtrapolator(toGradient(tle), getFrame()).
249             getBaseInitialState().
250             getPVCoordinates();
251 
252         // gather the derivatives of each Cartesian coordinate into a row
253         final RealMatrix jacobian = MatrixUtils.createRealMatrix(DEFAULT_STATE_DIMENSION, DEFAULT_STATE_DIMENSION);
254         jacobian.setRow(0, pv.getPosition().getX().getGradient());
255         jacobian.setRow(1, pv.getPosition().getY().getGradient());
256         jacobian.setRow(2, pv.getPosition().getZ().getGradient());
257         jacobian.setRow(3, pv.getVelocity().getX().getGradient());
258         jacobian.setRow(4, pv.getVelocity().getY().getGradient());
259         jacobian.setRow(5, pv.getVelocity().getZ().getGradient());
260 
261         return jacobian;
262 
263     }
264 
265     /** Convert a TLE into one whose orbital elements are gradient variables.
266      * @param tle TLE to convert
267      * @return converted TLE, whose orbital elements carry their own derivatives
268      */
269     // FIXME: should this one be in TLE class instead ?
270     private static FieldTLE<Gradient> toGradient(final TLE tle) {
271         final GradientField field = GradientField.getField(DEFAULT_STATE_DIMENSION);
272         return new FieldTLE<>(tle.getSatelliteNumber(), tle.getClassification(),
273                               tle.getLaunchYear(), tle.getLaunchNumber(), tle.getLaunchPiece(),
274                               tle.getEphemerisType(), tle.getElementNumber(),
275                               new FieldAbsoluteDate<>(field, tle.getDate()),
276                               Gradient.variable(DEFAULT_STATE_DIMENSION, 0, tle.getMeanMotion()),
277                               Gradient.constant(DEFAULT_STATE_DIMENSION, tle.getMeanMotionFirstDerivative()),
278                               Gradient.constant(DEFAULT_STATE_DIMENSION, tle.getMeanMotionSecondDerivative()),
279                               Gradient.variable(DEFAULT_STATE_DIMENSION, 1, tle.getE()),
280                               Gradient.variable(DEFAULT_STATE_DIMENSION, 2, tle.getI()),
281                               Gradient.variable(DEFAULT_STATE_DIMENSION, 3, tle.getPeriapsisArgument()),
282                               Gradient.variable(DEFAULT_STATE_DIMENSION, 4, tle.getRaan()),
283                               Gradient.variable(DEFAULT_STATE_DIMENSION, 5, tle.getMeanAnomaly()),
284                               tle.getRevolutionNumberAtEpoch(),
285                               Gradient.constant(DEFAULT_STATE_DIMENSION, tle.getBStar()),
286                               tle.getUtc());
287     }
288 
289     /** Get the current B-star value.
290      * @return current B-star value
291      */
292     protected double getBStar() {
293         return nonKeplerianDrivers.getDrivers().getFirst().getValue();
294     }
295 
296     /**
297      * Generate a TLE from a given spacecraft state and a template TLE.
298      * <p>
299      * The template TLE is only used to get identifiers like satellite
300      * number, launch year, etc.
301      * In other words, the keplerian elements contained in the generated
302      * TLE are based on the provided state and not the template TLE.
303      * </p>
304      * @param state spacecraft state
305      * @param newTemplateTLE template TLE
306      * @return a TLE corresponding to the given state
307      */
308     public TLE generate(final SpacecraftState state, final TLE newTemplateTLE) {
309         final KeplerianOrbit mean =
310             (KeplerianOrbit) OrbitParamsType.KEPLERIAN.convertType(converter.convertToMean(state.getOrbit()));
311         return TleGenerationUtil.newTLE(mean, newTemplateTLE);
312     }
313 
314     /**
315      * Generate a TLE from a given spacecraft state and a template TLE.
316      * <p>
317      * The template TLE is only used to get identifiers like satellite
318      * number, launch year, etc.
319      * In other words, the keplerian elements contained in the generated
320      * TLE are based on the provided state and not the template TLE.
321      * </p>
322      * @param <T> type of the elements
323      * @param state spacecraft state
324      * @param newTemplateTLE template TLE
325      * @return a TLE corresponding to the given state
326      */
327     public <T extends CalculusFieldElement<T>> FieldTLE<T> generate(final FieldSpacecraftState<T> state,
328                                                                     final FieldTLE<T> newTemplateTLE) {
329         final FieldKeplerianOrbit<T> mean =
330             (FieldKeplerianOrbit<T>) OrbitParamsType.KEPLERIAN.convertType(converter.convertToMean(state.getOrbit()));
331         return TleGenerationUtil.newTLE(mean, newTemplateTLE);
332     }
333 
334     /** {@inheritDoc} */
335     @Override
336     public TleGenerationAlgorithm clone() {
337 
338         final TleGenerationAlgorithm clone = (TleGenerationAlgorithm) super.clone();
339 
340         // de-couple b-star driver
341         final ParameterDriversList newDrivers = new ParameterDriversList();
342         final ParameterDriver driver = nonKeplerianDrivers.getDrivers().getFirst();
343         newDrivers.add(new ParameterDriver(driver.getName(), driver.getValue(), driver.getScale(),
344                                            driver.getMinValue(), driver.getMaxValue(),
345                                            driver.getValidity()));
346         clone.nonKeplerianDrivers = newDrivers;
347 
348         return clone;
349 
350     }
351 
352 }