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