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