1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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.attitudes.Attitude;
27 import org.orekit.errors.OrekitException;
28 import org.orekit.frames.Frame;
29 import org.orekit.orbits.CartesianOrbit;
30 import org.orekit.orbits.OrbitType;
31 import org.orekit.orbits.PositionAngleType;
32 import org.orekit.propagation.MatricesHarvester;
33 import org.orekit.propagation.SpacecraftState;
34 import org.orekit.propagation.ToleranceProvider;
35 import org.orekit.propagation.analytical.tle.generation.FixedPointTleGenerationAlgorithm;
36 import org.orekit.propagation.analytical.tle.generation.TleGenerationAlgorithm;
37 import org.orekit.time.AbsoluteDate;
38
39 public class TLEStateTransitionMatrixTest {
40
41
42 private TLE tleGPS;
43 private TLE tleSPOT;
44
45 @BeforeEach
46 public void setUp() {
47 Utils.setDataRoot("regular-data");
48
49
50 String line1GPS = "1 11783U 80032A 03300.87313441 .00000062 00000-0 10000-3 0 6416";
51 String line2GPS = "2 11783 62.0472 164.2367 0320924 39.0039 323.3716 2.03455768173530";
52 tleGPS = new TLE(line1GPS, line2GPS);
53
54
55 String line1SPOT = "1 22823U 93061A 03339.49496229 .00000173 00000-0 10336-3 0 133";
56 String line2SPOT = "2 22823 98.4132 359.2998 0017888 100.4310 259.8872 14.18403464527664";
57 tleSPOT = new TLE(line1SPOT, line2SPOT);
58 }
59
60 @Test
61 public void testPropagationSGP4() {
62 doTestStateJacobian(8.35e-10, tleSPOT);
63 }
64
65 @Test
66 public void testPropagationSDP4() {
67 doTestStateJacobian(2.53e-9, tleGPS);
68 }
69
70 @Test
71 public void testNullStmName() {
72 Assertions.assertThrows(OrekitException.class, () -> {
73 TLEPropagator propagator = TLEPropagator.selectExtrapolator(tleSPOT);
74 propagator.setupMatricesComputation(null, null, null);
75 });
76 }
77
78 private void doTestStateJacobian(double tolerance, TLE tle) {
79
80
81 TLEPropagator propagator = TLEPropagator.selectExtrapolator(tle);
82 final SpacecraftState initialState = propagator.getInitialState();
83 final double[] stateVector = new double[6];
84 OrbitType.CARTESIAN.mapOrbitToArray(initialState.getOrbit(), PositionAngleType.MEAN, stateVector, null);
85 final AbsoluteDate target = initialState.getDate().shiftedBy(initialState.getOrbit().getKeplerianPeriod());
86 MatricesHarvester harvester = propagator.setupMatricesComputation("stm", null, null);
87 RealMatrix dYdY0 = harvester.getStateTransitionMatrix(initialState);
88 Assertions.assertNull(dYdY0);
89 final SpacecraftState finalState = propagator.propagate(target);
90 dYdY0 = harvester.getStateTransitionMatrix(finalState);
91
92
93 TleGenerationAlgorithm algorithm = new FixedPointTleGenerationAlgorithm();
94
95
96 double[][] dYdY0Ref = new double[6][6];
97 TLEPropagator propagator2;
98 double[] steps = ToleranceProvider.getDefaultToleranceProvider(10.).getTolerances(initialState.getOrbit(), OrbitType.CARTESIAN)[0];
99 for (int i = 0; i < 6; ++i) {
100 propagator2 = TLEPropagator.selectExtrapolator(algorithm.generate(shiftState(initialState, OrbitType.CARTESIAN, -4 * steps[i], i), tle));
101 SpacecraftState sM4h = propagator2.propagate(target);
102 propagator2 = TLEPropagator.selectExtrapolator(algorithm.generate(shiftState(initialState, OrbitType.CARTESIAN, -3 * steps[i], i), tle));
103 SpacecraftState sM3h = propagator2.propagate(target);
104 propagator2 = TLEPropagator.selectExtrapolator(algorithm.generate(shiftState(initialState, OrbitType.CARTESIAN, -2 * steps[i], i), tle));
105 SpacecraftState sM2h = propagator2.propagate(target);
106 propagator2 = TLEPropagator.selectExtrapolator(algorithm.generate(shiftState(initialState, OrbitType.CARTESIAN, -1 * steps[i], i), tle));
107 SpacecraftState sM1h = propagator2.propagate(target);
108 propagator2 = TLEPropagator.selectExtrapolator(algorithm.generate(shiftState(initialState, OrbitType.CARTESIAN, +1 * steps[i], i), tle));
109 SpacecraftState sP1h = propagator2.propagate(target);
110 propagator2 = TLEPropagator.selectExtrapolator(algorithm.generate(shiftState(initialState, OrbitType.CARTESIAN, +2 * steps[i], i), tle));
111 SpacecraftState sP2h = propagator2.propagate(target);
112 propagator2 = TLEPropagator.selectExtrapolator(algorithm.generate(shiftState(initialState, OrbitType.CARTESIAN, +3 * steps[i], i), tle));
113 SpacecraftState sP3h = propagator2.propagate(target);
114 propagator2 = TLEPropagator.selectExtrapolator(algorithm.generate(shiftState(initialState, OrbitType.CARTESIAN, +4 * steps[i], i), tle));
115 SpacecraftState sP4h = propagator2.propagate(target);
116 fillJacobianColumn(dYdY0Ref, i, OrbitType.CARTESIAN, steps[i],
117 sM4h, sM3h, sM2h, sM1h, sP1h, sP2h, sP3h, sP4h);
118 }
119
120 for (int i = 0; i < 6; ++i) {
121 for (int j = 0; j < 6; ++j) {
122 if (stateVector[i] != 0) {
123 double error = FastMath.abs((dYdY0.getEntry(i, j) - dYdY0Ref[i][j]) / stateVector[i]) * steps[j];
124 Assertions.assertEquals(0, error, tolerance);
125 }
126 }
127 }
128 }
129
130 private void fillJacobianColumn(double[][] jacobian, int column,
131 OrbitType orbitType, double h,
132 SpacecraftState sM4h, SpacecraftState sM3h,
133 SpacecraftState sM2h, SpacecraftState sM1h,
134 SpacecraftState sP1h, SpacecraftState sP2h,
135 SpacecraftState sP3h, SpacecraftState sP4h) {
136 double[] aM4h = stateToArray(sM4h, orbitType)[0];
137 double[] aM3h = stateToArray(sM3h, orbitType)[0];
138 double[] aM2h = stateToArray(sM2h, orbitType)[0];
139 double[] aM1h = stateToArray(sM1h, orbitType)[0];
140 double[] aP1h = stateToArray(sP1h, orbitType)[0];
141 double[] aP2h = stateToArray(sP2h, orbitType)[0];
142 double[] aP3h = stateToArray(sP3h, orbitType)[0];
143 double[] aP4h = stateToArray(sP4h, orbitType)[0];
144 for (int i = 0; i < jacobian.length; ++i) {
145 jacobian[i][column] = ( -3 * (aP4h[i] - aM4h[i]) +
146 32 * (aP3h[i] - aM3h[i]) -
147 168 * (aP2h[i] - aM2h[i]) +
148 672 * (aP1h[i] - aM1h[i])) / (840 * h);
149 }
150 }
151
152 private SpacecraftState shiftState(SpacecraftState state, OrbitType orbitType,
153 double delta, int column) {
154
155 double[][] array = stateToArray(state, orbitType);
156 array[0][column] += delta;
157
158 return arrayToState(array, state.getFrame(), state.getDate(),
159 state.getOrbit().getMu(), state.getAttitude());
160
161 }
162
163
164
165 private double[][] stateToArray(SpacecraftState state, OrbitType orbitType) {
166 double[][] array = new double[2][6];
167
168 orbitType.mapOrbitToArray(state.getOrbit(), PositionAngleType.MEAN, array[0], array[1]);
169 return array;
170 }
171
172
173 private SpacecraftState arrayToState(double[][] array,
174 Frame frame, AbsoluteDate date, double mu,
175 Attitude attitude) {
176 CartesianOrbit orbit = (CartesianOrbit) OrbitType.CARTESIAN.mapArrayToOrbit(array[0], array[1], PositionAngleType.MEAN, date, mu, frame);
177 return new SpacecraftState(orbit, attitude);
178 }
179
180 }