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