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.CalculusFieldElement;
20  import org.hipparchus.Field;
21  import org.hipparchus.geometry.euclidean.threed.Rotation;
22  import org.hipparchus.geometry.euclidean.threed.RotationConvention;
23  import org.hipparchus.geometry.euclidean.threed.RotationOrder;
24  import org.hipparchus.geometry.euclidean.threed.Vector3D;
25  import org.hipparchus.random.RandomGenerator;
26  import org.hipparchus.random.Well19937a;
27  import org.hipparchus.util.Binary64Field;
28  import org.hipparchus.util.FastMath;
29  import org.hipparchus.util.MathUtils;
30  import org.junit.jupiter.api.AfterEach;
31  import org.junit.jupiter.api.Assertions;
32  import org.junit.jupiter.api.BeforeEach;
33  import org.junit.jupiter.api.Test;
34  import org.orekit.Utils;
35  import org.orekit.bodies.OneAxisEllipsoid;
36  import org.orekit.errors.OrekitException;
37  import org.orekit.errors.OrekitMessages;
38  import org.orekit.frames.Frame;
39  import org.orekit.frames.FramesFactory;
40  import org.orekit.frames.LOFType;
41  import org.orekit.orbits.CircularOrbit;
42  import org.orekit.orbits.FieldOrbit;
43  import org.orekit.orbits.Orbit;
44  import org.orekit.orbits.PositionAngleType;
45  import org.orekit.propagation.FieldSpacecraftState;
46  import org.orekit.propagation.Propagator;
47  import org.orekit.propagation.SpacecraftState;
48  import org.orekit.propagation.analytical.KeplerianPropagator;
49  import org.orekit.time.AbsoluteDate;
50  import org.orekit.time.DateComponents;
51  import org.orekit.time.FieldAbsoluteDate;
52  import org.orekit.time.TimeComponents;
53  import org.orekit.time.TimeScalesFactory;
54  import org.orekit.utils.AngularDerivativesFilter;
55  import org.orekit.utils.Constants;
56  import org.orekit.utils.IERSConventions;
57  import org.orekit.utils.PVCoordinates;
58  import org.orekit.utils.TimeStampedAngularCoordinates;
59  
60  import java.util.ArrayList;
61  import java.util.Arrays;
62  import java.util.List;
63  
64  
65  public class TabulatedLofOffsetTest {
66  
67      // Computation date
68      private AbsoluteDate date;
69  
70      // Body mu
71      private double mu;
72  
73      // Reference frame = ITRF
74      private Frame itrf;
75  
76      // Earth shape
77      OneAxisEllipsoid earth;
78  
79      //  Satellite position
80      CircularOrbit orbit;
81      PVCoordinates pvSatEME2000;
82  
83      @Test
84      void testConstantOffset() {
85  
86          RandomGenerator random = new Well19937a(0x1199d4bb8f53d2b6l);
87          for (LOFType type : LOFType.values()) {
88              for (int i = 0; i < 100; ++i) {
89                  double a1 = random.nextDouble() * MathUtils.TWO_PI;
90                  double a2 = random.nextDouble() * MathUtils.TWO_PI;
91                  double a3 = random.nextDouble() * MathUtils.TWO_PI;
92                  LofOffset law          = new LofOffset(orbit.getFrame(), type, RotationOrder.XYZ, a1, a2, a3);
93                  Rotation  offsetAtt    = law.getAttitude(orbit, orbit.getDate(), orbit.getFrame()).getRotation();
94                  LofOffset aligned      = new LofOffset(orbit.getFrame(), type);
95                  Rotation  alignedAtt   = aligned.getAttitude(orbit, orbit.getDate(), orbit.getFrame()).getRotation();
96                  Rotation  offsetProper = offsetAtt.compose(alignedAtt.revert(), RotationConvention.VECTOR_OPERATOR);
97                  TabulatedLofOffset tabulated =
98                          new TabulatedLofOffset(orbit.getFrame(), type,
99                                                 Arrays.asList(new TimeStampedAngularCoordinates(orbit.getDate().shiftedBy(-10),
100                                                                                                offsetProper,
101                                                                                                Vector3D.ZERO,
102                                                                                                Vector3D.ZERO),
103                                                              new TimeStampedAngularCoordinates(orbit.getDate().shiftedBy(0),
104                                                                                                offsetProper,
105                                                                                                Vector3D.ZERO,
106                                                                                                Vector3D.ZERO),
107                                                              new TimeStampedAngularCoordinates(orbit.getDate().shiftedBy(+10),
108                                                                                                offsetProper,
109                                                                                                Vector3D.ZERO,
110                                                                                                Vector3D.ZERO)),
111                                                2, AngularDerivativesFilter.USE_R);
112                 Rotation rebuilt = tabulated.getAttitude(orbit, orbit.getDate(), orbit.getFrame()).getRotation();
113                 Assertions.assertEquals(0.0, Rotation.distance(offsetAtt, rebuilt), 1.48e-15);
114                 Assertions.assertEquals(3, tabulated.getTable().size());
115             }
116         }
117     }
118 
119     @Test
120     void testYawCompensation() {
121 
122         // create a sample from Yaw compensation law
123         final LOFType type = LOFType.VNC;
124         final List<TimeStampedAngularCoordinates> sample = new ArrayList<TimeStampedAngularCoordinates>();
125         final AttitudeProvider yawCompensLaw =
126                 new YawCompensation(orbit.getFrame(), new NadirPointing(orbit.getFrame(), earth));
127         final Propagator originalPropagator = new KeplerianPropagator(orbit);
128         originalPropagator.setAttitudeProvider(yawCompensLaw);
129         originalPropagator.setStepHandler(1.0, currentState -> {
130                 Rotation  offsetAtt    = currentState.getAttitude().getRotation();
131                 LofOffset aligned      = new LofOffset(currentState.getFrame(), type);
132                 Rotation  alignedAtt   = aligned.getAttitude(currentState.getOrbit(), currentState.getDate(),
133                                                              currentState.getFrame()).getRotation();
134                 Rotation  offsetProper = offsetAtt.compose(alignedAtt.revert(), RotationConvention.VECTOR_OPERATOR);
135                 sample.add(new TimeStampedAngularCoordinates(currentState.getDate(),
136                                                              offsetProper, Vector3D.ZERO, Vector3D.ZERO));
137             });
138         final AbsoluteDate endDate = orbit.getDate().shiftedBy(2000);
139         originalPropagator.propagate(endDate);
140         originalPropagator.clearStepHandlers();
141 
142         // use the sample and compare it to original
143         final BoundedAttitudeProvider tabulated = new TabulatedLofOffset(orbit.getFrame(), type, sample,
144                                                                          6, AngularDerivativesFilter.USE_RR);
145         Assertions.assertEquals(0., orbit.getDate().durationFrom(tabulated.getMinDate()), Double.MIN_VALUE);
146         Assertions.assertEquals(0., endDate.durationFrom(tabulated.getMaxDate()), Double.MIN_VALUE);
147         final Propagator rebuildingPropagator = new KeplerianPropagator(orbit);
148         rebuildingPropagator.setAttitudeProvider(tabulated);
149         rebuildingPropagator.setStepHandler(0.3, currentState -> {
150                 final SpacecraftState rebuilt = originalPropagator.propagate(currentState.getDate());
151                 final Rotation r1 = currentState.getAttitude().getRotation();
152                 final Rotation r2 = rebuilt.getAttitude().getRotation();
153                 Assertions.assertEquals(0.0, Rotation.distance(r1, r2), 7.0e-6);
154                 checkField(Binary64Field.getInstance(), tabulated,
155                            currentState.getOrbit(), currentState.getDate(), currentState.getFrame());
156             });
157         rebuildingPropagator.propagate(orbit.getDate().shiftedBy(50), orbit.getDate().shiftedBy(1950));
158 
159     }
160 
161     private <T extends CalculusFieldElement<T>> void checkField(final Field<T> field, final AttitudeProvider provider,
162                                                                 final Orbit orbit, final AbsoluteDate date,
163                                                                 final Frame frame)
164         {
165         Attitude attitudeD = provider.getAttitude(orbit, date, frame);
166         final FieldOrbit<T> orbitF = new FieldSpacecraftState<>(field, new SpacecraftState(orbit)).getOrbit();
167         final FieldAbsoluteDate<T> dateF = new FieldAbsoluteDate<>(field, date);
168         FieldAttitude<T> attitudeF = provider.getAttitude(orbitF, dateF, frame);
169         Assertions.assertEquals(0.0, Rotation.distance(attitudeD.getRotation(), attitudeF.getRotation().toRotation()), 1.0e-15);
170         Assertions.assertEquals(0.0, Vector3D.distance(attitudeD.getSpin(), attitudeF.getSpin().toVector3D()), 1.0e-15);
171         Assertions.assertEquals(0.0, Vector3D.distance(attitudeD.getRotationAcceleration(), attitudeF.getRotationAcceleration().toVector3D()), 1.0e-15);
172     }
173 
174     @Test
175     void testNonPseudoInertialFrame() {
176         final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
177         try {
178             new TabulatedLofOffset(itrf, LOFType.QSW,
179                                    Arrays.asList(new TimeStampedAngularCoordinates(orbit.getDate().shiftedBy(-10),
180                                                                                    Rotation.IDENTITY, Vector3D.ZERO, Vector3D.ZERO),
181                                                  new TimeStampedAngularCoordinates(orbit.getDate().shiftedBy(0),
182                                                                                    Rotation.IDENTITY, Vector3D.ZERO, Vector3D.ZERO),
183                                                  new TimeStampedAngularCoordinates(orbit.getDate().shiftedBy(+10),
184                                                                                    Rotation.IDENTITY, Vector3D.ZERO, Vector3D.ZERO)),
185                                    2, AngularDerivativesFilter.USE_R);
186         } catch (OrekitException oe) {
187             Assertions.assertEquals(OrekitMessages.NON_PSEUDO_INERTIAL_FRAME, oe.getSpecifier());
188             Assertions.assertEquals(itrf.getName(), oe.getParts()[0]);
189         }
190     }
191 
192     @BeforeEach
193     public void setUp() {
194         try {
195 
196             Utils.setDataRoot("regular-data");
197 
198             // Computation date
199             date = new AbsoluteDate(new DateComponents(2008, 04, 07),
200                                     TimeComponents.H00,
201                                     TimeScalesFactory.getUTC());
202 
203             // Body mu
204             mu = 3.9860047e14;
205 
206             // Reference frame = ITRF
207             itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
208 
209             // Elliptic earth shape
210             earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
211                                          Constants.WGS84_EARTH_FLATTENING,
212                                          itrf);
213 
214             //  Satellite position
215             orbit =
216                 new CircularOrbit(7178000.0, 0.5e-8, -0.5e-8, FastMath.toRadians(50.), FastMath.toRadians(150.),
217                                        FastMath.toRadians(5.300), PositionAngleType.MEAN,
218                                        FramesFactory.getEME2000(), date, mu);
219             pvSatEME2000 = orbit.getPVCoordinates();
220 
221 
222         } catch (OrekitException oe) {
223             Assertions.fail(oe.getMessage());
224         }
225 
226     }
227 
228     @AfterEach
229     public void tearDown() {
230         date = null;
231         itrf = null;
232         earth = null;
233     }
234 
235 }
236