1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.propagation.numerical;
18
19 import java.io.FileNotFoundException;
20 import java.io.IOException;
21 import java.text.ParseException;
22 import java.util.ArrayList;
23 import java.util.Collections;
24 import java.util.List;
25 import java.util.concurrent.ExecutionException;
26 import java.util.concurrent.atomic.AtomicInteger;
27 import java.util.function.Consumer;
28 import java.util.stream.IntStream;
29 import java.util.stream.Stream;
30
31 import org.hamcrest.CoreMatchers;
32 import org.hamcrest.MatcherAssert;
33 import org.hipparchus.CalculusFieldElement;
34 import org.hipparchus.Field;
35 import org.hipparchus.exception.LocalizedCoreFormats;
36 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
37 import org.hipparchus.geometry.euclidean.threed.Rotation;
38 import org.hipparchus.geometry.euclidean.threed.Vector3D;
39 import org.hipparchus.ode.ODEIntegrator;
40 import org.hipparchus.ode.events.Action;
41 import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
42 import org.hipparchus.ode.nonstiff.ClassicalRungeKuttaIntegrator;
43 import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
44 import org.hipparchus.ode.nonstiff.LutherIntegrator;
45 import org.hipparchus.util.FastMath;
46 import org.junit.jupiter.api.AfterEach;
47 import org.junit.jupiter.api.Assertions;
48 import org.junit.jupiter.api.BeforeEach;
49 import org.junit.jupiter.api.Test;
50 import org.mockito.Mockito;
51 import org.orekit.OrekitMatchers;
52 import org.orekit.Utils;
53 import org.orekit.attitudes.*;
54 import org.orekit.bodies.CelestialBodyFactory;
55 import org.orekit.bodies.OneAxisEllipsoid;
56 import org.orekit.errors.OrekitException;
57 import org.orekit.errors.OrekitMessages;
58 import org.orekit.forces.ForceModel;
59 import org.orekit.forces.drag.DragForce;
60 import org.orekit.forces.drag.IsotropicDrag;
61 import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
62 import org.orekit.forces.gravity.ThirdBodyAttraction;
63 import org.orekit.forces.gravity.potential.GRGSFormatReader;
64 import org.orekit.forces.gravity.potential.GravityFieldFactory;
65 import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
66 import org.orekit.forces.gravity.potential.SHMFormatReader;
67 import org.orekit.forces.maneuvers.ImpulseManeuver;
68 import org.orekit.forces.radiation.IsotropicRadiationSingleCoefficient;
69 import org.orekit.forces.radiation.SolarRadiationPressure;
70 import org.orekit.frames.Frame;
71 import org.orekit.frames.FramesFactory;
72 import org.orekit.frames.LOFType;
73 import org.orekit.models.earth.atmosphere.DTM2000;
74 import org.orekit.models.earth.atmosphere.data.MarshallSolarActivityFutureEstimation;
75 import org.orekit.orbits.CartesianOrbit;
76 import org.orekit.orbits.EquinoctialOrbit;
77 import org.orekit.orbits.KeplerianOrbit;
78 import org.orekit.orbits.Orbit;
79 import org.orekit.orbits.OrbitType;
80 import org.orekit.orbits.PositionAngleType;
81 import org.orekit.propagation.AdditionalStateProvider;
82 import org.orekit.propagation.BoundedPropagator;
83 import org.orekit.propagation.EphemerisGenerator;
84 import org.orekit.propagation.FieldSpacecraftState;
85 import org.orekit.propagation.PropagationType;
86 import org.orekit.propagation.SpacecraftState;
87 import org.orekit.propagation.conversion.DormandPrince853IntegratorBuilder;
88 import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
89 import org.orekit.propagation.events.AbstractDetector;
90 import org.orekit.propagation.events.AdaptableInterval;
91 import org.orekit.propagation.events.ApsideDetector;
92 import org.orekit.propagation.events.DateDetector;
93 import org.orekit.propagation.events.EventDetector;
94 import org.orekit.propagation.events.FieldEventDetector;
95 import org.orekit.propagation.events.handlers.ContinueOnEvent;
96 import org.orekit.propagation.events.handlers.EventHandler;
97 import org.orekit.propagation.events.handlers.RecordAndContinue;
98 import org.orekit.propagation.events.handlers.StopOnEvent;
99 import org.orekit.propagation.integration.AbstractIntegratedPropagator;
100 import org.orekit.propagation.integration.AdditionalDerivativesProvider;
101 import org.orekit.propagation.integration.CombinedDerivatives;
102 import org.orekit.propagation.sampling.OrekitFixedStepHandler;
103 import org.orekit.propagation.sampling.OrekitStepHandler;
104 import org.orekit.propagation.sampling.OrekitStepInterpolator;
105 import org.orekit.time.AbsoluteDate;
106 import org.orekit.time.DateComponents;
107 import org.orekit.time.FieldAbsoluteDate;
108 import org.orekit.time.TimeComponents;
109 import org.orekit.time.TimeScale;
110 import org.orekit.time.TimeScalesFactory;
111 import org.orekit.utils.AngularCoordinates;
112 import org.orekit.utils.Constants;
113 import org.orekit.utils.FieldPVCoordinatesProvider;
114 import org.orekit.utils.IERSConventions;
115 import org.orekit.utils.PVCoordinates;
116 import org.orekit.utils.PVCoordinatesProvider;
117 import org.orekit.utils.ParameterDriver;
118 import org.orekit.utils.TimeStampedPVCoordinates;
119
120 class NumericalPropagatorTest {
121
122 private double mu;
123 private AbsoluteDate initDate;
124 private SpacecraftState initialState;
125 private NumericalPropagator propagator;
126
127 @Test
128 void testIssue1032() {
129 Assertions.assertEquals(PropagationType.OSCULATING, propagator.getPropagationType());
130 }
131
132 @Test
133 void testEventsWithTimeRangePropagation() {
134 final AtomicInteger counter = new AtomicInteger(0);
135 final double dt = 60.0;
136 final EventDetector singleDetector = new DateDetector(initDate.shiftedBy(dt / 2)).
137 withHandler((state, detector, increasing) -> {
138 counter.incrementAndGet();
139 return Action.CONTINUE;
140 });
141 ForceModel force = new ForceModelAdapter() {
142 @Override
143 public Stream<EventDetector> getEventDetectors() {
144 return Stream.of(singleDetector);
145 }
146 };
147
148
149 propagator.setOrbitType(OrbitType.CARTESIAN);
150 propagator.addForceModel(force);
151 AbsoluteDate target = initDate.shiftedBy(dt);
152 propagator.propagate(initDate.shiftedBy(1.0), target);
153
154 Assertions.assertEquals(1, counter.intValue());
155
156 }
157
158 @Test
159 void testForceModelInitialized() {
160
161
162 SpacecraftState[] actualState = new SpacecraftState[1];
163 AbsoluteDate[] actualDate = new AbsoluteDate[1];
164 ForceModel force = new ForceModelAdapter() {
165 @Override
166 public void init(SpacecraftState initialState, AbsoluteDate target) {
167 actualState[0] = initialState;
168 actualDate[0] = target;
169 }
170 };
171
172
173 propagator.setOrbitType(OrbitType.CARTESIAN);
174 propagator.addForceModel(force);
175 AbsoluteDate target = initDate.shiftedBy(60);
176 propagator.propagate(target);
177
178
179 MatcherAssert.assertThat(actualDate[0], CoreMatchers.is(target));
180 MatcherAssert.assertThat(actualState[0].getDate().durationFrom(initDate),
181 CoreMatchers.is(0.0));
182 MatcherAssert.assertThat(actualState[0].getPVCoordinates(),
183 OrekitMatchers.pvIs(initialState.getPVCoordinates()));
184 }
185
186 @Test
187 void testEphemerisModeWithHandler() {
188
189 AbsoluteDate end = initDate.shiftedBy(90 * 60);
190
191
192 final List<SpacecraftState> states = new ArrayList<>();
193 EphemerisGenerator generator = propagator.getEphemerisGenerator();
194 propagator.setStepHandler(interpolator -> states.add(interpolator.getCurrentState()));
195 propagator.propagate(end);
196 final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
197
198
199 Assertions.assertTrue(states.size() > 10);
200 for (SpacecraftState state : states) {
201 PVCoordinates actual =
202 ephemeris.propagate(state.getDate()).getPVCoordinates();
203 MatcherAssert.assertThat(actual, OrekitMatchers.pvIs(state.getPVCoordinates()));
204 }
205 }
206
207
208 @Test
209 void testEventAtEndOfEphemeris() {
210
211
212 AbsoluteDate end = initDate.shiftedBy(100)
213 .shiftedBy(3 * FastMath.ulp(100.0) / 4);
214 final EphemerisGenerator generator = propagator.getEphemerisGenerator();
215 propagator.propagate(end);
216 BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
217 CountingHandler handler = new CountingHandler();
218 DateDetector detector = new DateDetector(end).
219 withMaxCheck(10).
220 withThreshold(1e-9).
221 withMinGap(10).
222 withHandler(handler);
223
224 ephemeris.addEventDetector(detector);
225
226
227
228 SpacecraftState actual = ephemeris.propagate(end);
229
230
231 Assertions.assertEquals(actual.getDate().durationFrom(end), 0.0, 0.0);
232 Assertions.assertEquals(1, handler.eventCount);
233 }
234
235
236 @Test
237 void testEventAtBeginningOfEphemeris() {
238
239
240 AbsoluteDate end = initDate.shiftedBy(100)
241 .shiftedBy(3 * FastMath.ulp(100.0) / 4);
242 final EphemerisGenerator generator = propagator.getEphemerisGenerator();
243 propagator.propagate(end);
244 BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
245 CountingHandler handler = new CountingHandler();
246
247
248 AbsoluteDate eventDate = initDate.shiftedBy(FastMath.ulp(100.0) / 10.0);
249 DateDetector detector = new DateDetector(eventDate).
250 withMaxCheck(10).
251 withThreshold(1e-9).
252 withMinGap(10).
253 withHandler(handler);
254
255 ephemeris.addEventDetector(detector);
256
257
258
259 Assertions.assertEquals(ephemeris.propagate(end).getDate().durationFrom(end), 0.0, 0.0);
260
261 Assertions.assertEquals(ephemeris.propagate(initDate).getDate().durationFrom(initDate), 0.0, 0.0);
262 Assertions.assertEquals(2, handler.eventCount);
263 }
264
265
266 private static class CountingHandler implements EventHandler {
267
268
269
270
271
272 private int eventCount = 0;
273
274 @Override
275 public Action eventOccurred(SpacecraftState s,
276 EventDetector detector,
277 boolean increasing) {
278 eventCount++;
279 return Action.CONTINUE;
280 }
281
282 @Override
283 public SpacecraftState resetState(EventDetector detector,
284 SpacecraftState oldState) {
285 return null;
286 }
287 }
288
289
290
291
292
293 @Test
294 void testCloseEventDates() {
295
296 DateDetector d1 = new DateDetector(initDate.shiftedBy(15)).
297 withMaxCheck(10).
298 withThreshold(1).
299 withHandler(new ContinueOnEvent());
300 DateDetector d2 = new DateDetector(initDate.shiftedBy(15.5)).
301 withMaxCheck(10).
302 withThreshold(1).
303 withHandler(new ContinueOnEvent());
304 propagator.addEventDetector(d1);
305 propagator.addEventDetector(d2);
306
307
308 AbsoluteDate end = initDate.shiftedBy(30);
309 SpacecraftState actual = propagator.propagate(end);
310
311
312 Assertions.assertEquals(actual.getDate().durationFrom(end), 0.0, 0.0);
313 }
314
315 @Test
316 void testEphemerisDates() {
317
318 TimeScale tai = TimeScalesFactory.getTAI();
319 AbsoluteDate initialDate = new AbsoluteDate("2015-07-01", tai);
320 AbsoluteDate startDate = new AbsoluteDate("2015-07-03", tai).shiftedBy(-0.1);
321 AbsoluteDate endDate = new AbsoluteDate("2015-07-04", tai);
322 Frame eci = FramesFactory.getGCRF();
323 KeplerianOrbit orbit = new KeplerianOrbit(
324 600e3 + Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 0, 0, 0, 0, 0,
325 PositionAngleType.TRUE, eci, initialDate, mu);
326 OrbitType type = OrbitType.CARTESIAN;
327 double[][] tol = NumericalPropagator.tolerances(1e-3, orbit, type);
328 NumericalPropagator prop = new NumericalPropagator(
329 new DormandPrince853Integrator(0.1, 500, tol[0], tol[1]));
330 prop.setOrbitType(type);
331 prop.resetInitialState(new SpacecraftState(new CartesianOrbit(orbit)));
332
333
334 final EphemerisGenerator generator = prop.getEphemerisGenerator();
335 prop.propagate(startDate, endDate);
336 BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
337
338
339 TimeStampedPVCoordinates actualPV = ephemeris.getPVCoordinates(startDate, eci);
340 TimeStampedPVCoordinates expectedPV = orbit.getPVCoordinates(startDate, eci);
341 MatcherAssert.assertThat(actualPV.getPosition(),
342 OrekitMatchers.vectorCloseTo(expectedPV.getPosition(), 1.0));
343 MatcherAssert.assertThat(actualPV.getVelocity(),
344 OrekitMatchers.vectorCloseTo(expectedPV.getVelocity(), 1.0));
345 MatcherAssert.assertThat(ephemeris.getMinDate().durationFrom(startDate),
346 OrekitMatchers.closeTo(0, 0));
347 MatcherAssert.assertThat(ephemeris.getMaxDate().durationFrom(endDate),
348 OrekitMatchers.closeTo(0, 0));
349
350 AbsoluteDate date = endDate.shiftedBy(-0.11);
351 Assertions.assertEquals(
352 ephemeris.propagate(date).getDate().durationFrom(date), 0, 0);
353 }
354
355 @Test
356 void testEphemerisDatesBackward() {
357
358 TimeScale tai = TimeScalesFactory.getTAI();
359 AbsoluteDate initialDate = new AbsoluteDate("2015-07-05", tai);
360 AbsoluteDate startDate = new AbsoluteDate("2015-07-03", tai).shiftedBy(-0.1);
361 AbsoluteDate endDate = new AbsoluteDate("2015-07-04", tai);
362 Frame eci = FramesFactory.getGCRF();
363 KeplerianOrbit orbit = new KeplerianOrbit(
364 600e3 + Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 0, 0, 0, 0, 0,
365 PositionAngleType.TRUE, eci, initialDate, mu);
366 OrbitType type = OrbitType.CARTESIAN;
367 double[][] tol = NumericalPropagator.tolerances(1e-3, orbit, type);
368 NumericalPropagator prop = new NumericalPropagator(
369 new DormandPrince853Integrator(0.1, 500, tol[0], tol[1]));
370 prop.setOrbitType(type);
371 prop.resetInitialState(new SpacecraftState(new CartesianOrbit(orbit)));
372
373
374 final EphemerisGenerator generator = prop.getEphemerisGenerator();
375 prop.propagate(endDate, startDate);
376 BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
377
378
379 TimeStampedPVCoordinates actualPV = ephemeris.getPVCoordinates(startDate, eci);
380 TimeStampedPVCoordinates expectedPV = orbit.getPVCoordinates(startDate, eci);
381 MatcherAssert.assertThat(actualPV.getPosition(),
382 OrekitMatchers.vectorCloseTo(expectedPV.getPosition(), 1.0));
383 MatcherAssert.assertThat(actualPV.getVelocity(),
384 OrekitMatchers.vectorCloseTo(expectedPV.getVelocity(), 1.0));
385 MatcherAssert.assertThat(ephemeris.getMinDate().durationFrom(startDate),
386 OrekitMatchers.closeTo(0, 0));
387 MatcherAssert.assertThat(ephemeris.getMaxDate().durationFrom(endDate),
388 OrekitMatchers.closeTo(0, 0));
389
390 AbsoluteDate date = endDate.shiftedBy(-0.11);
391 Assertions.assertEquals(
392 ephemeris.propagate(date).getDate().durationFrom(date), 0, 0);
393 }
394
395 @Test
396 void testNoExtrapolation() {
397
398
399 final SpacecraftState finalState = propagator.propagate(initDate);
400
401
402 final Vector3D initialPosition = initialState.getPosition();
403 final Vector3D initialVelocity = initialState.getPVCoordinates().getVelocity();
404
405
406 final Vector3D finalPosition = finalState.getPosition();
407 final Vector3D finalVelocity = finalState.getPVCoordinates().getVelocity();
408
409
410 Assertions.assertEquals(initialPosition.getX(), finalPosition.getX(), 1.0e-10);
411 Assertions.assertEquals(initialPosition.getY(), finalPosition.getY(), 1.0e-10);
412 Assertions.assertEquals(initialPosition.getZ(), finalPosition.getZ(), 1.0e-10);
413 Assertions.assertEquals(initialVelocity.getX(), finalVelocity.getX(), 1.0e-10);
414 Assertions.assertEquals(initialVelocity.getY(), finalVelocity.getY(), 1.0e-10);
415 Assertions.assertEquals(initialVelocity.getZ(), finalVelocity.getZ(), 1.0e-10);
416
417 }
418
419 @Test
420 void testNotInitialised1() {
421 Assertions.assertThrows(OrekitException.class, () -> {
422 final AbstractIntegratedPropagator notInitialised =
423 new NumericalPropagator(new ClassicalRungeKuttaIntegrator(10.0));
424 notInitialised.propagate(AbsoluteDate.J2000_EPOCH);
425 });
426 }
427
428 @Test
429 void testNotInitialised2() {
430 Assertions.assertThrows(OrekitException.class, () -> {
431 final AbstractIntegratedPropagator notInitialised =
432 new NumericalPropagator(new ClassicalRungeKuttaIntegrator(10.0));
433 notInitialised.propagate(AbsoluteDate.J2000_EPOCH, AbsoluteDate.J2000_EPOCH.shiftedBy(3600));
434 });
435 }
436
437 @Test
438 void testKepler() {
439
440
441 final double dt = 3200;
442 final SpacecraftState finalState =
443 propagator.propagate(initDate.shiftedBy(-60), initDate.shiftedBy(dt));
444
445
446 final double n = FastMath.sqrt(initialState.getMu() / initialState.getA()) / initialState.getA();
447 Assertions.assertEquals(initialState.getA(), finalState.getA(), 1.0e-10);
448 Assertions.assertEquals(initialState.getEquinoctialEx(), finalState.getEquinoctialEx(), 1.0e-10);
449 Assertions.assertEquals(initialState.getEquinoctialEy(), finalState.getEquinoctialEy(), 1.0e-10);
450 Assertions.assertEquals(initialState.getHx(), finalState.getHx(), 1.0e-10);
451 Assertions.assertEquals(initialState.getHy(), finalState.getHy(), 1.0e-10);
452 Assertions.assertEquals(initialState.getLM() + n * dt, finalState.getLM(), 2.0e-9);
453
454 }
455
456 @Test
457 void testCartesian() {
458
459
460 final double dt = 3200;
461 propagator.setOrbitType(OrbitType.CARTESIAN);
462 final PVCoordinates finalState =
463 propagator.propagate(initDate.shiftedBy(dt)).getPVCoordinates();
464 final Vector3D pFin = finalState.getPosition();
465 final Vector3D vFin = finalState.getVelocity();
466
467
468 final PVCoordinates reference = initialState.shiftedBy(dt).getPVCoordinates();
469 final Vector3D pRef = reference.getPosition();
470 final Vector3D vRef = reference.getVelocity();
471 Assertions.assertEquals(0, pRef.subtract(pFin).getNorm(), 2e-4);
472 Assertions.assertEquals(0, vRef.subtract(vFin).getNorm(), 7e-8);
473
474 }
475
476 @Test
477 void testPropagationTypesElliptical() throws ParseException, IOException {
478
479 AbsoluteDate initDate = new AbsoluteDate();
480 SpacecraftState initialState;
481 final Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6);
482 final Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0);
483 initDate = AbsoluteDate.J2000_EPOCH;
484
485 final Orbit orbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
486 FramesFactory.getEME2000(), initDate, mu);
487 initialState = new SpacecraftState(orbit);
488 OrbitType type = OrbitType.EQUINOCTIAL;
489 double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit, type);
490 AdaptiveStepsizeIntegrator integrator =
491 new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
492 integrator.setInitialStepSize(60);
493 propagator = new NumericalPropagator(integrator);
494 propagator.setOrbitType(type);
495 propagator.setInitialState(initialState);
496
497 ForceModel gravityField =
498 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true),
499 GravityFieldFactory.getNormalizedProvider(5, 5));
500 propagator.addForceModel(gravityField);
501
502
503 final PVCoordinates pv = initialState.getPVCoordinates();
504 final double dP = 0.001;
505 final double dV = initialState.getMu() * dP /
506 (pv.getPosition().getNormSq() * pv.getVelocity().getNorm());
507
508 final PVCoordinates pvcM = propagateInType(initialState, dP, OrbitType.CARTESIAN, PositionAngleType.MEAN);
509 final PVCoordinates pviM = propagateInType(initialState, dP, OrbitType.CIRCULAR, PositionAngleType.MEAN);
510 final PVCoordinates pveM = propagateInType(initialState, dP, OrbitType.EQUINOCTIAL, PositionAngleType.MEAN);
511 final PVCoordinates pvkM = propagateInType(initialState, dP, OrbitType.KEPLERIAN, PositionAngleType.MEAN);
512
513 final PVCoordinates pvcE = propagateInType(initialState, dP, OrbitType.CARTESIAN, PositionAngleType.ECCENTRIC);
514 final PVCoordinates pviE = propagateInType(initialState, dP, OrbitType.CIRCULAR, PositionAngleType.ECCENTRIC);
515 final PVCoordinates pveE = propagateInType(initialState, dP, OrbitType.EQUINOCTIAL, PositionAngleType.ECCENTRIC);
516 final PVCoordinates pvkE = propagateInType(initialState, dP, OrbitType.KEPLERIAN, PositionAngleType.ECCENTRIC);
517
518 final PVCoordinates pvcT = propagateInType(initialState, dP, OrbitType.CARTESIAN, PositionAngleType.TRUE);
519 final PVCoordinates pviT = propagateInType(initialState, dP, OrbitType.CIRCULAR, PositionAngleType.TRUE);
520 final PVCoordinates pveT = propagateInType(initialState, dP, OrbitType.EQUINOCTIAL, PositionAngleType.TRUE);
521 final PVCoordinates pvkT = propagateInType(initialState, dP, OrbitType.KEPLERIAN, PositionAngleType.TRUE);
522
523 Assertions.assertEquals(0, pvcM.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 3.0);
524 Assertions.assertEquals(0, pvcM.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 2.0);
525 Assertions.assertEquals(0, pviM.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 0.6);
526 Assertions.assertEquals(0, pviM.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 0.4);
527 Assertions.assertEquals(0, pvkM.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 0.5);
528 Assertions.assertEquals(0, pvkM.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 0.3);
529 Assertions.assertEquals(0, pveM.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 0.2);
530 Assertions.assertEquals(0, pveM.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 0.2);
531
532 Assertions.assertEquals(0, pvcE.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 3.0);
533 Assertions.assertEquals(0, pvcE.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 2.0);
534 Assertions.assertEquals(0, pviE.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 0.03);
535 Assertions.assertEquals(0, pviE.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 0.04);
536 Assertions.assertEquals(0, pvkE.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 0.4);
537 Assertions.assertEquals(0, pvkE.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 0.3);
538 Assertions.assertEquals(0, pveE.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 0.2);
539 Assertions.assertEquals(0, pveE.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 0.07);
540
541 Assertions.assertEquals(0, pvcT.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 3.0);
542 Assertions.assertEquals(0, pvcT.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 2.0);
543 Assertions.assertEquals(0, pviT.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 0.3);
544 Assertions.assertEquals(0, pviT.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 0.2);
545 Assertions.assertEquals(0, pvkT.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 0.4);
546 Assertions.assertEquals(0, pvkT.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 0.2);
547
548 }
549
550 @Test
551 void testPropagationTypesHyperbolic() throws ParseException, IOException {
552
553 SpacecraftState state =
554 new SpacecraftState(new KeplerianOrbit(-10000000.0, 2.5, 0.3, 0, 0, 0.0,
555 PositionAngleType.TRUE,
556 FramesFactory.getEME2000(), initDate,
557 mu));
558
559 ForceModel gravityField =
560 new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true),
561 GravityFieldFactory.getNormalizedProvider(5, 5));
562 propagator.addForceModel(gravityField);
563
564
565 final PVCoordinates pv = state.getPVCoordinates();
566 final double dP = 0.001;
567 final double dV = state.getMu() * dP /
568 (pv.getPosition().getNormSq() * pv.getVelocity().getNorm());
569
570 final PVCoordinates pvcM = propagateInType(state, dP, OrbitType.CARTESIAN, PositionAngleType.MEAN);
571 final PVCoordinates pvkM = propagateInType(state, dP, OrbitType.KEPLERIAN, PositionAngleType.MEAN);
572
573 final PVCoordinates pvcE = propagateInType(state, dP, OrbitType.CARTESIAN, PositionAngleType.ECCENTRIC);
574 final PVCoordinates pvkE = propagateInType(state, dP, OrbitType.KEPLERIAN, PositionAngleType.ECCENTRIC);
575
576 final PVCoordinates pvcT = propagateInType(state, dP, OrbitType.CARTESIAN, PositionAngleType.TRUE);
577 final PVCoordinates pvkT = propagateInType(state, dP, OrbitType.KEPLERIAN, PositionAngleType.TRUE);
578
579 Assertions.assertEquals(0, pvcM.getPosition().subtract(pvkT.getPosition()).getNorm() / dP, 0.3);
580 Assertions.assertEquals(0, pvcM.getVelocity().subtract(pvkT.getVelocity()).getNorm() / dV, 0.4);
581 Assertions.assertEquals(0, pvkM.getPosition().subtract(pvkT.getPosition()).getNorm() / dP, 0.2);
582 Assertions.assertEquals(0, pvkM.getVelocity().subtract(pvkT.getVelocity()).getNorm() / dV, 0.3);
583
584 Assertions.assertEquals(0, pvcE.getPosition().subtract(pvkT.getPosition()).getNorm() / dP, 0.3);
585 Assertions.assertEquals(0, pvcE.getVelocity().subtract(pvkT.getVelocity()).getNorm() / dV, 0.4);
586 Assertions.assertEquals(0, pvkE.getPosition().subtract(pvkT.getPosition()).getNorm() / dP, 0.009);
587 Assertions.assertEquals(0, pvkE.getVelocity().subtract(pvkT.getVelocity()).getNorm() / dV, 0.006);
588
589 Assertions.assertEquals(0, pvcT.getPosition().subtract(pvkT.getPosition()).getNorm() / dP, 0.3);
590 Assertions.assertEquals(0, pvcT.getVelocity().subtract(pvkT.getVelocity()).getNorm() / dV, 0.4);
591
592 }
593
594 private PVCoordinates propagateInType(SpacecraftState state, double dP,
595 OrbitType type, PositionAngleType angle)
596 {
597
598 final double dt = 3200;
599 final double minStep = 0.001;
600 final double maxStep = 1000;
601
602 double[][] tol = NumericalPropagator.tolerances(dP, state.getOrbit(), type);
603 AdaptiveStepsizeIntegrator integrator =
604 new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
605 NumericalPropagator newPropagator = new NumericalPropagator(integrator);
606 newPropagator.setOrbitType(type);
607 newPropagator.setPositionAngleType(angle);
608 newPropagator.setInitialState(state);
609 for (ForceModel force: propagator.getAllForceModels()) {
610 newPropagator.addForceModel(force);
611 }
612 return newPropagator.propagate(state.getDate().shiftedBy(dt)).getPVCoordinates();
613
614 }
615
616 @Test
617 void testException() {
618 Assertions.assertThrows(OrekitException.class, () -> {
619 propagator.setStepHandler(new OrekitStepHandler() {
620 private int countDown = 3;
621 private AbsoluteDate previousCall = null;
622 public void init(SpacecraftState s0, AbsoluteDate t) {
623 }
624 public void handleStep(OrekitStepInterpolator interpolator) {
625 if (previousCall != null) {
626 Assertions.assertTrue(interpolator.getCurrentState().getDate().compareTo(previousCall) < 0);
627 }
628 if (--countDown == 0) {
629 throw new OrekitException(LocalizedCoreFormats.SIMPLE_MESSAGE, "dummy error");
630 }
631 }
632 });
633 propagator.propagate(initDate.shiftedBy(-3600));
634 });
635 }
636
637 @Test
638 void testStopEvent() {
639 final AbsoluteDate stopDate = initDate.shiftedBy(1000);
640 CheckingHandler checking = new CheckingHandler(Action.STOP);
641 propagator.addEventDetector(new DateDetector(stopDate).withHandler(checking));
642 Assertions.assertEquals(1, propagator.getEventsDetectors().size());
643 checking.assertEvent(false);
644 final SpacecraftState finalState = propagator.propagate(initDate.shiftedBy(3200));
645 checking.assertEvent(true);
646 Assertions.assertEquals(0, finalState.getDate().durationFrom(stopDate), 1.0e-10);
647 propagator.clearEventsDetectors();
648 Assertions.assertEquals(0, propagator.getEventsDetectors().size());
649 Assertions.assertTrue(checking.isFinished);
650 }
651
652 @Test
653 void testResetStateEvent() {
654 final AbsoluteDate resetDate = initDate.shiftedBy(1000);
655 CheckingHandler checking = new CheckingHandler(Action.RESET_STATE) {
656 public SpacecraftState resetState(EventDetector detector, SpacecraftState oldState) {
657 return new SpacecraftState(oldState.getOrbit(), oldState.getAttitude(), oldState.getMass() - 200.0);
658 }
659 };
660 propagator.addEventDetector(new DateDetector(resetDate).withHandler(checking));
661 checking.assertEvent(false);
662 final SpacecraftState finalState = propagator.propagate(initDate.shiftedBy(3200));
663 checking.assertEvent(true);
664 Assertions.assertEquals(initialState.getMass() - 200, finalState.getMass(), 1.0e-10);
665 }
666
667 @Test
668 void testResetDerivativesEvent() {
669 final AbsoluteDate resetDate = initDate.shiftedBy(1000);
670 CheckingHandler checking = new CheckingHandler(Action.RESET_DERIVATIVES);
671 propagator.addEventDetector(new DateDetector(resetDate).withHandler(checking));
672 final double dt = 3200;
673 checking.assertEvent(false);
674 Assertions.assertEquals(0.0, propagator.getInitialState().getDate().durationFrom(initDate), 1.0e-10);
675 propagator.setResetAtEnd(true);
676 final SpacecraftState finalState =
677 propagator.propagate(initDate.shiftedBy(dt));
678 Assertions.assertEquals(dt, propagator.getInitialState().getDate().durationFrom(initDate), 1.0e-10);
679 checking.assertEvent(true);
680 final double n = FastMath.sqrt(initialState.getMu() / initialState.getA()) / initialState.getA();
681 Assertions.assertEquals(initialState.getA(), finalState.getA(), 1.0e-10);
682 Assertions.assertEquals(initialState.getEquinoctialEx(), finalState.getEquinoctialEx(), 1.0e-10);
683 Assertions.assertEquals(initialState.getEquinoctialEy(), finalState.getEquinoctialEy(), 1.0e-10);
684 Assertions.assertEquals(initialState.getHx(), finalState.getHx(), 1.0e-10);
685 Assertions.assertEquals(initialState.getHy(), finalState.getHy(), 1.0e-10);
686 Assertions.assertEquals(initialState.getLM() + n * dt, finalState.getLM(), 6.0e-10);
687 }
688
689 @Test
690 void testContinueEvent() {
691 final AbsoluteDate resetDate = initDate.shiftedBy(1000);
692 CheckingHandler checking = new CheckingHandler(Action.CONTINUE);
693 propagator.addEventDetector(new DateDetector(resetDate).withHandler(checking));
694 final double dt = 3200;
695 checking.assertEvent(false);
696 Assertions.assertEquals(0.0, propagator.getInitialState().getDate().durationFrom(initDate), 1.0e-10);
697 propagator.setResetAtEnd(false);
698 final SpacecraftState finalState =
699 propagator.propagate(initDate.shiftedBy(dt));
700 Assertions.assertEquals(0.0, propagator.getInitialState().getDate().durationFrom(initDate), 1.0e-10);
701 checking.assertEvent(true);
702 final double n = FastMath.sqrt(initialState.getMu() / initialState.getA()) / initialState.getA();
703 Assertions.assertEquals(initialState.getA(), finalState.getA(), 1.0e-10);
704 Assertions.assertEquals(initialState.getEquinoctialEx(), finalState.getEquinoctialEx(), 1.0e-10);
705 Assertions.assertEquals(initialState.getEquinoctialEy(), finalState.getEquinoctialEy(), 1.0e-10);
706 Assertions.assertEquals(initialState.getHx(), finalState.getHx(), 1.0e-10);
707 Assertions.assertEquals(initialState.getHy(), finalState.getHy(), 1.0e-10);
708 Assertions.assertEquals(initialState.getLM() + n * dt, finalState.getLM(), 6.0e-10);
709 }
710
711 @Test
712 void testAdditionalStateEvent() {
713 propagator.addAdditionalDerivativesProvider(new AdditionalDerivativesProvider() {
714
715 public String getName() {
716 return "linear";
717 }
718
719 public int getDimension() {
720 return 1;
721 }
722
723 public CombinedDerivatives combinedDerivatives(SpacecraftState s) {
724 return new CombinedDerivatives(new double[] { 1.0 }, null);
725 }
726
727 });
728 try {
729 propagator.addAdditionalDerivativesProvider(new AdditionalDerivativesProvider() {
730
731 public String getName() {
732 return "linear";
733 }
734
735 public int getDimension() {
736 return 1;
737 }
738
739 public CombinedDerivatives combinedDerivatives(SpacecraftState s) {
740 return new CombinedDerivatives(new double[] { 1.0 }, null);
741 }
742
743 });
744 Assertions.fail("an exception should have been thrown");
745 } catch (OrekitException oe) {
746 Assertions.assertEquals(oe.getSpecifier(), OrekitMessages.ADDITIONAL_STATE_NAME_ALREADY_IN_USE);
747 }
748 try {
749 propagator.addAdditionalDerivativesProvider(new AdditionalDerivativesProvider() {
750 public String getName() {
751 return "linear";
752 }
753
754 public int getDimension() {
755 return 1;
756 }
757
758 public CombinedDerivatives combinedDerivatives(SpacecraftState s) {
759 return new CombinedDerivatives(new double[1], null);
760 }
761 });
762 Assertions.fail("an exception should have been thrown");
763 } catch (OrekitException oe) {
764 Assertions.assertEquals(oe.getSpecifier(), OrekitMessages.ADDITIONAL_STATE_NAME_ALREADY_IN_USE);
765 }
766 propagator.addAdditionalStateProvider(new AdditionalStateProvider() {
767 public String getName() {
768 return "constant";
769 }
770
771 public double[] getAdditionalState(SpacecraftState state) {
772 return new double[] { 1.0 };
773 }
774 });
775 Assertions.assertTrue(propagator.isAdditionalStateManaged("linear"));
776 Assertions.assertTrue(propagator.isAdditionalStateManaged("constant"));
777 Assertions.assertFalse(propagator.isAdditionalStateManaged("non-managed"));
778 Assertions.assertEquals(2, propagator.getManagedAdditionalStates().length);
779 propagator.setInitialState(propagator.getInitialState().addAdditionalState("linear", 1.5));
780
781 CheckingHandler checking = new CheckingHandler(Action.STOP);
782 propagator.addEventDetector(new AdditionalStateLinearDetector(10.0, 1.0e-8).withHandler(checking));
783
784 final double dt = 3200;
785 checking.assertEvent(false);
786 final SpacecraftState finalState =
787 propagator.propagate(initDate.shiftedBy(dt));
788 checking.assertEvent(true);
789 Assertions.assertEquals(3.0, finalState.getAdditionalState("linear")[0], 1.0e-8);
790 Assertions.assertEquals(1.5, finalState.getDate().durationFrom(initDate), 1.0e-8);
791
792 }
793
794 private static class AdditionalStateLinearDetector extends AbstractDetector<AdditionalStateLinearDetector> {
795
796 public AdditionalStateLinearDetector(double maxCheck, double threshold) {
797 this(AdaptableInterval.of(maxCheck), threshold, DEFAULT_MAX_ITER, new StopOnEvent());
798 }
799
800 private AdditionalStateLinearDetector(AdaptableInterval maxCheck, double threshold, int maxIter, EventHandler handler) {
801 super(maxCheck, threshold, maxIter, handler);
802 }
803
804 protected AdditionalStateLinearDetector create(final AdaptableInterval newMaxCheck, final double newThreshold,
805 final int newMaxIter, final EventHandler newHandler) {
806 return new AdditionalStateLinearDetector(newMaxCheck, newThreshold, newMaxIter, newHandler);
807 }
808
809 public double g(SpacecraftState s) {
810 return s.getAdditionalState("linear")[0] - 3.0;
811 }
812
813 }
814
815 @Test
816 void testResetAdditionalStateEvent() {
817 propagator.addAdditionalDerivativesProvider(new AdditionalDerivativesProvider() {
818
819 public String getName() {
820 return "linear";
821 }
822
823 public int getDimension() {
824 return 1;
825 }
826
827 public CombinedDerivatives combinedDerivatives(SpacecraftState s) {
828 return new CombinedDerivatives(new double[] { 1.0 }, null);
829 }
830 });
831 propagator.setInitialState(propagator.getInitialState().addAdditionalState("linear", 1.5));
832
833 CheckingHandler checking = new CheckingHandler(Action.RESET_STATE) {
834 public SpacecraftState resetState(EventDetector detector, SpacecraftState oldState)
835 {
836 return oldState.addAdditionalState("linear", oldState.getAdditionalState("linear")[0] * 2);
837 }
838 };
839
840 propagator.addEventDetector(new AdditionalStateLinearDetector(10.0, 1.0e-8).withHandler(checking));
841
842 final double dt = 3200;
843 checking.assertEvent(false);
844 final SpacecraftState finalState = propagator.propagate(initDate.shiftedBy(dt));
845 checking.assertEvent(true);
846 Assertions.assertEquals(dt + 4.5, finalState.getAdditionalState("linear")[0], 1.0e-8);
847 Assertions.assertEquals(dt, finalState.getDate().durationFrom(initDate), 1.0e-8);
848
849 }
850
851 @Test
852 void testEventDetectionBug() throws IOException, ParseException {
853
854 TimeScale utc = TimeScalesFactory.getUTC();
855 AbsoluteDate initialDate = new AbsoluteDate(2005, 1, 1, 0, 0, 0.0, utc);
856 double duration = 100000.;
857 AbsoluteDate endDate = new AbsoluteDate(initialDate, duration);
858
859
860 Frame EME2000 = FramesFactory.getEME2000();
861
862
863
864 double a = 35786000. + 6378137.0;
865 double e = 0.70;
866 double rApogee = a*(1+e);
867 double vApogee = FastMath.sqrt(mu*(1-e)/(a*(1+e)));
868 Orbit geo = new CartesianOrbit(new PVCoordinates(new Vector3D(rApogee, 0., 0.),
869 new Vector3D(0., vApogee, 0.)), EME2000,
870 initialDate, mu);
871
872
873 duration = geo.getKeplerianPeriod();
874 endDate = new AbsoluteDate(initialDate, duration);
875
876
877 final double minStep = 0.001;
878 final double maxStep = 1000;
879 final double initStep = 60;
880 final double[] absTolerance = {
881 0.001, 1.0e-9, 1.0e-9, 1.0e-6, 1.0e-6, 1.0e-6, 0.001};
882 final double[] relTolerance = {
883 1.0e-7, 1.0e-4, 1.0e-4, 1.0e-7, 1.0e-7, 1.0e-7, 1.0e-7};
884
885 AdaptiveStepsizeIntegrator integrator =
886 new DormandPrince853Integrator(minStep, maxStep, absTolerance, relTolerance);
887 integrator.setInitialStepSize(initStep);
888
889
890 propagator = new NumericalPropagator(integrator);
891 double mass = 1000.;
892 SpacecraftState initialState = new SpacecraftState(geo, mass);
893 propagator.setInitialState(initialState);
894 propagator.setOrbitType(OrbitType.CARTESIAN);
895
896
897
898 ApsideDetector event1 = new ApsideDetector(geo);
899 propagator.addEventDetector(event1);
900
901
902 propagator.clearStepHandlers();
903
904
905 SpacecraftState finalState = propagator.propagate(endDate);
906
907
908 Assertions.assertTrue(endDate.durationFrom(finalState.getDate()) > 40000.0);
909 }
910
911 @Test
912 void testEphemerisGenerationIssue14() throws IOException {
913
914
915 final double dt = 3200;
916 propagator.getInitialState();
917
918 propagator.setOrbitType(OrbitType.CARTESIAN);
919 final EphemerisGenerator generator = propagator.getEphemerisGenerator();
920 propagator.propagate(initDate.shiftedBy(dt));
921 final BoundedPropagator ephemeris1 = generator.getGeneratedEphemeris();
922 Assertions.assertEquals(initDate, ephemeris1.getMinDate());
923 Assertions.assertEquals(initDate.shiftedBy(dt), ephemeris1.getMaxDate());
924
925 propagator.getPVCoordinates(initDate.shiftedBy( 2 * dt), FramesFactory.getEME2000());
926 propagator.getPVCoordinates(initDate.shiftedBy(-2 * dt), FramesFactory.getEME2000());
927
928
929 Assertions.assertEquals(initDate, ephemeris1.getMinDate());
930 Assertions.assertEquals(initDate.shiftedBy(dt), ephemeris1.getMaxDate());
931
932 final BoundedPropagator ephemeris2 = generator.getGeneratedEphemeris();
933 Assertions.assertEquals(initDate.shiftedBy(-2 * dt), ephemeris2.getMinDate());
934 Assertions.assertEquals(initDate.shiftedBy( 2 * dt), ephemeris2.getMaxDate());
935
936
937 Assertions.assertEquals(initDate, ephemeris1.getMinDate());
938 Assertions.assertEquals(initDate.shiftedBy(dt), ephemeris1.getMaxDate());
939
940 }
941
942 @Test
943 void testEphemerisAdditionalState() throws IOException {
944
945
946 final double dt = -3200;
947 final double rate = 2.0;
948
949 propagator.addAdditionalStateProvider(new AdditionalStateProvider() {
950 public String getName() {
951 return "squaredA";
952 }
953 public double[] getAdditionalState(SpacecraftState state) {
954 return new double[] { state.getA() * state.getA() };
955 }
956 });
957 propagator.addAdditionalDerivativesProvider(new AdditionalDerivativesProvider() {
958 public String getName() {
959 return "extra";
960 }
961 public int getDimension() {
962 return 1;
963 }
964 public CombinedDerivatives combinedDerivatives(SpacecraftState s) {
965 return new CombinedDerivatives(new double[] { rate }, null);
966 }
967 });
968 propagator.setInitialState(propagator.getInitialState().addAdditionalState("extra", 1.5));
969
970 propagator.setOrbitType(OrbitType.CARTESIAN);
971 final EphemerisGenerator generator = propagator.getEphemerisGenerator();
972 propagator.propagate(initDate.shiftedBy(dt));
973 final BoundedPropagator ephemeris1 = generator.getGeneratedEphemeris();
974 Assertions.assertEquals(initDate.shiftedBy(dt), ephemeris1.getMinDate());
975 Assertions.assertEquals(initDate, ephemeris1.getMaxDate());
976 try {
977 ephemeris1.propagate(ephemeris1.getMinDate().shiftedBy(-10.0));
978 Assertions.fail("an exception should have been thrown");
979 } catch (OrekitException pe) {
980 Assertions.assertEquals(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE_BEFORE, pe.getSpecifier());
981 }
982 try {
983 ephemeris1.propagate(ephemeris1.getMaxDate().shiftedBy(+10.0));
984 Assertions.fail("an exception should have been thrown");
985 } catch (OrekitException pe) {
986 Assertions.assertEquals(OrekitMessages.OUT_OF_RANGE_EPHEMERIDES_DATE_AFTER, pe.getSpecifier());
987 }
988
989 double shift = -60;
990 SpacecraftState s = ephemeris1.propagate(initDate.shiftedBy(shift));
991 Assertions.assertEquals(2, s.getAdditionalStatesValues().size());
992 Assertions.assertTrue(s.hasAdditionalState("squaredA"));
993 Assertions.assertTrue(s.hasAdditionalState("extra"));
994 Assertions.assertEquals(s.getA() * s.getA(), s.getAdditionalState("squaredA")[0], 1.0e-10);
995 Assertions.assertEquals(1.5 + shift * rate, s.getAdditionalState("extra")[0], 1.0e-10);
996
997 try {
998 ephemeris1.resetInitialState(s);
999 Assertions.fail("an exception should have been thrown");
1000 } catch (OrekitException oe) {
1001 Assertions.assertEquals(OrekitMessages.NON_RESETABLE_STATE, oe.getSpecifier());
1002 }
1003
1004 }
1005
1006 @Test
1007 void testIssue157() {
1008 try {
1009 Orbit orbit = new KeplerianOrbit(13378000, 0.05, 0, 0, FastMath.PI, 0, PositionAngleType.MEAN,
1010 FramesFactory.getTOD(false),
1011 new AbsoluteDate(2003, 5, 6, TimeScalesFactory.getUTC()),
1012 Constants.EIGEN5C_EARTH_MU);
1013 NumericalPropagator.tolerances(1.0, orbit, OrbitType.KEPLERIAN);
1014 Assertions.fail("an exception should have been thrown");
1015 } catch (OrekitException pe) {
1016 Assertions.assertEquals(OrekitMessages.SINGULAR_JACOBIAN_FOR_ORBIT_TYPE, pe.getSpecifier());
1017 }
1018 }
1019
1020 @Test
1021 void testIssue704() {
1022
1023
1024 final Orbit orbit = initialState.getOrbit();
1025 final PVCoordinates pv = orbit.getPVCoordinates();
1026
1027
1028 final double dP = 10.0;
1029
1030
1031 final double r2 = pv.getPosition().getNormSq();
1032 final double v = pv.getVelocity().getNorm();
1033 final double dV = orbit.getMu() * dP / (v * r2);
1034
1035
1036 final double[][] tolCart1 = NumericalPropagator.tolerances(dP, orbit, OrbitType.CARTESIAN);
1037 final double[][] tolCart2 = NumericalPropagator.tolerances(dP, dV, orbit, OrbitType.CARTESIAN);
1038 for (int i = 0; i < tolCart1.length; i++) {
1039 Assertions.assertArrayEquals(tolCart1[i], tolCart2[i], Double.MIN_VALUE);
1040 }
1041
1042
1043 final double[][] tolKep1 = NumericalPropagator.tolerances(dP, orbit, OrbitType.KEPLERIAN);
1044 final double[][] tolKep2 = NumericalPropagator.tolerances(dP, dV, orbit, OrbitType.KEPLERIAN);
1045 for (int i = 0; i < tolCart1.length; i++) {
1046 Assertions.assertArrayEquals(tolKep1[i], tolKep2[i], Double.MIN_VALUE);
1047 }
1048
1049 }
1050
1051 private static class CheckingHandler implements EventHandler {
1052
1053 private final Action actionOnEvent;
1054 private boolean gotHere;
1055 private boolean isFinished = false;
1056
1057 public CheckingHandler(final Action actionOnEvent) {
1058 this.actionOnEvent = actionOnEvent;
1059 this.gotHere = false;
1060 }
1061
1062 public void assertEvent(boolean expected) {
1063 Assertions.assertEquals(expected, gotHere);
1064 }
1065
1066 public Action eventOccurred(SpacecraftState s, EventDetector detector, boolean increasing) {
1067 gotHere = true;
1068 return actionOnEvent;
1069 }
1070
1071 @Override
1072 public void finish(SpacecraftState finalState, EventDetector detector) {
1073 isFinished = true;
1074 }
1075 }
1076
1077 @Test
1078 void testParallelismIssue258() throws InterruptedException, ExecutionException, FileNotFoundException {
1079
1080 Utils.setDataRoot("regular-data:atmosphere:potential/grgs-format");
1081 GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
1082 final double mu = GravityFieldFactory.getNormalizedProvider(2, 2).getMu();
1083
1084
1085 final double a = 24396159;
1086 final double e = 0.72831215;
1087 final double i = FastMath.toRadians(7);
1088 final double omega = FastMath.toRadians(180);
1089 final double raan = FastMath.toRadians(261);
1090 final double lM = 0;
1091 final Frame inertialFrame = FramesFactory.getEME2000();
1092 final TimeScale utc = TimeScalesFactory.getUTC();
1093 final AbsoluteDate initialDate = new AbsoluteDate(2003, 1, 1, 00, 00, 00.000, utc);
1094 final Orbit initialOrbit = new CartesianOrbit( new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngleType.MEAN,
1095 inertialFrame, initialDate, mu));
1096 final SpacecraftState initialState = new SpacecraftState(initialOrbit, 1000);
1097
1098
1099 final List<SpacecraftState> states = new ArrayList<SpacecraftState>();
1100 final NumericalPropagator propagator = createPropagator(initialState, OrbitType.CARTESIAN, PositionAngleType.TRUE);
1101 final double samplingStep = 10000.0;
1102 propagator.setStepHandler(samplingStep, state -> states.add(state));
1103 propagator.propagate(initialDate.shiftedBy(5 * samplingStep));
1104
1105
1106 final double[][] referenceErrors = new double[states.size() - 1][];
1107 for (int startIndex = 0; startIndex < states.size() - 1; ++startIndex) {
1108 referenceErrors[startIndex] = recomputeFollowing(startIndex, states);
1109 }
1110
1111 final Consumer<SpacecraftState> checker = point -> {
1112 try {
1113 final int startIndex = states.indexOf(point);
1114 double[] errors = recomputeFollowing(startIndex, states);
1115 for (int k = 0; k < errors.length; ++k) {
1116 Assertions.assertEquals(
1117 referenceErrors[startIndex][k], errors[k],
1118 1.0e-9,startIndex + " to " + (startIndex + k + 1));
1119 }
1120 } catch (OrekitException oe) {
1121 Assertions.fail(oe.getLocalizedMessage());
1122 }
1123 };
1124
1125
1126 states.stream().forEach(checker);
1127
1128
1129 states.parallelStream().forEach(checker);
1130
1131 }
1132
1133 @Test
1134 void testShiftKeplerianEllipticTrueWithoutDerivatives() {
1135 doTestShift(createEllipticOrbit(), OrbitType.KEPLERIAN, PositionAngleType.TRUE, false,
1136 18.1, 72.0, 437.3, 1601.1, 3141.8);
1137 }
1138
1139 @Test
1140 void testShiftKeplerianEllipticTrueWithDerivatives() {
1141 doTestShift(createEllipticOrbit(), OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
1142 1.14, 9.1, 140.3, 1066.7, 3306.9);
1143 }
1144
1145 @Test
1146 void testShiftKeplerianEllipticEccentricWithoutDerivatives() {
1147 doTestShift(createEllipticOrbit(), OrbitType.KEPLERIAN, PositionAngleType.ECCENTRIC, false,
1148 18.1, 72.0, 437.3, 1601.1, 3141.8);
1149 }
1150
1151 @Test
1152 void testShiftKeplerianEllipticEcentricWithDerivatives() {
1153 doTestShift(createEllipticOrbit(), OrbitType.KEPLERIAN, PositionAngleType.ECCENTRIC, true,
1154 1.14, 9.1, 140.3, 1066.7, 3306.9);
1155 }
1156
1157 @Test
1158 void testShiftKeplerianEllipticMeanWithoutDerivatives() {
1159 doTestShift(createEllipticOrbit(), OrbitType.KEPLERIAN, PositionAngleType.MEAN, false,
1160 18.1, 72.0, 437.3, 1601.1, 3141.8);
1161 }
1162
1163 @Test
1164 void testShiftKeplerianEllipticMeanWithDerivatives() {
1165 doTestShift(createEllipticOrbit(), OrbitType.KEPLERIAN, PositionAngleType.MEAN, true,
1166 1.14, 9.1, 140.3, 1066.7, 3306.9);
1167 }
1168
1169 @Test
1170 void testShiftKeplerianHyperbolicTrueWithoutDerivatives() {
1171 doTestShift(createHyperbolicOrbit(), OrbitType.KEPLERIAN, PositionAngleType.TRUE, false,
1172 0.484, 1.94, 12.1, 48.3, 108.5);
1173 }
1174
1175 @Test
1176 void testShiftKeplerianHyperbolicTrueWithDerivatives() {
1177 doTestShift(createHyperbolicOrbit(), OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
1178 1.38e-4, 1.10e-3, 1.72e-2, 1.37e-1, 4.62e-1);
1179 }
1180
1181 @Test
1182 void testShiftKeplerianHyperbolicEccentricWithoutDerivatives() {
1183 doTestShift(createHyperbolicOrbit(), OrbitType.KEPLERIAN, PositionAngleType.ECCENTRIC, false,
1184 0.484, 1.94, 12.1, 48.3, 108.5);
1185 }
1186
1187 @Test
1188 void testShiftKeplerianHyperbolicEcentricWithDerivatives() {
1189 doTestShift(createHyperbolicOrbit(), OrbitType.KEPLERIAN, PositionAngleType.ECCENTRIC, true,
1190 1.38e-4, 1.10e-3, 1.72e-2, 1.37e-1, 4.62e-1);
1191 }
1192
1193 @Test
1194 void testShiftKeplerianHyperbolicMeanWithoutDerivatives() {
1195 doTestShift(createHyperbolicOrbit(), OrbitType.KEPLERIAN, PositionAngleType.MEAN, false,
1196 0.484, 1.94, 12.1, 48.3, 108.5);
1197 }
1198
1199 @Test
1200 void testShiftKeplerianHyperbolicMeanWithDerivatives() {
1201 doTestShift(createHyperbolicOrbit(), OrbitType.KEPLERIAN, PositionAngleType.MEAN, true,
1202 1.38e-4, 1.10e-3, 1.72e-2, 1.37e-1, 4.62e-1);
1203 }
1204
1205 @Test
1206 void testShiftCartesianEllipticTrueWithoutDerivatives() {
1207 doTestShift(createEllipticOrbit(), OrbitType.CARTESIAN, PositionAngleType.TRUE, false,
1208 18.1, 72.0, 437.3, 1601.1, 3141.8);
1209 }
1210
1211 @Test
1212 void testShiftCartesianEllipticTrueWithDerivatives() {
1213 doTestShift(createEllipticOrbit(), OrbitType.CARTESIAN, PositionAngleType.TRUE, true,
1214 1.14, 9.1, 140.3, 1066.7, 3306.9);
1215 }
1216
1217 @Test
1218 void testShiftCartesianEllipticEccentricWithoutDerivatives() {
1219 doTestShift(createEllipticOrbit(), OrbitType.CARTESIAN, PositionAngleType.ECCENTRIC, false,
1220 18.1, 72.0, 437.3, 1601.1, 3141.8);
1221 }
1222
1223 @Test
1224 void testShiftCartesianEllipticEcentricWithDerivatives() {
1225 doTestShift(createEllipticOrbit(), OrbitType.CARTESIAN, PositionAngleType.ECCENTRIC, true,
1226 1.14, 9.1, 140.3, 1066.7, 3306.9);
1227 }
1228
1229 @Test
1230 void testShiftCartesianEllipticMeanWithoutDerivatives() {
1231 doTestShift(createEllipticOrbit(), OrbitType.CARTESIAN, PositionAngleType.MEAN, false,
1232 18.1, 72.0, 437.3, 1601.1, 3141.8);
1233 }
1234
1235 @Test
1236 void testShiftCartesianEllipticMeanWithDerivatives() {
1237 doTestShift(createEllipticOrbit(), OrbitType.CARTESIAN, PositionAngleType.MEAN, true,
1238 1.14, 9.1, 140.3, 1066.7, 3306.9);
1239 }
1240
1241 @Test
1242 void testShiftCartesianHyperbolicTrueWithoutDerivatives() {
1243 doTestShift(createHyperbolicOrbit(), OrbitType.CARTESIAN, PositionAngleType.TRUE, false,
1244 0.48, 1.93, 12.1, 48.3, 108.5);
1245 }
1246
1247 @Test
1248 void testShiftCartesianHyperbolicTrueWithDerivatives() {
1249 doTestShift(createHyperbolicOrbit(), OrbitType.CARTESIAN, PositionAngleType.TRUE, true,
1250 1.38e-4, 1.10e-3, 1.72e-2, 1.37e-1, 4.62e-1);
1251 }
1252
1253 @Test
1254 void testShiftCartesianHyperbolicEccentricWithoutDerivatives() {
1255 doTestShift(createHyperbolicOrbit(), OrbitType.CARTESIAN, PositionAngleType.ECCENTRIC, false,
1256 0.48, 1.93, 12.1, 48.3, 108.5);
1257 }
1258
1259 @Test
1260 void testShiftCartesianHyperbolicEcentricWithDerivatives() {
1261 doTestShift(createHyperbolicOrbit(), OrbitType.CARTESIAN, PositionAngleType.ECCENTRIC, true,
1262 1.38e-4, 1.10e-3, 1.72e-2, 1.37e-1, 4.62e-1);
1263 }
1264
1265 @Test
1266 void testShiftCartesianHyperbolicMeanWithoutDerivatives() {
1267 doTestShift(createHyperbolicOrbit(), OrbitType.CARTESIAN, PositionAngleType.MEAN, false,
1268 0.48, 1.93, 12.1, 48.3, 108.5);
1269 }
1270
1271 @Test
1272 void testShiftCartesianHyperbolicMeanWithDerivatives() {
1273 doTestShift(createHyperbolicOrbit(), OrbitType.CARTESIAN, PositionAngleType.MEAN, true,
1274 1.38e-4, 1.10e-3, 1.72e-2, 1.37e-1, 4.62e-1);
1275 }
1276
1277 @Test
1278 void testShiftCircularTrueWithoutDerivatives() {
1279 doTestShift(createEllipticOrbit(), OrbitType.CIRCULAR, PositionAngleType.TRUE, false,
1280 18.1, 72.0, 437.3, 1601.1, 3141.8);
1281 }
1282
1283 @Test
1284 void testShiftCircularTrueWithDerivatives() {
1285 doTestShift(createEllipticOrbit(), OrbitType.CIRCULAR, PositionAngleType.TRUE, true,
1286 1.14, 9.1, 140.3, 1066.7, 3306.9);
1287 }
1288
1289 @Test
1290 void testShiftCircularEccentricWithoutDerivatives() {
1291 doTestShift(createEllipticOrbit(), OrbitType.CIRCULAR, PositionAngleType.ECCENTRIC, false,
1292 18.1, 72.0, 437.3, 1601.1, 3141.8);
1293 }
1294
1295 @Test
1296 void testShiftCircularEcentricWithDerivatives() {
1297 doTestShift(createEllipticOrbit(), OrbitType.CIRCULAR, PositionAngleType.ECCENTRIC, true,
1298 1.14, 9.1, 140.3, 1066.7, 3306.9);
1299 }
1300
1301 @Test
1302 void testShiftCircularMeanWithoutDerivatives() {
1303 doTestShift(createEllipticOrbit(), OrbitType.CIRCULAR, PositionAngleType.MEAN, false,
1304 18.1, 72.0, 437.3, 1601.1, 3141.8);
1305 }
1306
1307 @Test
1308 void testShiftCircularMeanWithDerivatives() {
1309 doTestShift(createEllipticOrbit(), OrbitType.CIRCULAR, PositionAngleType.MEAN, true,
1310 1.14, 9.1, 140.3, 1066.7, 3306.9);
1311 }
1312
1313 @Test
1314 void testShiftEquinoctialTrueWithoutDerivatives() {
1315 doTestShift(createEllipticOrbit(), OrbitType.EQUINOCTIAL, PositionAngleType.TRUE, false,
1316 18.1, 72.0, 437.3, 1601.1, 3141.8);
1317 }
1318
1319 @Test
1320 void testShiftEquinoctialTrueWithDerivatives() {
1321 doTestShift(createEllipticOrbit(), OrbitType.EQUINOCTIAL, PositionAngleType.TRUE, true,
1322 1.14, 9.1, 140.3, 1066.7, 3306.9);
1323 }
1324
1325 @Test
1326 void testShiftEquinoctialEccentricWithoutDerivatives() {
1327 doTestShift(createEllipticOrbit(), OrbitType.EQUINOCTIAL, PositionAngleType.ECCENTRIC, false,
1328 18.1, 72.0, 437.3, 1601.1, 3141.8);
1329 }
1330
1331 @Test
1332 void testShiftEquinoctialEcentricWithDerivatives() {
1333 doTestShift(createEllipticOrbit(), OrbitType.EQUINOCTIAL, PositionAngleType.ECCENTRIC, true,
1334 1.14, 9.1, 140.3, 1066.7, 3306.9);
1335 }
1336
1337 @Test
1338 void testShiftEquinoctialMeanWithoutDerivatives() {
1339 doTestShift(createEllipticOrbit(), OrbitType.EQUINOCTIAL, PositionAngleType.MEAN, false,
1340 18.1, 72.0, 437.3, 1601.1, 3141.8);
1341 }
1342
1343 @Test
1344 void testShiftEquinoctialMeanWithDerivatives() {
1345 doTestShift(createEllipticOrbit(), OrbitType.EQUINOCTIAL, PositionAngleType.MEAN, true,
1346 1.14, 9.1, 140.3, 1066.7, 3306.9);
1347 }
1348
1349 private static void doTestShift(final CartesianOrbit orbit, final OrbitType orbitType,
1350 final PositionAngleType angleType, final boolean withDerivatives,
1351 final double error60s, final double error120s,
1352 final double error300s, final double error600s,
1353 final double error900s)
1354 {
1355
1356 Utils.setDataRoot("regular-data:atmosphere:potential/grgs-format");
1357 GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
1358 final NumericalPropagator np = createPropagator(new SpacecraftState(orbit), orbitType, angleType);
1359
1360
1361
1362 final AbsoluteDate reference = orbit.getDate().shiftedBy(60.0);
1363 final ShiftChecker checker = new ShiftChecker(withDerivatives, orbitType, angleType,
1364 error60s,
1365 error120s, error300s,
1366 error600s, error900s);
1367 np.addEventDetector(new DateDetector(reference,
1368 reference.shiftedBy( 60.0),
1369 reference.shiftedBy(120.0),
1370 reference.shiftedBy(300.0),
1371 reference.shiftedBy(600.0),
1372 reference.shiftedBy(900.0)).
1373 withMaxCheck(30.0).
1374 withThreshold(1.0e-9).
1375 withHandler(checker));
1376 np.propagate(reference.shiftedBy(1000.0));
1377 }
1378
1379 private static class ShiftChecker implements EventHandler {
1380
1381 private final boolean withDerivatives;
1382 private final OrbitType orbitType;
1383 private final PositionAngleType angleType;
1384 private final double error60s;
1385 private final double error120s;
1386 private final double error300s;
1387 private final double error600s;
1388 private final double error900s;
1389 private SpacecraftState referenceState;
1390
1391 ShiftChecker(final boolean withDerivatives, final OrbitType orbitType,
1392 final PositionAngleType angleType, final double error60s,
1393 final double error120s, final double error300s,
1394 final double error600s, final double error900s) {
1395 this.withDerivatives = withDerivatives;
1396 this.orbitType = orbitType;
1397 this.angleType = angleType;
1398 this.error60s = error60s;
1399 this.error120s = error120s;
1400 this.error300s = error300s;
1401 this.error600s = error600s;
1402 this.error900s = error900s;
1403 this.referenceState = null;
1404 }
1405
1406 @Override
1407 public Action eventOccurred(final SpacecraftState s, final EventDetector detector,
1408 final boolean increasing) {
1409 if (referenceState == null) {
1410
1411 if (withDerivatives) {
1412 referenceState = s;
1413 } else {
1414
1415 final double[] stateVector = new double[6];
1416 final Orbit o = s.getOrbit();
1417 orbitType.mapOrbitToArray(o, angleType, stateVector, null);
1418 final Orbit fixedOrbit = orbitType.mapArrayToOrbit(stateVector, null, angleType,
1419 o.getDate(), o.getMu(), o.getFrame());
1420 referenceState = new SpacecraftState(fixedOrbit, s.getAttitude(), s.getMass());
1421 }
1422 } else {
1423
1424 final double dt = s.getDate().durationFrom(referenceState.getDate());
1425 final SpacecraftState shifted = referenceState.shiftedBy(dt);
1426 final double error = Vector3D.distance(shifted.getPosition(),
1427 s.getPosition());
1428 switch ((int) FastMath.rint(dt)) {
1429 case 60 :
1430 Assertions.assertEquals(error60s, error, 0.01 * error60s);
1431 break;
1432 case 120 :
1433 Assertions.assertEquals(error120s, error, 0.01 * error120s);
1434 break;
1435 case 300 :
1436 Assertions.assertEquals(error300s, error, 0.01 * error300s);
1437 break;
1438 case 600 :
1439 Assertions.assertEquals(error600s, error, 0.01 * error600s);
1440 break;
1441 case 900 :
1442 Assertions.assertEquals(error900s, error, 0.01 * error900s);
1443 break;
1444 default :
1445
1446 Assertions.fail("no error set for dt = " + dt);
1447 break;
1448 }
1449 }
1450 return Action.CONTINUE;
1451 }
1452
1453 }
1454
1455
1456
1457
1458
1459
1460
1461 @Test
1462 void testEventAndStepHandlerDeactivationIssue449() {
1463
1464
1465 RecordAndContinue recordAndContinue = new RecordAndContinue();
1466 DateDetector dateDetector = new DateDetector(initDate.shiftedBy(10.),
1467 initDate.shiftedBy(15.),
1468 initDate.shiftedBy(20.)).
1469 withMaxCheck(1).
1470 withThreshold(1e-1).
1471 withHandler(recordAndContinue);
1472
1473 propagator.addEventDetector(dateDetector);
1474
1475 final AbsoluteDate startDate = initDate.shiftedBy(30.);
1476 final AbsoluteDate finalDate = initDate.shiftedBy(40.);
1477
1478 final DateRecorderHandler dateRecorderHandler = new DateRecorderHandler(startDate, finalDate);
1479 propagator.setStepHandler(1.0, dateRecorderHandler);
1480
1481
1482 propagator.propagate(startDate, finalDate);
1483
1484
1485
1486 Assertions.assertEquals(0, recordAndContinue.getEvents().size());
1487
1488
1489 Assertions.assertEquals(0, dateRecorderHandler.handledDatesOutOfInterval.size());
1490 }
1491
1492 @Test
1493 void testResetStateForward() {
1494 final Frame eme2000 = FramesFactory.getEME2000();
1495 final AbsoluteDate date = new AbsoluteDate(new DateComponents(2008, 6, 23),
1496 new TimeComponents(14, 0, 0),
1497 TimeScalesFactory.getUTC());
1498 final Orbit orbit = new KeplerianOrbit(8000000.0, 0.01, 0.87, 2.44, 0.21, -1.05, PositionAngleType.MEAN,
1499 eme2000,
1500 date, Constants.EIGEN5C_EARTH_MU);
1501 final NumericalPropagator propagator =
1502 new NumericalPropagator(new LutherIntegrator(300.0),
1503 new LofOffset(eme2000, LOFType.LVLH));
1504 propagator.resetInitialState(new SpacecraftState(orbit));
1505
1506
1507 final AbsoluteDate maneuverDate = date.shiftedBy(1000.0);
1508 propagator.addEventDetector(new ImpulseManeuver(new DateDetector(maneuverDate),
1509 new Vector3D(0.0, 0.0, -100.0),
1510 350.0));
1511
1512 final Vector3D initialNormal = orbit.getPVCoordinates().getMomentum();
1513 propagator.setStepHandler(60.0, state -> {
1514 final Vector3D currentNormal = state.getPVCoordinates().getMomentum();
1515 if (state.getDate().isBefore(maneuverDate)) {
1516 Assertions.assertEquals(0.000, Vector3D.angle(initialNormal, currentNormal), 1.0e-3);
1517 } else {
1518 Assertions.assertEquals(0.014, Vector3D.angle(initialNormal, currentNormal), 1.0e-3);
1519 }
1520 });
1521
1522 propagator.propagate(orbit.getDate().shiftedBy(1500.0));
1523
1524 }
1525
1526 @Test
1527 void testResetStateBackward() {
1528 final Frame eme2000 = FramesFactory.getEME2000();
1529 final AbsoluteDate date = new AbsoluteDate(new DateComponents(2008, 6, 23),
1530 new TimeComponents(14, 0, 0),
1531 TimeScalesFactory.getUTC());
1532 final Orbit orbit = new KeplerianOrbit(8000000.0, 0.01, 0.87, 2.44, 0.21, -1.05, PositionAngleType.MEAN,
1533 eme2000,
1534 date, Constants.EIGEN5C_EARTH_MU);
1535 final NumericalPropagator propagator =
1536 new NumericalPropagator(new LutherIntegrator(300.0),
1537 new LofOffset(eme2000, LOFType.LVLH));
1538 propagator.resetInitialState(new SpacecraftState(orbit));
1539
1540
1541 final AbsoluteDate maneuverDate = date.shiftedBy(-1000.0);
1542 propagator.addEventDetector(new ImpulseManeuver(new DateDetector(maneuverDate),
1543 new Vector3D(0.0, 0.0, -100.0),
1544 350.0));
1545
1546 final Vector3D initialNormal = orbit.getPVCoordinates().getMomentum();
1547 propagator.setStepHandler(60.0, state -> {
1548 final Vector3D currentNormal = state.getPVCoordinates().getMomentum();
1549 if (state.getDate().isAfter(maneuverDate)) {
1550 Assertions.assertEquals(0.000, Vector3D.angle(initialNormal, currentNormal), 1.0e-3);
1551 } else {
1552 Assertions.assertEquals(0.014, Vector3D.angle(initialNormal, currentNormal), 1.0e-3);
1553 }
1554 });
1555
1556 propagator.propagate(orbit.getDate().shiftedBy(-1500.0));
1557
1558 }
1559
1560 @Test
1561 void testAdditionalDerivatives() {
1562
1563 final Frame eme2000 = FramesFactory.getEME2000();
1564 final AbsoluteDate date = new AbsoluteDate(new DateComponents(2008, 6, 23),
1565 new TimeComponents(14, 0, 0),
1566 TimeScalesFactory.getUTC());
1567 final Orbit initialOrbit = new KeplerianOrbit(8000000.0, 0.01, 0.87, 2.44, 0.21, -1.05, PositionAngleType.MEAN,
1568 eme2000,
1569 date, Constants.EIGEN5C_EARTH_MU);
1570 NumericalPropagatorBuilder builder = new NumericalPropagatorBuilder(initialOrbit, new DormandPrince853IntegratorBuilder(0.02, 0.2, 1., 0.0001), PositionAngleType.TRUE, 10);
1571 NumericalPropagator propagator = (NumericalPropagator) builder.buildPropagator();
1572
1573 IntStream.
1574 range(0, 2).
1575 mapToObj(i -> new EmptyDerivativeProvider("test_provider_" + i, new double[] { 10 * i, 20 * i })).
1576 forEach(provider -> addDerivativeProvider(propagator, provider));
1577
1578 EphemerisGenerator generator = propagator.getEphemerisGenerator();
1579 propagator.propagate(initialOrbit.getDate().shiftedBy(600));
1580 BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
1581 final SpacecraftState finalState = ephemeris.propagate(initialOrbit.getDate().shiftedBy(300));
1582 Assertions.assertEquals(2, finalState.getAdditionalStatesValues().size());
1583 Assertions.assertEquals(2, finalState.getAdditionalState("test_provider_0").length);
1584 Assertions.assertEquals(0.0, finalState.getAdditionalState("test_provider_0")[0], 1.0e-15);
1585 Assertions.assertEquals(0.0, finalState.getAdditionalState("test_provider_0")[1], 1.0e-15);
1586 Assertions.assertEquals(2, finalState.getAdditionalState("test_provider_1").length);
1587 Assertions.assertEquals(10.0, finalState.getAdditionalState("test_provider_1")[0], 1.0e-15);
1588 Assertions.assertEquals(20.0, finalState.getAdditionalState("test_provider_1")[1], 1.0e-15);
1589 }
1590
1591 private void addDerivativeProvider(NumericalPropagator propagator, EmptyDerivativeProvider provider) {
1592 SpacecraftState initialState = propagator.getInitialState();
1593 propagator.setInitialState(initialState.addAdditionalState(provider.getName(), provider.getInitialState()));
1594 propagator.addAdditionalDerivativesProvider(provider);
1595 }
1596
1597 private static class EmptyDerivativeProvider implements AdditionalDerivativesProvider {
1598
1599 private final String name;
1600 private final double[] state;
1601
1602 public EmptyDerivativeProvider(String name, double[] state) {
1603 this.name = name;
1604 this.state = state;
1605 }
1606
1607 @Override
1608 public CombinedDerivatives combinedDerivatives(SpacecraftState s) {
1609 return new CombinedDerivatives(new double[getDimension()], null);
1610 }
1611
1612 @Override
1613 public String getName() {
1614 return name;
1615 }
1616
1617 @Override
1618 public int getDimension() {
1619 return state.length;
1620 }
1621
1622 public double[] getInitialState() {
1623 return state;
1624 }
1625 }
1626
1627 @Test
1628 void testInfinitePropagation() {
1629
1630 Utils.setDataRoot("regular-data:atmosphere:potential/grgs-format");
1631 GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
1632
1633 final NumericalPropagator propag = createPropagator(initialState, OrbitType.KEPLERIAN, PositionAngleType.TRUE);
1634
1635
1636 final double convergenceThreshold = 1e-9;
1637 propag.addEventDetector(new DateDetector(initialState.getDate().shiftedBy(60)).
1638 withMaxCheck(1e10).
1639 withThreshold(convergenceThreshold));
1640
1641
1642 final SpacecraftState finalState = propag.propagate(AbsoluteDate.FUTURE_INFINITY);
1643
1644
1645 Assertions.assertEquals(60, finalState.getDate().durationFrom(initialState.getDate()), convergenceThreshold);
1646
1647 }
1648
1649 @Test
1650 void getIntegratorNameTest() {
1651
1652 final String expectedName = "Name";
1653 final ODEIntegrator mockedIntegrator = Mockito.mock(ODEIntegrator.class);
1654 Mockito.when(mockedIntegrator.getName()).thenReturn(expectedName);
1655
1656 final NumericalPropagator numericalPropagator = new NumericalPropagator(mockedIntegrator);
1657 final String actualName = numericalPropagator.getIntegratorName();
1658
1659 Assertions.assertEquals(expectedName, actualName);
1660 }
1661
1662 @Test
1663 void testIssue1395() {
1664
1665 final Orbit initialOrbit = createEllipticOrbit();
1666 final ClassicalRungeKuttaIntegrator rungeKuttaIntegrator = new ClassicalRungeKuttaIntegrator(10);
1667 final NumericalPropagator numericalPropagator = new NumericalPropagator(rungeKuttaIntegrator);
1668 final SpacecraftState state = new SpacecraftState(initialOrbit);
1669 final String name = "test";
1670 numericalPropagator.setInitialState(state.addAdditionalState(name, 0.));
1671 numericalPropagator.addAdditionalDerivativesProvider(mockDerivativeProvider(name));
1672 numericalPropagator.addForceModel(createForceModelBasedOnAdditionalState(name));
1673
1674 final AbsoluteDate epoch = initialOrbit.getDate();
1675 final SpacecraftState propagateState = Assertions.assertDoesNotThrow(() ->
1676 numericalPropagator.propagate(epoch.shiftedBy(10.)));
1677 Assertions.assertNotEquals(epoch, propagateState.getDate());
1678 }
1679
1680 private static AdditionalDerivativesProvider mockDerivativeProvider(final String name) {
1681 final AdditionalDerivativesProvider mockedProvider = Mockito.mock(AdditionalDerivativesProvider.class);
1682 final int dimension = 1;
1683 Mockito.when(mockedProvider.getDimension()).thenReturn(dimension);
1684 Mockito.when(mockedProvider.getName()).thenReturn(name);
1685 final double[] yDot = new double[dimension];
1686 final CombinedDerivatives combinedDerivatives = new CombinedDerivatives(yDot, null);
1687 Mockito.when(mockedProvider.combinedDerivatives(Mockito.any(SpacecraftState.class)))
1688 .thenReturn(combinedDerivatives);
1689 return mockedProvider;
1690 }
1691
1692 private static ForceModel createForceModelBasedOnAdditionalState(final String name) {
1693 return new ForceModel() {
1694
1695 @Override
1696 public void init(SpacecraftState initialState, AbsoluteDate target) {
1697 ForceModel.super.init(initialState, target);
1698 initialState.getAdditionalState(name);
1699 }
1700
1701 @Override
1702 public boolean dependsOnPositionOnly() {
1703 return false;
1704 }
1705
1706 @Override
1707 public Vector3D acceleration(SpacecraftState s, double[] parameters) {
1708 return Vector3D.ZERO;
1709 }
1710
1711 @Override
1712 public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(FieldSpacecraftState<T> s, T[] parameters) {
1713 return null;
1714 }
1715
1716 @Override
1717 public List<ParameterDriver> getParametersDrivers() {
1718 return new ArrayList<>();
1719 }
1720 };
1721 }
1722
1723
1724
1725
1726 private static class DateRecorderHandler implements OrekitFixedStepHandler {
1727
1728
1729 private final AbsoluteDate startDate;
1730
1731
1732 private final AbsoluteDate finalDate;
1733
1734
1735 public final List<AbsoluteDate> handledDatesOutOfInterval;
1736
1737 DateRecorderHandler(final AbsoluteDate startDate, final AbsoluteDate finalDate) {
1738 this.startDate = startDate;
1739 this.finalDate = finalDate;
1740 this.handledDatesOutOfInterval = new ArrayList<>();
1741 }
1742
1743 @Override
1744 public void handleStep(SpacecraftState currentState)
1745 {
1746 final AbsoluteDate date = currentState.getDate();
1747 if (date.compareTo(startDate) < 0 || date.compareTo(finalDate) > 0) {
1748 handledDatesOutOfInterval.add(currentState.getDate());
1749 }
1750 }
1751 }
1752
1753
1754
1755
1756
1757
1758
1759 private static double[] recomputeFollowing(final int startIndex, List<SpacecraftState> allPoints)
1760 {
1761 SpacecraftState startState = allPoints.get(startIndex);
1762 NumericalPropagator innerPropagator = createPropagator(startState, OrbitType.CARTESIAN, PositionAngleType.TRUE);
1763 double[] errors = new double[allPoints.size() - startIndex - 1];
1764 for (int endIndex = startIndex + 1; endIndex < allPoints.size(); ++endIndex) {
1765 final TimeStampedPVCoordinates reference = allPoints.get(endIndex).getPVCoordinates();
1766 final TimeStampedPVCoordinates recomputed = innerPropagator.propagate(reference.getDate()).getPVCoordinates();
1767 errors[endIndex - startIndex - 1] = Vector3D.distance(recomputed.getPosition(), reference.getPosition());
1768 }
1769 return errors;
1770 }
1771
1772 private synchronized static NumericalPropagator createPropagator(SpacecraftState spacecraftState,
1773 OrbitType orbitType, PositionAngleType angleType)
1774 {
1775
1776 final double minStep = 0.001;
1777 final double maxStep = 120.0;
1778 final double positionTolerance = 0.1;
1779 final int degree = 20;
1780 final int order = 20;
1781 final double spacecraftArea = 1.0;
1782 final double spacecraftDragCoefficient = 2.0;
1783 final double spacecraftReflectionCoefficient = 2.0;
1784
1785
1786 final double[][] tol = NumericalPropagator.tolerances(positionTolerance, spacecraftState.getOrbit(), orbitType);
1787 final ODEIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
1788 final NumericalPropagator np = new NumericalPropagator(integrator);
1789 np.setOrbitType(orbitType);
1790 np.setPositionAngleType(angleType);
1791 np.setInitialState(spacecraftState);
1792
1793
1794 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
1795 Constants.WGS84_EARTH_FLATTENING,
1796 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
1797 final NormalizedSphericalHarmonicsProvider harmonicsGravityProvider = GravityFieldFactory.getNormalizedProvider(degree, order);
1798 np.addForceModel(new HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), harmonicsGravityProvider));
1799
1800
1801 np.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
1802 np.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
1803
1804
1805 MarshallSolarActivityFutureEstimation msafe =
1806 new MarshallSolarActivityFutureEstimation("Jan2000F10-edited-data\\.txt",
1807 MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
1808 DTM2000 atmosphere = new DTM2000(msafe, CelestialBodyFactory.getSun(), earth);
1809 np.addForceModel(new DragForce(atmosphere, new IsotropicDrag(spacecraftArea, spacecraftDragCoefficient)));
1810
1811
1812 np.addForceModel(new SolarRadiationPressure(CelestialBodyFactory.getSun(), earth,
1813 new IsotropicRadiationSingleCoefficient(spacecraftArea, spacecraftReflectionCoefficient)));
1814
1815 return np;
1816
1817 }
1818
1819 private CartesianOrbit createEllipticOrbit() {
1820 final AbsoluteDate date = new AbsoluteDate("2003-05-01T00:00:20.000", TimeScalesFactory.getUTC());
1821 final Vector3D position = new Vector3D(6896874.444705, 1956581.072644, -147476.245054);
1822 final Vector3D velocity = new Vector3D(169.816407662, -1126.783301861, -7332.745712770);
1823 final TimeStampedPVCoordinates pv = new TimeStampedPVCoordinates(date, position, velocity);
1824 final Frame frame = FramesFactory.getEME2000();
1825 final double mu = Constants.EIGEN5C_EARTH_MU;
1826 return new CartesianOrbit(pv, frame, mu);
1827 }
1828
1829 private CartesianOrbit createHyperbolicOrbit() {
1830 final AbsoluteDate date = new AbsoluteDate("2003-05-01T00:00:20.000", TimeScalesFactory.getUTC());
1831 final Vector3D position = new Vector3D(224267911.905821, 290251613.109399, 45534292.777492);
1832 final Vector3D velocity = new Vector3D(-1494.068165293, 1124.771027677, 526.915286134);
1833 final TimeStampedPVCoordinates pv = new TimeStampedPVCoordinates(date, position, velocity);
1834 final Frame frame = FramesFactory.getEME2000();
1835 final double mu = Constants.EIGEN5C_EARTH_MU;
1836 return new CartesianOrbit(pv, frame, mu);
1837 }
1838
1839 @Test
1840 void testFinalAttitudeWithForcesNeedingRates() {
1841 final ForceModel dummyForceDependingOnRates = new ForceModel() {
1842 @Override
1843 public boolean dependsOnAttitudeRate() {
1844 return true;
1845 }
1846
1847 @Override
1848 public boolean dependsOnPositionOnly() {
1849 return false;
1850 }
1851
1852 @Override
1853 public Vector3D acceleration(SpacecraftState s, double[] parameters) {
1854 return Vector3D.ZERO;
1855 }
1856
1857 @Override
1858 public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(FieldSpacecraftState<T> s, T[] parameters) {
1859 return null;
1860 }
1861
1862 @Override
1863 public List<ParameterDriver> getParametersDrivers() {
1864 return Collections.emptyList();
1865 }
1866 };
1867 final List<ForceModel> forceModels = new ArrayList<>();
1868 forceModels.add(dummyForceDependingOnRates);
1869 testTemplateFinalAttitudeWithoutForcesNeedingRates(forceModels);
1870 }
1871
1872 @Test
1873 void testFinalAttitudeWithoutForcesNeedingRates() {
1874 testTemplateFinalAttitudeWithoutForcesNeedingRates(new ArrayList<>());
1875 }
1876
1877 private void testTemplateFinalAttitudeWithoutForcesNeedingRates(final List<ForceModel> forceModels) {
1878
1879 final AttitudeProvider attitudeProvider = createAttitudeProviderWithNonZeroRates();
1880 propagator.setAttitudeProvider(attitudeProvider);
1881 for (final ForceModel forceModel : forceModels) {
1882 propagator.addForceModel(forceModel);
1883 }
1884
1885 final SpacecraftState state = propagator.propagate(propagator.getInitialState().getDate().shiftedBy(10.));
1886
1887 final Attitude attitude = state.getAttitude();
1888 Assertions.assertNotEquals(Vector3D.ZERO, attitude.getSpin());
1889 Assertions.assertNotEquals(Vector3D.ZERO, attitude.getRotationAcceleration());
1890 }
1891
1892 private AttitudeProvider createAttitudeProviderWithNonZeroRates() {
1893 return new AttitudeProvider() {
1894 @Override
1895 public Attitude getAttitude(PVCoordinatesProvider pvProv, AbsoluteDate date, Frame frame) {
1896 return createAttitudeWithNonZeroRates(date, frame);
1897 }
1898
1899 @Override
1900 public <T extends CalculusFieldElement<T>> FieldAttitude<T> getAttitude(FieldPVCoordinatesProvider<T> pvProv, FieldAbsoluteDate<T> date, Frame frame) {
1901 return null;
1902 }
1903 };
1904 }
1905
1906 private Attitude createAttitudeWithNonZeroRates(final AbsoluteDate date, final Frame frame) {
1907 final AngularCoordinates angularCoordinates = new AngularCoordinates(Rotation.IDENTITY,
1908 Vector3D.PLUS_K, Vector3D.MINUS_I);
1909 return new Attitude(date, frame, angularCoordinates);
1910 }
1911
1912 @BeforeEach
1913 public void setUp() {
1914 Utils.setDataRoot("regular-data:potential/shm-format");
1915 GravityFieldFactory.addPotentialCoefficientsReader(new SHMFormatReader("^eigen_cg03c_coef$", false));
1916 mu = GravityFieldFactory.getUnnormalizedProvider(0, 0).getMu();
1917 final Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6);
1918 final Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0);
1919 initDate = AbsoluteDate.J2000_EPOCH;
1920 final Orbit orbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
1921 FramesFactory.getEME2000(), initDate, mu);
1922 initialState = new SpacecraftState(orbit);
1923 double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit, OrbitType.EQUINOCTIAL);
1924 AdaptiveStepsizeIntegrator integrator =
1925 new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
1926 integrator.setInitialStepSize(60);
1927 propagator = new NumericalPropagator(integrator);
1928 propagator.setInitialState(initialState);
1929 }
1930
1931 @AfterEach
1932 public void tearDown() {
1933 initDate = null;
1934 initialState = null;
1935 propagator = null;
1936 }
1937
1938
1939
1940
1941
1942 private static class ForceModelAdapter implements ForceModel {
1943
1944 @Override
1945 public boolean dependsOnPositionOnly() {
1946 return false;
1947 }
1948
1949 @Override
1950 public boolean isSupported(String name) {
1951 return false;
1952 }
1953
1954 @Override
1955 public void addContribution(SpacecraftState s, TimeDerivativesEquations adder) {
1956 }
1957
1958 @Override
1959 public <T extends CalculusFieldElement<T>> void
1960 addContribution(FieldSpacecraftState<T> s,
1961 FieldTimeDerivativesEquations<T> adder) {
1962 }
1963
1964
1965 @Override
1966 public Vector3D acceleration(final SpacecraftState s, final double[] parameters)
1967 {
1968 return Vector3D.ZERO;
1969 }
1970
1971
1972 @Override
1973 public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(final FieldSpacecraftState<T> s,
1974 final T[] parameters)
1975 {
1976 return FieldVector3D.getZero(s.getDate().getField());
1977 }
1978
1979 @Override
1980 public Stream<EventDetector> getEventDetectors() {
1981 return Stream.empty();
1982 }
1983
1984 @Override
1985 public <T extends CalculusFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventDetectors(final Field<T> field) {
1986 return Stream.empty();
1987 }
1988
1989 @Override
1990 public List<ParameterDriver> getParametersDrivers() {
1991 return Collections.emptyList();
1992 }
1993 }
1994
1995 }
1996