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