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