1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18 package org.orekit.propagation.analytical;
19
20 import org.hipparchus.geometry.euclidean.threed.Vector3D;
21 import org.hipparchus.linear.RealMatrix;
22 import org.hipparchus.util.FastMath;
23 import org.junit.jupiter.api.AfterEach;
24 import org.junit.jupiter.api.Assertions;
25 import org.junit.jupiter.api.BeforeEach;
26 import org.junit.jupiter.api.Test;
27 import org.orekit.Utils;
28 import org.orekit.attitudes.Attitude;
29 import org.orekit.errors.OrekitException;
30 import org.orekit.forces.gravity.potential.GravityFieldFactory;
31 import org.orekit.forces.gravity.potential.TideSystem;
32 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
33 import org.orekit.frames.Frame;
34 import org.orekit.frames.FramesFactory;
35 import org.orekit.orbits.CartesianOrbit;
36 import org.orekit.orbits.Orbit;
37 import org.orekit.orbits.OrbitType;
38 import org.orekit.orbits.PositionAngleType;
39 import org.orekit.propagation.MatricesHarvester;
40 import org.orekit.propagation.SpacecraftState;
41 import org.orekit.propagation.ToleranceProvider;
42 import org.orekit.time.AbsoluteDate;
43 import org.orekit.utils.PVCoordinates;
44
45 public class EcksteinHechlerStateTransitionMatrixTest {
46
47
48 @Test
49 public void testNullStmName() {
50 Assertions.assertThrows(OrekitException.class, () -> {
51
52
53
54
55 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
56 Vector3D position = new Vector3D(3220103., 69623., 6149822.);
57 Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
58
59 Orbit initialOrbit = new CartesianOrbit(new PVCoordinates(position, velocity),
60 FramesFactory.getEME2000(), initDate, provider.getMu());
61
62 EcksteinHechlerPropagator propagator = new EcksteinHechlerPropagator(initialOrbit, provider);
63 propagator.setupMatricesComputation(null, null, null);
64 });
65 }
66
67 @Test
68 public void testStateJacobian() {
69
70
71
72
73
74 AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
75 Vector3D position = new Vector3D(3220103., 69623., 6149822.);
76 Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
77 Orbit initialOrbit = new CartesianOrbit(new PVCoordinates(position, velocity),
78 FramesFactory.getEME2000(), initDate, provider.getMu());
79
80
81 EcksteinHechlerPropagator propagator = new EcksteinHechlerPropagator(initialOrbit, provider);
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 final SpacecraftState finalState = propagator.propagate(target);
88 RealMatrix dYdY0 = harvester.getStateTransitionMatrix(finalState);
89 Assertions.assertEquals(OrbitType.CARTESIAN, harvester.getOrbitType());
90 Assertions.assertEquals(PositionAngleType.MEAN, harvester.getPositionAngleType());
91
92
93 double[][] dYdY0Ref = new double[6][6];
94 EcksteinHechlerPropagator propagator2;
95 double[] steps = ToleranceProvider.getDefaultToleranceProvider(10.).getTolerances(initialState.getOrbit(), OrbitType.CARTESIAN)[0];
96 for (int i = 0; i < 6; ++i) {
97 propagator2 = new EcksteinHechlerPropagator(shiftState(initialState, OrbitType.CARTESIAN, -4 * steps[i], i).getOrbit(), provider);
98 SpacecraftState sM4h = propagator2.propagate(target);
99 propagator2 = new EcksteinHechlerPropagator(shiftState(initialState, OrbitType.CARTESIAN, -3 * steps[i], i).getOrbit(), provider);
100 SpacecraftState sM3h = propagator2.propagate(target);
101 propagator2 = new EcksteinHechlerPropagator(shiftState(initialState, OrbitType.CARTESIAN, -2 * steps[i], i).getOrbit(), provider);
102 SpacecraftState sM2h = propagator2.propagate(target);
103 propagator2 = new EcksteinHechlerPropagator(shiftState(initialState, OrbitType.CARTESIAN, -1 * steps[i], i).getOrbit(), provider);
104 SpacecraftState sM1h = propagator2.propagate(target);
105 propagator2 = new EcksteinHechlerPropagator(shiftState(initialState, OrbitType.CARTESIAN, +1 * steps[i], i).getOrbit(), provider);
106 SpacecraftState sP1h = propagator2.propagate(target);
107 propagator2 = new EcksteinHechlerPropagator(shiftState(initialState, OrbitType.CARTESIAN, +2 * steps[i], i).getOrbit(), provider);
108 SpacecraftState sP2h = propagator2.propagate(target);
109 propagator2 = new EcksteinHechlerPropagator(shiftState(initialState, OrbitType.CARTESIAN, +3 * steps[i], i).getOrbit(), provider);
110 SpacecraftState sP3h = propagator2.propagate(target);
111 propagator2 = new EcksteinHechlerPropagator(shiftState(initialState, OrbitType.CARTESIAN, +4 * steps[i], i).getOrbit(), provider);
112 SpacecraftState sP4h = propagator2.propagate(target);
113 fillJacobianColumn(dYdY0Ref, i, OrbitType.CARTESIAN, steps[i],
114 sM4h, sM3h, sM2h, sM1h, sP1h, sP2h, sP3h, sP4h);
115 }
116
117
118 for (int i = 0; i < 6; ++i) {
119 for (int j = 0; j < 6; ++j) {
120 if (stateVector[i] != 0) {
121 double error = FastMath.abs((dYdY0.getEntry(i, j) - dYdY0Ref[i][j]) / stateVector[i]) * steps[j];
122 Assertions.assertEquals(0, error, 1.1e-13);
123 }
124 }
125 }
126
127 }
128
129 private void fillJacobianColumn(double[][] jacobian, int column,
130 OrbitType orbitType, double h,
131 SpacecraftState sM4h, SpacecraftState sM3h,
132 SpacecraftState sM2h, SpacecraftState sM1h,
133 SpacecraftState sP1h, SpacecraftState sP2h,
134 SpacecraftState sP3h, SpacecraftState sP4h) {
135 double[] aM4h = stateToArray(sM4h, orbitType)[0];
136 double[] aM3h = stateToArray(sM3h, orbitType)[0];
137 double[] aM2h = stateToArray(sM2h, orbitType)[0];
138 double[] aM1h = stateToArray(sM1h, orbitType)[0];
139 double[] aP1h = stateToArray(sP1h, orbitType)[0];
140 double[] aP2h = stateToArray(sP2h, orbitType)[0];
141 double[] aP3h = stateToArray(sP3h, orbitType)[0];
142 double[] aP4h = stateToArray(sP4h, orbitType)[0];
143 for (int i = 0; i < jacobian.length; ++i) {
144 jacobian[i][column] = ( -3 * (aP4h[i] - aM4h[i]) +
145 32 * (aP3h[i] - aM3h[i]) -
146 168 * (aP2h[i] - aM2h[i]) +
147 672 * (aP1h[i] - aM1h[i])) / (840 * h);
148 }
149 }
150
151 private SpacecraftState shiftState(SpacecraftState state, OrbitType orbitType,
152 double delta, int column) {
153
154 double[][] array = stateToArray(state, orbitType);
155 array[0][column] += delta;
156
157 return arrayToState(array, state.getFrame(), state.getDate(),
158 state.getOrbit().getMu(), state.getAttitude());
159
160 }
161
162 private double[][] stateToArray(SpacecraftState state, OrbitType orbitType) {
163 double[][] array = new double[2][6];
164
165 orbitType.mapOrbitToArray(state.getOrbit(), PositionAngleType.MEAN, array[0], array[1]);
166 return array;
167 }
168
169 private SpacecraftState arrayToState(double[][] array,
170 Frame frame, AbsoluteDate date, double mu,
171 Attitude attitude) {
172 CartesianOrbit orbit = (CartesianOrbit) OrbitType.CARTESIAN.mapArrayToOrbit(array[0], array[1], PositionAngleType.MEAN, date, mu, frame);
173 return new SpacecraftState(orbit, attitude);
174 }
175
176 @BeforeEach
177 public void setUp() {
178 Utils.setDataRoot("regular-data");
179 double mu = 3.9860047e14;
180 double ae = 6.378137e6;
181 double[][] cnm = new double[][] {
182 { 0 }, { 0 }, { -1.08263e-3 }, { 2.54e-6 }, { 1.62e-6 }, { 2.3e-7 }, { -5.5e-7 }
183 };
184 double[][] snm = new double[][] {
185 { 0 }, { 0 }, { 0 }, { 0 }, { 0 }, { 0 }, { 0 }
186 };
187 provider = GravityFieldFactory.getUnnormalizedProvider(ae, mu, TideSystem.UNKNOWN, cnm, snm);
188 }
189
190 @AfterEach
191 public void tearDown() {
192 provider = null;
193 }
194
195 private UnnormalizedSphericalHarmonicsProvider provider;
196
197 }