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