1   /* Copyright 2002-2022 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 java.lang.reflect.InvocationTargetException;
20  import java.lang.reflect.Method;
21  import java.util.ArrayList;
22  import java.util.Arrays;
23  import java.util.List;
24  
25  import org.hipparchus.geometry.euclidean.threed.Rotation;
26  import org.hipparchus.geometry.euclidean.threed.Vector3D;
27  import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
28  import org.hipparchus.util.FastMath;
29  import org.junit.Assert;
30  import org.junit.Before;
31  import org.junit.Test;
32  import org.orekit.Utils;
33  import org.orekit.attitudes.AttitudeProvider;
34  import org.orekit.attitudes.AttitudesSequence;
35  import org.orekit.attitudes.CelestialBodyPointed;
36  import org.orekit.attitudes.LofOffset;
37  import org.orekit.bodies.CelestialBodyFactory;
38  import org.orekit.errors.OrekitException;
39  import org.orekit.errors.OrekitMessages;
40  import org.orekit.errors.TimeStampedCacheException;
41  import org.orekit.frames.Frame;
42  import org.orekit.frames.FramesFactory;
43  import org.orekit.frames.LOFType;
44  import org.orekit.orbits.KeplerianOrbit;
45  import org.orekit.orbits.Orbit;
46  import org.orekit.orbits.PositionAngle;
47  import org.orekit.propagation.AdditionalStateProvider;
48  import org.orekit.propagation.Propagator;
49  import org.orekit.propagation.SpacecraftState;
50  import org.orekit.propagation.events.DateDetector;
51  import org.orekit.propagation.events.handlers.ContinueOnEvent;
52  import org.orekit.propagation.numerical.NumericalPropagator;
53  import org.orekit.time.AbsoluteDate;
54  import org.orekit.time.DateComponents;
55  import org.orekit.time.TimeComponents;
56  import org.orekit.time.TimeScalesFactory;
57  import org.orekit.utils.AbsolutePVCoordinates;
58  import org.orekit.utils.AngularDerivativesFilter;
59  import org.orekit.utils.Constants;
60  import org.orekit.utils.PVCoordinates;
61  import org.orekit.utils.TimeStampedPVCoordinates;
62  
63  public class EphemerisTest {
64  
65      private AbsoluteDate initDate;
66      private AbsoluteDate finalDate;
67      private Frame        inertialFrame;
68      private Propagator   propagator;
69  
70      @Test
71      public void testAttitudeOverride() throws IllegalArgumentException, OrekitException {
72          final double positionTolerance = 1e-6;
73          final double velocityTolerance = 1e-5;
74          final double attitudeTolerance = 1e-6;
75  
76          int numberOfInterals = 1440;
77          double deltaT = finalDate.durationFrom(initDate)/((double)numberOfInterals);
78  
79          propagator.setAttitudeProvider(new LofOffset(inertialFrame, LOFType.LVLH_CCSDS));
80  
81          List<SpacecraftState> states = new ArrayList<SpacecraftState>(numberOfInterals + 1);
82          for (int j = 0; j<= numberOfInterals; j++) {
83              states.add(propagator.propagate(initDate.shiftedBy((j * deltaT))));
84          }
85  
86          int numInterpolationPoints = 2;
87          Ephemeris ephemPropagator = new Ephemeris(states, numInterpolationPoints);
88          Assert.assertEquals(0, ephemPropagator.getManagedAdditionalStates().length);
89  
90          //First test that we got position, velocity and attitude nailed
91          int numberEphemTestIntervals = 2880;
92          deltaT = finalDate.durationFrom(initDate)/((double)numberEphemTestIntervals);
93          for (int j = 0; j <= numberEphemTestIntervals; j++) {
94              AbsoluteDate currentDate = initDate.shiftedBy(j * deltaT);
95              SpacecraftState ephemState = ephemPropagator.propagate(currentDate);
96              SpacecraftState keplerState = propagator.propagate(currentDate);
97              double positionDelta = calculatePositionDelta(ephemState, keplerState);
98              double velocityDelta = calculateVelocityDelta(ephemState, keplerState);
99              double attitudeDelta = calculateAttitudeDelta(ephemState, keplerState);
100             Assert.assertEquals("LVLH_CCSDS Unmatched Position at: " + currentDate, 0.0, positionDelta, positionTolerance);
101             Assert.assertEquals("LVLH_CCSDS Unmatched Velocity at: " + currentDate, 0.0, velocityDelta, velocityTolerance);
102             Assert.assertEquals("LVLH_CCSDS Unmatched Attitude at: " + currentDate, 0.0, attitudeDelta, attitudeTolerance);
103         }
104 
105         //Now force an override on the attitude and check it against a Keplerian propagator
106         //setup identically to the first but with a different attitude
107         //If override isn't working this will fail.
108         propagator = new KeplerianPropagator(propagator.getInitialState().getOrbit());
109         propagator.setAttitudeProvider(new LofOffset(inertialFrame, LOFType.QSW));
110 
111         ephemPropagator.setAttitudeProvider(new LofOffset(inertialFrame, LOFType.QSW));
112         for (int j = 0; j <= numberEphemTestIntervals; j++) {
113             AbsoluteDate currentDate = initDate.shiftedBy(j * deltaT);
114             SpacecraftState ephemState = ephemPropagator.propagate(currentDate);
115             SpacecraftState keplerState = propagator.propagate(currentDate);
116             double positionDelta = calculatePositionDelta(ephemState, keplerState);
117             double velocityDelta = calculateVelocityDelta(ephemState, keplerState);
118             double attitudeDelta = calculateAttitudeDelta(ephemState, keplerState);
119             Assert.assertEquals("QSW Unmatched Position at: " + currentDate, 0.0, positionDelta, positionTolerance);
120             Assert.assertEquals("QSW Unmatched Velocity at: " + currentDate, 0.0, velocityDelta, velocityTolerance);
121             Assert.assertEquals("QSW Unmatched Attitude at: " + currentDate, 0.0, attitudeDelta, attitudeTolerance);
122         }
123 
124     }
125     
126     @Test
127     public void testAttitudeSequenceTransition() {
128     	        
129         // Initialize the orbit
130     	final AbsoluteDate initialDate = new AbsoluteDate(2003, 01, 01, 0, 0, 00.000, TimeScalesFactory.getUTC());
131         final Vector3D position  = new Vector3D(-39098981.4866597, -15784239.3610601, 78908.2289853595);
132         final Vector3D velocity  = new Vector3D(1151.00321021175, -2851.14864755189, -2.02133248357321);
133         final Orbit initialOrbit = new KeplerianOrbit(new PVCoordinates(position, velocity),
134                                                       FramesFactory.getGCRF(), initialDate,
135                                                       Constants.WGS84_EARTH_MU);
136         final SpacecraftState initialState = new SpacecraftState(initialOrbit);        
137 
138         // Define attitude laws
139         AttitudeProvider before = new CelestialBodyPointed(FramesFactory.getICRF(), CelestialBodyFactory.getSun(), Vector3D.PLUS_K, Vector3D.PLUS_I, Vector3D.PLUS_K);
140         AttitudeProvider after = new CelestialBodyPointed(FramesFactory.getICRF(), CelestialBodyFactory.getEarth(), Vector3D.PLUS_K, Vector3D.PLUS_I, Vector3D.PLUS_K);
141 
142         // Define attitude sequence
143         AbsoluteDate switchDate = initialDate.shiftedBy(86400.0);
144         double transitionTime = 600;
145         DateDetector switchDetector = new DateDetector(switchDate).withHandler(new ContinueOnEvent<DateDetector>());
146 
147         AttitudesSequence attitudeSequence = new AttitudesSequence();
148         attitudeSequence.resetActiveProvider(before);
149         attitudeSequence.addSwitchingCondition(before, after, switchDetector, true, false, transitionTime, AngularDerivativesFilter.USE_RR, null);
150 
151         NumericalPropagator propagator = new NumericalPropagator(new DormandPrince853Integrator(0.1, 500, 1e-9, 1e-9));
152         propagator.setInitialState(initialState);
153         
154         // Propagate and build ephemeris
155         final List<SpacecraftState> propagatedStates = new ArrayList<>();
156 
157         propagator.setStepHandler(60, currentState -> propagatedStates.add(currentState));
158         propagator.propagate(initialDate.shiftedBy(2*86400.0));
159         final Ephemeris ephemeris = new Ephemeris(propagatedStates, 8);
160 
161         // Add attitude switch event to ephemeris
162         ephemeris.setAttitudeProvider(attitudeSequence);
163         attitudeSequence.registerSwitchEvents(ephemeris);
164 
165         // Propagate with a step during the transition
166         AbsoluteDate endDate = initialDate.shiftedBy(2*86400.0);
167         SpacecraftState stateBefore = ephemeris.getInitialState();
168         ephemeris.propagate(switchDate.shiftedBy(transitionTime/2));
169         SpacecraftState stateAfter = ephemeris.propagate(endDate);
170         
171         
172         // Check that the attitudes are correct
173         Assert.assertEquals(before.getAttitude(stateBefore.getOrbit(), stateBefore.getDate(), stateBefore.getFrame()).getRotation().getQ0(),
174         		stateBefore.getAttitude().getRotation().getQ0(),
175         		1.0E-16);
176         Assert.assertEquals(before.getAttitude(stateBefore.getOrbit(), stateBefore.getDate(), stateBefore.getFrame()).getRotation().getQ1(),
177         		stateBefore.getAttitude().getRotation().getQ1(),
178         		1.0E-16);
179         Assert.assertEquals(before.getAttitude(stateBefore.getOrbit(), stateBefore.getDate(), stateBefore.getFrame()).getRotation().getQ2(),
180         		stateBefore.getAttitude().getRotation().getQ2(),
181         		1.0E-16);
182         Assert.assertEquals(before.getAttitude(stateBefore.getOrbit(), stateBefore.getDate(), stateBefore.getFrame()).getRotation().getQ3(),
183         		stateBefore.getAttitude().getRotation().getQ3(),
184         		1.0E-16);
185 
186         Assert.assertEquals(after.getAttitude(stateAfter.getOrbit(), stateAfter.getDate(), stateAfter.getFrame()).getRotation().getQ0(),
187         		stateAfter.getAttitude().getRotation().getQ0(),
188         		1.0E-16);
189         Assert.assertEquals(after.getAttitude(stateAfter.getOrbit(), stateAfter.getDate(), stateAfter.getFrame()).getRotation().getQ1(),
190         		stateAfter.getAttitude().getRotation().getQ1(),
191         		1.0E-16);
192         Assert.assertEquals(after.getAttitude(stateAfter.getOrbit(), stateAfter.getDate(), stateAfter.getFrame()).getRotation().getQ2(),
193         		stateAfter.getAttitude().getRotation().getQ2(),
194         		1.0E-16);
195         Assert.assertEquals(after.getAttitude(stateAfter.getOrbit(), stateAfter.getDate(), stateAfter.getFrame()).getRotation().getQ3(),
196         		stateAfter.getAttitude().getRotation().getQ3(),
197         		1.0E-16);
198     }
199 
200     @Test
201     public void testNonResettableState() {
202         try {
203             propagator.setAttitudeProvider(new LofOffset(inertialFrame, LOFType.LVLH_CCSDS));
204 
205             List<SpacecraftState> states = new ArrayList<SpacecraftState>();
206             for (double dt = 0; dt >= -1200; dt -= 60.0) {
207                 states.add(propagator.propagate(initDate.shiftedBy(dt)));
208             }
209 
210             new Ephemeris(states, 2).resetInitialState(propagator.getInitialState());
211             Assert.fail("an exception should have been thrown");
212         } catch (OrekitException oe) {
213             Assert.assertEquals(OrekitMessages.NON_RESETABLE_STATE, oe.getSpecifier());
214         }
215     }
216 
217     @Test
218     public void testAdditionalStates() {
219         final String name1  = "dt0";
220         final String name2  = "dt1";
221         propagator.setAttitudeProvider(new LofOffset(inertialFrame, LOFType.LVLH_CCSDS));
222 
223         List<SpacecraftState> states = new ArrayList<SpacecraftState>();
224         for (double dt = 0; dt >= -1200; dt -= 60.0) {
225             final SpacecraftState original = propagator.propagate(initDate.shiftedBy(dt));
226             final SpacecraftState expanded = original.addAdditionalState(name2, original.getDate().durationFrom(finalDate));
227             states.add(expanded);
228         }
229 
230         final Propagator ephem = new Ephemeris(states, 2);
231         ephem.addAdditionalStateProvider(new AdditionalStateProvider() {
232             public String getName() {
233                 return name1;
234             }
235             public double[] getAdditionalState(SpacecraftState state) {
236                 return new double[] { state.getDate().durationFrom(initDate) };
237             }
238         });
239 
240         final String[] additional = ephem.getManagedAdditionalStates();
241         Arrays.sort(additional);
242         Assert.assertEquals(2, additional.length);
243         Assert.assertEquals(name1, ephem.getManagedAdditionalStates()[0]);
244         Assert.assertEquals(name2, ephem.getManagedAdditionalStates()[1]);
245         Assert.assertTrue(ephem.isAdditionalStateManaged(name1));
246         Assert.assertTrue(ephem.isAdditionalStateManaged(name2));
247         Assert.assertFalse(ephem.isAdditionalStateManaged("not managed"));
248 
249         SpacecraftState s = ephem.propagate(initDate.shiftedBy(-270.0));
250         Assert.assertEquals(-270.0,   s.getAdditionalState(name1)[0], 1.0e-15);
251         Assert.assertEquals(-86670.0, s.getAdditionalState(name2)[0], 1.0e-15);
252 
253     }
254 
255     @Test
256     public void testProtectedMethods()
257         throws SecurityException, NoSuchMethodException,
258                InvocationTargetException, IllegalAccessException {
259         propagator.setAttitudeProvider(new LofOffset(inertialFrame, LOFType.LVLH_CCSDS));
260 
261         List<SpacecraftState> states = new ArrayList<SpacecraftState>();
262         for (double dt = 0; dt >= -1200; dt -= 60.0) {
263             final SpacecraftState original = propagator.propagate(initDate.shiftedBy(dt));
264             final SpacecraftState modified = new SpacecraftState(original.getOrbit(),
265                                                                  original.getAttitude(),
266                                                                  original.getMass() - 0.0625 * dt);
267             states.add(modified);
268         }
269 
270         final Propagator ephem = new Ephemeris(states, 2);
271         Method propagateOrbit = Ephemeris.class.getDeclaredMethod("propagateOrbit", AbsoluteDate.class);
272         propagateOrbit.setAccessible(true);
273         Method getMass        = Ephemeris.class.getDeclaredMethod("getMass", AbsoluteDate.class);
274         getMass.setAccessible(true);
275 
276         SpacecraftState s = ephem.propagate(initDate.shiftedBy(-270.0));
277         Orbit  o = (Orbit) propagateOrbit.invoke(ephem, s.getDate());
278         double m = ((Double) getMass.invoke(ephem, s.getDate())).doubleValue();
279         Assert.assertEquals(0.0,
280                             Vector3D.distance(s.getPVCoordinates().getPosition(),
281                                               o.getPVCoordinates().getPosition()),
282                             1.0e-15);
283         Assert.assertEquals(s.getMass(), m, 1.0e-15);
284 
285     }
286 
287     @Test
288     public void testExtrapolation() {
289         double dt = finalDate.durationFrom(initDate);
290         double timeStep = dt / 20.0;
291         List<SpacecraftState> states = new ArrayList<SpacecraftState>();
292 
293         for(double t = 0 ; t <= dt; t+=timeStep) {
294             states.add(propagator.propagate(initDate.shiftedBy(t)));
295         }
296 
297         final int interpolationPoints = 5;
298         Ephemeris ephemeris = new Ephemeris(states, interpolationPoints);
299         Assert.assertEquals(finalDate, ephemeris.getMaxDate());
300 
301         double tolerance = ephemeris.getExtrapolationThreshold();
302 
303         ephemeris.propagate(ephemeris.getMinDate());
304         ephemeris.propagate(ephemeris.getMaxDate());
305         ephemeris.propagate(ephemeris.getMinDate().shiftedBy(-tolerance / 2.0));
306         ephemeris.propagate(ephemeris.getMaxDate().shiftedBy(tolerance / 2.0));
307 
308         try {
309             ephemeris.propagate(ephemeris.getMinDate().shiftedBy(-2.0 * tolerance));
310             Assert.fail("an exception should have been thrown");
311         } catch (TimeStampedCacheException e) {
312             //supposed to fail since out of bounds
313         }
314 
315         try {
316             ephemeris.propagate(ephemeris.getMaxDate().shiftedBy(2.0 * tolerance));
317             Assert.fail("an exception should have been thrown");
318         } catch (TimeStampedCacheException e) {
319             //supposed to fail since out of bounds
320         }
321     }
322 
323     @Test
324     public void testIssue662() {
325 
326         // Input parameters
327         int numberOfInterals = 1440;
328         double deltaT = finalDate.durationFrom(initDate)/((double)numberOfInterals);
329 
330         // Build the list of spacecraft states
331         String additionalName  = "testValue";
332         double additionalValue = 1.0;
333         List<SpacecraftState> states = new ArrayList<SpacecraftState>(numberOfInterals + 1);
334         for (int j = 0; j<= numberOfInterals; j++) {
335             states.add(propagator.propagate(initDate.shiftedBy((j * deltaT))).addAdditionalState(additionalName, additionalValue));
336         }
337 
338         // Build the epemeris propagator
339         Ephemeris ephemPropagator = new Ephemeris(states, 2);
340 
341         // State before adding an attitude provider
342         SpacecraftState stateBefore = ephemPropagator.propagate(ephemPropagator.getMaxDate().shiftedBy(-60.0));
343         Assert.assertEquals(1,               stateBefore.getAdditionalState(additionalName).length);
344         Assert.assertEquals(additionalValue, stateBefore.getAdditionalState(additionalName)[0], Double.MIN_VALUE);
345 
346         // Set an attitude provider
347         ephemPropagator.setAttitudeProvider(new LofOffset(inertialFrame, LOFType.LVLH_CCSDS));
348 
349         // State after adding an attitude provider
350         SpacecraftState stateAfter = ephemPropagator.propagate(ephemPropagator.getMaxDate().shiftedBy(-60.0));
351         Assert.assertEquals(1,               stateAfter.getAdditionalState(additionalName).length);
352         Assert.assertEquals(additionalValue, stateAfter.getAdditionalState(additionalName)[0], Double.MIN_VALUE);
353 
354     }
355 
356     @Test
357     public void testIssue680() {
358 
359         // Initial PV coordinates
360         AbsolutePVCoordinates initPV = new AbsolutePVCoordinates(inertialFrame,
361                                                                  new TimeStampedPVCoordinates(initDate,
362                                                                                               new PVCoordinates(new Vector3D(-29536113.0, 30329259.0, -100125.0),
363                                                                                                                 new Vector3D(-2194.0, -2141.0, -8.0))));
364         // Input parameters
365         int numberOfInterals = 1440;
366         double deltaT = finalDate.durationFrom(initDate)/((double)numberOfInterals);
367 
368         // Build the list of spacecraft states
369         List<SpacecraftState> states = new ArrayList<SpacecraftState>(numberOfInterals + 1);
370         for (int j = 0; j<= numberOfInterals; j++) {
371             states.add(new SpacecraftState(initPV).shiftedBy(j * deltaT));
372         }
373 
374         // Build the epemeris propagator
375         Ephemeris ephemPropagator = new Ephemeris(states, 2);
376 
377         // Get initial state without attitude provider
378         SpacecraftState withoutAttitudeProvider = ephemPropagator.getInitialState();
379         Assert.assertEquals(0.0, Vector3D.distance(withoutAttitudeProvider.getAbsPVA().getPosition(), initPV.getPosition()), 1.0e-10);
380         Assert.assertEquals(0.0, Vector3D.distance(withoutAttitudeProvider.getAbsPVA().getVelocity(), initPV.getVelocity()), 1.0e-10);
381 
382         // Set an attitude provider
383         ephemPropagator.setAttitudeProvider(new LofOffset(inertialFrame, LOFType.LVLH_CCSDS));
384 
385         // Get initial state with attitude provider
386         SpacecraftState withAttitudeProvider = ephemPropagator.getInitialState();
387         Assert.assertEquals(0.0, Vector3D.distance(withAttitudeProvider.getAbsPVA().getPosition(), initPV.getPosition()), 1.0e-10);
388         Assert.assertEquals(0.0, Vector3D.distance(withAttitudeProvider.getAbsPVA().getVelocity(), initPV.getVelocity()), 1.0e-10);
389 
390     }
391 
392     @Before
393     public void setUp() throws IllegalArgumentException, OrekitException {
394         Utils.setDataRoot("regular-data");
395 
396         initDate = new AbsoluteDate(new DateComponents(2004, 01, 01),
397                 TimeComponents.H00,
398                 TimeScalesFactory.getUTC());
399 
400         finalDate = new AbsoluteDate(new DateComponents(2004, 01, 02),
401                  TimeComponents.H00,
402                  TimeScalesFactory.getUTC());
403 
404         double a = 7187990.1979844316;
405         double e = 0.5e-4;
406         double i = 1.7105407051081795;
407         double omega = 1.9674147913622104;
408         double OMEGA = FastMath.toRadians(261);
409         double lv = 0;
410         double mu  = 3.9860047e14;
411         inertialFrame = FramesFactory.getEME2000();
412 
413         Orbit initialState = new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngle.TRUE,
414                                             inertialFrame, initDate, mu);
415         propagator = new KeplerianPropagator(initialState);
416 
417     }
418 
419     private double calculatePositionDelta(SpacecraftState state1, SpacecraftState state2) {
420         return Vector3D.distance(state1.getPVCoordinates().getPosition(), state2.getPVCoordinates().getPosition());
421     }
422 
423     private double calculateVelocityDelta(SpacecraftState state1, SpacecraftState state2) {
424         return Vector3D.distance(state1.getPVCoordinates().getVelocity(), state2.getPVCoordinates().getVelocity());
425     }
426 
427     private double calculateAttitudeDelta(SpacecraftState state1, SpacecraftState state2) {
428         return Rotation.distance(state1.getAttitude().getRotation(), state2.getAttitude().getRotation());
429     }
430 
431 }