1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.orbits;
18
19 import java.util.ArrayList;
20 import java.util.List;
21
22 import org.hipparchus.Field;
23 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
24 import org.hipparchus.geometry.euclidean.threed.Vector3D;
25 import org.hipparchus.util.Binary64;
26 import org.hipparchus.util.Binary64Field;
27 import org.hipparchus.util.FastMath;
28 import org.hipparchus.util.MathArrays;
29 import org.junit.jupiter.api.Assertions;
30 import org.junit.jupiter.api.BeforeEach;
31 import org.junit.jupiter.api.DisplayName;
32 import org.junit.jupiter.api.Test;
33 import org.orekit.TestUtils;
34 import org.orekit.Utils;
35 import org.orekit.errors.OrekitIllegalArgumentException;
36 import org.orekit.frames.Frame;
37 import org.orekit.frames.FramesFactory;
38 import org.orekit.propagation.analytical.FieldEcksteinHechlerPropagator;
39 import org.orekit.time.AbstractTimeInterpolator;
40 import org.orekit.time.FieldAbsoluteDate;
41 import org.orekit.time.FieldTimeInterpolator;
42 import org.orekit.utils.CartesianDerivativesFilter;
43 import org.orekit.utils.FieldPVCoordinates;
44
45 class FieldOrbitHermiteInterpolatorTest {
46
47 private final Field<Binary64> field = Binary64Field.getInstance();
48
49 @BeforeEach
50 public void setUp() {
51 Utils.setDataRoot("regular-data");
52 }
53
54 @Test
55 public void testCartesianInterpolationWithDerivatives() {
56
57 doTestCartesianInterpolation(true, CartesianDerivativesFilter.USE_P,
58 394, 0.1968, 3.21, 0.021630,
59 2474, 2707.6418, 6.6, 26.28);
60
61
62 doTestCartesianInterpolation(true, CartesianDerivativesFilter.USE_PV,
63 394, 1.085E-8, 3.21, 2.41518E-10,
64 2474, 0.06249, 6.6, 0.001170);
65
66
67
68 doTestCartesianInterpolation(true, CartesianDerivativesFilter.USE_PVA,
69 394, 1.92e-8, 3.21, 1.15e-9,
70 2474, 5227, 6.55, 142);
71 }
72
73 @Test
74 public void testCartesianInterpolationWithoutDerivatives() {
75
76 doTestCartesianInterpolation(false, CartesianDerivativesFilter.USE_P,
77 394, 0.1968, 3.21, 0.02163,
78 2474, 2707.6419, 6.55, 26.2826);
79
80
81 doTestCartesianInterpolation(false, CartesianDerivativesFilter.USE_PV,
82 394, 1.0857E-8, 3.21, 2.4151E-10,
83 2474, 0.06249, 6.55, 0.001170);
84
85
86
87 doTestCartesianInterpolation(false, CartesianDerivativesFilter.USE_PVA,
88 394, 2.61, 3.21, 0.154,
89 2474, 2.28e12, 6.55, 6.22e10);
90 }
91
92 @Test
93 public void testCircularInterpolationWithDerivatives() {
94 doTestCircularInterpolation(true,
95 397, 2.53e-8,
96 610, 4.72e-6,
97 4870, 110);
98 }
99
100 @Test
101 public void testCircularInterpolationWithoutDerivatives() {
102 doTestCircularInterpolation(false,
103 397, 0.0372,
104 610.0, 1.23,
105 4870, 8869);
106 }
107
108 @Test
109 public void testEquinoctialInterpolationWithDerivatives() {
110 doTestEquinoctialInterpolation(true,
111 397, 1.28e-8,
112 610, 4.35e-6,
113 4870, 114);
114 }
115
116 @Test
117 public void testEquinoctialInterpolationWithoutDerivatives() {
118 doTestEquinoctialInterpolation(false,
119 397, 0.0372,
120 610.0, 1.23,
121 4879, 8871);
122 }
123
124 @Test
125 public void testKeplerianInterpolationWithDerivatives() {
126 doTestKeplerianInterpolation(true,
127 397, 4.01, 4.75e-4, 1.28e-7,
128 2159, 1.05e7, 1.19e-3, 0.773);
129 }
130
131 @Test
132 public void testKeplerianInterpolationWithoutDerivatives() {
133 doTestKeplerianInterpolation(false,
134 397, 62.0, 4.75e-4, 2.87e-6,
135 2159, 79365, 1.19e-3, 3.89e-3);
136 }
137
138 private void doTestCartesianInterpolation(boolean useDerivatives, CartesianDerivativesFilter pvaFilter,
139 double shiftPositionErrorWithin, double interpolationPositionErrorWithin,
140 double shiftVelocityErrorWithin, double interpolationVelocityErrorWithin,
141 double shiftPositionErrorFarPast, double interpolationPositionErrorFarPast,
142 double shiftVelocityErrorFarPast, double interpolationVelocityErrorFarPast) {
143
144 Binary64 zero = field.getZero();
145 final Binary64 ehMu = zero.add(3.9860047e14);
146
147 final double ae = 6.378137e6;
148 final double c20 = -1.08263e-3;
149 final double c30 = 2.54e-6;
150 final double c40 = 1.62e-6;
151 final double c50 = 2.3e-7;
152 final double c60 = -5.5e-7;
153
154 final FieldAbsoluteDate<Binary64> date = FieldAbsoluteDate.getJ2000Epoch(field).shiftedBy(584.);
155 final FieldVector3D<Binary64> position =
156 new FieldVector3D<>(zero.add(3220103.), zero.add(69623.), zero.add(6449822.));
157 final FieldVector3D<Binary64> velocity =
158 new FieldVector3D<>(zero.add(6414.7), zero.add(-2006.), zero.add(-3180.));
159 final FieldCartesianOrbit<Binary64> initialOrbit =
160 new FieldCartesianOrbit<>(new FieldPVCoordinates<>(position, velocity),
161 FramesFactory.getEME2000(), date, ehMu);
162
163 FieldEcksteinHechlerPropagator<Binary64> propagator =
164 new FieldEcksteinHechlerPropagator<>(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
165
166
167 List<FieldOrbit<Binary64>> sample = new ArrayList<>();
168 for (Binary64 dt = zero; dt.getReal() < 251.0; dt = dt.add(60.0)) {
169 FieldOrbit<Binary64> orbit = propagator.propagate(date.shiftedBy(dt)).getOrbit();
170 if (!useDerivatives) {
171
172 Binary64[] stateVector = MathArrays.buildArray(field, 6);
173 orbit.getType().mapOrbitToArray(orbit, PositionAngleType.TRUE, stateVector, null);
174 orbit = orbit.getType().mapArrayToOrbit(stateVector, null, PositionAngleType.TRUE,
175 orbit.getDate(), orbit.getMu(), orbit.getFrame());
176 }
177 sample.add(orbit);
178 }
179
180
181 final FieldTimeInterpolator<FieldOrbit<Binary64>, Binary64> interpolator =
182 new FieldOrbitHermiteInterpolator<>(sample.size(), 409, initialOrbit.getFrame(), pvaFilter);
183
184
185
186
187
188 double maxShiftPError = 0;
189 double maxInterpolationPError = 0;
190 double maxShiftVError = 0;
191 double maxInterpolationVError = 0;
192 for (Binary64 dt = zero; dt.getReal() < 240.0; dt = dt.add(1.0)) {
193 FieldAbsoluteDate<Binary64> t = initialOrbit.getDate().shiftedBy(dt);
194 FieldPVCoordinates<Binary64> propagated = propagator.propagate(t).getPVCoordinates();
195 FieldPVCoordinates<Binary64> shiftError = new FieldPVCoordinates<>(propagated,
196 initialOrbit.shiftedBy(dt.getReal())
197 .getPVCoordinates());
198 FieldPVCoordinates<Binary64> interpolationError = new FieldPVCoordinates<>(propagated,
199 interpolator.interpolate(t, sample)
200 .getPVCoordinates());
201 maxShiftPError = FastMath.max(maxShiftPError,
202 shiftError.getPosition().getNorm().getReal());
203 maxInterpolationPError = FastMath.max(maxInterpolationPError,
204 interpolationError.getPosition().getNorm().getReal());
205 maxShiftVError = FastMath.max(maxShiftVError,
206 shiftError.getVelocity().getNorm().getReal());
207 maxInterpolationVError = FastMath.max(maxInterpolationVError,
208 interpolationError.getVelocity().getNorm().getReal());
209 }
210 Assertions.assertEquals(shiftPositionErrorWithin, maxShiftPError, 0.01 * shiftPositionErrorWithin);
211 Assertions.assertEquals(interpolationPositionErrorWithin, maxInterpolationPError,
212 0.01 * interpolationPositionErrorWithin);
213 Assertions.assertEquals(shiftVelocityErrorWithin, maxShiftVError, 0.01 * shiftVelocityErrorWithin);
214 Assertions.assertEquals(interpolationVelocityErrorWithin, maxInterpolationVError,
215 0.01 * interpolationVelocityErrorWithin);
216
217
218 maxShiftPError = 0;
219 maxInterpolationPError = 0;
220 maxShiftVError = 0;
221 maxInterpolationVError = 0;
222 for (Binary64 dt = zero.add(500.0); dt.getReal() < 650.0; dt = dt.add(1.0)) {
223 FieldAbsoluteDate<Binary64> t = initialOrbit.getDate().shiftedBy(dt);
224 FieldPVCoordinates<Binary64> propagated = propagator.propagate(t).getPVCoordinates();
225 FieldPVCoordinates<Binary64> shiftError = new FieldPVCoordinates<>(propagated,
226 initialOrbit.shiftedBy(dt)
227 .getPVCoordinates());
228 FieldPVCoordinates<Binary64> interpolationError = new FieldPVCoordinates<>(propagated,
229 interpolator.interpolate(t, sample)
230 .getPVCoordinates());
231 maxShiftPError = FastMath.max(maxShiftPError,
232 shiftError.getPosition().getNorm().getReal());
233 maxInterpolationPError = FastMath.max(maxInterpolationPError,
234 interpolationError.getPosition().getNorm().getReal());
235 maxShiftVError = FastMath.max(maxShiftVError,
236 shiftError.getVelocity().getNorm().getReal());
237 maxInterpolationVError = FastMath.max(maxInterpolationVError,
238 interpolationError.getVelocity().getNorm().getReal());
239 }
240 Assertions.assertEquals(shiftPositionErrorFarPast, maxShiftPError, 0.01 * shiftPositionErrorFarPast);
241 Assertions.assertEquals(interpolationPositionErrorFarPast, maxInterpolationPError,
242 0.01 * interpolationPositionErrorFarPast);
243 Assertions.assertEquals(shiftVelocityErrorFarPast, maxShiftVError, 0.01 * shiftVelocityErrorFarPast);
244 Assertions.assertEquals(interpolationVelocityErrorFarPast, maxInterpolationVError,
245 0.01 * interpolationVelocityErrorFarPast);
246
247 }
248
249 private void doTestCircularInterpolation(boolean useDerivatives,
250 double shiftErrorWithin, double interpolationErrorWithin,
251 double shiftErrorSlightlyPast, double interpolationErrorSlightlyPast,
252 double shiftErrorFarPast, double interpolationErrorFarPast) {
253
254 Binary64 zero = field.getZero();
255 FieldAbsoluteDate<Binary64> date = new FieldAbsoluteDate<>(field);
256
257 final Binary64 ehMu = zero.add(3.9860047e14);
258 final double ae = 6.378137e6;
259 final double c20 = -1.08263e-3;
260 final double c30 = 2.54e-6;
261 final double c40 = 1.62e-6;
262 final double c50 = 2.3e-7;
263 final double c60 = -5.5e-7;
264
265 date = date.shiftedBy(584.);
266 final FieldVector3D<Binary64> position =
267 new FieldVector3D<>(zero.add(3220103.), zero.add(69623.), zero.add(6449822.));
268 final FieldVector3D<Binary64> velocity = new FieldVector3D<>(zero.add(6414.7), zero.add(-2006.), zero.add(-3180.));
269 final FieldCircularOrbit<Binary64> initialOrbit =
270 new FieldCircularOrbit<>(new FieldPVCoordinates<>(position, velocity),
271 FramesFactory.getEME2000(), date, ehMu);
272
273 FieldEcksteinHechlerPropagator<Binary64> propagator =
274 new FieldEcksteinHechlerPropagator<>(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
275
276
277 List<FieldOrbit<Binary64>> sample = new ArrayList<>();
278 for (Binary64 dt = zero; dt.getReal() < 300.0; dt = dt.add(60.0)) {
279 FieldOrbit<Binary64> orbit = OrbitType.CIRCULAR.convertType(propagator.propagate(date.shiftedBy(dt)).getOrbit());
280 if (!useDerivatives) {
281
282 Binary64[] stateVector = MathArrays.buildArray(field, 6);
283 orbit.getType().mapOrbitToArray(orbit, PositionAngleType.TRUE, stateVector, null);
284 orbit = orbit.getType().mapArrayToOrbit(stateVector, null, PositionAngleType.TRUE,
285 orbit.getDate(), orbit.getMu(), orbit.getFrame());
286 }
287 sample.add(orbit);
288 }
289
290
291 final FieldTimeInterpolator<FieldOrbit<Binary64>, Binary64> interpolator =
292 new FieldOrbitHermiteInterpolator<>(sample.size(), 759, initialOrbit.getFrame(),
293 CartesianDerivativesFilter.USE_PVA);
294
295
296 double maxShiftError = 0.0;
297 double maxInterpolationError = 0.0;
298 for (Binary64 dt = zero; dt.getReal() < 241.0; dt = dt.add(1.0)) {
299 FieldAbsoluteDate<Binary64> Binary64 = initialOrbit.getDate().shiftedBy(dt);
300 FieldVector3D<Binary64> shifted = initialOrbit.shiftedBy(dt.getReal()).getPosition();
301 FieldVector3D<Binary64> interpolated = interpolator.interpolate(Binary64, sample).getPosition();
302 FieldVector3D<Binary64> propagated = propagator.propagate(Binary64).getPosition();
303 maxShiftError = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm().getReal());
304 maxInterpolationError =
305 FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm().getReal());
306 }
307 Assertions.assertEquals(shiftErrorWithin, maxShiftError, 0.01 * shiftErrorWithin);
308 Assertions.assertEquals(interpolationErrorWithin, maxInterpolationError, 0.01 * interpolationErrorWithin);
309
310
311 maxShiftError = 0;
312 maxInterpolationError = 0;
313 for (Binary64 dt = zero.add(240); dt.getReal() < 300.0; dt = dt.add(1.0)) {
314 FieldAbsoluteDate<Binary64> Binary64 = initialOrbit.getDate().shiftedBy(dt);
315 FieldVector3D<Binary64> shifted = initialOrbit.shiftedBy(dt).getPosition();
316 FieldVector3D<Binary64> interpolated = interpolator.interpolate(Binary64, sample).getPosition();
317 FieldVector3D<Binary64> propagated = propagator.propagate(Binary64).getPosition();
318 maxShiftError = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm().getReal());
319 maxInterpolationError =
320 FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm().getReal());
321 }
322 Assertions.assertEquals(shiftErrorSlightlyPast, maxShiftError, 0.01 * shiftErrorSlightlyPast);
323 Assertions.assertEquals(interpolationErrorSlightlyPast, maxInterpolationError,
324 0.01 * interpolationErrorSlightlyPast);
325
326
327
328 maxShiftError = 0;
329 maxInterpolationError = 0;
330 for (Binary64 dt = zero.add(300); dt.getReal() < 1000; dt = dt.add(1.0)) {
331 FieldAbsoluteDate<Binary64> t = initialOrbit.getDate().shiftedBy(dt);
332 FieldVector3D<Binary64> shifted = initialOrbit.shiftedBy(dt).getPosition();
333 FieldVector3D<Binary64> interpolated = interpolator.interpolate(t, sample).getPosition();
334 FieldVector3D<Binary64> propagated = propagator.propagate(t).getPosition();
335 maxShiftError = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm().getReal());
336 maxInterpolationError =
337 FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm().getReal());
338 }
339 Assertions.assertEquals(shiftErrorFarPast, maxShiftError, 0.01 * shiftErrorFarPast);
340 Assertions.assertEquals(interpolationErrorFarPast, maxInterpolationError, 0.01 * interpolationErrorFarPast);
341
342 }
343
344 private void doTestEquinoctialInterpolation(boolean useDerivatives,
345 double shiftErrorWithin, double interpolationErrorWithin,
346 double shiftErrorSlightlyPast, double interpolationErrorSlightlyPast,
347 double shiftErrorFarPast, double interpolationErrorFarPast) {
348
349 Binary64 zero = field.getZero();
350 FieldAbsoluteDate<Binary64> date = new FieldAbsoluteDate<>(field);
351
352 final double ehMu = 3.9860047e14;
353 final double ae = 6.378137e6;
354 final double c20 = -1.08263e-3;
355 final double c30 = 2.54e-6;
356 final double c40 = 1.62e-6;
357 final double c50 = 2.3e-7;
358 final double c60 = -5.5e-7;
359
360 date = date.shiftedBy(584.);
361 final FieldVector3D<Binary64> position =
362 new FieldVector3D<>(zero.add(3220103.), zero.add(69623.), zero.add(6449822.));
363 final FieldVector3D<Binary64> velocity = new FieldVector3D<>(zero.add(6414.7), zero.add(-2006.), zero.add(-3180.));
364 final FieldEquinoctialOrbit<Binary64> initialOrbit =
365 new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
366 FramesFactory.getEME2000(), date, zero.add(ehMu));
367
368 FieldEcksteinHechlerPropagator<Binary64> propagator =
369 new FieldEcksteinHechlerPropagator<>(initialOrbit, ae, zero.add(ehMu), c20, c30, c40, c50, c60);
370
371
372 List<FieldOrbit<Binary64>> sample = new ArrayList<>();
373 for (double dt = 0; dt < 300.0; dt += 60.0) {
374 FieldOrbit<Binary64> orbit =
375 OrbitType.EQUINOCTIAL.convertType(propagator.propagate(date.shiftedBy(dt)).getOrbit());
376 if (!useDerivatives) {
377
378 Binary64[] stateVector = MathArrays.buildArray(field, 6);
379 orbit.getType().mapOrbitToArray(orbit, PositionAngleType.TRUE, stateVector, null);
380 orbit = orbit.getType().mapArrayToOrbit(stateVector, null, PositionAngleType.TRUE,
381 orbit.getDate(), orbit.getMu(), orbit.getFrame());
382 }
383 sample.add(orbit);
384 }
385
386
387 final FieldTimeInterpolator<FieldOrbit<Binary64>, Binary64> interpolator =
388 new FieldOrbitHermiteInterpolator<>(sample.size(), 759, initialOrbit.getFrame(),
389 CartesianDerivativesFilter.USE_PVA);
390
391
392 double maxShiftError = 0;
393 double maxInterpolationError = 0;
394 for (Binary64 dt = zero; dt.getReal() < 241.0; dt = dt.add(1.0)) {
395 FieldAbsoluteDate<Binary64> t = initialOrbit.getDate().shiftedBy(dt);
396 FieldVector3D<Binary64> shifted = initialOrbit.shiftedBy(dt.getReal()).getPosition();
397 FieldVector3D<Binary64> interpolated = interpolator.interpolate(t, sample).getPosition();
398 FieldVector3D<Binary64> propagated = propagator.propagate(t).getPosition();
399 maxShiftError = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm().getReal());
400 maxInterpolationError =
401 FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm().getReal());
402 }
403 Assertions.assertEquals(shiftErrorWithin, maxShiftError, 0.01 * shiftErrorWithin);
404 Assertions.assertEquals(interpolationErrorWithin, maxInterpolationError, 0.01 * interpolationErrorWithin);
405
406
407 maxShiftError = 0;
408 maxInterpolationError = 0;
409 for (Binary64 dt = zero.add(240); dt.getReal() < 300.0; dt = dt.add(1.0)) {
410 FieldAbsoluteDate<Binary64> t = initialOrbit.getDate().shiftedBy(dt);
411 FieldVector3D<Binary64> shifted = initialOrbit.shiftedBy(dt).getPosition();
412 FieldVector3D<Binary64> interpolated = interpolator.interpolate(t, sample).getPosition();
413 FieldVector3D<Binary64> propagated = propagator.propagate(t).getPosition();
414 maxShiftError = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm().getReal());
415 maxInterpolationError =
416 FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm().getReal());
417 }
418 Assertions.assertEquals(shiftErrorSlightlyPast, maxShiftError, 0.01 * shiftErrorSlightlyPast);
419 Assertions.assertEquals(interpolationErrorSlightlyPast, maxInterpolationError,
420 0.01 * interpolationErrorSlightlyPast);
421
422
423
424 maxShiftError = 0;
425 maxInterpolationError = 0;
426 for (Binary64 dt = zero.add(300); dt.getReal() < 1000; dt = dt.add(1.0)) {
427 FieldAbsoluteDate<Binary64> t = initialOrbit.getDate().shiftedBy(dt);
428 FieldVector3D<Binary64> shifted = initialOrbit.shiftedBy(dt).getPosition();
429 FieldVector3D<Binary64> interpolated = interpolator.interpolate(t, sample).getPosition();
430 FieldVector3D<Binary64> propagated = propagator.propagate(t).getPosition();
431 maxShiftError = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm().getReal());
432 maxInterpolationError =
433 FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm().getReal());
434 }
435 Assertions.assertEquals(shiftErrorFarPast, maxShiftError, 0.01 * shiftErrorFarPast);
436 Assertions.assertEquals(interpolationErrorFarPast, maxInterpolationError, 0.01 * interpolationErrorFarPast);
437
438 }
439
440 private void doTestKeplerianInterpolation(boolean useDerivatives,
441 double shiftPositionErrorWithin, double interpolationPositionErrorWithin,
442 double shiftEccentricityErrorWithin,
443 double interpolationEccentricityErrorWithin,
444 double shiftPositionErrorSlightlyPast,
445 double interpolationPositionErrorSlightlyPast,
446 double shiftEccentricityErrorSlightlyPast,
447 double interpolationEccentricityErrorSlightlyPast) {
448
449 final Binary64 zero = field.getZero();
450 final Binary64 ehMu = zero.add(3.9860047e14);
451 final double ae = 6.378137e6;
452 final double c20 = -1.08263e-3;
453 final double c30 = 2.54e-6;
454 final double c40 = 1.62e-6;
455 final double c50 = 2.3e-7;
456 final double c60 = -5.5e-7;
457
458 final FieldAbsoluteDate<Binary64> date = FieldAbsoluteDate.getJ2000Epoch(field).shiftedBy(584.);
459 final FieldVector3D<Binary64> position =
460 new FieldVector3D<>(field.getZero().add(3220103.), field.getZero().add(69623.),
461 field.getZero().add(6449822.));
462 final FieldVector3D<Binary64> velocity =
463 new FieldVector3D<>(field.getZero().add(6414.7), field.getZero().add(-2006.), field.getZero().add(-3180.));
464 final FieldKeplerianOrbit<Binary64> initialOrbit =
465 new FieldKeplerianOrbit<>(new FieldPVCoordinates<>(position, velocity),
466 FramesFactory.getEME2000(), date, ehMu);
467
468 FieldEcksteinHechlerPropagator<Binary64> propagator =
469 new FieldEcksteinHechlerPropagator<>(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
470
471
472 List<FieldOrbit<Binary64>> sample = new ArrayList<>();
473 for (double dt = 0; dt < 300.0; dt += 60.0) {
474 FieldOrbit<Binary64> orbit =
475 OrbitType.KEPLERIAN.convertType(propagator.propagate(date.shiftedBy(dt)).getOrbit());
476 if (!useDerivatives) {
477
478 Binary64[] stateVector = MathArrays.buildArray(field, 6);
479 orbit.getType().mapOrbitToArray(orbit, PositionAngleType.TRUE, stateVector, null);
480 orbit = orbit.getType().mapArrayToOrbit(stateVector, null, PositionAngleType.TRUE,
481 orbit.getDate(), orbit.getMu(), orbit.getFrame());
482 }
483 sample.add(orbit);
484 }
485
486
487 final FieldTimeInterpolator<FieldOrbit<Binary64>, Binary64> interpolator =
488 new FieldOrbitHermiteInterpolator<>(sample.size(), 359, initialOrbit.getFrame(),
489 CartesianDerivativesFilter.USE_PVA);
490
491
492
493
494 double maxShiftPositionError = 0;
495 double maxInterpolationPositionError = 0;
496 double maxShiftEccentricityError = 0;
497 double maxInterpolationEccentricityError = 0;
498 for (double dt = 0; dt < 241.0; dt += 1.0) {
499 FieldAbsoluteDate<Binary64> t = initialOrbit.getDate().shiftedBy(dt);
500 FieldVector3D<Binary64> shiftedP = initialOrbit.shiftedBy(dt).getPosition();
501 FieldVector3D<Binary64> interpolatedP = interpolator.interpolate(t, sample).getPosition();
502 FieldVector3D<Binary64> propagatedP = propagator.propagate(t).getPosition();
503 Binary64 shiftedE = initialOrbit.shiftedBy(zero.add(dt)).getE();
504 Binary64 interpolatedE = interpolator.interpolate(t, sample).getE();
505 Binary64 propagatedE = propagator.propagate(t).getOrbit().getE();
506 maxShiftPositionError =
507 FastMath.max(maxShiftPositionError, shiftedP.subtract(propagatedP).getNorm().getReal());
508 maxInterpolationPositionError =
509 FastMath.max(maxInterpolationPositionError, interpolatedP.subtract(propagatedP).getNorm().getReal());
510 maxShiftEccentricityError =
511 FastMath.max(maxShiftEccentricityError, shiftedE.subtract(propagatedE).abs().getReal());
512 maxInterpolationEccentricityError =
513 FastMath.max(maxInterpolationEccentricityError, interpolatedE.subtract(propagatedE).abs().getReal());
514 }
515 Assertions.assertEquals(shiftPositionErrorWithin, maxShiftPositionError, 0.01 * shiftPositionErrorWithin);
516 Assertions.assertEquals(interpolationPositionErrorWithin, maxInterpolationPositionError,
517 0.01 * interpolationPositionErrorWithin);
518 Assertions.assertEquals(shiftEccentricityErrorWithin, maxShiftEccentricityError,
519 0.01 * shiftEccentricityErrorWithin);
520 Assertions.assertEquals(interpolationEccentricityErrorWithin, maxInterpolationEccentricityError,
521 0.01 * interpolationEccentricityErrorWithin);
522
523
524
525
526
527 maxShiftPositionError = 0;
528 maxInterpolationPositionError = 0;
529 maxShiftEccentricityError = 0;
530 maxInterpolationEccentricityError = 0;
531 for (double dt = 240; dt < 600; dt += 1.0) {
532 FieldAbsoluteDate<Binary64> t = initialOrbit.getDate().shiftedBy(dt);
533 FieldVector3D<Binary64> shiftedP = initialOrbit.shiftedBy(zero.add(dt)).getPosition();
534 FieldVector3D<Binary64> interpolatedP = interpolator.interpolate(t, sample).getPosition();
535 FieldVector3D<Binary64> propagatedP = propagator.propagate(t).getPosition();
536 Binary64 shiftedE = initialOrbit.shiftedBy(zero.add(dt)).getE();
537 Binary64 interpolatedE = interpolator.interpolate(t, sample).getE();
538 Binary64 propagatedE = propagator.propagate(t).getOrbit().getE();
539 maxShiftPositionError =
540 FastMath.max(maxShiftPositionError, shiftedP.subtract(propagatedP).getNorm().getReal());
541 maxInterpolationPositionError =
542 FastMath.max(maxInterpolationPositionError, interpolatedP.subtract(propagatedP).getNorm().getReal());
543 maxShiftEccentricityError =
544 FastMath.max(maxShiftEccentricityError, shiftedE.subtract(propagatedE).abs().getReal());
545 maxInterpolationEccentricityError =
546 FastMath.max(maxInterpolationEccentricityError, interpolatedE.subtract(propagatedE).abs().getReal());
547 }
548 Assertions.assertEquals(shiftPositionErrorSlightlyPast, maxShiftPositionError,
549 0.01 * shiftPositionErrorSlightlyPast);
550 Assertions.assertEquals(interpolationPositionErrorSlightlyPast, maxInterpolationPositionError,
551 0.01 * interpolationPositionErrorSlightlyPast);
552 Assertions.assertEquals(shiftEccentricityErrorSlightlyPast, maxShiftEccentricityError,
553 0.01 * shiftEccentricityErrorSlightlyPast);
554 Assertions.assertEquals(interpolationEccentricityErrorSlightlyPast, maxInterpolationEccentricityError,
555 0.01 * interpolationEccentricityErrorSlightlyPast);
556
557 }
558
559 @Test
560 @DisplayName("test default constructor")
561 void testDefaultConstructor() {
562
563 final int interpolationPoints = 2;
564
565
566 final FieldOrbitHermiteInterpolator<Binary64> interpolator =
567 new FieldOrbitHermiteInterpolator<>(FramesFactory.getGCRF());
568
569
570 Assertions.assertEquals(AbstractTimeInterpolator.DEFAULT_EXTRAPOLATION_THRESHOLD_SEC,
571 interpolator.getExtrapolationThreshold());
572 Assertions.assertEquals(interpolationPoints,
573 interpolator.getNbInterpolationPoints());
574 Assertions.assertEquals(CartesianDerivativesFilter.USE_PVA, interpolator.getPVAFilter());
575 }
576
577 @Test
578 public void testErrorThrownWhenUsingOrbitWithFrameMismatch() {
579
580
581 final Field<Binary64> field = Binary64Field.getInstance();
582 final FieldOrbit<Binary64> orbit1 = TestUtils.getFakeFieldOrbit();
583 final FieldOrbit<Binary64> orbit2 = new FieldCartesianOrbit<>(orbit1.shiftedBy(1).getPVCoordinates(),
584 FramesFactory.getEME2000(),
585 orbit1.getMu());
586
587 final List<FieldOrbit<Binary64>> sample = new ArrayList<>();
588 sample.add(orbit1);
589 sample.add(orbit2);
590
591
592 final FieldOrbitHermiteInterpolator<Binary64> interpolator =
593 new FieldOrbitHermiteInterpolator<>(FramesFactory.getGCRF());
594
595
596 Assertions.assertThrows(OrekitIllegalArgumentException.class,
597 () -> interpolator.interpolate(new FieldAbsoluteDate<>(field), sample));
598 }
599
600
601
602
603
604
605 @Test
606 void testOutputFrameCartesian() {
607
608
609 final FieldOrbit<Binary64> orbit1 = TestUtils.getFakeFieldOrbit();
610 final FieldOrbit<Binary64> orbit2 = orbit1.shiftedBy(1);
611
612 final List<FieldOrbit<Binary64>> samples = new ArrayList<>();
613 samples.add(orbit1);
614 samples.add(orbit2);
615
616
617 final Frame outputFrame = FramesFactory.getEME2000();
618
619
620 final FieldOrbitHermiteInterpolator<Binary64> interpolator = new FieldOrbitHermiteInterpolator<>(outputFrame);
621
622
623 final FieldAbsoluteDate<Binary64> interpolationDate = orbit1.getDate();
624
625
626 final FieldOrbit<Binary64> interpolatedOrbit = interpolator.interpolate(interpolationDate, samples);
627
628
629
630 final Vector3D originalPosition = orbit1.getPosition().toVector3D();
631 final Vector3D actualPosition = interpolatedOrbit.getPosition(orbit1.getFrame()).toVector3D();
632 Assertions.assertEquals(originalPosition.getX(), actualPosition.getX(), 1.0e-14);
633 Assertions.assertEquals(originalPosition.getY(), actualPosition.getY(), 1.0e-14);
634 Assertions.assertEquals(originalPosition.getZ(), actualPosition.getZ(), 1.0e-14);
635
636
637 Assertions.assertInstanceOf(FieldCartesianOrbit.class, interpolatedOrbit);
638
639
640 Assertions.assertEquals(outputFrame, interpolatedOrbit.getFrame());
641 }
642
643
644
645
646
647
648 @Test
649 void testOutputFrameCommon() {
650
651
652 final FieldOrbit<Binary64> orbit1 = new FieldCircularOrbit<>(TestUtils.getTestFieldOrbit());
653 final FieldOrbit<Binary64> orbit2 = new FieldCircularOrbit<>(orbit1.shiftedBy(1));
654
655 final List<FieldOrbit<Binary64>> samples = new ArrayList<>();
656 samples.add(orbit1);
657 samples.add(orbit2);
658
659
660 final Frame outputFrame = FramesFactory.getEME2000();
661
662
663 final FieldOrbitHermiteInterpolator<Binary64> interpolator = new FieldOrbitHermiteInterpolator<>(outputFrame);
664
665
666 final FieldAbsoluteDate<Binary64> interpolationDate = orbit1.getDate();
667
668
669 final FieldOrbit<Binary64> interpolatedOrbit = interpolator.interpolate(interpolationDate, samples);
670
671
672
673 final Vector3D originalPosition = orbit1.getPosition().toVector3D();
674 final Vector3D actualPosition = interpolatedOrbit.getPosition(orbit1.getFrame()).toVector3D();
675 Assertions.assertEquals(originalPosition.getX(), actualPosition.getX(), 1.0e-8);
676 Assertions.assertEquals(originalPosition.getY(), actualPosition.getY(), 1.0e-8);
677 Assertions.assertEquals(originalPosition.getZ(), actualPosition.getZ(), 1.0e-8);
678
679
680 Assertions.assertEquals(outputFrame, interpolatedOrbit.getFrame());
681 }
682 }