1   /* Copyright 2002-2025 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS licenses this file to You under the Apache License, Version 2.0
6    * (the "License"); you may not use this file except in compliance with
7    * the License.  You may obtain a copy of the License at
8    *
9    *   http://www.apache.org/licenses/LICENSE-2.0
10   *
11   * Unless required by applicable law or agreed to in writing, software
12   * distributed under the License is distributed on an "AS IS" BASIS,
13   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14   * See the License for the specific language governing permissions and
15   * limitations under the License.
16   */
17  package org.orekit.orbits;
18  
19  import java.util.ArrayList;
20  import java.util.List;
21  
22  import org.hipparchus.geometry.euclidean.threed.Vector3D;
23  import org.hipparchus.util.FastMath;
24  import org.junit.jupiter.api.Assertions;
25  import org.junit.jupiter.api.BeforeEach;
26  import org.junit.jupiter.api.DisplayName;
27  import org.junit.jupiter.api.Test;
28  import org.orekit.Utils;
29  import org.orekit.frames.FramesFactory;
30  import org.orekit.propagation.analytical.EcksteinHechlerPropagator;
31  import org.orekit.time.AbsoluteDate;
32  import org.orekit.time.AbstractTimeInterpolator;
33  import org.orekit.time.TimeInterpolator;
34  import org.orekit.utils.CartesianDerivativesFilter;
35  import org.orekit.utils.PVCoordinates;
36  
37  class OrbitHermiteInterpolatorTest {
38  
39      @BeforeEach
40      public void setUp() {
41          Utils.setDataRoot("regular-data");
42      }
43  
44      @Test
45      public void testCartesianInterpolationWithDerivatives() {
46          doTestCartesianInterpolation(true, CartesianDerivativesFilter.USE_P,
47                                       394, 0.1968, 3.21, 0.021630,
48                                       2474, 2707.6418, 6.6, 26.28);
49  
50          doTestCartesianInterpolation(true, CartesianDerivativesFilter.USE_PV,
51                                       394, 8.2392E-9, 3.21, 1.3236E-10,
52                                       2474, 0.07474, 6.6, 0.001450);
53  
54          // Solution with PVA less precise than with PV only as the interpolating polynomial begins to oscillate heavily
55          // outside the interpolating interval
56          doTestCartesianInterpolation(true, CartesianDerivativesFilter.USE_PVA,
57                                       394, 9.60e-9, 3.21, 5.39e-10,
58                                       2474, 2614, 6.55, 71);
59  
60      }
61  
62      @Test
63      public void testCartesianInterpolationWithoutDerivatives() {
64          doTestCartesianInterpolation(false, CartesianDerivativesFilter.USE_P,
65                                       394, 0.1968, 3.21, 0.02163,
66                                       2474, 2707.6419, 6.55, 26.2826);
67  
68          doTestCartesianInterpolation(false, CartesianDerivativesFilter.USE_PV,
69                                       394, 8.2392E-9, 3.21, 1.3236E-10,
70                                       2474, 0.07474, 6.55, 0.001450);
71  
72          // Interpolation without derivatives is very wrong in PVA as we force first and second derivatives to be 0 i.e. we
73          // give false information to the interpolator
74          doTestCartesianInterpolation(false, CartesianDerivativesFilter.USE_PVA,
75                                       394, 2.61, 3.21, 0.154,
76                                       2474, 2.28e12, 6.55, 6.22e10);
77      }
78  
79      @Test
80      public void testCircularInterpolationWithDerivatives() {
81          doTestCircularInterpolation(true,
82                                      397, 2.26e-8,
83                                      610, 3.99e-6,
84                                      4870, 113.8);
85      }
86  
87      @Test
88      public void testCircularInterpolationWithoutDerivatives() {
89          doTestCircularInterpolation(false,
90                                      397, 0.0372,
91                                      610.0, 1.23,
92                                      4870, 8869);
93      }
94  
95      @Test
96      public void testEquinoctialInterpolationWithDerivatives() {
97          doTestEquinoctialInterpolation(true,
98                                         397, 9.84e-9,
99                                         610, 4.28e-6,
100                                        4870, 115);
101     }
102 
103     @Test
104     public void testEquinoctialInterpolationWithoutDerivatives() {
105         doTestEquinoctialInterpolation(false,
106                                        397, 0.0372,
107                                        610.0, 1.23,
108                                        4879, 8871);
109     }
110 
111     @Test
112     public void testKeplerianInterpolationWithDerivatives() {
113         doTestKeplerianInterpolation(true,
114                                      397, 4.01, 4.75e-4, 1.28e-7,
115                                      2159, 1.05e7, 1.19e-3, 0.773);
116     }
117 
118     @Test
119     public void testKeplerianInterpolationWithoutDerivatives() {
120         doTestKeplerianInterpolation(false,
121                                      397, 62.0, 4.75e-4, 2.87e-6,
122                                      2159, 79365, 1.19e-3, 3.89e-3);
123     }
124 
125     private void doTestCartesianInterpolation(boolean useDerivatives, CartesianDerivativesFilter pvaFilter,
126                                               double shiftPositionErrorWithin, double interpolationPositionErrorWithin,
127                                               double shiftVelocityErrorWithin, double interpolationVelocityErrorWithin,
128                                               double shiftPositionErrorFarPast, double interpolationPositionErrorFarPast,
129                                               double shiftVelocityErrorFarPast, double interpolationVelocityErrorFarPast) {
130 
131         final double ehMu = 3.9860047e14;
132         final double ae   = 6.378137e6;
133         final double c20  = -1.08263e-3;
134         final double c30  = 2.54e-6;
135         final double c40  = 1.62e-6;
136         final double c50  = 2.3e-7;
137         final double c60  = -5.5e-7;
138 
139         final AbsoluteDate date     = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
140         final Vector3D     position = new Vector3D(3220103., 69623., 6449822.);
141         final Vector3D     velocity = new Vector3D(6414.7, -2006., -3180.);
142         final CartesianOrbit initialOrbit = new CartesianOrbit(new PVCoordinates(position, velocity),
143                                                                FramesFactory.getEME2000(), date, ehMu);
144 
145         EcksteinHechlerPropagator propagator =
146                 new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
147 
148         // set up a 5 points sample
149         List<Orbit> sample = new ArrayList<>();
150         for (double dt = 0; dt < 251.0; dt += 60) {
151             Orbit orbit = propagator.propagate(date.shiftedBy(dt)).getOrbit();
152             if (!useDerivatives) {
153                 // remove derivatives
154                 double[] stateVector = new double[6];
155                 orbit.getType().mapOrbitToArray(orbit, PositionAngleType.TRUE, stateVector, null);
156                 orbit = orbit.getType().mapArrayToOrbit(stateVector, null, PositionAngleType.TRUE,
157                                                         orbit.getDate(), orbit.getMu(), orbit.getFrame());
158             }
159             sample.add(orbit);
160         }
161 
162         // Create interpolator
163         final TimeInterpolator<Orbit> interpolator =
164                 new OrbitHermiteInterpolator(sample.size(), 409, initialOrbit.getFrame(), pvaFilter);
165 
166         // well inside the sample, interpolation should be much better than Keplerian shift
167         // this is because we take the full non-Keplerian acceleration into account in
168         // the Cartesian parameters, which in this case is preserved by the
169         // Eckstein-Hechler propagator
170         double maxShiftPError         = 0;
171         double maxInterpolationPError = 0;
172         double maxShiftVError         = 0;
173         double maxInterpolationVError = 0;
174 
175         for (double dt = 0; dt < 240.0; dt += 1.0) {
176             AbsoluteDate  t            = date.shiftedBy(dt);
177             PVCoordinates propagated   = propagator.propagate(t).getPVCoordinates();
178             PVCoordinates interpolated = interpolator.interpolate(t, sample).getPVCoordinates();
179 
180             PVCoordinates shiftError = new PVCoordinates(propagated,
181                                                          initialOrbit.shiftedBy(dt).getPVCoordinates());
182             PVCoordinates interpolationError = new PVCoordinates(propagated, interpolated);
183 
184             maxShiftPError         = FastMath.max(maxShiftPError,
185                                                   shiftError.getPosition().getNorm());
186             maxInterpolationPError = FastMath.max(maxInterpolationPError,
187                                                   interpolationError.getPosition().getNorm());
188             maxShiftVError         = FastMath.max(maxShiftVError,
189                                                   shiftError.getVelocity().getNorm());
190             maxInterpolationVError = FastMath.max(maxInterpolationVError,
191                                                   interpolationError.getVelocity().getNorm());
192         }
193 
194         Assertions.assertEquals(shiftPositionErrorWithin, maxShiftPError, 0.01 * shiftPositionErrorWithin);
195         Assertions.assertEquals(interpolationPositionErrorWithin, maxInterpolationPError,
196                                 0.01 * interpolationPositionErrorWithin);
197         Assertions.assertEquals(shiftVelocityErrorWithin, maxShiftVError, 0.01 * shiftVelocityErrorWithin);
198         Assertions.assertEquals(interpolationVelocityErrorWithin, maxInterpolationVError,
199                                 0.01 * interpolationVelocityErrorWithin);
200 
201         // if we go far past sample end, interpolation becomes worse than Keplerian shift
202         maxShiftPError         = 0;
203         maxInterpolationPError = 0;
204         maxShiftVError         = 0;
205         maxInterpolationVError = 0;
206         for (double dt = 500.0; dt < 650.0; dt += 1.0) {
207             AbsoluteDate  t            = initialOrbit.getDate().shiftedBy(dt);
208             PVCoordinates propagated   = propagator.propagate(t).getPVCoordinates();
209             PVCoordinates interpolated = interpolator.interpolate(t, sample).getPVCoordinates();
210 
211             PVCoordinates shiftError = new PVCoordinates(propagated,
212                                                          initialOrbit.shiftedBy(dt).getPVCoordinates());
213             PVCoordinates interpolationError = new PVCoordinates(propagated, interpolated);
214 
215             maxShiftPError         = FastMath.max(maxShiftPError,
216                                                   shiftError.getPosition().getNorm());
217             maxInterpolationPError = FastMath.max(maxInterpolationPError,
218                                                   interpolationError.getPosition().getNorm());
219             maxShiftVError         = FastMath.max(maxShiftVError,
220                                                   shiftError.getVelocity().getNorm());
221             maxInterpolationVError = FastMath.max(maxInterpolationVError,
222                                                   interpolationError.getVelocity().getNorm());
223 
224         }
225 
226         Assertions.assertEquals(shiftPositionErrorFarPast, maxShiftPError, 0.01 * shiftPositionErrorFarPast);
227         Assertions.assertEquals(interpolationPositionErrorFarPast, maxInterpolationPError,
228                                 0.01 * interpolationPositionErrorFarPast);
229         Assertions.assertEquals(shiftVelocityErrorFarPast, maxShiftVError, 0.01 * shiftVelocityErrorFarPast);
230         Assertions.assertEquals(interpolationVelocityErrorFarPast, maxInterpolationVError,
231                                 0.01 * interpolationVelocityErrorFarPast);
232     }
233 
234     private void doTestCircularInterpolation(boolean useDerivatives,
235                                              double shiftErrorWithin, double interpolationErrorWithin,
236                                              double shiftErrorSlightlyPast, double interpolationErrorSlightlyPast,
237                                              double shiftErrorFarPast, double interpolationErrorFarPast) {
238 
239         final double ehMu = 3.9860047e14;
240         final double ae   = 6.378137e6;
241         final double c20  = -1.08263e-3;
242         final double c30  = 2.54e-6;
243         final double c40  = 1.62e-6;
244         final double c50  = 2.3e-7;
245         final double c60  = -5.5e-7;
246 
247         final AbsoluteDate date     = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
248         final Vector3D     position = new Vector3D(3220103., 69623., 6449822.);
249         final Vector3D     velocity = new Vector3D(6414.7, -2006., -3180.);
250         final CircularOrbit initialOrbit = new CircularOrbit(new PVCoordinates(position, velocity),
251                                                              FramesFactory.getEME2000(), date, ehMu);
252 
253         EcksteinHechlerPropagator propagator =
254                 new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
255 
256         // set up a 5 points sample
257         List<Orbit> sample = new ArrayList<>();
258         for (double dt = 0; dt < 300.0; dt += 60.0) {
259             Orbit orbit = OrbitType.CIRCULAR.convertType(propagator.propagate(date.shiftedBy(dt)).getOrbit());
260             if (!useDerivatives) {
261                 // remove derivatives
262                 double[] stateVector = new double[6];
263                 orbit.getType().mapOrbitToArray(orbit, PositionAngleType.TRUE, stateVector, null);
264                 orbit = orbit.getType().mapArrayToOrbit(stateVector, null, PositionAngleType.TRUE,
265                                                         orbit.getDate(), orbit.getMu(), orbit.getFrame());
266             }
267             sample.add(orbit);
268         }
269 
270         // Create interpolator
271         final TimeInterpolator<Orbit> interpolator =
272                 new OrbitHermiteInterpolator(sample.size(), 759, initialOrbit.getFrame(),
273                                              CartesianDerivativesFilter.USE_PVA);
274 
275         // well inside the sample, interpolation should be much better than Keplerian shift
276         double maxShiftError         = 0;
277         double maxInterpolationError = 0;
278         for (double dt = 0; dt < 241.0; dt += 1.0) {
279             AbsoluteDate t            = initialOrbit.getDate().shiftedBy(dt);
280             Vector3D     shifted      = initialOrbit.shiftedBy(dt).getPosition();
281             Vector3D     interpolated = interpolator.interpolate(t, sample).getPosition();
282             Vector3D     propagated   = propagator.propagate(t).getPosition();
283             maxShiftError         = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm());
284             maxInterpolationError = FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm());
285         }
286         Assertions.assertEquals(shiftErrorWithin, maxShiftError, 0.01 * shiftErrorWithin);
287         Assertions.assertEquals(interpolationErrorWithin, maxInterpolationError, 0.01 * interpolationErrorWithin);
288 
289         // slightly past sample end, interpolation should quickly increase, but remain reasonable
290         maxShiftError         = 0;
291         maxInterpolationError = 0;
292         for (double dt = 240; dt < 300.0; dt += 1.0) {
293             AbsoluteDate t            = initialOrbit.getDate().shiftedBy(dt);
294             Vector3D     shifted      = initialOrbit.shiftedBy(dt).getPosition();
295             Vector3D     interpolated = interpolator.interpolate(t, sample).getPosition();
296             Vector3D     propagated   = propagator.propagate(t).getPosition();
297             maxShiftError         = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm());
298             maxInterpolationError = FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm());
299         }
300         Assertions.assertEquals(shiftErrorSlightlyPast, maxShiftError, 0.01 * shiftErrorSlightlyPast);
301         Assertions.assertEquals(interpolationErrorSlightlyPast, maxInterpolationError,
302                                 0.01 * interpolationErrorSlightlyPast);
303 
304         // far past sample end, interpolation should become really wrong
305         maxShiftError         = 0;
306         maxInterpolationError = 0;
307         for (double dt = 300; dt < 1000; dt += 1.0) {
308             AbsoluteDate t            = initialOrbit.getDate().shiftedBy(dt);
309             Vector3D     shifted      = initialOrbit.shiftedBy(dt).getPosition();
310             Vector3D     interpolated = interpolator.interpolate(t, sample).getPosition();
311             Vector3D     propagated   = propagator.propagate(t).getPosition();
312             maxShiftError         = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm());
313             maxInterpolationError = FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm());
314         }
315         Assertions.assertEquals(shiftErrorFarPast, maxShiftError, 0.01 * shiftErrorFarPast);
316         Assertions.assertEquals(interpolationErrorFarPast, maxInterpolationError, 0.01 * interpolationErrorFarPast);
317 
318     }
319 
320     private void doTestEquinoctialInterpolation(boolean useDerivatives,
321                                                 double shiftErrorWithin, double interpolationErrorWithin,
322                                                 double shiftErrorSlightlyPast, double interpolationErrorSlightlyPast,
323                                                 double shiftErrorFarPast, double interpolationErrorFarPast) {
324 
325         final double ehMu = 3.9860047e14;
326         final double ae   = 6.378137e6;
327         final double c20  = -1.08263e-3;
328         final double c30  = 2.54e-6;
329         final double c40  = 1.62e-6;
330         final double c50  = 2.3e-7;
331         final double c60  = -5.5e-7;
332 
333         final AbsoluteDate date     = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
334         final Vector3D     position = new Vector3D(3220103., 69623., 6449822.);
335         final Vector3D     velocity = new Vector3D(6414.7, -2006., -3180.);
336         final EquinoctialOrbit initialOrbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
337                                                                    FramesFactory.getEME2000(), date, ehMu);
338 
339         EcksteinHechlerPropagator propagator =
340                 new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
341 
342         // set up a 5 points sample
343         List<Orbit> sample = new ArrayList<>();
344         for (double dt = 0; dt < 300.0; dt += 60.0) {
345             Orbit orbit = OrbitType.EQUINOCTIAL.convertType(propagator.propagate(date.shiftedBy(dt)).getOrbit());
346             if (!useDerivatives) {
347                 // remove derivatives
348                 double[] stateVector = new double[6];
349                 orbit.getType().mapOrbitToArray(orbit, PositionAngleType.TRUE, stateVector, null);
350                 orbit = orbit.getType().mapArrayToOrbit(stateVector, null, PositionAngleType.TRUE,
351                                                         orbit.getDate(), orbit.getMu(), orbit.getFrame());
352             }
353             sample.add(orbit);
354         }
355 
356         // Create interpolator
357         final TimeInterpolator<Orbit> interpolator =
358                 new OrbitHermiteInterpolator(sample.size(), 759, initialOrbit.getFrame(),
359                                              CartesianDerivativesFilter.USE_PVA);
360 
361         // well inside the sample, interpolation should be much better than Keplerian shift
362         double maxShiftError         = 0;
363         double maxInterpolationError = 0;
364         for (double dt = 0; dt < 241.0; dt += 1.0) {
365             AbsoluteDate t            = initialOrbit.getDate().shiftedBy(dt);
366             Vector3D     shifted      = initialOrbit.shiftedBy(dt).getPosition();
367             Vector3D     interpolated = interpolator.interpolate(t, sample).getPosition();
368             Vector3D     propagated   = propagator.propagate(t).getPosition();
369             maxShiftError         = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm());
370             maxInterpolationError = FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm());
371         }
372         Assertions.assertEquals(shiftErrorWithin, maxShiftError, 0.01 * shiftErrorWithin);
373         Assertions.assertEquals(interpolationErrorWithin, maxInterpolationError, 0.01 * interpolationErrorWithin);
374 
375         // slightly past sample end, interpolation should quickly increase, but remain reasonable
376         maxShiftError         = 0;
377         maxInterpolationError = 0;
378         for (double dt = 240; dt < 300.0; dt += 1.0) {
379             AbsoluteDate t            = initialOrbit.getDate().shiftedBy(dt);
380             Vector3D     shifted      = initialOrbit.shiftedBy(dt).getPosition();
381             Vector3D     interpolated = interpolator.interpolate(t, sample).getPosition();
382             Vector3D     propagated   = propagator.propagate(t).getPosition();
383             maxShiftError         = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm());
384             maxInterpolationError = FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm());
385         }
386         Assertions.assertEquals(shiftErrorSlightlyPast, maxShiftError, 0.01 * shiftErrorSlightlyPast);
387         Assertions.assertEquals(interpolationErrorSlightlyPast, maxInterpolationError,
388                                 0.01 * interpolationErrorSlightlyPast);
389 
390         // far past sample end, interpolation should become really wrong
391         // (in this test case, break even occurs at around 863 seconds, with a 3.9 km error)
392         maxShiftError         = 0;
393         maxInterpolationError = 0;
394         for (double dt = 300; dt < 1000; dt += 1.0) {
395             AbsoluteDate t            = initialOrbit.getDate().shiftedBy(dt);
396             Vector3D     shifted      = initialOrbit.shiftedBy(dt).getPosition();
397             Vector3D     interpolated = interpolator.interpolate(t, sample).getPosition();
398             Vector3D     propagated   = propagator.propagate(t).getPosition();
399             maxShiftError         = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm());
400             maxInterpolationError = FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm());
401         }
402         Assertions.assertEquals(shiftErrorFarPast, maxShiftError, 0.01 * shiftErrorFarPast);
403         Assertions.assertEquals(interpolationErrorFarPast, maxInterpolationError, 0.01 * interpolationErrorFarPast);
404 
405     }
406 
407     private void doTestKeplerianInterpolation(boolean useDerivatives,
408                                               double shiftPositionErrorWithin, double interpolationPositionErrorWithin,
409                                               double shiftEccentricityErrorWithin,
410                                               double interpolationEccentricityErrorWithin,
411                                               double shiftPositionErrorSlightlyPast,
412                                               double interpolationPositionErrorSlightlyPast,
413                                               double shiftEccentricityErrorSlightlyPast,
414                                               double interpolationEccentricityErrorSlightlyPast) {
415 
416         final double ehMu = 3.9860047e14;
417         final double ae   = 6.378137e6;
418         final double c20  = -1.08263e-3;
419         final double c30  = 2.54e-6;
420         final double c40  = 1.62e-6;
421         final double c50  = 2.3e-7;
422         final double c60  = -5.5e-7;
423 
424         final AbsoluteDate date     = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
425         final Vector3D     position = new Vector3D(3220103., 69623., 6449822.);
426         final Vector3D     velocity = new Vector3D(6414.7, -2006., -3180.);
427         final KeplerianOrbit initialOrbit = new KeplerianOrbit(new PVCoordinates(position, velocity),
428                                                                FramesFactory.getEME2000(), date, ehMu);
429 
430         EcksteinHechlerPropagator propagator =
431                 new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
432 
433         // set up a 5 points sample
434         List<Orbit> sample = new ArrayList<>();
435         for (double dt = 0; dt < 300.0; dt += 60.0) {
436             Orbit orbit = OrbitType.KEPLERIAN.convertType(propagator.propagate(date.shiftedBy(dt)).getOrbit());
437             if (!useDerivatives) {
438                 // remove derivatives
439                 double[] stateVector = new double[6];
440                 orbit.getType().mapOrbitToArray(orbit, PositionAngleType.TRUE, stateVector, null);
441                 orbit = orbit.getType().mapArrayToOrbit(stateVector, null, PositionAngleType.TRUE,
442                                                         orbit.getDate(), orbit.getMu(), orbit.getFrame());
443             }
444             sample.add(orbit);
445         }
446 
447         // Create interpolator
448         final TimeInterpolator<Orbit> interpolator =
449                 new OrbitHermiteInterpolator(sample.size(), 359, initialOrbit.getFrame(),
450                                              CartesianDerivativesFilter.USE_PVA);
451 
452         // well inside the sample, interpolation should be slightly better than Keplerian shift
453         // the relative bad behaviour here is due to eccentricity, which cannot be
454         // accurately interpolated with a polynomial in this case
455         double maxShiftPositionError             = 0;
456         double maxInterpolationPositionError     = 0;
457         double maxShiftEccentricityError         = 0;
458         double maxInterpolationEccentricityError = 0;
459         for (double dt = 0; dt < 241.0; dt += 1.0) {
460             AbsoluteDate t             = initialOrbit.getDate().shiftedBy(dt);
461             Vector3D     shiftedP      = initialOrbit.shiftedBy(dt).getPosition();
462             Vector3D     interpolatedP = interpolator.interpolate(t, sample).getPosition();
463             Vector3D     propagatedP   = propagator.propagate(t).getPosition();
464             double       shiftedE      = initialOrbit.shiftedBy(dt).getE();
465             double       interpolatedE = interpolator.interpolate(t, sample).getE();
466             double       propagatedE   = propagator.propagate(t).getOrbit().getE();
467             maxShiftPositionError             =
468                     FastMath.max(maxShiftPositionError, shiftedP.subtract(propagatedP).getNorm());
469             maxInterpolationPositionError     =
470                     FastMath.max(maxInterpolationPositionError, interpolatedP.subtract(propagatedP).getNorm());
471             maxShiftEccentricityError         =
472                     FastMath.max(maxShiftEccentricityError, FastMath.abs(shiftedE - propagatedE));
473             maxInterpolationEccentricityError =
474                     FastMath.max(maxInterpolationEccentricityError, FastMath.abs(interpolatedE - propagatedE));
475         }
476         Assertions.assertEquals(shiftPositionErrorWithin, maxShiftPositionError, 0.01 * shiftPositionErrorWithin);
477         Assertions.assertEquals(interpolationPositionErrorWithin, maxInterpolationPositionError,
478                                 0.01 * interpolationPositionErrorWithin);
479         Assertions.assertEquals(shiftEccentricityErrorWithin, maxShiftEccentricityError,
480                                 0.01 * shiftEccentricityErrorWithin);
481         Assertions.assertEquals(interpolationEccentricityErrorWithin, maxInterpolationEccentricityError,
482                                 0.01 * interpolationEccentricityErrorWithin);
483 
484         // slightly past sample end, bad eccentricity interpolation shows up
485         // (in this case, interpolated eccentricity exceeds 1.0 between 1900
486         // and 1910s, while semi-major axis remains positive, so this is not
487         // even a proper hyperbolic orbit...)
488         maxShiftPositionError             = 0;
489         maxInterpolationPositionError     = 0;
490         maxShiftEccentricityError         = 0;
491         maxInterpolationEccentricityError = 0;
492         for (double dt = 240; dt < 600; dt += 1.0) {
493             AbsoluteDate t             = initialOrbit.getDate().shiftedBy(dt);
494             Vector3D     shiftedP      = initialOrbit.shiftedBy(dt).getPosition();
495             Vector3D     interpolatedP = interpolator.interpolate(t, sample).getPosition();
496             Vector3D     propagatedP   = propagator.propagate(t).getPosition();
497             double       shiftedE      = initialOrbit.shiftedBy(dt).getE();
498             double       interpolatedE = interpolator.interpolate(t, sample).getE();
499             double       propagatedE   = propagator.propagate(t).getOrbit().getE();
500             maxShiftPositionError             =
501                     FastMath.max(maxShiftPositionError, shiftedP.subtract(propagatedP).getNorm());
502             maxInterpolationPositionError     =
503                     FastMath.max(maxInterpolationPositionError, interpolatedP.subtract(propagatedP).getNorm());
504             maxShiftEccentricityError         =
505                     FastMath.max(maxShiftEccentricityError, FastMath.abs(shiftedE - propagatedE));
506             maxInterpolationEccentricityError =
507                     FastMath.max(maxInterpolationEccentricityError, FastMath.abs(interpolatedE - propagatedE));
508         }
509         Assertions.assertEquals(shiftPositionErrorSlightlyPast, maxShiftPositionError,
510                                 0.01 * shiftPositionErrorSlightlyPast);
511         Assertions.assertEquals(interpolationPositionErrorSlightlyPast, maxInterpolationPositionError,
512                                 0.01 * interpolationPositionErrorSlightlyPast);
513         Assertions.assertEquals(shiftEccentricityErrorSlightlyPast, maxShiftEccentricityError,
514                                 0.01 * shiftEccentricityErrorSlightlyPast);
515         Assertions.assertEquals(interpolationEccentricityErrorSlightlyPast, maxInterpolationEccentricityError,
516                                 0.01 * interpolationEccentricityErrorSlightlyPast);
517 
518     }
519 
520     @Test
521     @DisplayName("test default constructor")
522     void testDefaultConstructor() {
523         // Given
524         final int interpolationPoints = 2;
525 
526         // When
527         final OrbitHermiteInterpolator interpolator =
528                 new OrbitHermiteInterpolator(interpolationPoints, FramesFactory.getGCRF());
529 
530         // Then
531         Assertions.assertEquals(AbstractTimeInterpolator.DEFAULT_EXTRAPOLATION_THRESHOLD_SEC,
532                                 interpolator.getExtrapolationThreshold());
533         Assertions.assertEquals(interpolationPoints,
534                                 interpolator.getNbInterpolationPoints());
535         Assertions.assertEquals(CartesianDerivativesFilter.USE_PVA, interpolator.getPVAFilter());
536     }
537 
538 }