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
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
1056 public CheckingHandler(final Action actionOnEvent) {
1057 this.actionOnEvent = actionOnEvent;
1058 this.gotHere = false;
1059 }
1060
1061 public void assertEvent(boolean expected) {
1062 Assertions.assertEquals(expected, gotHere);
1063 }
1064
1065 public Action eventOccurred(SpacecraftState s, EventDetector detector, boolean increasing) {
1066 gotHere = true;
1067 return actionOnEvent;
1068 }
1069
1070 }
1071
1072 @Test
1073 void testParallelismIssue258() throws InterruptedException, ExecutionException, FileNotFoundException {
1074
1075 Utils.setDataRoot("regular-data:atmosphere:potential/grgs-format");
1076 GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
1077 final double mu = GravityFieldFactory.getNormalizedProvider(2, 2).getMu();
1078
1079
1080 final double a = 24396159;
1081 final double e = 0.72831215;
1082 final double i = FastMath.toRadians(7);
1083 final double omega = FastMath.toRadians(180);
1084 final double raan = FastMath.toRadians(261);
1085 final double lM = 0;
1086 final Frame inertialFrame = FramesFactory.getEME2000();
1087 final TimeScale utc = TimeScalesFactory.getUTC();
1088 final AbsoluteDate initialDate = new AbsoluteDate(2003, 1, 1, 00, 00, 00.000, utc);
1089 final Orbit initialOrbit = new CartesianOrbit( new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngleType.MEAN,
1090 inertialFrame, initialDate, mu));
1091 final SpacecraftState initialState = new SpacecraftState(initialOrbit, 1000);
1092
1093
1094 final List<SpacecraftState> states = new ArrayList<SpacecraftState>();
1095 final NumericalPropagator propagator = createPropagator(initialState, OrbitType.CARTESIAN, PositionAngleType.TRUE);
1096 final double samplingStep = 10000.0;
1097 propagator.setStepHandler(samplingStep, state -> states.add(state));
1098 propagator.propagate(initialDate.shiftedBy(5 * samplingStep));
1099
1100
1101 final double[][] referenceErrors = new double[states.size() - 1][];
1102 for (int startIndex = 0; startIndex < states.size() - 1; ++startIndex) {
1103 referenceErrors[startIndex] = recomputeFollowing(startIndex, states);
1104 }
1105
1106 final Consumer<SpacecraftState> checker = point -> {
1107 try {
1108 final int startIndex = states.indexOf(point);
1109 double[] errors = recomputeFollowing(startIndex, states);
1110 for (int k = 0; k < errors.length; ++k) {
1111 Assertions.assertEquals(
1112 referenceErrors[startIndex][k], errors[k],
1113 1.0e-9,startIndex + " to " + (startIndex + k + 1));
1114 }
1115 } catch (OrekitException oe) {
1116 Assertions.fail(oe.getLocalizedMessage());
1117 }
1118 };
1119
1120
1121 states.stream().forEach(checker);
1122
1123
1124 states.parallelStream().forEach(checker);
1125
1126 }
1127
1128 @Test
1129 void testShiftKeplerianEllipticTrueWithoutDerivatives() {
1130 doTestShift(createEllipticOrbit(), OrbitType.KEPLERIAN, PositionAngleType.TRUE, false,
1131 18.1, 72.0, 437.3, 1601.1, 3141.8);
1132 }
1133
1134 @Test
1135 void testShiftKeplerianEllipticTrueWithDerivatives() {
1136 doTestShift(createEllipticOrbit(), OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
1137 1.14, 9.1, 140.3, 1066.7, 3306.9);
1138 }
1139
1140 @Test
1141 void testShiftKeplerianEllipticEccentricWithoutDerivatives() {
1142 doTestShift(createEllipticOrbit(), OrbitType.KEPLERIAN, PositionAngleType.ECCENTRIC, false,
1143 18.1, 72.0, 437.3, 1601.1, 3141.8);
1144 }
1145
1146 @Test
1147 void testShiftKeplerianEllipticEcentricWithDerivatives() {
1148 doTestShift(createEllipticOrbit(), OrbitType.KEPLERIAN, PositionAngleType.ECCENTRIC, true,
1149 1.14, 9.1, 140.3, 1066.7, 3306.9);
1150 }
1151
1152 @Test
1153 void testShiftKeplerianEllipticMeanWithoutDerivatives() {
1154 doTestShift(createEllipticOrbit(), OrbitType.KEPLERIAN, PositionAngleType.MEAN, false,
1155 18.1, 72.0, 437.3, 1601.1, 3141.8);
1156 }
1157
1158 @Test
1159 void testShiftKeplerianEllipticMeanWithDerivatives() {
1160 doTestShift(createEllipticOrbit(), OrbitType.KEPLERIAN, PositionAngleType.MEAN, true,
1161 1.14, 9.1, 140.3, 1066.7, 3306.9);
1162 }
1163
1164 @Test
1165 void testShiftKeplerianHyperbolicTrueWithoutDerivatives() {
1166 doTestShift(createHyperbolicOrbit(), OrbitType.KEPLERIAN, PositionAngleType.TRUE, false,
1167 0.484, 1.94, 12.1, 48.3, 108.5);
1168 }
1169
1170 @Test
1171 void testShiftKeplerianHyperbolicTrueWithDerivatives() {
1172 doTestShift(createHyperbolicOrbit(), OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
1173 1.38e-4, 1.10e-3, 1.72e-2, 1.37e-1, 4.62e-1);
1174 }
1175
1176 @Test
1177 void testShiftKeplerianHyperbolicEccentricWithoutDerivatives() {
1178 doTestShift(createHyperbolicOrbit(), OrbitType.KEPLERIAN, PositionAngleType.ECCENTRIC, false,
1179 0.484, 1.94, 12.1, 48.3, 108.5);
1180 }
1181
1182 @Test
1183 void testShiftKeplerianHyperbolicEcentricWithDerivatives() {
1184 doTestShift(createHyperbolicOrbit(), OrbitType.KEPLERIAN, PositionAngleType.ECCENTRIC, true,
1185 1.38e-4, 1.10e-3, 1.72e-2, 1.37e-1, 4.62e-1);
1186 }
1187
1188 @Test
1189 void testShiftKeplerianHyperbolicMeanWithoutDerivatives() {
1190 doTestShift(createHyperbolicOrbit(), OrbitType.KEPLERIAN, PositionAngleType.MEAN, false,
1191 0.484, 1.94, 12.1, 48.3, 108.5);
1192 }
1193
1194 @Test
1195 void testShiftKeplerianHyperbolicMeanWithDerivatives() {
1196 doTestShift(createHyperbolicOrbit(), OrbitType.KEPLERIAN, PositionAngleType.MEAN, true,
1197 1.38e-4, 1.10e-3, 1.72e-2, 1.37e-1, 4.62e-1);
1198 }
1199
1200 @Test
1201 void testShiftCartesianEllipticTrueWithoutDerivatives() {
1202 doTestShift(createEllipticOrbit(), OrbitType.CARTESIAN, PositionAngleType.TRUE, false,
1203 18.1, 72.0, 437.3, 1601.1, 3141.8);
1204 }
1205
1206 @Test
1207 void testShiftCartesianEllipticTrueWithDerivatives() {
1208 doTestShift(createEllipticOrbit(), OrbitType.CARTESIAN, PositionAngleType.TRUE, true,
1209 1.14, 9.1, 140.3, 1066.7, 3306.9);
1210 }
1211
1212 @Test
1213 void testShiftCartesianEllipticEccentricWithoutDerivatives() {
1214 doTestShift(createEllipticOrbit(), OrbitType.CARTESIAN, PositionAngleType.ECCENTRIC, false,
1215 18.1, 72.0, 437.3, 1601.1, 3141.8);
1216 }
1217
1218 @Test
1219 void testShiftCartesianEllipticEcentricWithDerivatives() {
1220 doTestShift(createEllipticOrbit(), OrbitType.CARTESIAN, PositionAngleType.ECCENTRIC, true,
1221 1.14, 9.1, 140.3, 1066.7, 3306.9);
1222 }
1223
1224 @Test
1225 void testShiftCartesianEllipticMeanWithoutDerivatives() {
1226 doTestShift(createEllipticOrbit(), OrbitType.CARTESIAN, PositionAngleType.MEAN, false,
1227 18.1, 72.0, 437.3, 1601.1, 3141.8);
1228 }
1229
1230 @Test
1231 void testShiftCartesianEllipticMeanWithDerivatives() {
1232 doTestShift(createEllipticOrbit(), OrbitType.CARTESIAN, PositionAngleType.MEAN, true,
1233 1.14, 9.1, 140.3, 1066.7, 3306.9);
1234 }
1235
1236 @Test
1237 void testShiftCartesianHyperbolicTrueWithoutDerivatives() {
1238 doTestShift(createHyperbolicOrbit(), OrbitType.CARTESIAN, PositionAngleType.TRUE, false,
1239 0.48, 1.93, 12.1, 48.3, 108.5);
1240 }
1241
1242 @Test
1243 void testShiftCartesianHyperbolicTrueWithDerivatives() {
1244 doTestShift(createHyperbolicOrbit(), OrbitType.CARTESIAN, PositionAngleType.TRUE, true,
1245 1.38e-4, 1.10e-3, 1.72e-2, 1.37e-1, 4.62e-1);
1246 }
1247
1248 @Test
1249 void testShiftCartesianHyperbolicEccentricWithoutDerivatives() {
1250 doTestShift(createHyperbolicOrbit(), OrbitType.CARTESIAN, PositionAngleType.ECCENTRIC, false,
1251 0.48, 1.93, 12.1, 48.3, 108.5);
1252 }
1253
1254 @Test
1255 void testShiftCartesianHyperbolicEcentricWithDerivatives() {
1256 doTestShift(createHyperbolicOrbit(), OrbitType.CARTESIAN, PositionAngleType.ECCENTRIC, true,
1257 1.38e-4, 1.10e-3, 1.72e-2, 1.37e-1, 4.62e-1);
1258 }
1259
1260 @Test
1261 void testShiftCartesianHyperbolicMeanWithoutDerivatives() {
1262 doTestShift(createHyperbolicOrbit(), OrbitType.CARTESIAN, PositionAngleType.MEAN, false,
1263 0.48, 1.93, 12.1, 48.3, 108.5);
1264 }
1265
1266 @Test
1267 void testShiftCartesianHyperbolicMeanWithDerivatives() {
1268 doTestShift(createHyperbolicOrbit(), OrbitType.CARTESIAN, PositionAngleType.MEAN, true,
1269 1.38e-4, 1.10e-3, 1.72e-2, 1.37e-1, 4.62e-1);
1270 }
1271
1272 @Test
1273 void testShiftCircularTrueWithoutDerivatives() {
1274 doTestShift(createEllipticOrbit(), OrbitType.CIRCULAR, PositionAngleType.TRUE, false,
1275 18.1, 72.0, 437.3, 1601.1, 3141.8);
1276 }
1277
1278 @Test
1279 void testShiftCircularTrueWithDerivatives() {
1280 doTestShift(createEllipticOrbit(), OrbitType.CIRCULAR, PositionAngleType.TRUE, true,
1281 1.14, 9.1, 140.3, 1066.7, 3306.9);
1282 }
1283
1284 @Test
1285 void testShiftCircularEccentricWithoutDerivatives() {
1286 doTestShift(createEllipticOrbit(), OrbitType.CIRCULAR, PositionAngleType.ECCENTRIC, false,
1287 18.1, 72.0, 437.3, 1601.1, 3141.8);
1288 }
1289
1290 @Test
1291 void testShiftCircularEcentricWithDerivatives() {
1292 doTestShift(createEllipticOrbit(), OrbitType.CIRCULAR, PositionAngleType.ECCENTRIC, true,
1293 1.14, 9.1, 140.3, 1066.7, 3306.9);
1294 }
1295
1296 @Test
1297 void testShiftCircularMeanWithoutDerivatives() {
1298 doTestShift(createEllipticOrbit(), OrbitType.CIRCULAR, PositionAngleType.MEAN, false,
1299 18.1, 72.0, 437.3, 1601.1, 3141.8);
1300 }
1301
1302 @Test
1303 void testShiftCircularMeanWithDerivatives() {
1304 doTestShift(createEllipticOrbit(), OrbitType.CIRCULAR, PositionAngleType.MEAN, true,
1305 1.14, 9.1, 140.3, 1066.7, 3306.9);
1306 }
1307
1308 @Test
1309 void testShiftEquinoctialTrueWithoutDerivatives() {
1310 doTestShift(createEllipticOrbit(), OrbitType.EQUINOCTIAL, PositionAngleType.TRUE, false,
1311 18.1, 72.0, 437.3, 1601.1, 3141.8);
1312 }
1313
1314 @Test
1315 void testShiftEquinoctialTrueWithDerivatives() {
1316 doTestShift(createEllipticOrbit(), OrbitType.EQUINOCTIAL, PositionAngleType.TRUE, true,
1317 1.14, 9.1, 140.3, 1066.7, 3306.9);
1318 }
1319
1320 @Test
1321 void testShiftEquinoctialEccentricWithoutDerivatives() {
1322 doTestShift(createEllipticOrbit(), OrbitType.EQUINOCTIAL, PositionAngleType.ECCENTRIC, false,
1323 18.1, 72.0, 437.3, 1601.1, 3141.8);
1324 }
1325
1326 @Test
1327 void testShiftEquinoctialEcentricWithDerivatives() {
1328 doTestShift(createEllipticOrbit(), OrbitType.EQUINOCTIAL, PositionAngleType.ECCENTRIC, true,
1329 1.14, 9.1, 140.3, 1066.7, 3306.9);
1330 }
1331
1332 @Test
1333 void testShiftEquinoctialMeanWithoutDerivatives() {
1334 doTestShift(createEllipticOrbit(), OrbitType.EQUINOCTIAL, PositionAngleType.MEAN, false,
1335 18.1, 72.0, 437.3, 1601.1, 3141.8);
1336 }
1337
1338 @Test
1339 void testShiftEquinoctialMeanWithDerivatives() {
1340 doTestShift(createEllipticOrbit(), OrbitType.EQUINOCTIAL, PositionAngleType.MEAN, true,
1341 1.14, 9.1, 140.3, 1066.7, 3306.9);
1342 }
1343
1344 private static void doTestShift(final CartesianOrbit orbit, final OrbitType orbitType,
1345 final PositionAngleType angleType, final boolean withDerivatives,
1346 final double error60s, final double error120s,
1347 final double error300s, final double error600s,
1348 final double error900s)
1349 {
1350
1351 Utils.setDataRoot("regular-data:atmosphere:potential/grgs-format");
1352 GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
1353 final NumericalPropagator np = createPropagator(new SpacecraftState(orbit), orbitType, angleType);
1354
1355
1356
1357 final AbsoluteDate reference = orbit.getDate().shiftedBy(60.0);
1358 final ShiftChecker checker = new ShiftChecker(withDerivatives, orbitType, angleType,
1359 error60s,
1360 error120s, error300s,
1361 error600s, error900s);
1362 np.addEventDetector(new DateDetector(reference,
1363 reference.shiftedBy( 60.0),
1364 reference.shiftedBy(120.0),
1365 reference.shiftedBy(300.0),
1366 reference.shiftedBy(600.0),
1367 reference.shiftedBy(900.0)).
1368 withMaxCheck(30.0).
1369 withThreshold(1.0e-9).
1370 withHandler(checker));
1371 np.propagate(reference.shiftedBy(1000.0));
1372 }
1373
1374 private static class ShiftChecker implements EventHandler {
1375
1376 private final boolean withDerivatives;
1377 private final OrbitType orbitType;
1378 private final PositionAngleType angleType;
1379 private final double error60s;
1380 private final double error120s;
1381 private final double error300s;
1382 private final double error600s;
1383 private final double error900s;
1384 private SpacecraftState referenceState;
1385
1386 ShiftChecker(final boolean withDerivatives, final OrbitType orbitType,
1387 final PositionAngleType angleType, final double error60s,
1388 final double error120s, final double error300s,
1389 final double error600s, final double error900s) {
1390 this.withDerivatives = withDerivatives;
1391 this.orbitType = orbitType;
1392 this.angleType = angleType;
1393 this.error60s = error60s;
1394 this.error120s = error120s;
1395 this.error300s = error300s;
1396 this.error600s = error600s;
1397 this.error900s = error900s;
1398 this.referenceState = null;
1399 }
1400
1401 @Override
1402 public Action eventOccurred(final SpacecraftState s, final EventDetector detector,
1403 final boolean increasing) {
1404 if (referenceState == null) {
1405
1406 if (withDerivatives) {
1407 referenceState = s;
1408 } else {
1409
1410 final double[] stateVector = new double[6];
1411 final Orbit o = s.getOrbit();
1412 orbitType.mapOrbitToArray(o, angleType, stateVector, null);
1413 final Orbit fixedOrbit = orbitType.mapArrayToOrbit(stateVector, null, angleType,
1414 o.getDate(), o.getMu(), o.getFrame());
1415 referenceState = new SpacecraftState(fixedOrbit, s.getAttitude(), s.getMass());
1416 }
1417 } else {
1418
1419 final double dt = s.getDate().durationFrom(referenceState.getDate());
1420 final SpacecraftState shifted = referenceState.shiftedBy(dt);
1421 final double error = Vector3D.distance(shifted.getPosition(),
1422 s.getPosition());
1423 switch ((int) FastMath.rint(dt)) {
1424 case 60 :
1425 Assertions.assertEquals(error60s, error, 0.01 * error60s);
1426 break;
1427 case 120 :
1428 Assertions.assertEquals(error120s, error, 0.01 * error120s);
1429 break;
1430 case 300 :
1431 Assertions.assertEquals(error300s, error, 0.01 * error300s);
1432 break;
1433 case 600 :
1434 Assertions.assertEquals(error600s, error, 0.01 * error600s);
1435 break;
1436 case 900 :
1437 Assertions.assertEquals(error900s, error, 0.01 * error900s);
1438 break;
1439 default :
1440
1441 Assertions.fail("no error set for dt = " + dt);
1442 break;
1443 }
1444 }
1445 return Action.CONTINUE;
1446 }
1447
1448 }
1449
1450
1451
1452
1453
1454
1455
1456 @Test
1457 void testEventAndStepHandlerDeactivationIssue449() {
1458
1459
1460 RecordAndContinue recordAndContinue = new RecordAndContinue();
1461 DateDetector dateDetector = new DateDetector(initDate.shiftedBy(10.),
1462 initDate.shiftedBy(15.),
1463 initDate.shiftedBy(20.)).
1464 withMaxCheck(1).
1465 withThreshold(1e-1).
1466 withHandler(recordAndContinue);
1467
1468 propagator.addEventDetector(dateDetector);
1469
1470 final AbsoluteDate startDate = initDate.shiftedBy(30.);
1471 final AbsoluteDate finalDate = initDate.shiftedBy(40.);
1472
1473 final DateRecorderHandler dateRecorderHandler = new DateRecorderHandler(startDate, finalDate);
1474 propagator.setStepHandler(1.0, dateRecorderHandler);
1475
1476
1477 propagator.propagate(startDate, finalDate);
1478
1479
1480
1481 Assertions.assertEquals(0, recordAndContinue.getEvents().size());
1482
1483
1484 Assertions.assertEquals(0, dateRecorderHandler.handledDatesOutOfInterval.size());
1485 }
1486
1487 @Test
1488 void testResetStateForward() {
1489 final Frame eme2000 = FramesFactory.getEME2000();
1490 final AbsoluteDate date = new AbsoluteDate(new DateComponents(2008, 6, 23),
1491 new TimeComponents(14, 0, 0),
1492 TimeScalesFactory.getUTC());
1493 final Orbit orbit = new KeplerianOrbit(8000000.0, 0.01, 0.87, 2.44, 0.21, -1.05, PositionAngleType.MEAN,
1494 eme2000,
1495 date, Constants.EIGEN5C_EARTH_MU);
1496 final NumericalPropagator propagator =
1497 new NumericalPropagator(new LutherIntegrator(300.0),
1498 new LofOffset(eme2000, LOFType.LVLH));
1499 propagator.resetInitialState(new SpacecraftState(orbit));
1500
1501
1502 final AbsoluteDate maneuverDate = date.shiftedBy(1000.0);
1503 propagator.addEventDetector(new ImpulseManeuver(new DateDetector(maneuverDate),
1504 new Vector3D(0.0, 0.0, -100.0),
1505 350.0));
1506
1507 final Vector3D initialNormal = orbit.getPVCoordinates().getMomentum();
1508 propagator.setStepHandler(60.0, state -> {
1509 final Vector3D currentNormal = state.getPVCoordinates().getMomentum();
1510 if (state.getDate().isBefore(maneuverDate)) {
1511 Assertions.assertEquals(0.000, Vector3D.angle(initialNormal, currentNormal), 1.0e-3);
1512 } else {
1513 Assertions.assertEquals(0.014, Vector3D.angle(initialNormal, currentNormal), 1.0e-3);
1514 }
1515 });
1516
1517 propagator.propagate(orbit.getDate().shiftedBy(1500.0));
1518
1519 }
1520
1521 @Test
1522 void testResetStateBackward() {
1523 final Frame eme2000 = FramesFactory.getEME2000();
1524 final AbsoluteDate date = new AbsoluteDate(new DateComponents(2008, 6, 23),
1525 new TimeComponents(14, 0, 0),
1526 TimeScalesFactory.getUTC());
1527 final Orbit orbit = new KeplerianOrbit(8000000.0, 0.01, 0.87, 2.44, 0.21, -1.05, PositionAngleType.MEAN,
1528 eme2000,
1529 date, Constants.EIGEN5C_EARTH_MU);
1530 final NumericalPropagator propagator =
1531 new NumericalPropagator(new LutherIntegrator(300.0),
1532 new LofOffset(eme2000, LOFType.LVLH));
1533 propagator.resetInitialState(new SpacecraftState(orbit));
1534
1535
1536 final AbsoluteDate maneuverDate = date.shiftedBy(-1000.0);
1537 propagator.addEventDetector(new ImpulseManeuver(new DateDetector(maneuverDate),
1538 new Vector3D(0.0, 0.0, -100.0),
1539 350.0));
1540
1541 final Vector3D initialNormal = orbit.getPVCoordinates().getMomentum();
1542 propagator.setStepHandler(60.0, state -> {
1543 final Vector3D currentNormal = state.getPVCoordinates().getMomentum();
1544 if (state.getDate().isAfter(maneuverDate)) {
1545 Assertions.assertEquals(0.000, Vector3D.angle(initialNormal, currentNormal), 1.0e-3);
1546 } else {
1547 Assertions.assertEquals(0.014, Vector3D.angle(initialNormal, currentNormal), 1.0e-3);
1548 }
1549 });
1550
1551 propagator.propagate(orbit.getDate().shiftedBy(-1500.0));
1552
1553 }
1554
1555 @Test
1556 void testAdditionalDerivatives() {
1557
1558 final Frame eme2000 = FramesFactory.getEME2000();
1559 final AbsoluteDate date = new AbsoluteDate(new DateComponents(2008, 6, 23),
1560 new TimeComponents(14, 0, 0),
1561 TimeScalesFactory.getUTC());
1562 final Orbit initialOrbit = new KeplerianOrbit(8000000.0, 0.01, 0.87, 2.44, 0.21, -1.05, PositionAngleType.MEAN,
1563 eme2000,
1564 date, Constants.EIGEN5C_EARTH_MU);
1565 NumericalPropagatorBuilder builder = new NumericalPropagatorBuilder(initialOrbit, new DormandPrince853IntegratorBuilder(0.02, 0.2, 1.), PositionAngleType.TRUE, 10);
1566 NumericalPropagator propagator = (NumericalPropagator) builder.buildPropagator();
1567
1568 IntStream.
1569 range(0, 2).
1570 mapToObj(i -> new EmptyDerivativeProvider("test_provider_" + i, new double[] { 10 * i, 20 * i })).
1571 forEach(provider -> addDerivativeProvider(propagator, provider));
1572
1573 EphemerisGenerator generator = propagator.getEphemerisGenerator();
1574 propagator.propagate(initialOrbit.getDate().shiftedBy(600));
1575 BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
1576 final SpacecraftState finalState = ephemeris.propagate(initialOrbit.getDate().shiftedBy(300));
1577 Assertions.assertEquals(2, finalState.getAdditionalStatesValues().size());
1578 Assertions.assertEquals(2, finalState.getAdditionalState("test_provider_0").length);
1579 Assertions.assertEquals(0.0, finalState.getAdditionalState("test_provider_0")[0], 1.0e-15);
1580 Assertions.assertEquals(0.0, finalState.getAdditionalState("test_provider_0")[1], 1.0e-15);
1581 Assertions.assertEquals(2, finalState.getAdditionalState("test_provider_1").length);
1582 Assertions.assertEquals(10.0, finalState.getAdditionalState("test_provider_1")[0], 1.0e-15);
1583 Assertions.assertEquals(20.0, finalState.getAdditionalState("test_provider_1")[1], 1.0e-15);
1584 }
1585
1586 private void addDerivativeProvider(NumericalPropagator propagator, EmptyDerivativeProvider provider) {
1587 SpacecraftState initialState = propagator.getInitialState();
1588 propagator.setInitialState(initialState.addAdditionalState(provider.getName(), provider.getInitialState()));
1589 propagator.addAdditionalDerivativesProvider(provider);
1590 }
1591
1592 private static class EmptyDerivativeProvider implements AdditionalDerivativesProvider {
1593
1594 private final String name;
1595 private final double[] state;
1596
1597 public EmptyDerivativeProvider(String name, double[] state) {
1598 this.name = name;
1599 this.state = state;
1600 }
1601
1602 @Override
1603 public CombinedDerivatives combinedDerivatives(SpacecraftState s) {
1604 return new CombinedDerivatives(new double[getDimension()], null);
1605 }
1606
1607 @Override
1608 public String getName() {
1609 return name;
1610 }
1611
1612 @Override
1613 public int getDimension() {
1614 return state.length;
1615 }
1616
1617 public double[] getInitialState() {
1618 return state;
1619 }
1620 }
1621
1622 @Test
1623 void testInfinitePropagation() {
1624
1625 Utils.setDataRoot("regular-data:atmosphere:potential/grgs-format");
1626 GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
1627
1628 final NumericalPropagator propag = createPropagator(initialState, OrbitType.KEPLERIAN, PositionAngleType.TRUE);
1629
1630
1631 final double convergenceThreshold = 1e-9;
1632 propag.addEventDetector(new DateDetector(initialState.getDate().shiftedBy(60)).
1633 withMaxCheck(1e10).
1634 withThreshold(convergenceThreshold));
1635
1636
1637 final SpacecraftState finalState = propag.propagate(AbsoluteDate.FUTURE_INFINITY);
1638
1639
1640 Assertions.assertEquals(60, finalState.getDate().durationFrom(initialState.getDate()), convergenceThreshold);
1641
1642 }
1643
1644 @Test
1645 void getIntegratorNameTest() {
1646
1647 final String expectedName = "Name";
1648 final ODEIntegrator mockedIntegrator = Mockito.mock(ODEIntegrator.class);
1649 Mockito.when(mockedIntegrator.getName()).thenReturn(expectedName);
1650
1651 final NumericalPropagator numericalPropagator = new NumericalPropagator(mockedIntegrator);
1652 final String actualName = numericalPropagator.getIntegratorName();
1653
1654 Assertions.assertEquals(expectedName, actualName);
1655 }
1656
1657 @Test
1658 void testIssue1395() {
1659
1660 final Orbit initialOrbit = createEllipticOrbit();
1661 final ClassicalRungeKuttaIntegrator rungeKuttaIntegrator = new ClassicalRungeKuttaIntegrator(10);
1662 final NumericalPropagator numericalPropagator = new NumericalPropagator(rungeKuttaIntegrator);
1663 final SpacecraftState state = new SpacecraftState(initialOrbit);
1664 final String name = "test";
1665 numericalPropagator.setInitialState(state.addAdditionalState(name, 0.));
1666 numericalPropagator.addAdditionalDerivativesProvider(mockDerivativeProvider(name));
1667 numericalPropagator.addForceModel(createForceModelBasedOnAdditionalState(name));
1668
1669 final AbsoluteDate epoch = initialOrbit.getDate();
1670 final SpacecraftState propagateState = Assertions.assertDoesNotThrow(() ->
1671 numericalPropagator.propagate(epoch.shiftedBy(10.)));
1672 Assertions.assertNotEquals(epoch, propagateState.getDate());
1673 }
1674
1675 private static AdditionalDerivativesProvider mockDerivativeProvider(final String name) {
1676 final AdditionalDerivativesProvider mockedProvider = Mockito.mock(AdditionalDerivativesProvider.class);
1677 final int dimension = 1;
1678 Mockito.when(mockedProvider.getDimension()).thenReturn(dimension);
1679 Mockito.when(mockedProvider.getName()).thenReturn(name);
1680 final double[] yDot = new double[dimension];
1681 final CombinedDerivatives combinedDerivatives = new CombinedDerivatives(yDot, null);
1682 Mockito.when(mockedProvider.combinedDerivatives(Mockito.any(SpacecraftState.class)))
1683 .thenReturn(combinedDerivatives);
1684 return mockedProvider;
1685 }
1686
1687 private static ForceModel createForceModelBasedOnAdditionalState(final String name) {
1688 return new ForceModel() {
1689
1690 @Override
1691 public void init(SpacecraftState initialState, AbsoluteDate target) {
1692 ForceModel.super.init(initialState, target);
1693 initialState.getAdditionalState(name);
1694 }
1695
1696 @Override
1697 public boolean dependsOnPositionOnly() {
1698 return false;
1699 }
1700
1701 @Override
1702 public Vector3D acceleration(SpacecraftState s, double[] parameters) {
1703 return Vector3D.ZERO;
1704 }
1705
1706 @Override
1707 public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(FieldSpacecraftState<T> s, T[] parameters) {
1708 return null;
1709 }
1710
1711 @Override
1712 public List<ParameterDriver> getParametersDrivers() {
1713 return new ArrayList<>();
1714 }
1715 };
1716 }
1717
1718
1719
1720
1721 private static class DateRecorderHandler implements OrekitFixedStepHandler {
1722
1723
1724 private final AbsoluteDate startDate;
1725
1726
1727 private final AbsoluteDate finalDate;
1728
1729
1730 public final List<AbsoluteDate> handledDatesOutOfInterval;
1731
1732 DateRecorderHandler(final AbsoluteDate startDate, final AbsoluteDate finalDate) {
1733 this.startDate = startDate;
1734 this.finalDate = finalDate;
1735 this.handledDatesOutOfInterval = new ArrayList<>();
1736 }
1737
1738 @Override
1739 public void handleStep(SpacecraftState currentState)
1740 {
1741 final AbsoluteDate date = currentState.getDate();
1742 if (date.compareTo(startDate) < 0 || date.compareTo(finalDate) > 0) {
1743 handledDatesOutOfInterval.add(currentState.getDate());
1744 }
1745 }
1746 }
1747
1748
1749
1750
1751
1752
1753
1754 private static double[] recomputeFollowing(final int startIndex, List<SpacecraftState> allPoints)
1755 {
1756 SpacecraftState startState = allPoints.get(startIndex);
1757 NumericalPropagator innerPropagator = createPropagator(startState, OrbitType.CARTESIAN, PositionAngleType.TRUE);
1758 double[] errors = new double[allPoints.size() - startIndex - 1];
1759 for (int endIndex = startIndex + 1; endIndex < allPoints.size(); ++endIndex) {
1760 final TimeStampedPVCoordinates reference = allPoints.get(endIndex).getPVCoordinates();
1761 final TimeStampedPVCoordinates recomputed = innerPropagator.propagate(reference.getDate()).getPVCoordinates();
1762 errors[endIndex - startIndex - 1] = Vector3D.distance(recomputed.getPosition(), reference.getPosition());
1763 }
1764 return errors;
1765 }
1766
1767 private synchronized static NumericalPropagator createPropagator(SpacecraftState spacecraftState,
1768 OrbitType orbitType, PositionAngleType angleType)
1769 {
1770
1771 final double minStep = 0.001;
1772 final double maxStep = 120.0;
1773 final double positionTolerance = 0.1;
1774 final int degree = 20;
1775 final int order = 20;
1776 final double spacecraftArea = 1.0;
1777 final double spacecraftDragCoefficient = 2.0;
1778 final double spacecraftReflectionCoefficient = 2.0;
1779
1780
1781 final double[][] tol = NumericalPropagator.tolerances(positionTolerance, spacecraftState.getOrbit(), orbitType);
1782 final ODEIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
1783 final NumericalPropagator np = new NumericalPropagator(integrator);
1784 np.setOrbitType(orbitType);
1785 np.setPositionAngleType(angleType);
1786 np.setInitialState(spacecraftState);
1787
1788
1789 final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
1790 Constants.WGS84_EARTH_FLATTENING,
1791 FramesFactory.getITRF(IERSConventions.IERS_2010, true));
1792 final NormalizedSphericalHarmonicsProvider harmonicsGravityProvider = GravityFieldFactory.getNormalizedProvider(degree, order);
1793 np.addForceModel(new HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), harmonicsGravityProvider));
1794
1795
1796 np.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
1797 np.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
1798
1799
1800 MarshallSolarActivityFutureEstimation msafe =
1801 new MarshallSolarActivityFutureEstimation("Jan2000F10-edited-data\\.txt",
1802 MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
1803 DTM2000 atmosphere = new DTM2000(msafe, CelestialBodyFactory.getSun(), earth);
1804 np.addForceModel(new DragForce(atmosphere, new IsotropicDrag(spacecraftArea, spacecraftDragCoefficient)));
1805
1806
1807 np.addForceModel(new SolarRadiationPressure(CelestialBodyFactory.getSun(), earth,
1808 new IsotropicRadiationSingleCoefficient(spacecraftArea, spacecraftReflectionCoefficient)));
1809
1810 return np;
1811
1812 }
1813
1814 private CartesianOrbit createEllipticOrbit() {
1815 final AbsoluteDate date = new AbsoluteDate("2003-05-01T00:00:20.000", TimeScalesFactory.getUTC());
1816 final Vector3D position = new Vector3D(6896874.444705, 1956581.072644, -147476.245054);
1817 final Vector3D velocity = new Vector3D(169.816407662, -1126.783301861, -7332.745712770);
1818 final TimeStampedPVCoordinates pv = new TimeStampedPVCoordinates(date, position, velocity);
1819 final Frame frame = FramesFactory.getEME2000();
1820 final double mu = Constants.EIGEN5C_EARTH_MU;
1821 return new CartesianOrbit(pv, frame, mu);
1822 }
1823
1824 private CartesianOrbit createHyperbolicOrbit() {
1825 final AbsoluteDate date = new AbsoluteDate("2003-05-01T00:00:20.000", TimeScalesFactory.getUTC());
1826 final Vector3D position = new Vector3D(224267911.905821, 290251613.109399, 45534292.777492);
1827 final Vector3D velocity = new Vector3D(-1494.068165293, 1124.771027677, 526.915286134);
1828 final TimeStampedPVCoordinates pv = new TimeStampedPVCoordinates(date, position, velocity);
1829 final Frame frame = FramesFactory.getEME2000();
1830 final double mu = Constants.EIGEN5C_EARTH_MU;
1831 return new CartesianOrbit(pv, frame, mu);
1832 }
1833
1834 @Test
1835 void testFinalAttitudeWithForcesNeedingRates() {
1836 final ForceModel dummyForceDependingOnRates = new ForceModel() {
1837 @Override
1838 public boolean dependsOnAttitudeRate() {
1839 return true;
1840 }
1841
1842 @Override
1843 public boolean dependsOnPositionOnly() {
1844 return false;
1845 }
1846
1847 @Override
1848 public Vector3D acceleration(SpacecraftState s, double[] parameters) {
1849 return Vector3D.ZERO;
1850 }
1851
1852 @Override
1853 public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(FieldSpacecraftState<T> s, T[] parameters) {
1854 return null;
1855 }
1856
1857 @Override
1858 public List<ParameterDriver> getParametersDrivers() {
1859 return Collections.emptyList();
1860 }
1861 };
1862 final List<ForceModel> forceModels = new ArrayList<>();
1863 forceModels.add(dummyForceDependingOnRates);
1864 testTemplateFinalAttitudeWithoutForcesNeedingRates(forceModels);
1865 }
1866
1867 @Test
1868 void testFinalAttitudeWithoutForcesNeedingRates() {
1869 testTemplateFinalAttitudeWithoutForcesNeedingRates(new ArrayList<>());
1870 }
1871
1872 private void testTemplateFinalAttitudeWithoutForcesNeedingRates(final List<ForceModel> forceModels) {
1873
1874 final AttitudeProvider attitudeProvider = createAttitudeProviderWithNonZeroRates();
1875 propagator.setAttitudeProvider(attitudeProvider);
1876 for (final ForceModel forceModel : forceModels) {
1877 propagator.addForceModel(forceModel);
1878 }
1879
1880 final SpacecraftState state = propagator.propagate(propagator.getInitialState().getDate().shiftedBy(10.));
1881
1882 final Attitude attitude = state.getAttitude();
1883 Assertions.assertNotEquals(Vector3D.ZERO, attitude.getSpin());
1884 Assertions.assertNotEquals(Vector3D.ZERO, attitude.getRotationAcceleration());
1885 }
1886
1887 private AttitudeProvider createAttitudeProviderWithNonZeroRates() {
1888 return new AttitudeProvider() {
1889 @Override
1890 public Attitude getAttitude(PVCoordinatesProvider pvProv, AbsoluteDate date, Frame frame) {
1891 return createAttitudeWithNonZeroRates(date, frame);
1892 }
1893
1894 @Override
1895 public <T extends CalculusFieldElement<T>> FieldAttitude<T> getAttitude(FieldPVCoordinatesProvider<T> pvProv, FieldAbsoluteDate<T> date, Frame frame) {
1896 return null;
1897 }
1898 };
1899 }
1900
1901 private Attitude createAttitudeWithNonZeroRates(final AbsoluteDate date, final Frame frame) {
1902 final AngularCoordinates angularCoordinates = new AngularCoordinates(Rotation.IDENTITY,
1903 Vector3D.PLUS_K, Vector3D.MINUS_I);
1904 return new Attitude(date, frame, angularCoordinates);
1905 }
1906
1907 @BeforeEach
1908 public void setUp() {
1909 Utils.setDataRoot("regular-data:potential/shm-format");
1910 GravityFieldFactory.addPotentialCoefficientsReader(new SHMFormatReader("^eigen_cg03c_coef$", false));
1911 mu = GravityFieldFactory.getUnnormalizedProvider(0, 0).getMu();
1912 final Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6);
1913 final Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0);
1914 initDate = AbsoluteDate.J2000_EPOCH;
1915 final Orbit orbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
1916 FramesFactory.getEME2000(), initDate, mu);
1917 initialState = new SpacecraftState(orbit);
1918 double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit, OrbitType.EQUINOCTIAL);
1919 AdaptiveStepsizeIntegrator integrator =
1920 new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
1921 integrator.setInitialStepSize(60);
1922 propagator = new NumericalPropagator(integrator);
1923 propagator.setInitialState(initialState);
1924 }
1925
1926 @AfterEach
1927 public void tearDown() {
1928 initDate = null;
1929 initialState = null;
1930 propagator = null;
1931 }
1932
1933
1934
1935
1936
1937 private static class ForceModelAdapter implements ForceModel {
1938
1939 @Override
1940 public boolean dependsOnPositionOnly() {
1941 return false;
1942 }
1943
1944 @Override
1945 public boolean isSupported(String name) {
1946 return false;
1947 }
1948
1949 @Override
1950 public void addContribution(SpacecraftState s, TimeDerivativesEquations adder) {
1951 }
1952
1953 @Override
1954 public <T extends CalculusFieldElement<T>> void
1955 addContribution(FieldSpacecraftState<T> s,
1956 FieldTimeDerivativesEquations<T> adder) {
1957 }
1958
1959
1960 @Override
1961 public Vector3D acceleration(final SpacecraftState s, final double[] parameters)
1962 {
1963 return Vector3D.ZERO;
1964 }
1965
1966
1967 @Override
1968 public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(final FieldSpacecraftState<T> s,
1969 final T[] parameters)
1970 {
1971 return FieldVector3D.getZero(s.getDate().getField());
1972 }
1973
1974 @Override
1975 public Stream<EventDetector> getEventDetectors() {
1976 return Stream.empty();
1977 }
1978
1979 @Override
1980 public <T extends CalculusFieldElement<T>> Stream<FieldEventDetector<T>> getFieldEventDetectors(final Field<T> field) {
1981 return Stream.empty();
1982 }
1983
1984 @Override
1985 public List<ParameterDriver> getParametersDrivers() {
1986 return Collections.emptyList();
1987 }
1988 }
1989
1990 }
1991