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