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