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