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