1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.propagation.semianalytical.dsst;
18
19 import java.io.FileNotFoundException;
20 import java.io.UnsupportedEncodingException;
21 import java.util.Collections;
22
23 import org.hipparchus.linear.RealMatrix;
24 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
25 import org.hipparchus.util.FastMath;
26 import org.junit.Assert;
27 import org.junit.Before;
28 import org.junit.Test;
29 import org.orekit.Utils;
30 import org.orekit.attitudes.Attitude;
31 import org.orekit.bodies.CelestialBodyFactory;
32 import org.orekit.errors.OrekitException;
33 import org.orekit.forces.gravity.potential.GravityFieldFactory;
34 import org.orekit.forces.gravity.potential.SHMFormatReader;
35 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
36 import org.orekit.forces.radiation.RadiationSensitive;
37 import org.orekit.frames.Frame;
38 import org.orekit.frames.FramesFactory;
39 import org.orekit.orbits.EquinoctialOrbit;
40 import org.orekit.orbits.KeplerianOrbit;
41 import org.orekit.orbits.Orbit;
42 import org.orekit.orbits.OrbitType;
43 import org.orekit.orbits.PositionAngle;
44 import org.orekit.propagation.MatricesHarvester;
45 import org.orekit.propagation.PropagationType;
46 import org.orekit.propagation.SpacecraftState;
47 import org.orekit.propagation.integration.AdditionalDerivativesProvider;
48 import org.orekit.propagation.numerical.NumericalPropagator;
49 import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
50 import org.orekit.propagation.semianalytical.dsst.forces.DSSTSolarRadiationPressure;
51 import org.orekit.propagation.semianalytical.dsst.forces.DSSTTesseral;
52 import org.orekit.propagation.semianalytical.dsst.forces.DSSTThirdBody;
53 import org.orekit.propagation.semianalytical.dsst.forces.DSSTZonal;
54 import org.orekit.time.AbsoluteDate;
55 import org.orekit.utils.Constants;
56
57
58 public class DSSTStateTransitionMatrixGeneratorTest {
59
60 @Before
61 public void setUp() {
62 Utils.setDataRoot("regular-data:potential/shm-format");
63 GravityFieldFactory.addPotentialCoefficientsReader(new SHMFormatReader("^eigen_cg03c_coef$", false));
64 }
65
66 @Test
67 public void testInterrupt() {
68
69 UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(5, 5);
70 double dt = 900;
71 double dP = 0.001;
72
73
74 DSSTPropagator propagator1 = setUpPropagator(PropagationType.MEAN, dP, provider);
75 AbsoluteDate t0 = propagator1.getInitialState().getDate();
76 propagator1.
77 getAllForceModels().
78 forEach(fm -> fm.
79 getParametersDrivers().
80 stream().
81 filter(d -> d.getName().equals(RadiationSensitive.REFLECTION_COEFFICIENT)).
82 forEach(d -> d.setSelected(true)));
83 final MatricesHarvester harvester1 = propagator1.setupMatricesComputation("stm", null, null);
84 initializeShortPeriod(harvester1, propagator1);
85 final SpacecraftState state1 = propagator1.propagate(t0.shiftedBy(dt));
86 final RealMatrix stm1 = harvester1.getStateTransitionMatrix(state1);
87 final RealMatrix jacobian1 = harvester1.getParametersJacobian(state1);
88
89
90 DSSTPropagator propagator2 = setUpPropagator(PropagationType.MEAN, dP, provider);
91
92
93 final DSSTStateTransitionMatrixGenerator dummyStmGenerator =
94 new DSSTStateTransitionMatrixGenerator("dummy-1",
95 Collections.emptyList(),
96 propagator2.getAttitudeProvider());
97 propagator2.addAdditionalDerivativesProvider(dummyStmGenerator);
98 propagator2.setInitialState(propagator2.getInitialState().addAdditionalState(dummyStmGenerator.getName(), new double[36]),
99 propagator2.getPropagationType());
100 propagator2.addAdditionalDerivativesProvider(new DSSTIntegrableJacobianColumnGenerator(dummyStmGenerator, "dummy-2"));
101 propagator2.setInitialState(propagator2.getInitialState().addAdditionalState("dummy-2", new double[6]),
102 propagator2.getPropagationType());
103 propagator2.addAdditionalDerivativesProvider(new AdditionalDerivativesProvider() {
104 public String getName() { return "dummy-3"; }
105 public int getDimension() { return 1; }
106 public double[] derivatives(SpacecraftState s) { return new double[1]; }
107 });
108 propagator2.setInitialState(propagator2.getInitialState().addAdditionalState("dummy-3", new double[1]),
109 propagator2.getPropagationType());
110 propagator2.
111 getAllForceModels().
112 forEach(fm -> fm.
113 getParametersDrivers().
114 stream().
115 filter(d -> d.getName().equals(RadiationSensitive.REFLECTION_COEFFICIENT)).
116 forEach(d -> d.setSelected(true)));
117 final MatricesHarvester harvester2 = propagator2.setupMatricesComputation("stm", null, null);
118 initializeShortPeriod(harvester2, propagator2);
119 final SpacecraftState intermediate = propagator2.propagate(t0.shiftedBy(dt / 2));
120 final RealMatrix stmI = harvester2.getStateTransitionMatrix(intermediate);
121 final RealMatrix jacobianI = harvester2.getParametersJacobian(intermediate);
122
123
124 Assert.assertEquals(0.158482, stmI.subtract(stm1).getNorm1() / stm1.getNorm1(), 1.0e-6);
125 Assert.assertEquals(0.499959, jacobianI.subtract(jacobian1).getNorm1() / jacobian1.getNorm1(), 1.0e-6);
126
127
128 final SpacecraftState state2 = propagator2.propagate(t0.shiftedBy(dt));
129 final RealMatrix stm2 = harvester2.getStateTransitionMatrix(state2);
130 final RealMatrix jacobian2 = harvester2.getParametersJacobian(state2);
131
132
133 Assert.assertEquals(0.0, stm2.subtract(stm1).getNorm1(), 5.0e-12 * stm1.getNorm1());
134 Assert.assertEquals(0.0, jacobian2.subtract(jacobian1).getNorm1(), 7.0e-10 * jacobian1.getNorm1());
135
136 }
137
138 @Test
139 public void testPropagationTypesElliptical() throws FileNotFoundException, UnsupportedEncodingException, OrekitException {
140 doTestPropagation(PropagationType.MEAN, 7.0e-16);
141 }
142
143 @Test
144 public void testPropagationTypesEllipticalWithShortPeriod() throws FileNotFoundException, UnsupportedEncodingException, OrekitException {
145 doTestPropagation(PropagationType.OSCULATING, 3.3e-4);
146 }
147
148 private void doTestPropagation(PropagationType type, double tolerance)
149 throws FileNotFoundException, UnsupportedEncodingException {
150
151 UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(5, 5);
152
153 double dt = 900;
154 double dP = 0.001;
155 final OrbitType orbitType = OrbitType.EQUINOCTIAL;
156
157
158 DSSTPropagator propagator = setUpPropagator(type, dP, provider);
159 propagator.setMu(provider.getMu());
160 final SpacecraftState initialState = propagator.getInitialState();
161 final double[] stateVector = new double[6];
162 OrbitType.EQUINOCTIAL.mapOrbitToArray(initialState.getOrbit(), PositionAngle.MEAN, stateVector, null);
163 final AbsoluteDate target = propagator.getInitialState().getDate().shiftedBy(dt);
164 PickUpHandler pickUp = new PickUpHandler(propagator, null, null, null);
165 propagator.getMultiplexer().add(pickUp);
166 propagator.propagate(target);
167
168
169 double[][] dYdY0Ref = new double[6][6];
170 DSSTPropagator propagator2 = setUpPropagator(type, dP, provider);
171 propagator2.setMu(provider.getMu());
172 double[] steps = NumericalPropagator.tolerances(1000000 * dP, initialState.getOrbit(), orbitType)[0];
173 for (int i = 0; i < 6; ++i) {
174 propagator2.setInitialState(shiftState(initialState, orbitType, -4 * steps[i], i), type);
175 SpacecraftState sM4h = propagator2.propagate(target);
176 propagator2.setInitialState(shiftState(initialState, orbitType, -3 * steps[i], i), type);
177 SpacecraftState sM3h = propagator2.propagate(target);
178 propagator2.setInitialState(shiftState(initialState, orbitType, -2 * steps[i], i), type);
179 SpacecraftState sM2h = propagator2.propagate(target);
180 propagator2.setInitialState(shiftState(initialState, orbitType, -1 * steps[i], i), type);
181 SpacecraftState sM1h = propagator2.propagate(target);
182 propagator2.setInitialState(shiftState(initialState, orbitType, 1 * steps[i], i), type);
183 SpacecraftState sP1h = propagator2.propagate(target);
184 propagator2.setInitialState(shiftState(initialState, orbitType, 2 * steps[i], i), type);
185 SpacecraftState sP2h = propagator2.propagate(target);
186 propagator2.setInitialState(shiftState(initialState, orbitType, 3 * steps[i], i), type);
187 SpacecraftState sP3h = propagator2.propagate(target);
188 propagator2.setInitialState(shiftState(initialState, orbitType, 4 * steps[i], i), type);
189 SpacecraftState sP4h = propagator2.propagate(target);
190 fillJacobianColumn(dYdY0Ref, i, orbitType, steps[i],
191 sM4h, sM3h, sM2h, sM1h, sP1h, sP2h, sP3h, sP4h);
192 }
193
194 for (int i = 0; i < 6; ++i) {
195 for (int j = 0; j < 6; ++j) {
196 if (stateVector[i] != 0) {
197 double error = FastMath.abs((pickUp.getStm().getEntry(i, j) - dYdY0Ref[i][j]) / stateVector[i]) * steps[j];
198 Assert.assertEquals(0, error, tolerance);
199 }
200 }
201 }
202 }
203
204 private void fillJacobianColumn(double[][] jacobian, int column,
205 OrbitType orbitType, double h,
206 SpacecraftState sM4h, SpacecraftState sM3h,
207 SpacecraftState sM2h, SpacecraftState sM1h,
208 SpacecraftState sP1h, SpacecraftState sP2h,
209 SpacecraftState sP3h, SpacecraftState sP4h) {
210 double[] aM4h = stateToArray(sM4h, orbitType)[0];
211 double[] aM3h = stateToArray(sM3h, orbitType)[0];
212 double[] aM2h = stateToArray(sM2h, orbitType)[0];
213 double[] aM1h = stateToArray(sM1h, orbitType)[0];
214 double[] aP1h = stateToArray(sP1h, orbitType)[0];
215 double[] aP2h = stateToArray(sP2h, orbitType)[0];
216 double[] aP3h = stateToArray(sP3h, orbitType)[0];
217 double[] aP4h = stateToArray(sP4h, orbitType)[0];
218 for (int i = 0; i < jacobian.length; ++i) {
219 jacobian[i][column] = ( -3 * (aP4h[i] - aM4h[i]) +
220 32 * (aP3h[i] - aM3h[i]) -
221 168 * (aP2h[i] - aM2h[i]) +
222 672 * (aP1h[i] - aM1h[i])) / (840 * h);
223 }
224 }
225
226 private SpacecraftState shiftState(SpacecraftState state, OrbitType orbitType,
227 double delta, int column) {
228
229 double[][] array = stateToArray(state, orbitType);
230 array[0][column] += delta;
231
232 return arrayToState(array, orbitType, state.getFrame(), state.getDate(),
233 state.getMu(), state.getAttitude());
234
235 }
236
237 private double[][] stateToArray(SpacecraftState state, OrbitType orbitType) {
238 double[][] array = new double[2][6];
239
240 orbitType.mapOrbitToArray(state.getOrbit(), PositionAngle.MEAN, array[0], array[1]);
241 return array;
242 }
243
244 private SpacecraftState arrayToState(double[][] array, OrbitType orbitType,
245 Frame frame, AbsoluteDate date, double mu,
246 Attitude attitude) {
247 EquinoctialOrbit orbit = (EquinoctialOrbit) orbitType.mapArrayToOrbit(array[0], array[1], PositionAngle.MEAN, date, mu, frame);
248 return new SpacecraftState(orbit, attitude);
249 }
250
251 private DSSTPropagator setUpPropagator(PropagationType type, double dP, UnnormalizedSphericalHarmonicsProvider provider) {
252
253 final double minStep = 6000.0;
254 final double maxStep = 86400.0;
255
256 Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
257
258 DSSTForceModel tesseral = new DSSTTesseral(earthFrame,
259 Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
260 4, 4, 4, 8, 4, 4, 2);
261
262 DSSTForceModel zonal = new DSSTZonal(provider, 4, 3, 9);
263 DSSTForceModel srp = new DSSTSolarRadiationPressure(1.2, 100., CelestialBodyFactory.getSun(),
264 Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
265 provider.getMu());
266
267 DSSTForceModel moon = new DSSTThirdBody(CelestialBodyFactory.getMoon(), provider.getMu());
268
269 Orbit initialOrbit =
270 new KeplerianOrbit(8000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.MEAN,
271 FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
272 provider.getMu());
273 final EquinoctialOrbit orbit = (EquinoctialOrbit) OrbitType.EQUINOCTIAL.convertType(initialOrbit);
274
275
276 double[][] tol = NumericalPropagator.tolerances(dP, orbit, OrbitType.EQUINOCTIAL);
277 DSSTPropagator propagator =
278 new DSSTPropagator(new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]), type);
279 propagator.addForceModel(srp);
280 propagator.addForceModel(tesseral);
281 propagator.addForceModel(zonal);
282 propagator.addForceModel(moon);
283
284 propagator.setInitialState(new SpacecraftState(orbit), type);
285 return propagator;
286 }
287
288 private void initializeShortPeriod(final MatricesHarvester harvester, final DSSTPropagator propagator) {
289
290 final SpacecraftState initial = propagator.initialIsOsculating() ?
291 DSSTPropagator.computeMeanState(propagator.getInitialState(), propagator.getAttitudeProvider(), propagator.getAllForceModels()) :
292 propagator.getInitialState();
293 ((DSSTHarvester) harvester).initializeFieldShortPeriodTerms(initial);
294 }
295
296
297
298
299
300 @Test
301 public void testIssue713() {
302 UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(5, 5);
303
304 double dP = 0.001;
305
306
307 DSSTPropagator propagatorMEAN = setUpPropagator(PropagationType.MEAN, dP, provider);
308 propagatorMEAN.setMu(provider.getMu());
309 SpacecraftState initialStateMEAN = propagatorMEAN.getInitialState();
310 DSSTHarvester harvesterMEAN = (DSSTHarvester) propagatorMEAN.setupMatricesComputation("stm", null, null);
311 propagatorMEAN.
312 getAllForceModels().
313 forEach(fm -> fm.
314 getParametersDrivers().
315 stream().
316 filter(d -> d.getName().equals(RadiationSensitive.REFLECTION_COEFFICIENT)).
317 forEach(d -> d.setSelected(true)));
318 harvesterMEAN.setReferenceState(initialStateMEAN);
319 SpacecraftState finalMEAN = propagatorMEAN.propagate(initialStateMEAN.getDate());
320 RealMatrix dYdY0MEAN = harvesterMEAN.getStateTransitionMatrix(finalMEAN);
321 for (int i = 0; i < 6; ++i) {
322 for (int j = 0; j < 6; ++j) {
323 Assert.assertEquals(i == j ? 1.0 : 0.0, dYdY0MEAN.getEntry(i, j), 1e-9);
324 }
325 }
326 RealMatrix dYdPMEAN = harvesterMEAN.getParametersJacobian(finalMEAN);
327 Assert.assertEquals(6, dYdPMEAN.getRowDimension());
328 Assert.assertEquals(1, dYdPMEAN.getColumnDimension());
329 for (int i = 0; i < 6; ++i) {
330 Assert.assertEquals(0.0, dYdPMEAN.getEntry(i, 0), 1e-9);
331 }
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365 }
366
367 }