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