1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.propagation;
18
19
20 import java.util.Arrays;
21 import java.util.List;
22 import java.util.concurrent.atomic.AtomicInteger;
23
24 import org.hipparchus.exception.LocalizedCoreFormats;
25 import org.hipparchus.geometry.euclidean.threed.Vector3D;
26 import org.hipparchus.ode.ODEIntegrator;
27 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
28 import org.hipparchus.util.FastMath;
29 import org.junit.After;
30 import org.junit.Assert;
31 import org.junit.Before;
32 import org.junit.Test;
33 import org.orekit.Utils;
34 import org.orekit.attitudes.AttitudeProvider;
35 import org.orekit.attitudes.BodyCenterPointing;
36 import org.orekit.bodies.OneAxisEllipsoid;
37 import org.orekit.errors.OrekitException;
38 import org.orekit.errors.OrekitMessages;
39 import org.orekit.forces.ForceModel;
40 import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
41 import org.orekit.forces.gravity.potential.GravityFieldFactory;
42 import org.orekit.forces.gravity.potential.ICGEMFormatReader;
43 import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
44 import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
45 import org.orekit.frames.FramesFactory;
46 import org.orekit.orbits.KeplerianOrbit;
47 import org.orekit.orbits.Orbit;
48 import org.orekit.orbits.OrbitType;
49 import org.orekit.orbits.PositionAngle;
50 import org.orekit.propagation.analytical.EcksteinHechlerPropagator;
51 import org.orekit.propagation.events.DateDetector;
52 import org.orekit.propagation.events.handlers.StopOnEvent;
53 import org.orekit.propagation.integration.AdditionalDerivativesProvider;
54 import org.orekit.propagation.numerical.NumericalPropagator;
55 import org.orekit.propagation.semianalytical.dsst.DSSTPropagator;
56 import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
57 import org.orekit.propagation.semianalytical.dsst.forces.DSSTZonal;
58 import org.orekit.time.AbsoluteDate;
59 import org.orekit.time.DateComponents;
60 import org.orekit.time.TimeComponents;
61 import org.orekit.time.TimeScalesFactory;
62 import org.orekit.utils.Constants;
63 import org.orekit.utils.IERSConventions;
64
65
66 public class PropagatorsParallelizerTest {
67
68 @Test
69 public void testIssue717() {
70
71
72 Utils.setDataRoot("regular-data:potential/icgem-format");
73 GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("^eigen-6s-truncated$", false));
74 UnnormalizedSphericalHarmonicsProvider gravity = GravityFieldFactory.getUnnormalizedProvider(8, 8);
75
76
77 Orbit orbit = new KeplerianOrbit(15000000.0, 0.125, 1.25,
78 0.250, 1.375, 0.0625, PositionAngle.MEAN,
79 FramesFactory.getEME2000(),
80 new AbsoluteDate(2000, 2, 24, 11, 35, 47.0, TimeScalesFactory.getUTC()),
81 gravity.getMu());
82
83
84
85 final double[][] tol = DSSTPropagator.tolerances(0.01, orbit);
86 final DSSTPropagator propagator = new DSSTPropagator(new DormandPrince853Integrator(0.01, 600.0, tol[0], tol[1]), PropagationType.OSCULATING);
87
88
89 final DSSTForceModel zonal = new DSSTZonal(gravity, 4, 3, 9);
90 propagator.addForceModel(zonal);
91
92 propagator.setInitialState(new SpacecraftState(orbit));
93
94
95 final double deltaT = 30.0;
96
97 final PropagatorsParallelizer parallelizer =
98 new PropagatorsParallelizer(Arrays.asList(propagator), interpolators -> {interpolators.get(0).getCurrentState().getDate();});
99
100 final SpacecraftState state = parallelizer.propagate(orbit.getDate().shiftedBy(deltaT).shiftedBy(+1.0), orbit.getDate().shiftedBy(-2.0 * deltaT).shiftedBy(-1.0)).get(0);
101
102
103 Assert.assertNotNull(state);
104
105 }
106
107 @Test
108 public void testNumericalNotInitialized() {
109
110 final AbsoluteDate startDate = orbit.getDate();
111 final AbsoluteDate endDate = startDate.shiftedBy(3600.0);
112 List<Propagator> propagators = Arrays.asList(buildEcksteinHechler(),
113 buildNotInitializedNumerical());
114
115 PropagatorsParallelizer parallelizer =
116 new PropagatorsParallelizer(propagators, interpolators -> Assert.fail("should not be called"));
117 try {
118 parallelizer.propagate(startDate, endDate);
119 Assert.fail("an exception should have been thrown");
120 } catch (OrekitException oe) {
121 Assert.assertEquals(OrekitMessages.INITIAL_STATE_NOT_SPECIFIED_FOR_ORBIT_PROPAGATION,
122 oe.getSpecifier());
123 }
124
125 }
126
127 @Test
128 public void testAnalyticalAndNumericalSameOrbit() {
129
130 final AbsoluteDate startDate = orbit.getDate();
131 final AbsoluteDate endDate = startDate.shiftedBy(3600.0);
132 List<Propagator> propagators = Arrays.asList(buildEcksteinHechler(),
133 buildNumerical());
134
135 PropagatorsParallelizer parallelizer =
136 new PropagatorsParallelizer(propagators,
137 interpolators -> {
138 Assert.assertEquals(2, interpolators.size());
139 AbsoluteDate aPrev = interpolators.get(0).getPreviousState().getDate();
140 AbsoluteDate aCurr = interpolators.get(0).getCurrentState().getDate();
141 AbsoluteDate nPrev = interpolators.get(1).getPreviousState().getDate();
142 AbsoluteDate nCurr = interpolators.get(1).getCurrentState().getDate();
143 Assert.assertEquals(0.0, aPrev.durationFrom(nPrev), 3.0e-13);
144 Assert.assertEquals(0.0, aCurr.durationFrom(nCurr), 3.0e-13);
145 Vector3D aPos = interpolators.get(0).getCurrentState().getPVCoordinates().getPosition();
146 Vector3D nPos = interpolators.get(1).getCurrentState().getPVCoordinates().getPosition();
147 Assert.assertTrue(Vector3D.distance(aPos, nPos) < 111.0);
148 });
149 List<SpacecraftState> results = parallelizer.propagate(startDate, endDate);
150
151 Assert.assertEquals(2, results.size());
152 for (final SpacecraftState state : results) {
153 Assert.assertEquals(0.0, state.getDate().durationFrom(endDate), 1.0e-15);
154 }
155
156 }
157
158 @Test
159 public void testVsAnalyticalMonoSat() {
160
161 final AbsoluteDate startDate = orbit.getDate();
162 final AbsoluteDate endDate = startDate.shiftedBy(3600.0);
163 Propagator mono = buildEcksteinHechler();
164 final EphemerisGenerator generator = mono.getEphemerisGenerator();
165 mono.propagate(startDate, endDate);
166 final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
167
168 List<Propagator> propagators = Arrays.asList(buildEcksteinHechler(),
169 buildNumerical());
170
171 PropagatorsParallelizer parallelizer =
172 new PropagatorsParallelizer(propagators,
173 interpolators -> {
174 AbsoluteDate aCurr = interpolators.get(0).getCurrentState().getDate();
175 Vector3D aPos = interpolators.get(0).getCurrentState().getPVCoordinates().getPosition();
176 Vector3D ePos = ephemeris.getPVCoordinates(aCurr, orbit.getFrame()).getPosition();
177 Assert.assertEquals(0, Vector3D.distance(ePos, aPos), 1.0e-15);
178 });
179 List<SpacecraftState> results = parallelizer.propagate(startDate, endDate);
180
181 Assert.assertEquals(2, parallelizer.getPropagators().size());
182 Assert.assertEquals(2, results.size());
183 for (final SpacecraftState state : results) {
184 Assert.assertEquals(0.0, state.getDate().durationFrom(endDate), 1.0e-15);
185 }
186
187 }
188
189 @Test
190 public void testVsNumericalMonoSat() {
191
192 final AbsoluteDate startDate = orbit.getDate();
193 final AbsoluteDate endDate = startDate.shiftedBy(3600.0);
194 Propagator mono = buildNumerical();
195 final EphemerisGenerator generator = mono.getEphemerisGenerator();
196 mono.propagate(startDate, endDate);
197 final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
198
199 List<Propagator> propagators = Arrays.asList(buildEcksteinHechler(),
200 buildNumerical());
201
202 PropagatorsParallelizer parallelizer =
203 new PropagatorsParallelizer(propagators,
204 interpolators -> {
205 AbsoluteDate nCurr = interpolators.get(1).getCurrentState().getDate();
206 Vector3D nPos = interpolators.get(1).getCurrentState().getPVCoordinates().getPosition();
207 Vector3D ePos = ephemeris.getPVCoordinates(nCurr, orbit.getFrame()).getPosition();
208 Assert.assertEquals(0, Vector3D.distance(ePos, nPos), 1.0e-15);
209 });
210 List<SpacecraftState> results = parallelizer.propagate(startDate, endDate);
211
212 Assert.assertEquals(2, results.size());
213 for (final SpacecraftState state : results) {
214 Assert.assertEquals(0.0, state.getDate().durationFrom(endDate), 1.0e-15);
215 }
216
217 }
218
219 @Test
220 public void testOrekitException() {
221 final AbsoluteDate startDate = orbit.getDate();
222 final AbsoluteDate endDate = startDate.shiftedBy(3600.0);
223 List<Propagator> propagators = Arrays.asList(buildEcksteinHechler(),
224 buildNumerical());
225 propagators.get(0).addEventDetector(new DateDetector(startDate.shiftedBy(900.0)).
226 withHandler((state, detector, increasing) -> {
227 throw new OrekitException(LocalizedCoreFormats.SIMPLE_MESSAGE,
228 "inTest");
229 }));
230 try {
231 new PropagatorsParallelizer(propagators, interpolators -> {}).propagate(startDate, endDate);
232 Assert.fail("an exception should have been thrown");
233 } catch (OrekitException oe) {
234 Assert.assertNull(oe.getCause());
235 Assert.assertEquals(LocalizedCoreFormats.SIMPLE_MESSAGE, oe.getSpecifier());
236 Assert.assertEquals("inTest", (String) oe.getParts()[0]);
237 }
238 }
239
240 @Test
241 public void testEarlyStop() {
242 final AbsoluteDate startDate = orbit.getDate();
243 final AbsoluteDate endDate = startDate.shiftedBy(3600.0);
244 List<Propagator> propagators = Arrays.asList(buildEcksteinHechler(),
245 buildNumerical());
246 propagators.get(0).addEventDetector(new DateDetector(startDate.shiftedBy(900.0)).
247 withHandler((state, detector, increasing) -> {
248 throw new RuntimeException("boo!");
249 }));
250 try {
251 new PropagatorsParallelizer(propagators, interpolators -> {}).propagate(startDate, endDate);
252 Assert.fail("an exception should have been thrown");
253 } catch (OrekitException oe) {
254 Assert.assertNotNull(oe.getCause());
255 Assert.assertTrue(oe.getCause() instanceof RuntimeException);
256 Assert.assertEquals(LocalizedCoreFormats.SIMPLE_MESSAGE, oe.getSpecifier());
257 Assert.assertTrue(((String) oe.getParts()[0]).endsWith("boo!"));
258 }
259 }
260
261 @Test
262 public void testStopOnEarlyEvent() {
263 final AbsoluteDate startDate = orbit.getDate();
264 final AbsoluteDate endDate = startDate.shiftedBy(3600.0);
265 final AbsoluteDate stopDate = startDate.shiftedBy(0.01);
266 List<Propagator> propagators = Arrays.asList(buildEcksteinHechler(),
267 buildNumerical());
268 propagators.get(0).addEventDetector(new DateDetector(stopDate).withHandler(new StopOnEvent<>()));
269 List<SpacecraftState> results = new PropagatorsParallelizer(propagators, interpolators -> {}).
270 propagate(startDate, endDate);
271 Assert.assertEquals(2, results.size());
272 Assert.assertEquals(0.0, results.get(0).getDate().durationFrom(stopDate), 1.0e-15);
273 Assert.assertEquals(0.0, results.get(1).getDate().durationFrom(stopDate), 1.0e-15);
274 }
275
276 @Test
277 public void testStopOnLateEvent() {
278 final AbsoluteDate startDate = orbit.getDate();
279 final AbsoluteDate endDate = startDate.shiftedBy(3600.0);
280 final AbsoluteDate stopDate = startDate.shiftedBy(900.0);
281 List<Propagator> propagators = Arrays.asList(buildEcksteinHechler(),
282 buildNumerical());
283 propagators.get(0).addEventDetector(new DateDetector(stopDate).withHandler(new StopOnEvent<>()));
284 List<SpacecraftState> results = new PropagatorsParallelizer(propagators, interpolators -> {}).
285 propagate(startDate, endDate);
286 Assert.assertEquals(2, results.size());
287 Assert.assertEquals(0.0, results.get(0).getDate().durationFrom(stopDate), 1.0e-15);
288 Assert.assertEquals(0.0, results.get(1).getDate().durationFrom(stopDate), 1.0e-15);
289 }
290
291 @Test
292 public void testInternalStepHandler() {
293 final AbsoluteDate startDate = orbit.getDate();
294 final AbsoluteDate endDate = startDate.shiftedBy(3600.0);
295 List<Propagator> propagators = Arrays.asList(buildEcksteinHechler(),
296 buildNumerical());
297 final AtomicInteger called0 = new AtomicInteger();
298 propagators.get(0).getMultiplexer().add(interpolator -> called0.set(11));
299 final AtomicInteger called1 = new AtomicInteger();
300 propagators.get(1).getMultiplexer().add(interpolator -> called1.set(22));
301 List<SpacecraftState> results = new PropagatorsParallelizer(propagators, interpolators -> {}).
302 propagate(startDate, endDate);
303 Assert.assertEquals(2, results.size());
304 Assert.assertEquals(0.0, results.get(0).getDate().durationFrom(endDate), 1.0e-15);
305 Assert.assertEquals(0.0, results.get(1).getDate().durationFrom(endDate), 1.0e-15);
306 Assert.assertEquals(11, called0.get());
307 Assert.assertEquals(22, called1.get());
308 }
309
310 @Test
311 public void testIntegrableGenerator() {
312 final AbsoluteDate startDate = orbit.getDate();
313 final AbsoluteDate endDate = startDate.shiftedBy(3600.0);
314 final NumericalPropagator p0 = buildNumerical();
315 final NumericalPropagator p1 = buildNumerical();
316 final String name = "generator";
317 final double base0 = 2.0e-3;
318 final double base1 = 2.5e-3;
319 p0.addAdditionalDerivativesProvider(new Exponential(name, base0));
320 p0.setInitialState(p0.getInitialState().addAdditionalState(name, 1.0));
321 p1.addAdditionalDerivativesProvider(new Exponential(name, base1));
322 p1.setInitialState(p1.getInitialState().addAdditionalState(name, 1.0));
323 List<SpacecraftState> results = new PropagatorsParallelizer(Arrays.asList(p0, p1), interpolators -> {}).
324 propagate(startDate, endDate);
325 double expected0 = FastMath.exp(base0 * endDate.durationFrom(startDate));
326 double expected1 = FastMath.exp(base1 * endDate.durationFrom(startDate));
327 Assert.assertEquals(expected0, results.get(0).getAdditionalState(name)[0], 6.0e-9 * expected0);
328 Assert.assertEquals(expected1, results.get(1).getAdditionalState(name)[0], 5.0e-8 * expected1);
329 }
330
331 private static class Exponential implements AdditionalDerivativesProvider {
332 final String name;
333 final double base;
334 Exponential(final String name, final double base) {
335 this.name = name;
336 this.base = base;
337 }
338 public String getName() {
339 return name;
340 }
341 public boolean yield(final SpacecraftState state) {
342 return !state.hasAdditionalState(name);
343 }
344 public double[] derivatives(final SpacecraftState state) {
345 return new double[] { base * state.getAdditionalState(name)[0] };
346 }
347 public int getDimension() {
348 return 1;
349 }
350 }
351
352 private EcksteinHechlerPropagator buildEcksteinHechler() {
353 return new EcksteinHechlerPropagator(orbit, attitudeLaw, mass, unnormalizedGravityField);
354 }
355
356 private NumericalPropagator buildNumerical() {
357 NumericalPropagator numericalPropagator = buildNotInitializedNumerical();
358 numericalPropagator.setInitialState(new SpacecraftState(orbit,
359 attitudeLaw.getAttitude(orbit,
360 orbit.getDate(),
361 orbit.getFrame()),
362 mass));
363 return numericalPropagator;
364 }
365
366 private NumericalPropagator buildNotInitializedNumerical() {
367 OrbitType type = OrbitType.CARTESIAN;
368 double minStep = 0.001;
369 double maxStep = 300;
370 double[][] tolerances = NumericalPropagator.tolerances(10.0, orbit, type);
371 ODEIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tolerances[0], tolerances[1]);
372 NumericalPropagator numericalPropagator = new NumericalPropagator(integrator);
373 ForceModel gravity = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true),
374 normalizedGravityField);
375 numericalPropagator.addForceModel(gravity);
376 return numericalPropagator;
377 }
378
379 @Before
380 public void setUp() {
381 try {
382 Utils.setDataRoot("regular-data:potential/icgem-format");
383 unnormalizedGravityField = GravityFieldFactory.getUnnormalizedProvider(6, 0);
384 normalizedGravityField = GravityFieldFactory.getNormalizedProvider(6, 0);
385
386 mass = 2500;
387 double a = 7187990.1979844316;
388 double e = 0.5e-4;
389 double i = 1.7105407051081795;
390 double omega = 1.9674147913622104;
391 double OMEGA = FastMath.toRadians(261);
392 double lv = 0;
393
394 AbsoluteDate date = new AbsoluteDate(new DateComponents(2004, 01, 01),
395 TimeComponents.H00,
396 TimeScalesFactory.getUTC());
397 orbit = new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngle.TRUE,
398 FramesFactory.getEME2000(), date, normalizedGravityField.getMu());
399 OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
400 Constants.WGS84_EARTH_FLATTENING,
401 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
402 attitudeLaw = new BodyCenterPointing(orbit.getFrame(), earth);
403
404 } catch (OrekitException oe) {
405 Assert.fail(oe.getLocalizedMessage());
406 }
407 }
408
409 @After
410 public void tearDown() {
411 mass = Double.NaN;
412 orbit = null;
413 attitudeLaw = null;
414 unnormalizedGravityField = null;
415 normalizedGravityField = null;
416 }
417
418 private double mass;
419 private Orbit orbit;
420 private AttitudeProvider attitudeLaw;
421 private UnnormalizedSphericalHarmonicsProvider unnormalizedGravityField;
422 private NormalizedSphericalHarmonicsProvider normalizedGravityField;
423
424 }