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.TestUtils;
29  import org.orekit.Utils;
30  import org.orekit.errors.OrekitIllegalArgumentException;
31  import org.orekit.frames.Frame;
32  import org.orekit.frames.FramesFactory;
33  import org.orekit.propagation.analytical.EcksteinHechlerPropagator;
34  import org.orekit.time.AbsoluteDate;
35  import org.orekit.time.AbstractTimeInterpolator;
36  import org.orekit.time.TimeInterpolator;
37  import org.orekit.utils.CartesianDerivativesFilter;
38  import org.orekit.utils.PVCoordinates;
39  
40  class OrbitHermiteInterpolatorTest {
41  
42      @BeforeEach
43      public void setUp() {
44          Utils.setDataRoot("regular-data");
45      }
46  
47      @Test
48      public void testCartesianInterpolationWithDerivatives() {
49          doTestCartesianInterpolation(true, CartesianDerivativesFilter.USE_P,
50                                       394, 0.1968, 3.21, 0.021630,
51                                       2474, 2707.6418, 6.6, 26.28);
52  
53          doTestCartesianInterpolation(true, CartesianDerivativesFilter.USE_PV,
54                                       394, 8.2392E-9, 3.21, 1.3236E-10,
55                                       2474, 0.07474, 6.6, 0.001450);
56  
57          // Solution with PVA less precise than with PV only as the interpolating polynomial begins to oscillate heavily
58          // outside the interpolating interval
59          doTestCartesianInterpolation(true, CartesianDerivativesFilter.USE_PVA,
60                                       394, 9.60e-9, 3.21, 5.39e-10,
61                                       2474, 2614, 6.55, 71);
62  
63      }
64  
65      @Test
66      public void testCartesianInterpolationWithoutDerivatives() {
67          doTestCartesianInterpolation(false, CartesianDerivativesFilter.USE_P,
68                                       394, 0.1968, 3.21, 0.02163,
69                                       2474, 2707.6419, 6.55, 26.2826);
70  
71          doTestCartesianInterpolation(false, CartesianDerivativesFilter.USE_PV,
72                                       394, 8.2392E-9, 3.21, 1.3236E-10,
73                                       2474, 0.07474, 6.55, 0.001450);
74  
75          // Interpolation without derivatives is very wrong in PVA as we force first and second derivatives to be 0 i.e. we
76          // give false information to the interpolator
77          doTestCartesianInterpolation(false, CartesianDerivativesFilter.USE_PVA,
78                                       394, 2.61, 3.21, 0.154,
79                                       2474, 2.28e12, 6.55, 6.22e10);
80      }
81  
82      @Test
83      public void testCircularInterpolationWithDerivatives() {
84          doTestCircularInterpolation(true,
85                                      397, 2.26e-8,
86                                      610, 3.99e-6,
87                                      4870, 113.8);
88      }
89  
90      @Test
91      public void testCircularInterpolationWithoutDerivatives() {
92          doTestCircularInterpolation(false,
93                                      397, 0.0372,
94                                      610.0, 1.23,
95                                      4870, 8869);
96      }
97  
98      @Test
99      public void testEquinoctialInterpolationWithDerivatives() {
100         doTestEquinoctialInterpolation(true,
101                                        397, 9.84e-9,
102                                        610, 4.28e-6,
103                                        4870, 115);
104     }
105 
106     @Test
107     public void testEquinoctialInterpolationWithoutDerivatives() {
108         doTestEquinoctialInterpolation(false,
109                                        397, 0.0372,
110                                        610.0, 1.23,
111                                        4879, 8871);
112     }
113 
114     @Test
115     public void testKeplerianInterpolationWithDerivatives() {
116         doTestKeplerianInterpolation(true,
117                                      397, 4.01, 4.75e-4, 1.28e-7,
118                                      2159, 1.05e7, 1.19e-3, 0.773);
119     }
120 
121     @Test
122     public void testKeplerianInterpolationWithoutDerivatives() {
123         doTestKeplerianInterpolation(false,
124                                      397, 62.0, 4.75e-4, 2.87e-6,
125                                      2159, 79365, 1.19e-3, 3.89e-3);
126     }
127 
128     private void doTestCartesianInterpolation(boolean useDerivatives, CartesianDerivativesFilter pvaFilter,
129                                               double shiftPositionErrorWithin, double interpolationPositionErrorWithin,
130                                               double shiftVelocityErrorWithin, double interpolationVelocityErrorWithin,
131                                               double shiftPositionErrorFarPast, double interpolationPositionErrorFarPast,
132                                               double shiftVelocityErrorFarPast, double interpolationVelocityErrorFarPast) {
133 
134         final double ehMu = 3.9860047e14;
135         final double ae   = 6.378137e6;
136         final double c20  = -1.08263e-3;
137         final double c30  = 2.54e-6;
138         final double c40  = 1.62e-6;
139         final double c50  = 2.3e-7;
140         final double c60  = -5.5e-7;
141 
142         final AbsoluteDate date     = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
143         final Vector3D     position = new Vector3D(3220103., 69623., 6449822.);
144         final Vector3D     velocity = new Vector3D(6414.7, -2006., -3180.);
145         final CartesianOrbit initialOrbit = new CartesianOrbit(new PVCoordinates(position, velocity),
146                                                                FramesFactory.getEME2000(), date, ehMu);
147 
148         EcksteinHechlerPropagator propagator =
149                 new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
150 
151         // set up a 5 points sample
152         List<Orbit> sample = new ArrayList<>();
153         for (double dt = 0; dt < 251.0; dt += 60) {
154             Orbit orbit = propagator.propagate(date.shiftedBy(dt)).getOrbit();
155             if (!useDerivatives) {
156                 // remove derivatives
157                 double[] stateVector = new double[6];
158                 orbit.getType().mapOrbitToArray(orbit, PositionAngleType.TRUE, stateVector, null);
159                 orbit = orbit.getType().mapArrayToOrbit(stateVector, null, PositionAngleType.TRUE,
160                                                         orbit.getDate(), orbit.getMu(), orbit.getFrame());
161             }
162             sample.add(orbit);
163         }
164 
165         // Create interpolator
166         final TimeInterpolator<Orbit> interpolator =
167                 new OrbitHermiteInterpolator(sample.size(), 409, initialOrbit.getFrame(), pvaFilter);
168 
169         // well inside the sample, interpolation should be much better than Keplerian shift
170         // this is because we take the full non-Keplerian acceleration into account in
171         // the Cartesian parameters, which in this case is preserved by the
172         // Eckstein-Hechler propagator
173         double maxShiftPError         = 0;
174         double maxInterpolationPError = 0;
175         double maxShiftVError         = 0;
176         double maxInterpolationVError = 0;
177 
178         for (double dt = 0; dt < 240.0; dt += 1.0) {
179             AbsoluteDate  t            = date.shiftedBy(dt);
180             PVCoordinates propagated   = propagator.propagate(t).getPVCoordinates();
181             PVCoordinates interpolated = interpolator.interpolate(t, sample).getPVCoordinates();
182 
183             PVCoordinates shiftError = new PVCoordinates(propagated,
184                                                          initialOrbit.shiftedBy(dt).getPVCoordinates());
185             PVCoordinates interpolationError = new PVCoordinates(propagated, interpolated);
186 
187             maxShiftPError         = FastMath.max(maxShiftPError,
188                                                   shiftError.getPosition().getNorm());
189             maxInterpolationPError = FastMath.max(maxInterpolationPError,
190                                                   interpolationError.getPosition().getNorm());
191             maxShiftVError         = FastMath.max(maxShiftVError,
192                                                   shiftError.getVelocity().getNorm());
193             maxInterpolationVError = FastMath.max(maxInterpolationVError,
194                                                   interpolationError.getVelocity().getNorm());
195         }
196 
197         Assertions.assertEquals(shiftPositionErrorWithin, maxShiftPError, 0.01 * shiftPositionErrorWithin);
198         Assertions.assertEquals(interpolationPositionErrorWithin, maxInterpolationPError,
199                                 0.01 * interpolationPositionErrorWithin);
200         Assertions.assertEquals(shiftVelocityErrorWithin, maxShiftVError, 0.01 * shiftVelocityErrorWithin);
201         Assertions.assertEquals(interpolationVelocityErrorWithin, maxInterpolationVError,
202                                 0.01 * interpolationVelocityErrorWithin);
203 
204         // if we go far past sample end, interpolation becomes worse than Keplerian shift
205         maxShiftPError         = 0;
206         maxInterpolationPError = 0;
207         maxShiftVError         = 0;
208         maxInterpolationVError = 0;
209         for (double dt = 500.0; dt < 650.0; dt += 1.0) {
210             AbsoluteDate  t            = initialOrbit.getDate().shiftedBy(dt);
211             PVCoordinates propagated   = propagator.propagate(t).getPVCoordinates();
212             PVCoordinates interpolated = interpolator.interpolate(t, sample).getPVCoordinates();
213 
214             PVCoordinates shiftError = new PVCoordinates(propagated,
215                                                          initialOrbit.shiftedBy(dt).getPVCoordinates());
216             PVCoordinates interpolationError = new PVCoordinates(propagated, interpolated);
217 
218             maxShiftPError         = FastMath.max(maxShiftPError,
219                                                   shiftError.getPosition().getNorm());
220             maxInterpolationPError = FastMath.max(maxInterpolationPError,
221                                                   interpolationError.getPosition().getNorm());
222             maxShiftVError         = FastMath.max(maxShiftVError,
223                                                   shiftError.getVelocity().getNorm());
224             maxInterpolationVError = FastMath.max(maxInterpolationVError,
225                                                   interpolationError.getVelocity().getNorm());
226 
227         }
228 
229         Assertions.assertEquals(shiftPositionErrorFarPast, maxShiftPError, 0.01 * shiftPositionErrorFarPast);
230         Assertions.assertEquals(interpolationPositionErrorFarPast, maxInterpolationPError,
231                                 0.01 * interpolationPositionErrorFarPast);
232         Assertions.assertEquals(shiftVelocityErrorFarPast, maxShiftVError, 0.01 * shiftVelocityErrorFarPast);
233         Assertions.assertEquals(interpolationVelocityErrorFarPast, maxInterpolationVError,
234                                 0.01 * interpolationVelocityErrorFarPast);
235     }
236 
237     private void doTestCircularInterpolation(boolean useDerivatives,
238                                              double shiftErrorWithin, double interpolationErrorWithin,
239                                              double shiftErrorSlightlyPast, double interpolationErrorSlightlyPast,
240                                              double shiftErrorFarPast, double interpolationErrorFarPast) {
241 
242         final double ehMu = 3.9860047e14;
243         final double ae   = 6.378137e6;
244         final double c20  = -1.08263e-3;
245         final double c30  = 2.54e-6;
246         final double c40  = 1.62e-6;
247         final double c50  = 2.3e-7;
248         final double c60  = -5.5e-7;
249 
250         final AbsoluteDate date     = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
251         final Vector3D     position = new Vector3D(3220103., 69623., 6449822.);
252         final Vector3D     velocity = new Vector3D(6414.7, -2006., -3180.);
253         final CircularOrbit initialOrbit = new CircularOrbit(new PVCoordinates(position, velocity),
254                                                              FramesFactory.getEME2000(), date, ehMu);
255 
256         EcksteinHechlerPropagator propagator =
257                 new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
258 
259         // set up a 5 points sample
260         List<Orbit> sample = new ArrayList<>();
261         for (double dt = 0; dt < 300.0; dt += 60.0) {
262             Orbit orbit = OrbitType.CIRCULAR.convertType(propagator.propagate(date.shiftedBy(dt)).getOrbit());
263             if (!useDerivatives) {
264                 // remove derivatives
265                 double[] stateVector = new double[6];
266                 orbit.getType().mapOrbitToArray(orbit, PositionAngleType.TRUE, stateVector, null);
267                 orbit = orbit.getType().mapArrayToOrbit(stateVector, null, PositionAngleType.TRUE,
268                                                         orbit.getDate(), orbit.getMu(), orbit.getFrame());
269             }
270             sample.add(orbit);
271         }
272 
273         // Create interpolator
274         final TimeInterpolator<Orbit> interpolator =
275                 new OrbitHermiteInterpolator(sample.size(), 759, initialOrbit.getFrame(),
276                                              CartesianDerivativesFilter.USE_PVA);
277 
278         // well inside the sample, interpolation should be much better than Keplerian shift
279         double maxShiftError         = 0;
280         double maxInterpolationError = 0;
281         for (double dt = 0; dt < 241.0; dt += 1.0) {
282             AbsoluteDate t            = initialOrbit.getDate().shiftedBy(dt);
283             Vector3D     shifted      = initialOrbit.shiftedBy(dt).getPosition();
284             Vector3D     interpolated = interpolator.interpolate(t, sample).getPosition();
285             Vector3D     propagated   = propagator.propagate(t).getPosition();
286             maxShiftError         = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm());
287             maxInterpolationError = FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm());
288         }
289         Assertions.assertEquals(shiftErrorWithin, maxShiftError, 0.01 * shiftErrorWithin);
290         Assertions.assertEquals(interpolationErrorWithin, maxInterpolationError, 0.01 * interpolationErrorWithin);
291 
292         // slightly past sample end, interpolation should quickly increase, but remain reasonable
293         maxShiftError         = 0;
294         maxInterpolationError = 0;
295         for (double dt = 240; dt < 300.0; dt += 1.0) {
296             AbsoluteDate t            = initialOrbit.getDate().shiftedBy(dt);
297             Vector3D     shifted      = initialOrbit.shiftedBy(dt).getPosition();
298             Vector3D     interpolated = interpolator.interpolate(t, sample).getPosition();
299             Vector3D     propagated   = propagator.propagate(t).getPosition();
300             maxShiftError         = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm());
301             maxInterpolationError = FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm());
302         }
303         Assertions.assertEquals(shiftErrorSlightlyPast, maxShiftError, 0.01 * shiftErrorSlightlyPast);
304         Assertions.assertEquals(interpolationErrorSlightlyPast, maxInterpolationError,
305                                 0.01 * interpolationErrorSlightlyPast);
306 
307         // far past sample end, interpolation should become really wrong
308         maxShiftError         = 0;
309         maxInterpolationError = 0;
310         for (double dt = 300; dt < 1000; dt += 1.0) {
311             AbsoluteDate t            = initialOrbit.getDate().shiftedBy(dt);
312             Vector3D     shifted      = initialOrbit.shiftedBy(dt).getPosition();
313             Vector3D     interpolated = interpolator.interpolate(t, sample).getPosition();
314             Vector3D     propagated   = propagator.propagate(t).getPosition();
315             maxShiftError         = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm());
316             maxInterpolationError = FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm());
317         }
318         Assertions.assertEquals(shiftErrorFarPast, maxShiftError, 0.01 * shiftErrorFarPast);
319         Assertions.assertEquals(interpolationErrorFarPast, maxInterpolationError, 0.01 * interpolationErrorFarPast);
320 
321     }
322 
323     private void doTestEquinoctialInterpolation(boolean useDerivatives,
324                                                 double shiftErrorWithin, double interpolationErrorWithin,
325                                                 double shiftErrorSlightlyPast, double interpolationErrorSlightlyPast,
326                                                 double shiftErrorFarPast, double interpolationErrorFarPast) {
327 
328         final double ehMu = 3.9860047e14;
329         final double ae   = 6.378137e6;
330         final double c20  = -1.08263e-3;
331         final double c30  = 2.54e-6;
332         final double c40  = 1.62e-6;
333         final double c50  = 2.3e-7;
334         final double c60  = -5.5e-7;
335 
336         final AbsoluteDate date     = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
337         final Vector3D     position = new Vector3D(3220103., 69623., 6449822.);
338         final Vector3D     velocity = new Vector3D(6414.7, -2006., -3180.);
339         final EquinoctialOrbit initialOrbit = new EquinoctialOrbit(new PVCoordinates(position, velocity),
340                                                                    FramesFactory.getEME2000(), date, ehMu);
341 
342         EcksteinHechlerPropagator propagator =
343                 new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
344 
345         // set up a 5 points sample
346         List<Orbit> sample = new ArrayList<>();
347         for (double dt = 0; dt < 300.0; dt += 60.0) {
348             Orbit orbit = OrbitType.EQUINOCTIAL.convertType(propagator.propagate(date.shiftedBy(dt)).getOrbit());
349             if (!useDerivatives) {
350                 // remove derivatives
351                 double[] stateVector = new double[6];
352                 orbit.getType().mapOrbitToArray(orbit, PositionAngleType.TRUE, stateVector, null);
353                 orbit = orbit.getType().mapArrayToOrbit(stateVector, null, PositionAngleType.TRUE,
354                                                         orbit.getDate(), orbit.getMu(), orbit.getFrame());
355             }
356             sample.add(orbit);
357         }
358 
359         // Create interpolator
360         final TimeInterpolator<Orbit> interpolator =
361                 new OrbitHermiteInterpolator(sample.size(), 759, initialOrbit.getFrame(),
362                                              CartesianDerivativesFilter.USE_PVA);
363 
364         // well inside the sample, interpolation should be much better than Keplerian shift
365         double maxShiftError         = 0;
366         double maxInterpolationError = 0;
367         for (double dt = 0; dt < 241.0; dt += 1.0) {
368             AbsoluteDate t            = initialOrbit.getDate().shiftedBy(dt);
369             Vector3D     shifted      = initialOrbit.shiftedBy(dt).getPosition();
370             Vector3D     interpolated = interpolator.interpolate(t, sample).getPosition();
371             Vector3D     propagated   = propagator.propagate(t).getPosition();
372             maxShiftError         = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm());
373             maxInterpolationError = FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm());
374         }
375         Assertions.assertEquals(shiftErrorWithin, maxShiftError, 0.01 * shiftErrorWithin);
376         Assertions.assertEquals(interpolationErrorWithin, maxInterpolationError, 0.01 * interpolationErrorWithin);
377 
378         // slightly past sample end, interpolation should quickly increase, but remain reasonable
379         maxShiftError         = 0;
380         maxInterpolationError = 0;
381         for (double dt = 240; dt < 300.0; dt += 1.0) {
382             AbsoluteDate t            = initialOrbit.getDate().shiftedBy(dt);
383             Vector3D     shifted      = initialOrbit.shiftedBy(dt).getPosition();
384             Vector3D     interpolated = interpolator.interpolate(t, sample).getPosition();
385             Vector3D     propagated   = propagator.propagate(t).getPosition();
386             maxShiftError         = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm());
387             maxInterpolationError = FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm());
388         }
389         Assertions.assertEquals(shiftErrorSlightlyPast, maxShiftError, 0.01 * shiftErrorSlightlyPast);
390         Assertions.assertEquals(interpolationErrorSlightlyPast, maxInterpolationError,
391                                 0.01 * interpolationErrorSlightlyPast);
392 
393         // far past sample end, interpolation should become really wrong
394         // (in this test case, break even occurs at around 863 seconds, with a 3.9 km error)
395         maxShiftError         = 0;
396         maxInterpolationError = 0;
397         for (double dt = 300; dt < 1000; dt += 1.0) {
398             AbsoluteDate t            = initialOrbit.getDate().shiftedBy(dt);
399             Vector3D     shifted      = initialOrbit.shiftedBy(dt).getPosition();
400             Vector3D     interpolated = interpolator.interpolate(t, sample).getPosition();
401             Vector3D     propagated   = propagator.propagate(t).getPosition();
402             maxShiftError         = FastMath.max(maxShiftError, shifted.subtract(propagated).getNorm());
403             maxInterpolationError = FastMath.max(maxInterpolationError, interpolated.subtract(propagated).getNorm());
404         }
405         Assertions.assertEquals(shiftErrorFarPast, maxShiftError, 0.01 * shiftErrorFarPast);
406         Assertions.assertEquals(interpolationErrorFarPast, maxInterpolationError, 0.01 * interpolationErrorFarPast);
407 
408     }
409 
410     private void doTestKeplerianInterpolation(boolean useDerivatives,
411                                               double shiftPositionErrorWithin, double interpolationPositionErrorWithin,
412                                               double shiftEccentricityErrorWithin,
413                                               double interpolationEccentricityErrorWithin,
414                                               double shiftPositionErrorSlightlyPast,
415                                               double interpolationPositionErrorSlightlyPast,
416                                               double shiftEccentricityErrorSlightlyPast,
417                                               double interpolationEccentricityErrorSlightlyPast) {
418 
419         final double ehMu = 3.9860047e14;
420         final double ae   = 6.378137e6;
421         final double c20  = -1.08263e-3;
422         final double c30  = 2.54e-6;
423         final double c40  = 1.62e-6;
424         final double c50  = 2.3e-7;
425         final double c60  = -5.5e-7;
426 
427         final AbsoluteDate date     = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
428         final Vector3D     position = new Vector3D(3220103., 69623., 6449822.);
429         final Vector3D     velocity = new Vector3D(6414.7, -2006., -3180.);
430         final KeplerianOrbit initialOrbit = new KeplerianOrbit(new PVCoordinates(position, velocity),
431                                                                FramesFactory.getEME2000(), date, ehMu);
432 
433         EcksteinHechlerPropagator propagator =
434                 new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
435 
436         // set up a 5 points sample
437         List<Orbit> sample = new ArrayList<>();
438         for (double dt = 0; dt < 300.0; dt += 60.0) {
439             Orbit orbit = OrbitType.KEPLERIAN.convertType(propagator.propagate(date.shiftedBy(dt)).getOrbit());
440             if (!useDerivatives) {
441                 // remove derivatives
442                 double[] stateVector = new double[6];
443                 orbit.getType().mapOrbitToArray(orbit, PositionAngleType.TRUE, stateVector, null);
444                 orbit = orbit.getType().mapArrayToOrbit(stateVector, null, PositionAngleType.TRUE,
445                                                         orbit.getDate(), orbit.getMu(), orbit.getFrame());
446             }
447             sample.add(orbit);
448         }
449 
450         // Create interpolator
451         final TimeInterpolator<Orbit> interpolator =
452                 new OrbitHermiteInterpolator(sample.size(), 359, initialOrbit.getFrame(),
453                                              CartesianDerivativesFilter.USE_PVA);
454 
455         // well inside the sample, interpolation should be slightly better than Keplerian shift
456         // the relative bad behaviour here is due to eccentricity, which cannot be
457         // accurately interpolated with a polynomial in this case
458         double maxShiftPositionError             = 0;
459         double maxInterpolationPositionError     = 0;
460         double maxShiftEccentricityError         = 0;
461         double maxInterpolationEccentricityError = 0;
462         for (double dt = 0; dt < 241.0; dt += 1.0) {
463             AbsoluteDate t             = initialOrbit.getDate().shiftedBy(dt);
464             Vector3D     shiftedP      = initialOrbit.shiftedBy(dt).getPosition();
465             Vector3D     interpolatedP = interpolator.interpolate(t, sample).getPosition();
466             Vector3D     propagatedP   = propagator.propagate(t).getPosition();
467             double       shiftedE      = initialOrbit.shiftedBy(dt).getE();
468             double       interpolatedE = interpolator.interpolate(t, sample).getE();
469             double       propagatedE   = propagator.propagate(t).getOrbit().getE();
470             maxShiftPositionError             =
471                     FastMath.max(maxShiftPositionError, shiftedP.subtract(propagatedP).getNorm());
472             maxInterpolationPositionError     =
473                     FastMath.max(maxInterpolationPositionError, interpolatedP.subtract(propagatedP).getNorm());
474             maxShiftEccentricityError         =
475                     FastMath.max(maxShiftEccentricityError, FastMath.abs(shiftedE - propagatedE));
476             maxInterpolationEccentricityError =
477                     FastMath.max(maxInterpolationEccentricityError, FastMath.abs(interpolatedE - propagatedE));
478         }
479         Assertions.assertEquals(shiftPositionErrorWithin, maxShiftPositionError, 0.01 * shiftPositionErrorWithin);
480         Assertions.assertEquals(interpolationPositionErrorWithin, maxInterpolationPositionError,
481                                 0.01 * interpolationPositionErrorWithin);
482         Assertions.assertEquals(shiftEccentricityErrorWithin, maxShiftEccentricityError,
483                                 0.01 * shiftEccentricityErrorWithin);
484         Assertions.assertEquals(interpolationEccentricityErrorWithin, maxInterpolationEccentricityError,
485                                 0.01 * interpolationEccentricityErrorWithin);
486 
487         // slightly past sample end, bad eccentricity interpolation shows up
488         // (in this case, interpolated eccentricity exceeds 1.0 between 1900
489         // and 1910s, while semi-major axis remains positive, so this is not
490         // even a proper hyperbolic orbit...)
491         maxShiftPositionError             = 0;
492         maxInterpolationPositionError     = 0;
493         maxShiftEccentricityError         = 0;
494         maxInterpolationEccentricityError = 0;
495         for (double dt = 240; dt < 600; dt += 1.0) {
496             AbsoluteDate t             = initialOrbit.getDate().shiftedBy(dt);
497             Vector3D     shiftedP      = initialOrbit.shiftedBy(dt).getPosition();
498             Vector3D     interpolatedP = interpolator.interpolate(t, sample).getPosition();
499             Vector3D     propagatedP   = propagator.propagate(t).getPosition();
500             double       shiftedE      = initialOrbit.shiftedBy(dt).getE();
501             double       interpolatedE = interpolator.interpolate(t, sample).getE();
502             double       propagatedE   = propagator.propagate(t).getOrbit().getE();
503             maxShiftPositionError             =
504                     FastMath.max(maxShiftPositionError, shiftedP.subtract(propagatedP).getNorm());
505             maxInterpolationPositionError     =
506                     FastMath.max(maxInterpolationPositionError, interpolatedP.subtract(propagatedP).getNorm());
507             maxShiftEccentricityError         =
508                     FastMath.max(maxShiftEccentricityError, FastMath.abs(shiftedE - propagatedE));
509             maxInterpolationEccentricityError =
510                     FastMath.max(maxInterpolationEccentricityError, FastMath.abs(interpolatedE - propagatedE));
511         }
512         Assertions.assertEquals(shiftPositionErrorSlightlyPast, maxShiftPositionError,
513                                 0.01 * shiftPositionErrorSlightlyPast);
514         Assertions.assertEquals(interpolationPositionErrorSlightlyPast, maxInterpolationPositionError,
515                                 0.01 * interpolationPositionErrorSlightlyPast);
516         Assertions.assertEquals(shiftEccentricityErrorSlightlyPast, maxShiftEccentricityError,
517                                 0.01 * shiftEccentricityErrorSlightlyPast);
518         Assertions.assertEquals(interpolationEccentricityErrorSlightlyPast, maxInterpolationEccentricityError,
519                                 0.01 * interpolationEccentricityErrorSlightlyPast);
520 
521     }
522 
523     @Test
524     @DisplayName("test default constructor")
525     void testDefaultConstructor() {
526         // Given
527         final int interpolationPoints = 2;
528 
529         // When
530         final OrbitHermiteInterpolator interpolator =
531                 new OrbitHermiteInterpolator(interpolationPoints, FramesFactory.getGCRF());
532 
533         // Then
534         Assertions.assertEquals(AbstractTimeInterpolator.DEFAULT_EXTRAPOLATION_THRESHOLD_SEC,
535                                 interpolator.getExtrapolationThreshold());
536         Assertions.assertEquals(interpolationPoints,
537                                 interpolator.getNbInterpolationPoints());
538         Assertions.assertEquals(CartesianDerivativesFilter.USE_PVA, interpolator.getPVAFilter());
539     }
540 
541    @Test
542    public void testErrorThrownWhenUsingOrbitWithFrameMismatch() {
543       // GIVEN
544       // Define sample with frame mismatch
545       final Orbit orbit1 = TestUtils.getFakeOrbit();
546       final Orbit orbit2 = new CartesianOrbit(orbit1.shiftedBy(1).getPVCoordinates(),
547                                               FramesFactory.getEME2000(),
548                                               orbit1.getMu());
549 
550       final List<Orbit> sample = new ArrayList<>();
551       sample.add(orbit1);
552       sample.add(orbit2);
553 
554       // Define interpolator
555       final OrbitHermiteInterpolator interpolator =
556               new OrbitHermiteInterpolator(FramesFactory.getGCRF());
557 
558       // WHEN & THEN
559       Assertions.assertThrows(OrekitIllegalArgumentException.class,
560                               () -> interpolator.interpolate(new AbsoluteDate(), sample));
561    }
562 
563     /**
564      * Test related to issue 1844 with Cartesian orbits.
565      *
566      * @see <a href="https://gitlab.orekit.org/orekit/orekit/-/issues/1844">Issue 1844</a>
567      */
568     @Test
569     void testOutputFrameCartesian() {
570         // GIVEN
571         // Define samples
572         final Orbit orbit1 = TestUtils.getFakeOrbit();
573         final Orbit orbit2 = orbit1.shiftedBy(1);
574 
575         final List<Orbit> samples = new ArrayList<>();
576         samples.add(orbit1);
577         samples.add(orbit2);
578 
579         // Define output frame
580         final Frame outputFrame = FramesFactory.getEME2000();
581 
582         // Define interpolator
583         final OrbitHermiteInterpolator interpolator = new OrbitHermiteInterpolator(outputFrame);
584 
585         // Define interpolation date
586         final AbsoluteDate interpolationDate = orbit1.getDate();
587 
588         // WHEN
589         final Orbit interpolatedOrbit = interpolator.interpolate(interpolationDate, samples);
590 
591         // THEN
592         // Assert against original orbit
593         final Vector3D originalPosition = orbit1.getPosition();
594         final Vector3D actualPosition   = interpolatedOrbit.getPosition(orbit1.getFrame());
595         Assertions.assertEquals(originalPosition.getX(), actualPosition.getX(), 1.0e-14);
596         Assertions.assertEquals(originalPosition.getY(), actualPosition.getY(), 1.0e-14);
597         Assertions.assertEquals(originalPosition.getZ(), actualPosition.getZ(), 1.0e-14);
598 
599         // Assert output type
600         Assertions.assertInstanceOf(CartesianOrbit.class, interpolatedOrbit);
601 
602         // Assert frame
603         Assertions.assertEquals(outputFrame, interpolatedOrbit.getFrame());
604     }
605 
606     /**
607      * Test related to issue 1844 with Common orbits (other than Cartesian).
608      *
609      * @see <a href="https://gitlab.orekit.org/orekit/orekit/-/issues/1844">Issue 1844</a>
610      */
611     @Test
612     void testOutputFrameCommon() {
613         // GIVEN
614         // Define samples
615         final Orbit orbit1 = new CircularOrbit(TestUtils.getTestOrbit());
616         final Orbit orbit2 = new CircularOrbit(orbit1.shiftedBy(1));
617 
618         final List<Orbit> samples = new ArrayList<>();
619         samples.add(orbit1);
620         samples.add(orbit2);
621 
622         // Define output frame
623         final Frame outputFrame = FramesFactory.getEME2000();
624 
625         // Define interpolator
626         final OrbitHermiteInterpolator interpolator = new OrbitHermiteInterpolator(outputFrame);
627 
628         // Define interpolation date
629         final AbsoluteDate interpolationDate = orbit1.getDate();
630 
631         // WHEN
632         final Orbit interpolatedOrbit = interpolator.interpolate(interpolationDate, samples);
633 
634         // THEN
635         // Assert against original orbit
636         final Vector3D originalPosition = orbit1.getPosition();
637         final Vector3D actualPosition   = interpolatedOrbit.getPosition(orbit1.getFrame());
638         Assertions.assertEquals(originalPosition.getX(), actualPosition.getX(), 1.0e-8);
639         Assertions.assertEquals(originalPosition.getY(), actualPosition.getY(), 1.0e-8);
640         Assertions.assertEquals(originalPosition.getZ(), actualPosition.getZ(), 1.0e-8);
641     
642         // Assert frame
643         Assertions.assertEquals(outputFrame, interpolatedOrbit.getFrame());
644     }
645 
646 }