1   /* Copyright 2002-2025 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  
18  package org.orekit.propagation.analytical.tle;
19  
20  import org.hipparchus.linear.RealMatrix;
21  import org.hipparchus.util.FastMath;
22  import org.junit.jupiter.api.Assertions;
23  import org.junit.jupiter.api.BeforeEach;
24  import org.junit.jupiter.api.Test;
25  import org.orekit.Utils;
26  import org.orekit.orbits.OrbitType;
27  import org.orekit.orbits.PositionAngleType;
28  import org.orekit.propagation.SpacecraftState;
29  import org.orekit.time.AbsoluteDate;
30  import org.orekit.utils.ParameterDriver;
31  import org.orekit.utils.ParameterDriversList;
32  
33  public class TLEParametersDerivativesTest {
34  
35      /** Spot 5 TLE. */
36      private TLE tleSPOT;
37  
38      @BeforeEach
39      public void setUp() {
40          Utils.setDataRoot("regular-data");
41          // SPOT TLE propagation will use SGP4
42          String line1SPOT = "1 22823U 93061A   03339.49496229  .00000173  00000-0  10336-3 0   133";
43          String line2SPOT = "2 22823  98.4132 359.2998 0017888 100.4310 259.8872 14.18403464527664";
44          tleSPOT = new TLE(line1SPOT, line2SPOT);
45      }
46  
47      @Test
48      public void testBStarEstimation() {
49          doTestParametersDerivatives(TLE.B_STAR, 5.16e-3, tleSPOT);
50      }
51  
52      @Test
53      public void testNoEstimatedParameters() {
54          // compute state Jacobian using PartialDerivatives
55          TLEPropagator propagator = TLEPropagator.selectExtrapolator(tleSPOT);
56          final SpacecraftState initialState = propagator.getInitialState();
57          final double[] stateVector = new double[6];
58          OrbitType.CARTESIAN.mapOrbitToArray(initialState.getOrbit(), PositionAngleType.MEAN, stateVector, null);
59          TLEHarvester harvester = (TLEHarvester) propagator.setupMatricesComputation("stm", null, null);
60          harvester.freezeColumnsNames();
61          RealMatrix dYdP = harvester.getParametersJacobian(initialState);
62          Assertions.assertNull(dYdP);
63      }
64  
65      private void doTestParametersDerivatives(String parameterName, double tolerance, TLE tle) {
66  
67          // compute state Jacobian using PartialDerivatives
68          ParameterDriversList bound = new ParameterDriversList();
69  
70          for (final ParameterDriver driver : tle.getParametersDrivers()) {
71              if (driver.getName().equals(parameterName)) {
72                  driver.setSelected(true);
73                  bound.add(driver);
74              } else {
75                  driver.setSelected(false);
76              }
77          }
78          // compute state Jacobian using PartialDerivatives
79          TLEPropagator propagator = TLEPropagator.selectExtrapolator(tle);
80          final SpacecraftState initialState = propagator.getInitialState();
81          final double[] stateVector = new double[6];
82          OrbitType.CARTESIAN.mapOrbitToArray(initialState.getOrbit(), PositionAngleType.MEAN, stateVector, null);
83          final AbsoluteDate target = initialState.getDate().shiftedBy(initialState.getOrbit().getKeplerianPeriod());
84          TLEHarvester harvester = (TLEHarvester) propagator.setupMatricesComputation("stm", null, null);
85          harvester.freezeColumnsNames();
86          RealMatrix dYdP = harvester.getParametersJacobian(initialState);
87          for (int i = 0; i < 6; ++i) {
88              for (int j = 0; j < 1; ++j) {
89                  Assertions.assertEquals(0.0, dYdP.getEntry(i, j), tolerance);
90              }
91          }
92          final SpacecraftState finalState = propagator.propagate(target);
93          dYdP = harvester.getParametersJacobian(finalState);
94  
95          // compute reference Jacobian using finite differences
96          OrbitType orbitType = OrbitType.CARTESIAN;
97          TLEPropagator propagator2 = TLEPropagator.selectExtrapolator(tle);
98          double[][] dYdPRef = new double[6][1];
99  
100         ParameterDriver selected = bound.getDrivers().get(0);
101         double p0 = selected.getReferenceValue();
102         double h  = selected.getScale();
103         selected.setValue(p0 - 4 * h);
104         propagator2 = TLEPropagator.selectExtrapolator(newTLE(tle, selected.getValue()));
105         SpacecraftState sM4h = propagator2.propagate(target);
106         selected.setValue(p0 - 3 * h);
107         propagator2 = TLEPropagator.selectExtrapolator(newTLE(tle, selected.getValue()));
108         SpacecraftState sM3h = propagator2.propagate(target);
109         selected.setValue(p0 - 2 * h);
110         propagator2 = TLEPropagator.selectExtrapolator(newTLE(tle, selected.getValue()));
111         SpacecraftState sM2h = propagator2.propagate(target);
112         selected.setValue(p0 - 1 * h);
113         propagator2 = TLEPropagator.selectExtrapolator(newTLE(tle, selected.getValue()));
114         SpacecraftState sM1h = propagator2.propagate(target);
115         selected.setValue(p0 + 1 * h);
116         propagator2 = TLEPropagator.selectExtrapolator(newTLE(tle, selected.getValue()));
117         SpacecraftState sP1h = propagator2.propagate(target);
118         selected.setValue(p0 + 2 * h);
119         propagator2 = TLEPropagator.selectExtrapolator(newTLE(tle, selected.getValue()));
120         SpacecraftState sP2h = propagator2.propagate(target);
121         selected.setValue(p0 + 3 * h);
122         propagator2 = TLEPropagator.selectExtrapolator(newTLE(tle, selected.getValue()));
123         SpacecraftState sP3h = propagator2.propagate(target);
124         selected.setValue(p0 + 4 * h);
125         propagator2 = TLEPropagator.selectExtrapolator(newTLE(tle, selected.getValue()));
126         SpacecraftState sP4h = propagator2.propagate(target);
127         fillJacobianColumn(dYdPRef, 0, orbitType, h,
128                            sM4h, sM3h, sM2h, sM1h, sP1h, sP2h, sP3h, sP4h);
129 
130         for (int i = 0; i < 6; ++i) {
131             Assertions.assertEquals(dYdPRef[i][0], dYdP.getEntry(i, 0), FastMath.abs(tolerance));
132         }
133 
134     }
135 
136     private void fillJacobianColumn(double[][] jacobian, int column,
137                                     OrbitType orbitType, double h,
138                                     SpacecraftState sM4h, SpacecraftState sM3h,
139                                     SpacecraftState sM2h, SpacecraftState sM1h,
140                                     SpacecraftState sP1h, SpacecraftState sP2h,
141                                     SpacecraftState sP3h, SpacecraftState sP4h) {
142         double[] aM4h = stateToArray(sM4h, orbitType)[0];
143         double[] aM3h = stateToArray(sM3h, orbitType)[0];
144         double[] aM2h = stateToArray(sM2h, orbitType)[0];
145         double[] aM1h = stateToArray(sM1h, orbitType)[0];
146         double[] aP1h = stateToArray(sP1h, orbitType)[0];
147         double[] aP2h = stateToArray(sP2h, orbitType)[0];
148         double[] aP3h = stateToArray(sP3h, orbitType)[0];
149         double[] aP4h = stateToArray(sP4h, orbitType)[0];
150         for (int i = 0; i < jacobian.length; ++i) {
151             jacobian[i][column] = ( -3 * (aP4h[i] - aM4h[i]) +
152                                     32 * (aP3h[i] - aM3h[i]) -
153                                    168 * (aP2h[i] - aM2h[i]) +
154                                    672 * (aP1h[i] - aM1h[i])) / (840 * h);
155         }
156     }
157 
158     private double[][] stateToArray(SpacecraftState state, OrbitType orbitType) {
159           double[][] array = new double[2][6];
160 
161           orbitType.mapOrbitToArray(state.getOrbit(), PositionAngleType.MEAN, array[0], array[1]);
162           return array;
163       }
164 
165     private TLE newTLE(final TLE template, final double newBStar) {
166         return new TLE(template.getSatelliteNumber(), template.getClassification(),
167                        template.getLaunchYear(), template.getLaunchNumber(), template.getLaunchPiece(),
168                        template.getEphemerisType(), template.getElementNumber(), template.getDate(),
169                        template.getMeanMotion(), template.getMeanMotionFirstDerivative(), template.getMeanMotionSecondDerivative(),
170                        template.getE(), template.getI(), template.getPerigeeArgument(),
171                        template.getRaan(), template.getMeanAnomaly(), template.getRevolutionNumberAtEpoch(),
172                        newBStar);
173     }
174 
175 }