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.FieldRotation;
22  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
23  import org.hipparchus.geometry.euclidean.threed.Rotation;
24  import org.hipparchus.geometry.euclidean.threed.Vector3D;
25  import org.hipparchus.util.Binary64;
26  import org.hipparchus.util.Binary64Field;
27  import org.hipparchus.util.FastMath;
28  import org.junit.jupiter.api.Assertions;
29  import org.junit.jupiter.api.BeforeEach;
30  import org.junit.jupiter.api.Test;
31  import org.orekit.Utils;
32  import org.orekit.annotation.DefaultDataContext;
33  import org.orekit.frames.Frame;
34  import org.orekit.frames.FramesFactory;
35  import org.orekit.orbits.*;
36  import org.orekit.propagation.FieldPropagator;
37  import org.orekit.propagation.FieldSpacecraftState;
38  import org.orekit.propagation.Propagator;
39  import org.orekit.propagation.SpacecraftState;
40  import org.orekit.propagation.analytical.FieldKeplerianPropagator;
41  import org.orekit.propagation.analytical.KeplerianPropagator;
42  import org.orekit.time.AbsoluteDate;
43  import org.orekit.time.DateComponents;
44  import org.orekit.time.FieldAbsoluteDate;
45  import org.orekit.time.TimeComponents;
46  import org.orekit.time.TimeScalesFactory;
47  import org.orekit.utils.AngularCoordinates;
48  import org.orekit.utils.FieldAngularCoordinates;
49  import org.orekit.utils.FieldPVCoordinates;
50  import org.orekit.utils.PVCoordinates;
51  
52  public class FixedRateTest {
53  
54      @Test
55      @DefaultDataContext
56      void testZeroRate() {
57          AbsoluteDate date = new AbsoluteDate(new DateComponents(2004, 3, 2),
58                                               new TimeComponents(13, 17, 7.865),
59                                               TimeScalesFactory.getUTC());
60          final Frame frame = FramesFactory.getEME2000();
61          FixedRate law = new FixedRate(new Attitude(date, frame,
62                                                     new Rotation(0.48, 0.64, 0.36, 0.48, false),
63                                                     Vector3D.ZERO, Vector3D.ZERO));
64          PVCoordinates pv =
65              new PVCoordinates(new Vector3D(28812595.32012577, 5948437.4640250085, 0),
66                                new Vector3D(0, 0, 3680.853673522056));
67          Orbit orbit = new KeplerianOrbit(pv, frame, date, 3.986004415e14);
68          Rotation attitude0 = law.getAttitude(orbit, date, frame).getRotation();
69          Assertions.assertEquals(0, Rotation.distance(attitude0, law.getReferenceAttitude().getRotation()), 1.0e-10);
70          Rotation attitude1 = law.getAttitude(orbit.shiftedBy(10.0), date.shiftedBy(10.0), frame).getRotation();
71          Assertions.assertEquals(0, Rotation.distance(attitude1, law.getReferenceAttitude().getRotation()), 1.0e-10);
72          Rotation attitude2 = law.getAttitude(orbit.shiftedBy(20.0), date.shiftedBy(20.0), frame).getRotation();
73          Assertions.assertEquals(0, Rotation.distance(attitude2, law.getReferenceAttitude().getRotation()), 1.0e-10);
74  
75      }
76  
77      @Test
78      @DefaultDataContext
79      void testNonZeroRate() {
80          final AbsoluteDate date = new AbsoluteDate(new DateComponents(2004, 3, 2),
81                                                     new TimeComponents(13, 17, 7.865),
82                                                     TimeScalesFactory.getUTC());
83          final double rate = 2 * FastMath.PI / (12 * 60);
84          final Frame frame = FramesFactory.getEME2000();
85          final Frame gcrf  = FramesFactory.getGCRF();
86          FixedRate law = new FixedRate(new Attitude(date, frame,
87                                                     new Rotation(0.48, 0.64, 0.36, 0.48, false),
88                                                     new Vector3D(rate, Vector3D.PLUS_K), Vector3D.ZERO));
89          final Rotation ref = law.getReferenceAttitude().getRotation().applyTo(gcrf.getTransformTo(frame, date).getRotation());
90          PVCoordinates pv =
91              new PVCoordinates(new Vector3D(28812595.32012577, 5948437.4640250085, 0),
92                                new Vector3D(0, 0, 3680.853673522056));
93          Orbit orbit = new KeplerianOrbit(pv, FramesFactory.getEME2000(), date, 3.986004415e14);
94          Rotation attitude0 = law.getAttitude(orbit, date, gcrf).getRotation();
95          Assertions.assertEquals(0, Rotation.distance(attitude0, ref), 1.0e-10);
96          Rotation attitude1 = law.getAttitude(orbit.shiftedBy(10.0), date.shiftedBy(10.0), gcrf).getRotation();
97          Assertions.assertEquals(10 * rate, Rotation.distance(attitude1, ref), 1.0e-10);
98          Rotation attitude2 = law.getAttitude(orbit.shiftedBy(-20.0), date.shiftedBy(-20.0), gcrf).getRotation();
99          Assertions.assertEquals(20 * rate, Rotation.distance(attitude2, ref), 1.0e-10);
100         Assertions.assertEquals(30 * rate, Rotation.distance(attitude2, attitude1), 1.0e-10);
101         Rotation attitude3 = law.getAttitude(orbit.shiftedBy(0.0), date, frame).getRotation();
102         Assertions.assertEquals(0, Rotation.distance(attitude3, law.getReferenceAttitude().getRotation()), 1.0e-10);
103 
104     }
105 
106     @Test
107     @DefaultDataContext
108     void testSpin() {
109 
110         AbsoluteDate date = new AbsoluteDate(new DateComponents(1970, 1, 1),
111                                              new TimeComponents(3, 25, 45.6789),
112                                              TimeScalesFactory.getUTC());
113 
114         final double rate = 2 * FastMath.PI / (12 * 60);
115         AttitudeProvider law =
116             new FixedRate(new Attitude(date, FramesFactory.getEME2000(),
117                                        new Rotation(0.48, 0.64, 0.36, 0.48, false),
118                                        new Vector3D(rate, Vector3D.PLUS_K),
119                                        Vector3D.ZERO));
120 
121         KeplerianOrbit orbit =
122             new KeplerianOrbit(7178000.0, 1.e-4, FastMath.toRadians(50.),
123                               FastMath.toRadians(10.), FastMath.toRadians(20.),
124                               FastMath.toRadians(30.), PositionAngleType.MEAN,
125                               FramesFactory.getEME2000(), date, 3.986004415e14);
126 
127         Propagator propagator = new KeplerianPropagator(orbit, law);
128 
129         double h = 0.01;
130         SpacecraftState sMinus = propagator.propagate(date.shiftedBy(-h));
131         SpacecraftState s0     = propagator.propagate(date);
132         SpacecraftState sPlus  = propagator.propagate(date.shiftedBy(h));
133 
134         // check spin is consistent with attitude evolution
135         double errorAngleMinus     = Rotation.distance(sMinus.shiftedBy(h).getAttitude().getRotation(),
136                                                        s0.getAttitude().getRotation());
137         double evolutionAngleMinus = Rotation.distance(sMinus.getAttitude().getRotation(),
138                                                        s0.getAttitude().getRotation());
139         Assertions.assertEquals(0.0, errorAngleMinus, 1.0e-6 * evolutionAngleMinus);
140         double errorAnglePlus      = Rotation.distance(s0.getAttitude().getRotation(),
141                                                        sPlus.shiftedBy(-h).getAttitude().getRotation());
142         double evolutionAnglePlus  = Rotation.distance(s0.getAttitude().getRotation(),
143                                                        sPlus.getAttitude().getRotation());
144         Assertions.assertEquals(0.0, errorAnglePlus, 1.0e-6 * evolutionAnglePlus);
145 
146         Vector3D spin0 = s0.getAttitude().getSpin();
147         Vector3D reference = AngularCoordinates.estimateRate(sMinus.getAttitude().getRotation(),
148                                                              sPlus.getAttitude().getRotation(),
149                                                              2 * h);
150         Assertions.assertEquals(0.0, spin0.subtract(reference).getNorm(), 1.0e-14);
151 
152     }
153 
154     @Test
155     @DefaultDataContext
156     void testZeroRateField() {
157         doTestZeroRate(Binary64Field.getInstance());
158     }
159 
160     @DefaultDataContext
161     private <T extends CalculusFieldElement<T>> void doTestZeroRate(final Field<T> field) {
162         final T zero = field.getZero();
163         FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field,
164                                                             new DateComponents(2004, 3, 2),
165                                                             new TimeComponents(13, 17, 7.865),
166                                                             TimeScalesFactory.getUTC());
167         final Frame frame = FramesFactory.getEME2000();
168         final Frame gcrf  = FramesFactory.getGCRF();
169         FixedRate law = new FixedRate(new Attitude(date.toAbsoluteDate(), frame,
170                                                    new Rotation(0.48, 0.64, 0.36, 0.48, false),
171                                                    Vector3D.ZERO, Vector3D.ZERO));
172         final Rotation ref = law.getReferenceAttitude().getRotation().applyTo(gcrf.getTransformTo(frame, date.toAbsoluteDate()).getRotation());
173         FieldPVCoordinates<T> pv =
174             new FieldPVCoordinates<>(field.getOne(),
175                                      new PVCoordinates(new Vector3D(28812595.32012577, 5948437.4640250085, 0),
176                                                        new Vector3D(0, 0, 3680.853673522056)));
177         FieldOrbit<T> orbit = new FieldKeplerianOrbit<>(pv, frame, date, zero.add(3.986004415e14));
178         FieldRotation<T> attitude0 = law.getAttitude(orbit, date, gcrf).getRotation();
179         Assertions.assertEquals(0, Rotation.distance(attitude0.toRotation(), ref), 1.0e-10);
180         FieldRotation<T> attitude1 = law.getAttitude(orbit.shiftedBy(zero.add(10.0)), date.shiftedBy(10.0), gcrf).getRotation();
181         Assertions.assertEquals(0, Rotation.distance(attitude1.toRotation(), ref), 1.0e-10);
182         FieldRotation<T> attitude2 = law.getAttitude(orbit.shiftedBy(zero.add(20.0)), date.shiftedBy(20.0), gcrf).getRotation();
183         Assertions.assertEquals(0, Rotation.distance(attitude2.toRotation(), ref), 1.0e-10);
184 
185     }
186 
187     @Test
188     @DefaultDataContext
189     void testNonZeroRateField() {
190         doTestNonZeroRate(Binary64Field.getInstance());
191     }
192 
193     @DefaultDataContext
194     private <T extends CalculusFieldElement<T>> void doTestNonZeroRate(final Field<T> field) {
195         final T zero = field.getZero();
196         FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field,
197                                                             new DateComponents(2004, 3, 2),
198                                                             new TimeComponents(13, 17, 7.865),
199                                                             TimeScalesFactory.getUTC());
200         final T rate = zero.add(2 * FastMath.PI / (12 * 60));
201         final Frame frame = FramesFactory.getEME2000();
202         FixedRate law = new FixedRate(new Attitude(date.toAbsoluteDate(), frame,
203                                                    new Rotation(0.48, 0.64, 0.36, 0.48, false),
204                                                    new Vector3D(rate.getReal(), Vector3D.PLUS_K), Vector3D.ZERO));
205         FieldPVCoordinates<T> pv =
206                         new FieldPVCoordinates<>(field.getOne(),
207                                                  new PVCoordinates(new Vector3D(28812595.32012577, 5948437.4640250085, 0),
208                                                                    new Vector3D(0, 0, 3680.853673522056)));
209         FieldOrbit<T> orbit = new FieldKeplerianOrbit<>(pv, FramesFactory.getEME2000(), date, zero.add(3.986004415e14));
210         FieldRotation<T> attitude0 = law.getAttitude(orbit, date, frame).getRotation();
211         Assertions.assertEquals(0, Rotation.distance(attitude0.toRotation(), law.getReferenceAttitude().getRotation()), 1.0e-10);
212         FieldRotation<T> attitude1 = law.getAttitude(orbit.shiftedBy(zero.add(10.0)), date.shiftedBy(10.0), frame).getRotation();
213         Assertions.assertEquals(10 * rate.getReal(), Rotation.distance(attitude1.toRotation(), law.getReferenceAttitude().getRotation()), 1.0e-10);
214         FieldRotation<T> attitude2 = law.getAttitude(orbit.shiftedBy(zero.add(-20.0)), date.shiftedBy(-20.0), frame).getRotation();
215         Assertions.assertEquals(20 * rate.getReal(), Rotation.distance(attitude2.toRotation(), law.getReferenceAttitude().getRotation()), 1.0e-10);
216         Assertions.assertEquals(30 * rate.getReal(), Rotation.distance(attitude2.toRotation(), attitude1.toRotation()), 1.0e-10);
217         FieldRotation<T> attitude3 = law.getAttitude(orbit.shiftedBy(zero.add(0.0)), date, frame).getRotation();
218         Assertions.assertEquals(0, Rotation.distance(attitude3.toRotation(), law.getReferenceAttitude().getRotation()), 1.0e-10);
219 
220     }
221 
222     @Test
223     @DefaultDataContext
224     void testSpinField() {
225         doTestSpin(Binary64Field.getInstance());
226     }
227 
228     @DefaultDataContext
229     private <T extends CalculusFieldElement<T>> void doTestSpin(final Field<T> field) {
230 
231         final T zero = field.getZero();
232         FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field,
233                                                             new DateComponents(1970, 1, 1),
234                                                             new TimeComponents(3, 25, 45.6789),
235                                                             TimeScalesFactory.getUTC());
236 
237         final T rate = zero.add(2 * FastMath.PI / (12 * 60));
238         AttitudeProvider law =
239                         new FixedRate(new Attitude(date.toAbsoluteDate(), FramesFactory.getEME2000(),
240                                                    new Rotation(0.48, 0.64, 0.36, 0.48, false),
241                                                    new Vector3D(rate.getReal(), Vector3D.PLUS_K),
242                                                    Vector3D.ZERO));
243 
244         FieldKeplerianOrbit<T> orbit =
245             new FieldKeplerianOrbit<>(zero.add(7178000.0),
246                                       zero.add(1.e-4),
247                                       zero.add(FastMath.toRadians(50.)),
248                                       zero.add(FastMath.toRadians(10.)),
249                                       zero.add(FastMath.toRadians(20.)),
250                                       zero.add(FastMath.toRadians(30.)), PositionAngleType.MEAN,
251                                       FramesFactory.getEME2000(), date, zero.add(3.986004415e14));
252 
253         FieldPropagator<T> propagator = new FieldKeplerianPropagator<>(orbit, law);
254 
255         T h = zero.add(0.01);
256         FieldSpacecraftState<T> sMinus = propagator.propagate(date.shiftedBy(h.negate()));
257         FieldSpacecraftState<T> s0     = propagator.propagate(date);
258         FieldSpacecraftState<T> sPlus  = propagator.propagate(date.shiftedBy(h));
259 
260         // check spin is consistent with attitude evolution
261         double errorAngleMinus     = FieldRotation.distance(sMinus.shiftedBy(h).getAttitude().getRotation(),
262                                                             s0.getAttitude().getRotation()).getReal();
263         double evolutionAngleMinus = FieldRotation.distance(sMinus.getAttitude().getRotation(),
264                                                             s0.getAttitude().getRotation()).getReal();
265         Assertions.assertEquals(0.0, errorAngleMinus, 1.0e-6 * evolutionAngleMinus);
266         double errorAnglePlus      = FieldRotation.distance(s0.getAttitude().getRotation(),
267                                                             sPlus.shiftedBy(h.negate()).getAttitude().getRotation()).getReal();
268         double evolutionAnglePlus  = FieldRotation.distance(s0.getAttitude().getRotation(),
269                                                             sPlus.getAttitude().getRotation()).getReal();
270         Assertions.assertEquals(0.0, errorAnglePlus, 1.0e-6 * evolutionAnglePlus);
271 
272         FieldVector3D<T> spin0 = s0.getAttitude().getSpin();
273         FieldVector3D<T> reference = FieldAngularCoordinates.estimateRate(sMinus.getAttitude().getRotation(),
274                                                                           sPlus.getAttitude().getRotation(),
275                                                                           h.multiply(2));
276         Assertions.assertEquals(0.0, spin0.subtract(reference).getNorm().getReal(), 1.0e-14);
277 
278     }
279 
280     @Test
281     void testGetAttitudeRotation() {
282         // GIVEN
283         final PVCoordinates pv =
284                 new PVCoordinates(new Vector3D(28812595.32012577, 5948437.4640250085, 0),
285                         new Vector3D(0, 0, 3680.853673522056));
286         final Orbit orbit = new CartesianOrbit(pv, FramesFactory.getGCRF(), AbsoluteDate.ARBITRARY_EPOCH, 3.986004415e14);
287         final Attitude attitude = new Attitude(orbit.getDate().shiftedBy(-100.), FramesFactory.getEME2000(),
288                 new Rotation(0.48, 0.64, 0.36, 0.48, false),
289                 Vector3D.ZERO, Vector3D.ZERO);
290         final FixedRate fixedRate = new FixedRate(attitude);
291         // WHEN
292         final Rotation actualRotation = fixedRate.getAttitudeRotation(orbit, orbit.getDate(), orbit.getFrame());
293         // THEN
294         final Rotation expectedRotation = fixedRate.getAttitude(orbit, orbit.getDate(), orbit.getFrame()).getRotation();
295         Assertions.assertEquals(0., Rotation.distance(expectedRotation, actualRotation), 2e-7);
296     }
297 
298     @Test
299     void testFieldGetAttitudeRotation() {
300         // GIVEN
301         final PVCoordinates pv =
302                 new PVCoordinates(new Vector3D(28812595.32012577, 5948437.4640250085, 0),
303                         new Vector3D(0, 0, 3680.853673522056));
304         final Orbit orbit = new CartesianOrbit(pv, FramesFactory.getGCRF(), AbsoluteDate.ARBITRARY_EPOCH, 3.986004415e14);
305         final FieldOrbit<Binary64> fieldOrbit = new FieldCartesianOrbit<>(Binary64Field.getInstance(), orbit);
306         final Attitude attitude = new Attitude(orbit.getDate().shiftedBy(-100.), FramesFactory.getEME2000(),
307                 new Rotation(0.48, 0.64, 0.36, 0.48, false),
308                 Vector3D.ZERO, Vector3D.ZERO);
309         final FixedRate fixedRate = new FixedRate(attitude);
310         // WHEN
311         final FieldRotation<Binary64> actualRotation = fixedRate.getAttitudeRotation(fieldOrbit, fieldOrbit.getDate(), orbit.getFrame());
312         // THEN
313         final FieldRotation<Binary64> expectedRotation = fixedRate.getAttitude(fieldOrbit, fieldOrbit.getDate(), orbit.getFrame()).getRotation();
314         Assertions.assertEquals(0., Rotation.distance(expectedRotation.toRotation(), actualRotation.toRotation()), 2e-7);
315     }
316 
317     @BeforeEach
318     public void setUp() {
319         Utils.setDataRoot("regular-data");
320     }
321 
322 }
323