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.attitudes;
18  
19  import org.hipparchus.geometry.euclidean.threed.Rotation;
20  import org.hipparchus.geometry.euclidean.threed.Vector3D;
21  import org.hipparchus.util.FastMath;
22  import org.junit.jupiter.api.Assertions;
23  import org.junit.jupiter.api.DisplayName;
24  import org.junit.jupiter.api.Test;
25  import org.orekit.Utils;
26  import org.orekit.annotation.DefaultDataContext;
27  import org.orekit.bodies.OneAxisEllipsoid;
28  import org.orekit.errors.OrekitIllegalArgumentException;
29  import org.orekit.errors.OrekitMessages;
30  import org.orekit.frames.Frame;
31  import org.orekit.frames.FramesFactory;
32  import org.orekit.orbits.CircularOrbit;
33  import org.orekit.propagation.analytical.EcksteinHechlerPropagator;
34  import org.orekit.time.AbsoluteDate;
35  import org.orekit.time.TimeInterpolator;
36  import org.orekit.utils.AngularDerivativesFilter;
37  import org.orekit.utils.Constants;
38  import org.orekit.utils.IERSConventions;
39  import org.orekit.utils.PVCoordinates;
40  import org.orekit.utils.TimeStampedAngularCoordinates;
41  import org.orekit.utils.TimeStampedAngularCoordinatesHermiteInterpolator;
42  
43  import java.util.ArrayList;
44  import java.util.List;
45  
46  import static org.mockito.Mockito.mock;
47  
48  class AttitudeInterpolatorTest {
49  
50      @Test
51      @DefaultDataContext
52      void testInterpolation() {
53  
54          // Given
55          Utils.setDataRoot("regular-data");
56          final double ehMu = 3.9860047e14;
57          final double ae   = 6.378137e6;
58          final double c20  = -1.08263e-3;
59          final double c30  = 2.54e-6;
60          final double c40  = 1.62e-6;
61          final double c50  = 2.3e-7;
62          final double c60  = -5.5e-7;
63  
64          final AbsoluteDate date     = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
65          final Vector3D     position = new Vector3D(3220103., 69623., 6449822.);
66          final Vector3D     velocity = new Vector3D(6414.7, -2006., -3180.);
67          final CircularOrbit initialOrbit = new CircularOrbit(new PVCoordinates(position, velocity),
68                                                               FramesFactory.getEME2000(), date, ehMu);
69  
70          EcksteinHechlerPropagator propagator =
71                  new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
72          OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
73                                                        Constants.WGS84_EARTH_FLATTENING,
74                                                        FramesFactory.getITRF(IERSConventions.IERS_2010, true));
75          propagator.setAttitudeProvider(new BodyCenterPointing(initialOrbit.getFrame(), earth));
76          final Attitude initialAttitude = propagator.propagate(initialOrbit.getDate()).getAttitude();
77  
78          // set up a 5 points sample
79          List<Attitude> sample = new ArrayList<>();
80          for (double dt = 0; dt < 251.0; dt += 60.0) {
81              sample.add(propagator.propagate(date.shiftedBy(dt)).getAttitude());
82          }
83  
84          // Create interpolator
85          final double extrapolationThreshold = 59;
86          final TimeInterpolator<TimeStampedAngularCoordinates> angularInterpolator =
87                  new TimeStampedAngularCoordinatesHermiteInterpolator(sample.size(), extrapolationThreshold,
88                                                                       AngularDerivativesFilter.USE_RR);
89  
90          final AttitudeInterpolator attitudeInterpolator =
91                  new AttitudeInterpolator(initialOrbit.getFrame(), angularInterpolator);
92  
93          // well inside the sample, interpolation should be better than quadratic shift
94          double maxShiftAngleError         = 0;
95          double maxInterpolationAngleError = 0;
96          double maxShiftRateError          = 0;
97          double maxInterpolationRateError  = 0;
98          for (double dt = 0; dt < 240.0; dt += 1.0) {
99              AbsoluteDate t          = initialOrbit.getDate().shiftedBy(dt);
100             Attitude     propagated = propagator.propagate(t).getAttitude();
101             double shiftAngleError = Rotation.distance(propagated.getRotation(),
102                                                        initialAttitude.shiftedBy(dt).getRotation());
103             double interpolationAngleError = Rotation.distance(propagated.getRotation(),
104                                                                attitudeInterpolator.interpolate(t, sample).getRotation());
105             double shiftRateError = Vector3D.distance(propagated.getSpin(),
106                                                       initialAttitude.shiftedBy(dt).getSpin());
107             double interpolationRateError = Vector3D.distance(propagated.getSpin(),
108                                                               attitudeInterpolator.interpolate(t, sample).getSpin());
109             maxShiftAngleError         = FastMath.max(maxShiftAngleError, shiftAngleError);
110             maxInterpolationAngleError = FastMath.max(maxInterpolationAngleError, interpolationAngleError);
111             maxShiftRateError          = FastMath.max(maxShiftRateError, shiftRateError);
112             maxInterpolationRateError  = FastMath.max(maxInterpolationRateError, interpolationRateError);
113         }
114         Assertions.assertTrue(maxShiftAngleError > 4.0e-6);
115         Assertions.assertTrue(maxInterpolationAngleError < 1.5e-13);
116         Assertions.assertTrue(maxShiftRateError > 6.0e-8);
117         Assertions.assertTrue(maxInterpolationRateError < 2.5e-14);
118 
119         // past sample end, interpolation error should increase, but still be far better than quadratic shift
120         maxShiftAngleError         = 0;
121         maxInterpolationAngleError = 0;
122         maxShiftRateError          = 0;
123         maxInterpolationRateError  = 0;
124         for (double dt = 250.0; dt < 300.0; dt += 1.0) {
125             AbsoluteDate t          = initialOrbit.getDate().shiftedBy(dt);
126             Attitude     propagated = propagator.propagate(t).getAttitude();
127             double shiftAngleError = Rotation.distance(propagated.getRotation(),
128                                                        initialAttitude.shiftedBy(dt).getRotation());
129             double interpolationAngleError = Rotation.distance(propagated.getRotation(),
130                                                                attitudeInterpolator.interpolate(t, sample).getRotation());
131             double shiftRateError = Vector3D.distance(propagated.getSpin(),
132                                                       initialAttitude.shiftedBy(dt).getSpin());
133             double interpolationRateError = Vector3D.distance(propagated.getSpin(),
134                                                               attitudeInterpolator.interpolate(t, sample).getSpin());
135             maxShiftAngleError         = FastMath.max(maxShiftAngleError, shiftAngleError);
136             maxInterpolationAngleError = FastMath.max(maxInterpolationAngleError, interpolationAngleError);
137             maxShiftRateError          = FastMath.max(maxShiftRateError, shiftRateError);
138             maxInterpolationRateError  = FastMath.max(maxInterpolationRateError, interpolationRateError);
139         }
140         Assertions.assertTrue(maxShiftAngleError > 9.0e-6);
141         Assertions.assertTrue(maxInterpolationAngleError < 6.0e-11);
142         Assertions.assertTrue(maxShiftRateError > 9.0e-8);
143         Assertions.assertTrue(maxInterpolationRateError < 4.0e-12);
144 
145         Assertions.assertEquals(initialOrbit.getFrame(), attitudeInterpolator.getReferenceFrame());
146         Assertions.assertEquals(angularInterpolator, attitudeInterpolator.getAngularInterpolator());
147 
148     }
149 
150     @Test
151     @DisplayName("Test constructor")
152     void testConstructor() {
153         // Given
154         final Frame frameMock = mock(Frame.class);
155 
156         @SuppressWarnings("unchecked")
157         final TimeInterpolator<TimeStampedAngularCoordinates> angularInterpolatorMock = mock(TimeInterpolator.class);
158 
159         // When
160         final AttitudeInterpolator attitudeInterpolator = new AttitudeInterpolator(frameMock, angularInterpolatorMock);
161 
162         // Then
163         Assertions.assertEquals(frameMock, attitudeInterpolator.getReferenceFrame());
164         Assertions.assertEquals(angularInterpolatorMock, attitudeInterpolator.getAngularInterpolator());
165     }
166 
167     @Test
168     @DisplayName("test error thrown when sample is too small")
169     @DefaultDataContext
170     void testErrorThrownWhenSampleIsTooSmall() {
171         // Given
172         final AbsoluteDate interpolationDate = new AbsoluteDate();
173 
174         final List<Attitude> attitudes = new ArrayList<>();
175         attitudes.add(mock(Attitude.class));
176 
177         final TimeInterpolator<TimeStampedAngularCoordinates> angularInterpolator =
178                 new TimeStampedAngularCoordinatesHermiteInterpolator(2, AngularDerivativesFilter.USE_R);
179 
180         final TimeInterpolator<Attitude> attitudeInterpolator =
181                 new AttitudeInterpolator(FramesFactory.getGCRF(), angularInterpolator);
182 
183         // When & Then
184         OrekitIllegalArgumentException thrown = Assertions.assertThrows(OrekitIllegalArgumentException.class,
185                                                                         () -> attitudeInterpolator.interpolate(interpolationDate, attitudes));
186 
187         Assertions.assertEquals(OrekitMessages.NOT_ENOUGH_CACHED_NEIGHBORS, thrown.getSpecifier());
188         Assertions.assertEquals(1, ((Integer) thrown.getParts()[0]).intValue());
189         Assertions.assertEquals(2, ((Integer) thrown.getParts()[1]).intValue());
190 
191     }
192 
193     /**
194      * Test related to issue 1878.
195      *
196      * @see <a href="https://gitlab.orekit.org/orekit/orekit/-/issues/1878">Issue 1878</a>
197      */
198     @Test
199     @SuppressWarnings("unchecked")
200     void testGetSubInterpolator() {
201 
202         // Define sub interpolator
203         final TimeInterpolator<TimeStampedAngularCoordinates> mockTimeInterpolator =
204                 mock(TimeInterpolator.class);
205 
206         // GIVEN
207         final AttitudeInterpolator attitudeInterpolator =
208                 new AttitudeInterpolator(FramesFactory.getGCRF(), mockTimeInterpolator);
209 
210         // WHEN
211         final List<TimeInterpolator<?>> actualSubInterpolators = attitudeInterpolator.getSubInterpolators();
212 
213         // THEN
214         Assertions.assertEquals(1, actualSubInterpolators.size());
215         Assertions.assertSame(mockTimeInterpolator, actualSubInterpolators.get(0));
216     }
217 
218 }