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