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.propagation.analytical;
18  
19  import org.hamcrest.MatcherAssert;
20  import org.hipparchus.analysis.polynomials.SmoothStepFactory;
21  import org.hipparchus.geometry.euclidean.threed.Rotation;
22  import org.hipparchus.geometry.euclidean.threed.Vector3D;
23  import org.hipparchus.linear.BlockRealMatrix;
24  import org.hipparchus.linear.MatrixUtils;
25  import org.hipparchus.linear.RealMatrix;
26  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
27  import org.hipparchus.util.FastMath;
28  import org.junit.jupiter.api.Assertions;
29  import org.junit.jupiter.api.DisplayName;
30  import org.junit.jupiter.api.Test;
31  import org.junit.jupiter.api.Timeout;
32  import org.mockito.Mockito;
33  import org.orekit.OrekitMatchers;
34  import org.orekit.TestUtils;
35  import org.orekit.Utils;
36  import org.orekit.attitudes.*;
37  import org.orekit.bodies.CelestialBodyFactory;
38  import org.orekit.bodies.GeodeticPoint;
39  import org.orekit.bodies.OneAxisEllipsoid;
40  import org.orekit.errors.*;
41  import org.orekit.frames.Frame;
42  import org.orekit.frames.FramesFactory;
43  import org.orekit.frames.LOFType;
44  import org.orekit.frames.TopocentricFrame;
45  import org.orekit.models.earth.EarthStandardAtmosphereRefraction;
46  import org.orekit.models.earth.ReferenceEllipsoid;
47  import org.orekit.orbits.*;
48  import org.orekit.propagation.*;
49  import org.orekit.propagation.events.DateDetector;
50  import org.orekit.propagation.events.ElevationDetector;
51  import org.orekit.propagation.events.handlers.ContinueOnEvent;
52  import org.orekit.propagation.events.handlers.RecordAndContinue;
53  import org.orekit.propagation.numerical.NumericalPropagator;
54  import org.orekit.time.*;
55  import org.orekit.utils.*;
56  
57  import java.lang.reflect.InvocationTargetException;
58  import java.lang.reflect.Method;
59  import java.util.ArrayList;
60  import java.util.Arrays;
61  import java.util.List;
62  import java.util.Optional;
63  import java.util.concurrent.TimeUnit;
64  import java.util.stream.Stream;
65  
66  public class EphemerisTest {
67  
68      private AbsoluteDate initDate;
69      private AbsoluteDate finalDate;
70      private Frame        inertialFrame;
71      private Propagator   propagator;
72      private Orbit        initialState;
73  
74      @Test
75      @DisplayName("test default Ephemeris constructor")
76      void testDefaultEphemerisConstructor() {
77          // Given
78          // Create spacecraft state sample
79          final Orbit           defaultOrbit = TestUtils.getDefaultOrbit(new AbsoluteDate());
80          final SpacecraftState firstState   = new SpacecraftState(defaultOrbit);
81          final SpacecraftState secondState  = firstState.shiftedBy(1);
82  
83          final List<SpacecraftState> states = new ArrayList<>();
84          states.add(firstState);
85          states.add(secondState);
86  
87          final Ephemeris defaultEphemeris = new Ephemeris(states, states.size());
88  
89          // When
90          final SpacecraftStateInterpolator defaultStateInterpolator =
91                  (SpacecraftStateInterpolator) defaultEphemeris.getStateInterpolator();
92  
93          // Then
94          final Frame                       inertialFrame              = defaultOrbit.getFrame();
95          final SpacecraftStateInterpolator referenceStateInterpolator = new SpacecraftStateInterpolator(inertialFrame);
96  
97          Assertions.assertEquals(referenceStateInterpolator.getExtrapolationThreshold(),
98                                  defaultStateInterpolator.getExtrapolationThreshold());
99          Assertions.assertEquals(referenceStateInterpolator.getOutputFrame(),
100                                 defaultStateInterpolator.getOutputFrame());
101         Assertions.assertEquals(referenceStateInterpolator.getAttitudeInterpolator().isPresent(),
102                                 defaultStateInterpolator.getAttitudeInterpolator().isPresent());
103         Assertions.assertEquals(referenceStateInterpolator.getAdditionalStateInterpolator().isPresent(),
104                                 defaultStateInterpolator.getAdditionalStateInterpolator().isPresent());
105 
106         Assertions.assertInstanceOf(OrbitHermiteInterpolator.class,
107                                     defaultStateInterpolator.getOrbitInterpolator().get());
108         Assertions.assertInstanceOf(TimeStampedDoubleHermiteInterpolator.class,
109                                     defaultStateInterpolator.getMassInterpolator().get());
110         Assertions.assertInstanceOf(AttitudeInterpolator.class,
111                                     defaultStateInterpolator.getAttitudeInterpolator().get());
112         Assertions.assertInstanceOf(TimeStampedDoubleHermiteInterpolator.class,
113                                     defaultStateInterpolator.getAdditionalStateInterpolator().get());
114 
115         Assertions.assertFalse(defaultEphemeris.getCovarianceInterpolator().isPresent());
116 
117     }
118 
119     @Test
120     @DisplayName("test error thrown when using an empty list to create an Ephemeris")
121     void testErrorThrownWhenUsingEmptyStateList() {
122         // Given
123         final SpacecraftState firstState  = Mockito.mock(SpacecraftState.class);
124         final SpacecraftState secondState = Mockito.mock(SpacecraftState.class);
125 
126         Mockito.when(firstState.isOrbitDefined()).thenReturn(true);
127         Mockito.when(secondState.isOrbitDefined()).thenReturn(true);
128 
129         final List<SpacecraftState> states = new ArrayList<>();
130 
131         // Create interpolator
132         final Frame inertialFrameMock = Mockito.mock(Frame.class);
133         Mockito.when(inertialFrameMock.isPseudoInertial()).thenReturn(true);
134 
135         final int nbInterpolationPoints = 3;
136         final TimeInterpolator<SpacecraftState> stateInterpolator =
137                 new SpacecraftStateInterpolator(nbInterpolationPoints, inertialFrameMock, inertialFrameMock);
138 
139         // When & Then
140         OrekitIllegalArgumentException thrown = Assertions.assertThrows(OrekitIllegalArgumentException.class,
141                                                                         () -> new Ephemeris(states, stateInterpolator));
142 
143         Assertions.assertEquals(OrekitMessages.NOT_ENOUGH_DATA, thrown.getSpecifier());
144         Assertions.assertEquals(0, ((Integer) thrown.getParts()[0]).intValue());
145 
146     }
147 
148     @Test
149     @DisplayName("test error thrown when using inconsistent states sample and state interpolator")
150     void testErrorThrownWhenUsingInconsistentStatesSampleAndStateInterpolator() {
151         // Given
152         final SpacecraftState firstState  = Mockito.mock(SpacecraftState.class);
153         final SpacecraftState secondState = Mockito.mock(SpacecraftState.class);
154 
155         Mockito.when(firstState.isOrbitDefined()).thenReturn(true);
156         Mockito.when(secondState.isOrbitDefined()).thenReturn(true);
157 
158         final List<SpacecraftState> states = new ArrayList<>();
159         states.add(firstState);
160         states.add(secondState);
161 
162         // Create interpolator
163         final Frame inertialFrameMock = Mockito.mock(Frame.class);
164         Mockito.when(inertialFrameMock.isPseudoInertial()).thenReturn(true);
165 
166         final int nbInterpolationPoints = 3;
167         final TimeInterpolator<SpacecraftState> stateInterpolator =
168                 new SpacecraftStateInterpolator(nbInterpolationPoints, inertialFrameMock, inertialFrameMock);
169 
170         // When & Then
171         OrekitIllegalArgumentException thrown = Assertions.assertThrows(OrekitIllegalArgumentException.class,
172                                                                         () -> new Ephemeris(states, stateInterpolator));
173 
174         Assertions.assertEquals(OrekitMessages.NOT_ENOUGH_DATA, thrown.getSpecifier());
175         Assertions.assertEquals(2, ((Integer) thrown.getParts()[0]).intValue());
176 
177     }
178 
179     @Test
180     @DisplayName("test error thrown when using inconsistent states and covariances")
181     void testErrorThrownWhenUsingInconsistentStatesAndCovariances() {
182         // Given
183 
184         // Create spacecraft state sample
185         final SpacecraftState firstState  = Mockito.mock(SpacecraftState.class);
186         final SpacecraftState secondState = Mockito.mock(SpacecraftState.class);
187 
188         Mockito.when(firstState.isOrbitDefined()).thenReturn(true);
189         Mockito.when(secondState.isOrbitDefined()).thenReturn(true);
190 
191         final List<SpacecraftState> states = new ArrayList<>();
192         states.add(firstState);
193         states.add(secondState);
194 
195         // Create covariance state sample
196         final StateCovariance firstCovariance  = Mockito.mock(StateCovariance.class);
197         final StateCovariance secondCovariance = Mockito.mock(StateCovariance.class);
198 
199         final List<StateCovariance> covariances = new ArrayList<>();
200         covariances.add(firstCovariance);
201         covariances.add(secondCovariance);
202 
203         // Simulate inconsistency between state and covariance
204         final AbsoluteDate referenceDate = new AbsoluteDate();
205         Mockito.when(firstState.getDate()).thenReturn(referenceDate);
206         Mockito.when(firstCovariance.getDate()).thenReturn(referenceDate.shiftedBy(1));
207 
208         // Create state interpolator
209         final Frame inertialFrameMock = Mockito.mock(Frame.class);
210         Mockito.when(inertialFrameMock.isPseudoInertial()).thenReturn(true);
211 
212         final TimeInterpolator<SpacecraftState> stateInterpolator = new SpacecraftStateInterpolator(inertialFrameMock);
213 
214         // Create covariance interpolator
215         @SuppressWarnings("unchecked")
216         final TimeInterpolator<TimeStampedPair<Orbit, StateCovariance>> covarianceInterpolator =
217                 Mockito.mock(TimeInterpolator.class);
218 
219         // When & Then
220         Exception thrown = Assertions.assertThrows(OrekitIllegalStateException.class,
221                                                    () -> new Ephemeris(states, stateInterpolator,
222                                                                        covariances, covarianceInterpolator));
223 
224         Assertions.assertEquals(
225                 "state date 2000-01-01T11:58:55.816Z does not match its covariance date 2000-01-01T11:58:56.816Z",
226                 thrown.getMessage());
227     }
228 
229     @Test
230     @DisplayName("test error thrown when using different states and covariances sample size ")
231     void testErrorThrownWhenUsingDifferentStatesAndCovariancesSampleSize() {
232         // Given
233 
234         // Create spacecraft state sample
235         final SpacecraftState firstState  = Mockito.mock(SpacecraftState.class);
236         final SpacecraftState secondState = Mockito.mock(SpacecraftState.class);
237 
238         Mockito.when(firstState.isOrbitDefined()).thenReturn(true);
239         Mockito.when(secondState.isOrbitDefined()).thenReturn(true);
240 
241         final List<SpacecraftState> states = new ArrayList<>();
242         states.add(firstState);
243         states.add(secondState);
244 
245         // Create covariance state sample
246         final StateCovariance firstCovariance = Mockito.mock(StateCovariance.class);
247 
248         final List<StateCovariance> covariances = new ArrayList<>();
249         covariances.add(firstCovariance);
250 
251         // Create state interpolator
252         final Frame inertialFrameMock = Mockito.mock(Frame.class);
253         Mockito.when(inertialFrameMock.isPseudoInertial()).thenReturn(true);
254 
255         final TimeInterpolator<SpacecraftState> stateInterpolator = new SpacecraftStateInterpolator(inertialFrameMock);
256 
257         // Create covariance interpolator
258         @SuppressWarnings("unchecked")
259         final TimeInterpolator<TimeStampedPair<Orbit, StateCovariance>> covarianceInterpolator =
260                 Mockito.mock(TimeInterpolator.class);
261 
262         // When & Then
263         Exception thrown = Assertions.assertThrows(OrekitIllegalArgumentException.class,
264                                                    () -> new Ephemeris(states, stateInterpolator,
265                                                                        covariances, covarianceInterpolator));
266 
267         Assertions.assertEquals(
268                 "inconsistent dimensions: 2 != 1", thrown.getMessage());
269     }
270 
271     @Test
272     @DisplayName("test error thrown when trying to reset intermediate state ")
273     void testErrorThrownWhenTryingToResetIntermediateState() {
274         // Given
275         // Create spacecraft state sample
276         final Orbit           defaultOrbit = TestUtils.getDefaultOrbit(new AbsoluteDate());
277         final SpacecraftState firstState   = new SpacecraftState(defaultOrbit);
278         final SpacecraftState secondState  = firstState.shiftedBy(1);
279 
280         final List<SpacecraftState> states = new ArrayList<>();
281         states.add(firstState);
282         states.add(secondState);
283 
284         // Create state interpolator
285         final Frame inertialFrame = FramesFactory.getEME2000();
286 
287         final TimeInterpolator<SpacecraftState> stateInterpolator = new SpacecraftStateInterpolator(inertialFrame);
288 
289         final Ephemeris ephemeris = new Ephemeris(states, stateInterpolator);
290 
291         // When & Then
292         Exception thrown =
293                 Assertions.assertThrows(OrekitException.class, () -> ephemeris.resetIntermediateState(firstState, true));
294 
295         Assertions.assertEquals("reset state not allowed", thrown.getMessage());
296     }
297 
298     @Test
299     @DisplayName("test that an empty optional is returned when getting covariance from a state only ephemeris")
300     void testEmptyCovarianceGetter() {
301         // GIVEN
302         final AbsoluteDate initialDate  = new AbsoluteDate();
303         final Orbit        initialOrbit = TestUtils.getDefaultOrbit(initialDate);
304 
305         // Setup propagator
306         final Orbit                 finalOrbit = initialOrbit.shiftedBy(1);
307         final List<SpacecraftState> states     = new ArrayList<>();
308         states.add(new SpacecraftState(initialOrbit));
309         states.add(new SpacecraftState(finalOrbit));
310 
311         final Ephemeris ephemeris = new Ephemeris(states, 2);
312 
313         // When
314         final Optional<StateCovariance> optionalCovariance = ephemeris.getCovariance(Mockito.mock(AbsoluteDate.class));
315 
316         // Then
317         Assertions.assertFalse(optionalCovariance.isPresent());
318 
319     }
320 
321     @Test
322     @SuppressWarnings("unchecked")
323     @DisplayName("test that the expected spacecraft state is returned when no attitude provider is given at construction")
324     void testExpectedStateReturnedWhenNoAttitudeProvider() {
325         // GIVEN
326         final AbsoluteDate initialDate  = new AbsoluteDate();
327         final Orbit        initialOrbit = TestUtils.getDefaultOrbit(initialDate);
328 
329         // Setup behaviour
330         final AbsoluteDate    propDate          = initialDate.shiftedBy(0.5);
331         final SpacecraftState mockExpectedState = Mockito.mock(SpacecraftState.class);
332 
333         // Setup propagator
334         final Orbit                 finalOrbit = initialOrbit.shiftedBy(1);
335         final List<SpacecraftState> states     = new ArrayList<>();
336 
337         final SpacecraftState initialState = new SpacecraftState(initialOrbit);
338         final SpacecraftState finalState   = new SpacecraftState(finalOrbit);
339 
340         states.add(initialState);
341         states.add(finalState);
342 
343         // Setup mock interpolator
344         final TimeInterpolator<SpacecraftState> mockInterpolator = Mockito.mock(TimeInterpolator.class);
345 
346         Mockito.when(mockInterpolator.getNbInterpolationPoints()).thenReturn(2);
347         Mockito.when(mockInterpolator.getExtrapolationThreshold()).thenReturn(0.001);
348 
349         Mockito.when(mockInterpolator.interpolate(Mockito.eq(initialDate), Mockito.any(Stream.class)))
350                .thenReturn(initialState);
351         Mockito.when(mockInterpolator.interpolate(Mockito.eq(propDate), Mockito.any(Stream.class)))
352                .thenReturn(mockExpectedState);
353 
354         final Ephemeris ephemeris = new Ephemeris(states, mockInterpolator, null);
355 
356         // When
357         final SpacecraftState propState = ephemeris.basicPropagate(propDate);
358 
359         // Then
360         Assertions.assertEquals(mockExpectedState, propState);
361 
362     }
363 
364     @Test
365     public void testAttitudeOverride() throws IllegalArgumentException, OrekitException {
366         setUp();
367 
368         final double positionTolerance = 1e-6;
369         final double velocityTolerance = 1e-5;
370         final double attitudeTolerance = 1e-6;
371 
372         int    numberOfIntervals = 1440;
373         double deltaT            = finalDate.durationFrom(initDate) / ((double) numberOfIntervals);
374 
375         propagator.setAttitudeProvider(new LofOffset(inertialFrame, LOFType.LVLH_CCSDS));
376 
377         List<SpacecraftState> states = new ArrayList<>(numberOfIntervals + 1);
378         for (int j = 0; j <= numberOfIntervals; j++) {
379             states.add(propagator.propagate(initDate.shiftedBy((j * deltaT))));
380         }
381 
382         // Create interpolator
383         final TimeInterpolator<SpacecraftState> interpolator = new SpacecraftStateInterpolator(inertialFrame);
384 
385         // Create ephemeris with attitude override
386         Ephemeris ephemPropagator = new Ephemeris(states, interpolator, new LofOffset(inertialFrame, LOFType.LVLH_CCSDS));
387         Assertions.assertEquals(0, ephemPropagator.getManagedAdditionalData().length);
388 
389         //First test that we got position, velocity and attitude nailed
390         int numberEphemTestIntervals = 2880;
391         deltaT = finalDate.durationFrom(initDate) / ((double) numberEphemTestIntervals);
392         doTestInterpolation(numberEphemTestIntervals, deltaT, ephemPropagator,
393                             positionTolerance, velocityTolerance, attitudeTolerance,
394                             "LVLH_CCSDS");
395 
396         //Now force an override on the attitude and check it against a Keplerian propagator
397         //setup identically to the first but with a different attitude
398         //If override isn't working this will fail.
399         propagator = new KeplerianPropagator(propagator.getInitialState().getOrbit());
400         propagator.setAttitudeProvider(new LofOffset(inertialFrame, LOFType.QSW));
401 
402         ephemPropagator.setAttitudeProvider(new LofOffset(inertialFrame, LOFType.QSW));
403         doTestInterpolation(numberEphemTestIntervals, deltaT, ephemPropagator,
404                             positionTolerance, velocityTolerance, attitudeTolerance,
405                             "QSW");
406     }
407 
408     private void doTestInterpolation(final int numberEphemTestIntervals, final double deltaT,
409                                      final Ephemeris ephemPropagator,
410                                      final double positionTolerance,
411                                      final double velocityTolerance,
412                                      final double attitudeTolerance,
413                                      final String errorMessage) {
414         for (int j = 0; j <= numberEphemTestIntervals; j++) {
415             AbsoluteDate    currentDate   = initDate.shiftedBy(j * deltaT);
416             SpacecraftState ephemState    = ephemPropagator.propagate(currentDate);
417             SpacecraftState keplerState   = propagator.propagate(currentDate);
418             double          positionDelta = calculatePositionDelta(ephemState, keplerState);
419             double          velocityDelta = calculateVelocityDelta(ephemState, keplerState);
420             double          attitudeDelta = calculateAttitudeDelta(ephemState, keplerState);
421             Assertions.assertEquals(0.0, positionDelta, positionTolerance, errorMessage + " Unmatched Position at: " + currentDate);
422             Assertions.assertEquals(0.0, velocityDelta, velocityTolerance, errorMessage + " Unmatched Velocity at: " + currentDate);
423             Assertions.assertEquals(0.0, attitudeDelta, attitudeTolerance, errorMessage + " Unmatched Attitude at: " + currentDate);
424         }
425     }
426 
427     @Test
428     public void testAttitudeSequenceTransition() {
429         setUp();
430 
431         // Initialize the orbit
432         final AbsoluteDate initialDate = new AbsoluteDate(2003, 01, 01, 0, 0, 00.000, TimeScalesFactory.getUTC());
433         final Vector3D     position    = new Vector3D(-39098981.4866597, -15784239.3610601, 78908.2289853595);
434         final Vector3D     velocity    = new Vector3D(1151.00321021175, -2851.14864755189, -2.02133248357321);
435         final Orbit initialOrbit = new KeplerianOrbit(new PVCoordinates(position, velocity),
436                                                       FramesFactory.getGCRF(), initialDate,
437                                                       Constants.WGS84_EARTH_MU);
438         final SpacecraftState initialState = new SpacecraftState(initialOrbit);
439 
440         // Define attitude laws
441         AttitudeProvider before =
442                 new CelestialBodyPointed(FramesFactory.getICRF(), CelestialBodyFactory.getSun(), Vector3D.PLUS_K,
443                                          Vector3D.PLUS_I, Vector3D.PLUS_K);
444         AttitudeProvider after =
445                 new CelestialBodyPointed(FramesFactory.getICRF(), CelestialBodyFactory.getEarth(), Vector3D.PLUS_K,
446                                          Vector3D.PLUS_I, Vector3D.PLUS_K);
447 
448         // Define attitude sequence
449         AbsoluteDate switchDate     = initialDate.shiftedBy(86400.0);
450         double       transitionTime = 600;
451         DateDetector switchDetector = new DateDetector(switchDate).withHandler(new ContinueOnEvent());
452 
453         AttitudesSequence attitudeSequence = new AttitudesSequence();
454         attitudeSequence.resetActiveProvider(before);
455         attitudeSequence.addSwitchingCondition(before, after, switchDetector, true, false, transitionTime,
456                                                AngularDerivativesFilter.USE_RR, null);
457 
458         NumericalPropagator propagator = new NumericalPropagator(new DormandPrince853Integrator(0.1, 500, 1e-9, 1e-9));
459         propagator.setInitialState(initialState);
460 
461         // Propagate and build ephemeris
462         final List<SpacecraftState> propagatedStates = new ArrayList<>();
463 
464         propagator.setStepHandler(60, propagatedStates::add);
465         propagator.propagate(initialDate.shiftedBy(2 * 86400.0));
466 
467         // Create interpolator
468         final int nbInterpolationPoints = 8;
469         final TimeInterpolator<SpacecraftState> interpolator =
470                 new SpacecraftStateInterpolator(nbInterpolationPoints, inertialFrame, inertialFrame);
471 
472         final Ephemeris ephemeris = new Ephemeris(propagatedStates, interpolator);
473 
474         // Add attitude switch event to ephemeris
475         ephemeris.setAttitudeProvider(attitudeSequence);
476 
477         // Propagate with a step during the transition
478         AbsoluteDate    endDate     = initialDate.shiftedBy(2 * 86400.0);
479         SpacecraftState stateBefore = ephemeris.getInitialState();
480         ephemeris.propagate(switchDate.shiftedBy(transitionTime / 2));
481         SpacecraftState stateAfter = ephemeris.propagate(endDate);
482 
483         // Check that the attitudes are correct
484         Assertions.assertEquals(
485                 before.getAttitude(stateBefore.getOrbit(), stateBefore.getDate(), stateBefore.getFrame()).getRotation()
486                       .getQ0(),
487                 stateBefore.getAttitude().getRotation().getQ0(),
488                 1.0E-16);
489         Assertions.assertEquals(
490                 before.getAttitude(stateBefore.getOrbit(), stateBefore.getDate(), stateBefore.getFrame()).getRotation()
491                       .getQ1(),
492                 stateBefore.getAttitude().getRotation().getQ1(),
493                 1.0E-16);
494         Assertions.assertEquals(
495                 before.getAttitude(stateBefore.getOrbit(), stateBefore.getDate(), stateBefore.getFrame()).getRotation()
496                       .getQ2(),
497                 stateBefore.getAttitude().getRotation().getQ2(),
498                 1.0E-16);
499         Assertions.assertEquals(
500                 before.getAttitude(stateBefore.getOrbit(), stateBefore.getDate(), stateBefore.getFrame()).getRotation()
501                       .getQ3(),
502                 stateBefore.getAttitude().getRotation().getQ3(),
503                 1.0E-16);
504 
505         Assertions.assertEquals(
506                 after.getAttitude(stateAfter.getOrbit(), stateAfter.getDate(), stateAfter.getFrame()).getRotation().getQ0(),
507                 stateAfter.getAttitude().getRotation().getQ0(),
508                 1.0E-16);
509         Assertions.assertEquals(
510                 after.getAttitude(stateAfter.getOrbit(), stateAfter.getDate(), stateAfter.getFrame()).getRotation().getQ1(),
511                 stateAfter.getAttitude().getRotation().getQ1(),
512                 1.0E-16);
513         Assertions.assertEquals(
514                 after.getAttitude(stateAfter.getOrbit(), stateAfter.getDate(), stateAfter.getFrame()).getRotation().getQ2(),
515                 stateAfter.getAttitude().getRotation().getQ2(),
516                 1.0E-16);
517         Assertions.assertEquals(
518                 after.getAttitude(stateAfter.getOrbit(), stateAfter.getDate(), stateAfter.getFrame()).getRotation().getQ3(),
519                 stateAfter.getAttitude().getRotation().getQ3(),
520                 1.0E-16);
521     }
522 
523     @Test
524     public void testNonResettableState() {
525         setUp();
526         try {
527             propagator.setAttitudeProvider(new LofOffset(inertialFrame, LOFType.LVLH_CCSDS));
528 
529             List<SpacecraftState> states = new ArrayList<>();
530             for (double dt = 0; dt >= -1200; dt -= 60.0) {
531                 states.add(propagator.propagate(initDate.shiftedBy(dt)));
532             }
533 
534             // Create interpolator
535             final TimeInterpolator<SpacecraftState> interpolator = new SpacecraftStateInterpolator(inertialFrame);
536 
537             new Ephemeris(states, interpolator).resetInitialState(propagator.getInitialState());
538             Assertions.fail("an exception should have been thrown");
539         }
540         catch (OrekitException oe) {
541             Assertions.assertEquals(OrekitMessages.NON_RESETABLE_STATE, oe.getSpecifier());
542         }
543     }
544 
545     @Test
546     public void testAdditionalStates() {
547         setUp();
548 
549         final String name1 = "dt0";
550         final String name2 = "dt1";
551         propagator.setAttitudeProvider(new LofOffset(inertialFrame, LOFType.LVLH_CCSDS));
552 
553         List<SpacecraftState> states = new ArrayList<>();
554         for (double dt = 0; dt >= -1200; dt -= 60.0) {
555             final SpacecraftState original = propagator.propagate(initDate.shiftedBy(dt));
556             final SpacecraftState expanded = original.addAdditionalData(name2, original.getDate().durationFrom(finalDate));
557             states.add(expanded);
558         }
559 
560         // Create interpolator
561         final TimeInterpolator<SpacecraftState> interpolator = new SpacecraftStateInterpolator(inertialFrame);
562 
563         final Propagator ephem = new Ephemeris(states, interpolator);
564         ephem.addAdditionalDataProvider(new AdditionalDataProvider<double[]>() {
565             public String getName() {
566                 return name1;
567             }
568 
569             public double[] getAdditionalData(SpacecraftState state) {
570                 return new double[] { state.getDate().durationFrom(initDate) };
571             }
572         });
573 
574         final String[] additional = ephem.getManagedAdditionalData();
575         Arrays.sort(additional);
576         Assertions.assertEquals(2, additional.length);
577         Assertions.assertEquals(name1, ephem.getManagedAdditionalData()[0]);
578         Assertions.assertEquals(name2, ephem.getManagedAdditionalData()[1]);
579         Assertions.assertTrue(ephem.isAdditionalDataManaged(name1));
580         Assertions.assertTrue(ephem.isAdditionalDataManaged(name2));
581         Assertions.assertFalse(ephem.isAdditionalDataManaged("not managed"));
582 
583         SpacecraftState s = ephem.propagate(initDate.shiftedBy(-270.0));
584         Assertions.assertEquals(-270.0, s.getAdditionalState(name1)[0], 1.0e-15);
585         Assertions.assertEquals(-86670.0, s.getAdditionalState(name2)[0], 1.0e-15);
586 
587     }
588 
589     @Test
590     public void testProtectedMethods()
591             throws SecurityException, NoSuchMethodException,
592             InvocationTargetException, IllegalAccessException {
593         setUp();
594 
595         propagator.setAttitudeProvider(new LofOffset(inertialFrame, LOFType.LVLH_CCSDS));
596 
597         List<SpacecraftState> states = new ArrayList<>();
598         for (double dt = 0; dt >= -1200; dt -= 60.0) {
599             final SpacecraftState original = propagator.propagate(initDate.shiftedBy(dt));
600             final SpacecraftState modified = new SpacecraftState(original.getOrbit(),
601                                                                  original.getAttitude(),
602                                                                  original.getMass() - 0.0625 * dt);
603             states.add(modified);
604         }
605 
606         // Create interpolator
607         final TimeInterpolator<SpacecraftState> interpolator = new SpacecraftStateInterpolator(inertialFrame);
608 
609         final Propagator ephem          = new Ephemeris(states, interpolator);
610         Method           propagateOrbit = Ephemeris.class.getDeclaredMethod("propagateOrbit", AbsoluteDate.class);
611         propagateOrbit.setAccessible(true);
612         Method getMass = Ephemeris.class.getDeclaredMethod("getMass", AbsoluteDate.class);
613         getMass.setAccessible(true);
614 
615         SpacecraftState s = ephem.propagate(initDate.shiftedBy(-270.0));
616         Orbit           o = (Orbit) propagateOrbit.invoke(ephem, s.getDate());
617         double          m = (Double) getMass.invoke(ephem, s.getDate());
618         Assertions.assertEquals(0.0,
619                                 Vector3D.distance(s.getPosition(),
620                                                   o.getPosition()),
621                                 1.0e-15);
622         Assertions.assertEquals(s.getMass(), m, 1.0e-15);
623 
624     }
625 
626     @Test
627     public void testGetCovariance() {
628 
629         // Given
630         setUp();
631 
632         double                dt          = finalDate.durationFrom(initDate);
633         double                timeStep    = dt / 20.0;
634         List<SpacecraftState> states      = new ArrayList<>();
635         List<StateCovariance> covariances = new ArrayList<>();
636 
637         final StateCovariance initialCovariance = getDefaultCovariance(initDate, inertialFrame);
638 
639         for (double t = 0; t <= dt; t += timeStep) {
640             final AbsoluteDate    currentDate       = initDate.shiftedBy(t);
641             final SpacecraftState currentState      = propagator.propagate(currentDate);
642             final StateCovariance currentCovariance = initialCovariance.shiftedBy(initialState, t);
643             states.add(currentState);
644             covariances.add(currentCovariance);
645         }
646 
647         // Create state interpolator
648         final int nbInterpolationPoints = 5;
649         final TimeInterpolator<SpacecraftState> stateInterpolator =
650                 new SpacecraftStateInterpolator(nbInterpolationPoints, inertialFrame, inertialFrame);
651 
652         // Create covariance interpolators
653         final SmoothStepFactory.SmoothStepFunction blendingFunction = SmoothStepFactory.getQuadratic();
654 
655         final OrbitBlender orbitBlender =
656                 new OrbitBlender(blendingFunction, new KeplerianPropagator(initialState), inertialFrame);
657 
658         final TimeInterpolator<TimeStampedPair<Orbit, StateCovariance>> covarianceBlender =
659                 new StateCovarianceBlender(blendingFunction, orbitBlender, inertialFrame, OrbitType.CARTESIAN,
660                                            PositionAngleType.MEAN);
661 
662         final TimeInterpolator<TimeStampedPair<Orbit, StateCovariance>> covarianceHermite =
663                 new StateCovarianceKeplerianHermiteInterpolator(2, orbitBlender, CartesianDerivativesFilter.USE_PVA,
664                                                                 inertialFrame, OrbitType.CARTESIAN, PositionAngleType.MEAN);
665 
666         Ephemeris ephemerisUsingBlender = new Ephemeris(states, stateInterpolator, covariances, covarianceBlender);
667         Ephemeris ephemerisUsingHermite = new Ephemeris(states, stateInterpolator, covariances, covarianceHermite);
668 
669         // When
670         final AbsoluteDate initialDateShiftedByHalfTimeStep = initDate.shiftedBy(timeStep * 0.5);
671         final AbsoluteDate initialDateShiftedByTimeStep     = initDate.shiftedBy(timeStep);
672 
673         final StateCovariance blendedCovarianceAtT0 = ephemerisUsingBlender.getCovariance(initDate).get();
674         final StateCovariance blendedCovarianceBetweenT0AndT1 =
675                 ephemerisUsingBlender.getCovariance(initialDateShiftedByHalfTimeStep).get();
676         final StateCovariance blendedCovarianceAtT1 =
677                 ephemerisUsingBlender.getCovariance(initialDateShiftedByTimeStep).get();
678 
679         final StateCovariance hermiteCovarianceAtT0 = ephemerisUsingHermite.getCovariance(initDate).get();
680         final StateCovariance hermiteCovarianceBetweenT0AndT1 =
681                 ephemerisUsingHermite.getCovariance(initialDateShiftedByHalfTimeStep).get();
682         final StateCovariance hermiteCovarianceAtT1 =
683                 ephemerisUsingHermite.getCovariance(initialDateShiftedByTimeStep).get();
684 
685         // Then
686         final RealMatrix initialMatrix = initialCovariance.getMatrix();
687         final RealMatrix initialMatrixShiftedByHalfTimeStep =
688                 initialCovariance.shiftedBy(initialState, timeStep * 0.5).getMatrix();
689         final RealMatrix initialMatrixShiftedByTimeStep =
690                 initialCovariance.shiftedBy(initialState, timeStep).getMatrix();
691         final RealMatrix zeroMatrix =
692                 MatrixUtils.createRealMatrix(StateCovariance.STATE_DIMENSION, StateCovariance.STATE_DIMENSION);
693 
694         StateCovarianceTest.compareCovariance(zeroMatrix,
695                                               blendedCovarianceAtT0.getMatrix().subtract(initialMatrix),
696                                               1e-11);
697         StateCovarianceTest.compareCovariance(zeroMatrix,
698                                               blendedCovarianceBetweenT0AndT1.getMatrix()
699                                                                              .subtract(initialMatrixShiftedByHalfTimeStep),
700                                               1e-9);
701         StateCovarianceTest.compareCovariance(zeroMatrix, blendedCovarianceAtT1.getMatrix()
702                                                                                .subtract(initialMatrixShiftedByTimeStep),
703                                               1e-9);
704 
705         StateCovarianceTest.compareCovariance(zeroMatrix,
706                                               hermiteCovarianceAtT0.getMatrix().subtract(initialMatrix),
707                                               1e-11);
708         StateCovarianceTest.compareCovariance(zeroMatrix,
709                                               hermiteCovarianceBetweenT0AndT1.getMatrix()
710                                                                              .subtract(initialMatrixShiftedByHalfTimeStep),
711                                               1e-9);
712         StateCovarianceTest.compareCovariance(zeroMatrix, hermiteCovarianceAtT1.getMatrix()
713                                                                                .subtract(initialMatrixShiftedByTimeStep),
714                                               1e-8);
715 
716     }
717 
718     private StateCovariance getDefaultCovariance(final AbsoluteDate date, final Frame frame) {
719 
720         final RealMatrix covarianceMatrix = new BlockRealMatrix(new double[][] { { 1, 0, 0, 0, 0, 0 },
721                                                                                  { 0, 1, 0, 0, 0, 0 },
722                                                                                  { 0, 0, 1, 0, 0, 0 },
723                                                                                  { 0, 0, 0, 1e-3, 0, 0 },
724                                                                                  { 0, 0, 0, 0, 1e-3, 0 },
725                                                                                  { 0, 0, 0, 0, 0, 1e-3 }, });
726 
727         return new StateCovariance(covarianceMatrix, date, frame, OrbitType.CARTESIAN, PositionAngleType.MEAN);
728     }
729 
730     @Test
731     public void testExtrapolation() {
732         setUp();
733 
734         double                dt       = finalDate.durationFrom(initDate);
735         double                timeStep = dt / 20.0;
736         List<SpacecraftState> states   = new ArrayList<>();
737 
738         for (double t = 0; t <= dt; t += timeStep) {
739             states.add(propagator.propagate(initDate.shiftedBy(t)));
740         }
741 
742         // Create interpolator
743         final int nbInterpolationPoints = 5;
744         final TimeInterpolator<SpacecraftState> interpolator =
745                 new SpacecraftStateInterpolator(nbInterpolationPoints, inertialFrame, inertialFrame);
746 
747         Ephemeris ephemeris = new Ephemeris(states, interpolator);
748         Assertions.assertEquals(finalDate, ephemeris.getMaxDate());
749 
750         double tolerance = interpolator.getExtrapolationThreshold();
751 
752         final TimeStampedPVCoordinates interpolatedPV =
753                 ephemeris.getPVCoordinates(initDate.shiftedBy(timeStep), inertialFrame);
754 
755         final TimeStampedPVCoordinates propagatedPV = states.get(1).getPVCoordinates();
756 
757         // Assert getPVCoordinates method
758         AbsolutePVCoordinatesTest.assertPV(propagatedPV, interpolatedPV, 1e-8);
759 
760         // Assert getFrame method
761         Assertions.assertEquals(inertialFrame, ephemeris.getFrame());
762 
763         ephemeris.propagate(ephemeris.getMinDate());
764         ephemeris.propagate(ephemeris.getMaxDate());
765         ephemeris.propagate(ephemeris.getMinDate().shiftedBy(-tolerance / 2.0));
766         ephemeris.propagate(ephemeris.getMaxDate().shiftedBy(tolerance / 2.0));
767 
768         try {
769             ephemeris.propagate(ephemeris.getMinDate().shiftedBy(-2.0 * tolerance));
770             Assertions.fail("an exception should have been thrown");
771         }
772         catch (TimeStampedCacheException e) {
773             //supposed to fail since out of bounds
774         }
775 
776         try {
777             ephemeris.propagate(ephemeris.getMaxDate().shiftedBy(2.0 * tolerance));
778             Assertions.fail("an exception should have been thrown");
779         }
780         catch (TimeStampedCacheException e) {
781             //supposed to fail since out of bounds
782         }
783     }
784 
785     @Test
786     public void testIssue662() {
787         setUp();
788 
789         // Input parameters
790         int    numberOfIntervals = 1440;
791         double deltaT            = finalDate.durationFrom(initDate) / ((double) numberOfIntervals);
792 
793         // Build the list of spacecraft states
794         String                additionalName  = "testValue";
795         double                additionalValue = 1.0;
796         List<SpacecraftState> states          = new ArrayList<>(numberOfIntervals + 1);
797         for (int j = 0; j <= numberOfIntervals; j++) {
798             states.add(propagator.propagate(initDate.shiftedBy((j * deltaT)))
799                                  .addAdditionalData(additionalName, additionalValue));
800         }
801 
802         // Build the ephemeris propagator
803         // Create interpolator
804         final TimeInterpolator<SpacecraftState> interpolator = new SpacecraftStateInterpolator(inertialFrame);
805 
806         Ephemeris ephemPropagator = new Ephemeris(states, interpolator);
807 
808         // State before adding an attitude provider
809         SpacecraftState stateBefore = ephemPropagator.propagate(ephemPropagator.getMaxDate().shiftedBy(-60.0));
810         Assertions.assertEquals(1, stateBefore.getAdditionalState(additionalName).length);
811         Assertions.assertEquals(additionalValue, stateBefore.getAdditionalState(additionalName)[0], Double.MIN_VALUE);
812 
813         // Set an attitude provider
814         ephemPropagator.setAttitudeProvider(new LofOffset(inertialFrame, LOFType.LVLH_CCSDS));
815 
816         // State after adding an attitude provider
817         SpacecraftState stateAfter = ephemPropagator.propagate(ephemPropagator.getMaxDate().shiftedBy(-60.0));
818         Assertions.assertEquals(1, stateAfter.getAdditionalState(additionalName).length);
819         Assertions.assertEquals(additionalValue, stateAfter.getAdditionalState(additionalName)[0], Double.MIN_VALUE);
820 
821     }
822 
823     @Test
824     public void testIssue680() {
825         setUp();
826 
827         // Initial PV coordinates
828         AbsolutePVCoordinates initPV = new AbsolutePVCoordinates(inertialFrame,
829                                                                  new TimeStampedPVCoordinates(initDate,
830                                                                                               new PVCoordinates(new Vector3D(
831                                                                                                       -29536113.0,
832                                                                                                       30329259.0, -100125.0),
833                                                                                                                 new Vector3D(
834                                                                                                                         -2194.0,
835                                                                                                                         -2141.0,
836                                                                                                                         -8.0))));
837         // Input parameters
838         int    numberOfIntervals = 1440;
839         double deltaT            = finalDate.durationFrom(initDate) / ((double) numberOfIntervals);
840 
841         // Build the list of spacecraft states
842         List<SpacecraftState> states = new ArrayList<>(numberOfIntervals + 1);
843         for (int j = 0; j <= numberOfIntervals; j++) {
844             states.add(new SpacecraftState(initPV).shiftedBy(j * deltaT));
845         }
846 
847         // Build the ephemeris propagator
848         // Create interpolator
849         final TimeInterpolator<SpacecraftState> interpolator = new SpacecraftStateInterpolator(inertialFrame);
850 
851         Ephemeris ephemPropagator = new Ephemeris(states, interpolator);
852 
853         // Get initial state without attitude provider
854         SpacecraftState withoutAttitudeProvider = ephemPropagator.getInitialState();
855         Assertions.assertEquals(0.0,
856                                 Vector3D.distance(withoutAttitudeProvider.getAbsPVA().getPosition(), initPV.getPosition()),
857                                 1.0e-10);
858         Assertions.assertEquals(0.0,
859                                 Vector3D.distance(withoutAttitudeProvider.getAbsPVA().getVelocity(), initPV.getVelocity()),
860                                 1.0e-10);
861 
862         // Set an attitude provider
863         ephemPropagator.setAttitudeProvider(new LofOffset(inertialFrame, LOFType.LVLH_CCSDS));
864 
865         // Get initial state with attitude provider
866         SpacecraftState withAttitudeProvider = ephemPropagator.getInitialState();
867         Assertions.assertEquals(0.0, Vector3D.distance(withAttitudeProvider.getAbsPVA().getPosition(), initPV.getPosition()),
868                                 1.0e-10);
869         Assertions.assertEquals(0.0, Vector3D.distance(withAttitudeProvider.getAbsPVA().getVelocity(), initPV.getVelocity()),
870                                 1.0e-10);
871 
872     }
873 
874     /** Error with specific propagators & additional state provider throwing a NullPointerException when propagating */
875     @Test
876     public void testIssue949() {
877         // GIVEN
878         final AbsoluteDate initialDate  = new AbsoluteDate();
879         final Orbit        initialOrbit = TestUtils.getDefaultOrbit(initialDate);
880 
881         // Setup propagator
882         final Orbit                 finalOrbit = initialOrbit.shiftedBy(1);
883         final List<SpacecraftState> states     = new ArrayList<>();
884         states.add(new SpacecraftState(initialOrbit));
885         states.add(new SpacecraftState(finalOrbit));
886 
887         final Ephemeris ephemeris = new Ephemeris(states, 2);
888 
889         // Setup additional data provider which use the initial state in its init method
890         final AdditionalDataProvider<double[]> additionalDataProvider = TestUtils.getAdditionalProviderWithInit();
891         ephemeris.addAdditionalDataProvider(additionalDataProvider);
892 
893         // WHEN & THEN
894         Assertions.assertDoesNotThrow(() -> ephemeris.propagate(ephemeris.getMaxDate()), "No error should have been thrown");
895 
896     }
897 
898     @Test
899     @Timeout(value = 5000, unit = TimeUnit.SECONDS)
900     void testIssue1277() {
901         // GIVEN
902         // Create orbit
903         double       mu    = Constants.IERS2010_EARTH_MU;
904         AbsoluteDate date  = new AbsoluteDate();
905         Frame        frame = FramesFactory.getEME2000();
906         Orbit orbit = new KeplerianOrbit(6378000 + 500000, 0.01, FastMath.toRadians(15), 0, 0, 0,
907                                          PositionAngleType.MEAN, frame,date, mu);
908         // Create propagator
909         Propagator propagator = new KeplerianPropagator(orbit);
910 
911         // Add step handler for saving states
912         List<SpacecraftState> states = new ArrayList<>();
913         propagator.getMultiplexer().add(10, states::add);
914 
915         // Run propagation over one week
916         propagator.propagate(date.shiftedBy(Constants.JULIAN_DAY * 7));
917 
918         // Create event detector
919         GeodeticPoint    station = new GeodeticPoint(0, 0, 0);
920         Frame            itrf    = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
921         OneAxisEllipsoid wgs84   = ReferenceEllipsoid.getWgs84(itrf);
922         TopocentricFrame topo    = new TopocentricFrame(wgs84, station, "station");
923 
924         RecordAndContinue handler = new RecordAndContinue();
925         ElevationDetector detector = new ElevationDetector(topo).withConstantElevation(FastMath.toRadians(5))
926                                                                 .withRefraction(new EarthStandardAtmosphereRefraction())
927                                                                 .withMaxCheck(60).withHandler(handler);
928 
929         // Create ephemeris from states
930         Ephemeris ephemeris = new Ephemeris(states, 4);
931 
932         // Add detector
933         ephemeris.addEventDetector(detector);
934 
935         // WHEN & THEN
936         ephemeris.propagate(ephemeris.getMinDate(), ephemeris.getMaxDate());
937     }
938 
939     /* Ephemeris propagate fails for interpolation point value of 1 */
940     @Test
941     void testMinInterpolationPoints() {
942         // GIVEN
943         final AbsoluteDate initialDate  = new AbsoluteDate();
944         final Orbit        initialOrbit = TestUtils.getDefaultOrbit(initialDate);
945 
946         // Setup propagator
947         final Orbit                 finalOrbit = initialOrbit.shiftedBy(1);
948         final List<SpacecraftState> states     = new ArrayList<>();
949         states.add(new SpacecraftState(initialOrbit));
950         states.add(new SpacecraftState(finalOrbit));
951 
952         final Ephemeris ephemeris = new Ephemeris(states, 1);
953 
954         // WHEN
955         final AbsoluteDate t = ephemeris.getMaxDate();
956         final SpacecraftState actual = ephemeris.propagate(t);
957 
958         // THEN
959         MatcherAssert.assertThat(actual.getDate(),
960                 OrekitMatchers.durationFrom(t, OrekitMatchers.closeTo(0, 0)));
961         MatcherAssert.assertThat(actual.getPVCoordinates(),
962                 OrekitMatchers.pvCloseTo(finalOrbit.getPVCoordinates(), 0.0));
963     }
964 
965     public void setUp() throws IllegalArgumentException, OrekitException {
966         Utils.setDataRoot("regular-data");
967 
968         initDate = new AbsoluteDate(new DateComponents(2004, 01, 01),
969                                     TimeComponents.H00,
970                                     TimeScalesFactory.getUTC());
971 
972         finalDate = new AbsoluteDate(new DateComponents(2004, 01, 02),
973                                      TimeComponents.H00,
974                                      TimeScalesFactory.getUTC());
975 
976         double a     = 7187990.1979844316;
977         double e     = 0.5e-4;
978         double i     = 1.7105407051081795;
979         double omega = 1.9674147913622104;
980         double OMEGA = FastMath.toRadians(261);
981         double lv    = 0;
982         double mu    = 3.9860047e14;
983         inertialFrame = FramesFactory.getEME2000();
984 
985         initialState = new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngleType.TRUE,
986                                           inertialFrame, initDate, mu);
987         propagator   = new KeplerianPropagator(initialState);
988 
989     }
990 
991     private double calculatePositionDelta(SpacecraftState state1, SpacecraftState state2) {
992         return Vector3D.distance(state1.getPosition(), state2.getPosition());
993     }
994 
995     private double calculateVelocityDelta(SpacecraftState state1, SpacecraftState state2) {
996         return Vector3D.distance(state1.getPVCoordinates().getVelocity(), state2.getPVCoordinates().getVelocity());
997     }
998 
999     private double calculateAttitudeDelta(SpacecraftState state1, SpacecraftState state2) {
1000         return Rotation.distance(state1.getAttitude().getRotation(), state2.getAttitude().getRotation());
1001     }
1002 
1003 }