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.propagation.analytical.tle;
18  
19  import org.hamcrest.MatcherAssert;
20  import org.hipparchus.CalculusFieldElement;
21  import org.hipparchus.Field;
22  import org.hipparchus.analysis.differentiation.Gradient;
23  import org.hipparchus.analysis.differentiation.GradientField;
24  import org.hipparchus.geometry.euclidean.threed.FieldLine;
25  import org.hipparchus.geometry.euclidean.threed.FieldRotation;
26  import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
27  import org.hipparchus.geometry.euclidean.threed.Vector3D;
28  import org.hipparchus.util.Binary64Field;
29  import org.hipparchus.util.FastMath;
30  import org.junit.jupiter.api.Assertions;
31  import org.junit.jupiter.api.BeforeEach;
32  import org.junit.jupiter.api.Test;
33  import org.orekit.OrekitMatchers;
34  import org.orekit.Utils;
35  import org.orekit.attitudes.BodyCenterPointing;
36  import org.orekit.bodies.OneAxisEllipsoid;
37  import org.orekit.frames.FieldStaticTransform;
38  import org.orekit.frames.Frame;
39  import org.orekit.frames.FramesFactory;
40  import org.orekit.propagation.FieldBoundedPropagator;
41  import org.orekit.propagation.FieldEphemerisGenerator;
42  import org.orekit.propagation.FieldSpacecraftState;
43  import org.orekit.propagation.Propagator;
44  import org.orekit.propagation.SpacecraftState;
45  import org.orekit.propagation.sampling.FieldOrekitFixedStepHandler;
46  import org.orekit.time.AbsoluteDate;
47  import org.orekit.time.FieldAbsoluteDate;
48  import org.orekit.utils.Constants;
49  import org.orekit.utils.FieldPVCoordinates;
50  import org.orekit.utils.IERSConventions;
51  import org.orekit.utils.PVCoordinates;
52  
53  
54  public class FieldTLEPropagatorTest {
55  
56      private double period;
57  
58      @Test
59      public void testsecondaryMode() {
60          doTestsecondaryMode(Binary64Field.getInstance());
61      }
62  
63      @Test
64      public void testEphemerisMode() {
65          doTestEphemerisMode(Binary64Field.getInstance());
66      }
67  
68      @Test
69      public void testBodyCenterInPointingDirection() {
70          doTestBodyCenterInPointingDirection(Binary64Field.getInstance());
71      }
72  
73      @Test
74      public void testComparisonWithNonField() {
75          doTestComparisonWithNonField(Binary64Field.getInstance());
76      }
77  
78      public <T extends CalculusFieldElement<T>> void doTestsecondaryMode(Field<T> field) {
79  
80          // setup a TLE for a GPS satellite
81          String line1 = "1 37753U 11036A   12090.13205652 -.00000006  00000-0  00000+0 0  2272";
82          String line2 = "2 37753  55.0032 176.5796 0004733  13.2285 346.8266  2.00565440  5153";
83          FieldTLE<T> tle = new FieldTLE<>(field, line1, line2);
84          
85          final T[] parameters = tle.getParameters(field);
86          FieldTLEPropagator<T> propagator = FieldTLEPropagator.selectExtrapolator(tle, parameters);
87          FieldAbsoluteDate<T> initDate = tle.getDate();
88          FieldSpacecraftState<T> initialState = propagator.getInitialState();
89  
90          // Simulate a full period of a GPS satellite
91          // -----------------------------------------
92          FieldSpacecraftState<T> finalState = propagator.propagate(initDate.shiftedBy(period));
93  
94          // Check results
95          Assertions.assertEquals(initialState.getOrbit().getA().getReal(), finalState.getOrbit().getA().getReal(), 1e-1);
96          Assertions.assertEquals(initialState.getOrbit().getEquinoctialEx().getReal(), finalState.getOrbit().getEquinoctialEx().getReal(), 1e-1);
97          Assertions.assertEquals(initialState.getOrbit().getEquinoctialEy().getReal(), finalState.getOrbit().getEquinoctialEy().getReal(), 1e-1);
98          Assertions.assertEquals(initialState.getOrbit().getHx().getReal(), finalState.getOrbit().getHx().getReal(), 1e-3);
99          Assertions.assertEquals(initialState.getOrbit().getHy().getReal(), finalState.getOrbit().getHy().getReal(), 1e-3);
100         Assertions.assertEquals(initialState.getOrbit().getLM().getReal(), finalState.getOrbit().getLM().getReal(), 1e-3);
101 
102     }
103 
104     public <T extends CalculusFieldElement<T>> void doTestEphemerisMode(Field<T> field) {
105 
106         // setup a TLE for a GPS satellite
107         String line1 = "1 37753U 11036A   12090.13205652 -.00000006  00000-0  00000+0 0  2272";
108         String line2 = "2 37753  55.0032 176.5796 0004733  13.2285 346.8266  2.00565440  5153";
109         FieldTLE<T> tle = new FieldTLE<>(field, line1, line2);
110         
111         final T[] parameters = tle.getParameters(field);
112         FieldTLEPropagator<T> propagator = FieldTLEPropagator.selectExtrapolator(tle, parameters);
113         final FieldEphemerisGenerator<T> generator = propagator.getEphemerisGenerator();
114 
115         FieldAbsoluteDate<T> initDate = tle.getDate();
116         FieldSpacecraftState<T> initialState = propagator.getInitialState();
117 
118         // Simulate a full period of a GPS satellite
119         // -----------------------------------------
120         FieldAbsoluteDate<T> endDate = initDate.shiftedBy(period);
121         propagator.propagate(endDate);
122 
123         // get the ephemeris
124         FieldBoundedPropagator<T> boundedProp = generator.getGeneratedEphemeris();
125 
126         // get the initial state from the ephemeris and check if it is the same as
127         // the initial state from the TLE
128         FieldSpacecraftState<T> boundedState = boundedProp.propagate(initDate);
129 
130         // Check results
131         Assertions.assertEquals(initialState.getOrbit().getA().getReal(), boundedState.getOrbit().getA().getReal(), 0.);
132         Assertions.assertEquals(initialState.getOrbit().getEquinoctialEx().getReal(), boundedState.getOrbit().getEquinoctialEx().getReal(), 0.);
133         Assertions.assertEquals(initialState.getOrbit().getEquinoctialEy().getReal(), boundedState.getOrbit().getEquinoctialEy().getReal(), 0.);
134         Assertions.assertEquals(initialState.getOrbit().getHx().getReal(), boundedState.getOrbit().getHx().getReal(), 0.);
135         Assertions.assertEquals(initialState.getOrbit().getHy().getReal(), boundedState.getOrbit().getHy().getReal(), 0.);
136         Assertions.assertEquals(initialState.getOrbit().getLM().getReal(), boundedState.getOrbit().getLM().getReal(), 1e-14);
137 
138         FieldSpacecraftState<T> finalState = boundedProp.propagate(endDate);
139 
140         // Check results
141         Assertions.assertEquals(initialState.getOrbit().getA().getReal(), finalState.getOrbit().getA().getReal(), 1e-1);
142         Assertions.assertEquals(initialState.getOrbit().getEquinoctialEx().getReal(), finalState.getOrbit().getEquinoctialEx().getReal(), 1e-1);
143         Assertions.assertEquals(initialState.getOrbit().getEquinoctialEy().getReal(), finalState.getOrbit().getEquinoctialEy().getReal(), 1e-1);
144         Assertions.assertEquals(initialState.getOrbit().getHx().getReal(), finalState.getOrbit().getHx().getReal(), 1e-3);
145         Assertions.assertEquals(initialState.getOrbit().getHy().getReal(), finalState.getOrbit().getHy().getReal(), 1e-3);
146         Assertions.assertEquals(initialState.getOrbit().getLM().getReal(), finalState.getOrbit().getLM().getReal(), 1e-3);
147 
148     }
149 
150     /** Test if body center belongs to the direction pointed by the satellite
151      */
152     public <T extends CalculusFieldElement<T>> void doTestBodyCenterInPointingDirection(Field<T> field) {
153 
154         // setup a TLE for a GPS satellite
155         String line1 = "1 37753U 11036A   12090.13205652 -.00000006  00000-0  00000+0 0  2272";
156         String line2 = "2 37753  55.0032 176.5796 0004733  13.2285 346.8266  2.00565440  5153";
157         FieldTLE<T> tle = new FieldTLE<>(field, line1, line2);
158 
159         // setup a 0 T element.
160         T T_zero = field.getZero();
161 
162         final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
163         final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
164                                                             Constants.WGS84_EARTH_FLATTENING,
165                                                             itrf);
166         FieldDistanceChecker<T> checker = new FieldDistanceChecker<T>(itrf);
167 
168         // with Earth pointing attitude, distance should be small
169         final T[] parameters = tle.getParameters(field);
170         FieldTLEPropagator<T> propagator =
171                 FieldTLEPropagator.selectExtrapolator(tle,
172                                                  new BodyCenterPointing(FramesFactory.getTEME(), earth),
173                                                  T_zero.add(Propagator.DEFAULT_MASS), parameters);
174         propagator.setStepHandler(T_zero.add(900.0), checker);
175         propagator.propagate(tle.getDate().shiftedBy(period));
176         Assertions.assertEquals(0.0, checker.getMaxDistance(), 2.0e-7);
177 
178         // with default attitude mode, distance should be large
179         propagator = FieldTLEPropagator.selectExtrapolator(tle, parameters);
180         propagator.setStepHandler(T_zero.add(900.0), checker);
181         propagator.propagate(tle.getDate().shiftedBy(period));
182         MatcherAssert.assertThat(checker.getMinDistance(),
183                 OrekitMatchers.greaterThan(1.5218e7));
184         Assertions.assertEquals(2.6572e7, checker.getMaxDistance(), 1000.0);
185 
186     }
187 
188     private static class FieldDistanceChecker<T extends CalculusFieldElement<T>> implements FieldOrekitFixedStepHandler<T> {
189 
190         private final Frame itrf;
191         private double minDistance;
192         private double maxDistance;
193 
194         public FieldDistanceChecker(Frame itrf) {
195             this.itrf = itrf;
196         }
197 
198         public double getMinDistance() {
199             return minDistance;
200         }
201 
202         public double getMaxDistance() {
203             return maxDistance;
204         }
205 
206         @Override
207         public void init(FieldSpacecraftState<T> s0, FieldAbsoluteDate<T> t, T step) {
208             minDistance = Double.POSITIVE_INFINITY;
209             maxDistance = Double.NEGATIVE_INFINITY;
210         }
211 
212         @Override
213         public void handleStep(FieldSpacecraftState<T> currentState) {
214             // Get satellite attitude rotation, i.e rotation from inertial frame to satellite frame
215             FieldRotation<T> rotSat = currentState.getAttitude().getRotation();
216 
217             // Transform Z axis from satellite frame to inertial frame
218             FieldVector3D<T> zSat = rotSat.applyInverseTo(Vector3D.PLUS_K);
219 
220             // Transform Z axis from inertial frame to ITRF
221             FieldStaticTransform<T> transform = currentState.getFrame().getStaticTransformTo(itrf, currentState.getDate());
222             FieldVector3D<T> zSatITRF = transform.transformVector(zSat);
223 
224             // Transform satellite position/velocity from inertial frame to ITRF
225             FieldVector3D<T> posSatITRF = transform.transformPosition(currentState.getPosition());
226 
227             // Line containing satellite point and following pointing direction
228             FieldLine<T> pointingLine = new FieldLine<T>(posSatITRF, posSatITRF.add(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, zSatITRF), 1.0e-10);
229 
230             double distance = pointingLine.distance(Vector3D.ZERO).getReal();
231             minDistance = FastMath.min(minDistance, distance);
232             maxDistance = FastMath.max(maxDistance, distance);
233         }
234 
235     }
236 
237     public <T extends CalculusFieldElement<T>> void doTestComparisonWithNonField(Field<T> field) {
238 
239         // propagation time.
240         final double propagtime = 10 * 60;
241 
242         // setup a TLE for a GPS satellite
243         String line1 = "1 37753U 11036A   12090.13205652 -.00000006  00000-0  00000+0 0  2272";
244         String line2 = "2 37753  55.0032 176.5796 0004733  13.2285 346.8266  2.00565440  5153";
245 
246         // build FieldTLE and TLE instances for GPS
247         FieldTLE<T> fieldtleGPS = new FieldTLE<>(field, line1, line2);
248         TLE tleGPS = new TLE(line1, line2);
249 
250         // setup a TLE for ISS
251         String line3 = "1 25544U 98067A   20162.14487814  .00001100  00000-0  27734-4 0  9997";
252         String line4 = "2 25544  51.6445  23.3222 0002345  38.1770 106.6280 15.49436440230920";
253 
254         // build FieldTLE and TLE instances for ISS
255         FieldTLE<T> fieldtleISS = new FieldTLE<>(field, line3, line4);
256         TLE tleISS = new TLE(line3, line4);
257 
258         // propagate Field GPS orbit
259         final T[] parametersGPS = fieldtleGPS.getParameters(field);
260         FieldTLEPropagator<T> fieldpropagator = FieldTLEPropagator.selectExtrapolator(fieldtleGPS, parametersGPS);
261         FieldAbsoluteDate<T> fieldinitDate = fieldtleGPS.getDate();
262         FieldAbsoluteDate<T> fieldendDate = fieldinitDate.shiftedBy(propagtime);
263         FieldPVCoordinates<T> fieldfinalGPS = fieldpropagator.getPVCoordinates(fieldendDate, parametersGPS);
264 
265         // propagate GPS orbit
266         TLEPropagator propagator = TLEPropagator.selectExtrapolator(tleGPS);
267         AbsoluteDate initDate = tleGPS.getDate();
268         AbsoluteDate endDate = initDate.shiftedBy(propagtime);
269         PVCoordinates finalGPS = propagator.getPVCoordinates(endDate);
270 
271         // propagate Field ISS orbit
272         final T[] parametersISS = fieldtleISS.getParameters(field);
273         fieldpropagator = FieldTLEPropagator.selectExtrapolator(fieldtleISS, parametersISS);
274         fieldinitDate = fieldtleISS.getDate();
275         fieldendDate = fieldinitDate.shiftedBy(propagtime);
276         FieldSpacecraftState<T> fieldfinalISS = fieldpropagator.propagate(fieldendDate);
277 
278         // propagate ISS orbit
279         propagator = TLEPropagator.selectExtrapolator(tleISS);
280         initDate = tleISS.getDate();
281         endDate = initDate.shiftedBy(propagtime);
282         SpacecraftState finalISS = propagator.propagate(endDate);
283 
284         // check
285         Assertions.assertEquals(0, Vector3D.distance(finalGPS.getPosition(), fieldfinalGPS.getPosition().toVector3D()), 3.8e-9);
286         Assertions.assertEquals(0, Vector3D.distance(finalGPS.getVelocity(), fieldfinalGPS.getVelocity().toVector3D()), 0.);
287 
288         Assertions.assertEquals(0, Vector3D.distance(finalISS.getPosition(), fieldfinalISS.getPVCoordinates().getPosition().toVector3D()), 0.);
289         Assertions.assertEquals(0, Vector3D.distance(finalISS.getPVCoordinates().getVelocity(), fieldfinalISS.getPVCoordinates().getVelocity().toVector3D()), 0.);
290 
291     }
292 
293     @Test
294     void testResetInitialState() {
295         // GIVEN
296         final FieldTLE<Gradient> tle = getGradientTLE();
297         final Field<Gradient> field = tle.getDate().getField();
298         final Gradient[] parameters = tle.getParameters(field);
299         final FieldTLEPropagator<Gradient> tlePropagator = FieldTLEPropagator.selectExtrapolator(tle, parameters);
300         final FieldSpacecraftState<Gradient> initialState = tlePropagator.getInitialState();
301         final Gradient unexpectedMass = initialState.getMass();
302         final Gradient expectedMass = unexpectedMass.multiply(2);
303         final FieldSpacecraftState<Gradient> newState = new FieldSpacecraftState<>(initialState.getOrbit(),
304                                                                                    initialState.getAttitude()).withMass(expectedMass);
305 
306         // WHEN
307         tlePropagator.resetInitialState(newState);
308 
309         // THEN
310         final FieldSpacecraftState<Gradient> actualState = tlePropagator.getInitialState();
311         Assertions.assertEquals(expectedMass.getReal(), tlePropagator.getMass(actualState.getDate()).getReal());
312         Assertions.assertEquals(expectedMass.getReal(), actualState.getMass().getReal());
313         Assertions.assertNotEquals(unexpectedMass.getReal(), actualState.getMass().getReal());
314     }
315 
316     private FieldTLE<Gradient> getGradientTLE() {
317         final String line1 = "1 37753U 11036A   12090.13205652 -.00000006  00000-0  00000+0 0  2272";
318         final String line2 = "2 37753  55.0032 176.5796 0004733  13.2285 346.8266  2.00565440  5153";
319         final GradientField field = GradientField.getField(1);
320         return new FieldTLE<>(field, line1, line2);
321     }
322 
323     @Test
324     void testResetIntermediateStateForward() {
325         testResetIntermediateStateTemplate(true);
326     }
327 
328     @Test
329     void testResetIntermediateStateBackward() {
330         testResetIntermediateStateTemplate(false);
331     }
332 
333     void testResetIntermediateStateTemplate(final boolean isForward) {
334         // GIVEN
335         final FieldTLE<Gradient> tle = getGradientTLE();
336         final Field<Gradient> field = tle.getDate().getField();
337         final Gradient[] parameters = tle.getParameters(field);
338         final FieldTLEPropagator<Gradient> tlePropagator = FieldTLEPropagator.selectExtrapolator(tle, parameters);
339         final double expectedMass = 2000.;
340         final FieldSpacecraftState<Gradient> propagatedState = tlePropagator.propagate(tle.getDate().shiftedBy(1));
341         final FieldSpacecraftState<Gradient> modifiedState = new FieldSpacecraftState<>(propagatedState.getOrbit())
342                                                                  .withMass(field.getZero().newInstance(expectedMass));
343 
344         // WHEN
345         tlePropagator.resetIntermediateState(modifiedState, isForward);
346 
347         // THEN
348         final double tinyTimeShift = (isForward) ? 1e-3 : -1e-3;
349         final double actualMass = tlePropagator.getMass(modifiedState.getDate().shiftedBy(tinyTimeShift)).getReal();
350         Assertions.assertEquals(expectedMass, actualMass);
351     }
352 
353     @BeforeEach
354     public void setUp() {
355         Utils.setDataRoot("regular-data");
356 
357         // the period of the GPS satellite
358         period = 717.97 * 60.0;
359     }
360 
361 }
362