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