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.gnss;
18  
19  import org.hipparchus.geometry.euclidean.threed.Vector3D;
20  import org.hipparchus.linear.RealMatrix;
21  import org.hipparchus.util.Binary64;
22  import org.hipparchus.util.Binary64Field;
23  import org.hipparchus.util.FastMath;
24  import org.hipparchus.util.MathUtils;
25  import org.hipparchus.util.Precision;
26  import org.junit.jupiter.api.Assertions;
27  import org.junit.jupiter.api.BeforeAll;
28  import org.junit.jupiter.api.Test;
29  import org.orekit.TestUtils;
30  import org.orekit.Utils;
31  import org.orekit.attitudes.AttitudeProvider;
32  import org.orekit.attitudes.FrameAlignedProvider;
33  import org.orekit.data.DataContext;
34  import org.orekit.frames.Frame;
35  import org.orekit.frames.FramesFactory;
36  import org.orekit.gnss.SEMParser;
37  import org.orekit.gnss.SatelliteSystem;
38  import org.orekit.orbits.OrbitType;
39  import org.orekit.propagation.AdditionalDataProvider;
40  import org.orekit.propagation.FieldSpacecraftState;
41  import org.orekit.propagation.MatricesHarvester;
42  import org.orekit.propagation.Propagator;
43  import org.orekit.propagation.SpacecraftState;
44  import org.orekit.propagation.analytical.gnss.data.CommonGnssData;
45  import org.orekit.propagation.analytical.gnss.data.FieldGPSAlmanac;
46  import org.orekit.propagation.analytical.gnss.data.GNSSOrbitalElements;
47  import org.orekit.propagation.analytical.gnss.data.GPSAlmanac;
48  import org.orekit.propagation.analytical.gnss.data.GPSLegacyNavigationMessage;
49  import org.orekit.propagation.analytical.tle.TLE;
50  import org.orekit.propagation.analytical.tle.TLEPropagator;
51  import org.orekit.time.AbsoluteDate;
52  import org.orekit.time.DateComponents;
53  import org.orekit.time.FieldAbsoluteDate;
54  import org.orekit.time.GNSSDate;
55  import org.orekit.time.TimeInterpolator;
56  import org.orekit.time.TimeScalesFactory;
57  import org.orekit.utils.CartesianDerivativesFilter;
58  import org.orekit.utils.Constants;
59  import org.orekit.utils.DoubleArrayDictionary;
60  import org.orekit.utils.FieldPVCoordinates;
61  import org.orekit.utils.IERSConventions;
62  import org.orekit.utils.PVCoordinates;
63  import org.orekit.utils.TimeStampedPVCoordinates;
64  import org.orekit.utils.TimeStampedPVCoordinatesHermiteInterpolator;
65  
66  import java.util.ArrayList;
67  import java.util.HashMap;
68  import java.util.List;
69  import java.util.Map;
70  
71  class GPSPropagatorTest {
72  
73      private static List<GPSAlmanac> almanacs;
74  
75      @BeforeAll
76      public static void setUpBeforeClass() {
77          Utils.setDataRoot("gnss");
78          GNSSDate.setRolloverReference(new DateComponents(DateComponents.GPS_EPOCH, 7 * 512));
79          // Get the parser to read a SEM file
80          SEMParser reader = new SEMParser(null);
81          // Reads the SEM file
82          reader.loadData();
83          // Gets the first SEM almanac
84          almanacs = reader.getAlmanacs();
85      }
86  
87      @Test
88      void testClockCorrections() {
89          final GNSSPropagator propagator = almanacs.get(0).getPropagator();
90          propagator.addAdditionalDataProvider(new ClockCorrectionsProvider(almanacs.get(0),
91                                                                            almanacs.get(0).getCycleDuration()));
92          // Propagate at the GPS date and one GPS cycle later
93          final AbsoluteDate date0 = almanacs.get(0).getDate();
94          double dtRelMin = 0;
95          double dtRelMax = 0;
96          for (double dt = 0; dt < 0.5 * Constants.JULIAN_DAY; dt += 1.0) {
97              SpacecraftState state = propagator.propagate(date0.shiftedBy(dt));
98              double[] corrections = state.getAdditionalState(ClockCorrectionsProvider.CLOCK_CORRECTIONS);
99              Assertions.assertEquals(3, corrections.length);
100             Assertions.assertEquals(1.33514404296875E-05, corrections[0], 1.0e-19);
101             dtRelMin = FastMath.min(dtRelMin, corrections[1]);
102             dtRelMax = FastMath.max(dtRelMax, corrections[1]);
103             Assertions.assertEquals(0.0, corrections[2], Precision.SAFE_MIN);
104         }
105         Assertions.assertEquals(0.0,        almanacs.get(0).getToc(), 1.0e-12);
106         Assertions.assertEquals(-1.1679e-8, dtRelMin, 1.0e-12);
107         Assertions.assertEquals(+1.1679e-8, dtRelMax, 1.0e-12);
108     }
109 
110     @Test
111     void testFieldClockCorrections() {
112         final FieldGPSAlmanac<Binary64> gpsAlmanac = almanacs.get(0).toField(Binary64Field.getInstance());
113         final FieldGnssPropagator<Binary64> propagator = gpsAlmanac.getPropagator();
114         propagator.addAdditionalDataProvider(new FieldClockCorrectionsProvider<>(gpsAlmanac,
115                                                                                   gpsAlmanac.getCycleDuration()));
116         // Propagate at the GPS date and one GPS cycle later
117         final FieldAbsoluteDate<Binary64> date0 = gpsAlmanac.getDate();
118         double dtRelMin = 0;
119         double dtRelMax = 0;
120         for (double dt = 0; dt < 0.5 * Constants.JULIAN_DAY; dt += 1.0) {
121             FieldSpacecraftState<Binary64> state = propagator.propagate(date0.shiftedBy(dt));
122             Binary64[] corrections = state.getAdditionalState(ClockCorrectionsProvider.CLOCK_CORRECTIONS);
123             Assertions.assertEquals(3, corrections.length);
124             Assertions.assertEquals(1.33514404296875E-05, corrections[0].getReal(), 1.0e-19);
125             dtRelMin = FastMath.min(dtRelMin, corrections[1].getReal());
126             dtRelMax = FastMath.max(dtRelMax, corrections[1].getReal());
127             Assertions.assertEquals(0.0, corrections[2].getReal(), Precision.SAFE_MIN);
128         }
129         Assertions.assertEquals(0.0,        gpsAlmanac.getToc().getReal(), 1.0e-12);
130         Assertions.assertEquals(-1.1679e-8, dtRelMin, 1.0e-12);
131         Assertions.assertEquals(+1.1679e-8, dtRelMax, 1.0e-12);
132     }
133 
134     @Test
135     void testGPSCycle() {
136         // Builds the GPSPropagator from the almanac
137         final GNSSPropagator propagator = almanacs.get(0).getPropagator(DataContext.getDefault().getFrames(),
138                 Utils.defaultLaw(), FramesFactory.getEME2000(), FramesFactory.getITRF(IERSConventions.IERS_2010, false), 1521.0);
139         // Propagate at the GPS date and one GPS cycle later
140         final AbsoluteDate date0 = almanacs.get(0).getDate();
141         final Vector3D p0 = propagator.propagateInEcef(date0).getPosition();
142         final double gpsCycleDuration = almanacs.get(0).getCycleDuration();
143         final AbsoluteDate date1 = date0.shiftedBy(gpsCycleDuration);
144         final Vector3D p1 = propagator.propagateInEcef(date1).getPosition();
145 
146         // Checks
147         Assertions.assertEquals(0., p0.distance(p1), 0.);
148     }
149 
150     @Test
151     void testFrames() {
152         // Builds the GPSPropagator from the almanac
153         final GNSSPropagator propagator = almanacs.get(0).getPropagator();
154         Assertions.assertEquals("EME2000", propagator.getFrame().getName());
155         Assertions.assertEquals(3.986005e14, almanacs.get(0).getMu(), 1.0e6);
156         // Defines some date
157         final AbsoluteDate date = new AbsoluteDate(2016, 3, 3, 12, 0, 0., TimeScalesFactory.getUTC());
158         // Get PVCoordinates at the date in the ECEF
159         final PVCoordinates pv0 = propagator.propagateInEcef(date);
160         // Get PVCoordinates at the date in the ECEF
161         final PVCoordinates pv1 = propagator.getPVCoordinates(date, propagator.getECEF());
162 
163         // Checks
164         Assertions.assertEquals(0., pv0.getPosition().distance(pv1.getPosition()), 3.3e-8);
165         Assertions.assertEquals(0., pv0.getVelocity().distance(pv1.getVelocity()), 3.9e-12);
166     }
167 
168     @Test
169     void testResetInitialState() {
170         final GNSSPropagator propagator = almanacs.get(0).getPropagator();
171         final SpacecraftState old = propagator.getInitialState();
172         propagator.resetInitialState(new SpacecraftState(old.getOrbit(), old.getAttitude(), old.getMass() + 1000));
173         Assertions.assertEquals(old.getMass() + 1000, propagator.getInitialState().getMass(), 1.0e-9);
174     }
175 
176     @Test
177     void testResetIntermediateState() {
178         GNSSPropagator propagator = new GNSSPropagatorBuilder(almanacs.get(0)).build();
179         final SpacecraftState old = propagator.getInitialState();
180         propagator.resetIntermediateState(new SpacecraftState(old.getOrbit(), old.getAttitude(), old.getMass() + 1000),
181                                           true);
182         Assertions.assertEquals(old.getMass() + 1000, propagator.getInitialState().getMass(), 1.0e-9);
183     }
184 
185     @Test
186     void testTLE() {
187 
188         List<GNSSPropagator> gpsPropagators = new ArrayList<>();
189         for (final GPSAlmanac almanac : almanacs) {
190             gpsPropagators.add(almanac.getPropagator());
191         }
192 
193         // the following map corresponds to the GPS constellation status in early 2016
194         final Map<Integer, TLE> prnToTLE = new HashMap<>();
195         prnToTLE.put( 1,
196                      new TLE("1 37753U 11036A   16059.51505483 -.00000016  00000-0  00000+0 0  9995",
197                              "2 37753  55.2230 119.7200 0049958  23.9363 306.3749  2.00566105 33828"));
198         prnToTLE.put( 2,
199                      new TLE("1 28474U 04045A   16059.68518942 -.00000018 +00000-0 +00000-0 0  9992",
200                              "2 28474 054.0048 117.2153 0156693 236.8152 092.4773 02.00556694082981"));
201         prnToTLE.put( 3,
202                      new TLE("1 40294U 14068A   16059.64183862 +.00000012 +00000-0 +00000-0 0  9990",
203                              "2 40294 054.9632 179.3690 0003634 223.4011 136.5596 02.00553679009768"));
204         prnToTLE.put( 4,
205                       new TLE("1 34661U 09014A   16059.51658765 -.00000072  00000-0  00000+0 0  9997",
206                               "2 34661  56.3471   0.6400 0080177  48.0430 129.2467  2.00572451 50844"));
207         prnToTLE.put( 5,
208                       new TLE("1 35752U 09043A   16059.39819238  .00000011  00000-0  00000+0 0  9993",
209                               "2 35752  54.2243 178.7652 0044753  24.8598  29.4422  2.00555538 47893"));
210         prnToTLE.put( 6,
211                       new TLE("1 39741U 14026A   16059.18044747 -.00000016  00000-0  00000+0 0  9990",
212                               "2 39741  55.2140 119.2493 0005660 259.2190 100.7882  2.00566982 13061"));
213         prnToTLE.put( 7,
214                       new TLE("1 32711U 08012A   16059.36304856 -.00000033 +00000-0 +00000-0 0  9998",
215                               "2 32711 055.4269 300.7399 0091867 207.6311 151.9340 02.00564257058321"));
216         prnToTLE.put( 8,
217                       new TLE("1 40730U 15033A   16059.44106931 -.00000026 +00000-0 +00000-0 0  9994",
218                               "2 40730 055.1388 059.0069 0020452 282.1769 077.6168 02.00566073004562"));
219         prnToTLE.put( 9,
220                       new TLE("1 40105U 14045A   16059.27451329  .00000045  00000-0  00000+0 0  9996",
221                               "2 40105  54.7529 238.9873 0004485 121.4766 238.5557  2.00568637 11512"));
222         prnToTLE.put(10,
223                      new TLE("1 41019U 15062A   16059.49433942  .00000013  00000-0  00000+0 0  9991",
224                              "2 41019  54.9785 179.1399 0012328 204.9013 155.0292  2.00561967  2382"));
225         prnToTLE.put(11,
226                      new TLE("1 25933U 99055A   16059.51073770 -.00000024  00000-0  00000+0 0  9997",
227                              "2 25933  51.3239  98.4815 0159812  86.1576 266.7718  2.00565163120122"));
228         prnToTLE.put(12,
229                      new TLE("1 29601U 06052A   16059.62966898 -.00000070 +00000-0 +00000-0 0  9994",
230                              "2 29601 056.7445 002.2755 0057667 037.0706 323.3313 02.00552237067968"));
231         prnToTLE.put(13,
232                      new TLE("1 24876U 97035A   16059.41696335  .00000046  00000-0  00000+0 0  9998",
233                              "2 24876  55.6966 245.8203 0044339 114.8899 245.5712  2.00562657136305"));
234         prnToTLE.put(14,
235                      new TLE("1 26605U 00071A   16059.56211888  .00000047  00000-0  00000+0 0  9997",
236                              "2 26605  55.2663 243.7251 0085518 248.7231  95.5323  2.00557009112094"));
237         prnToTLE.put(15,
238                      new TLE("1 32260U 07047A   16059.45678257 +.00000044 +00000-0 +00000-0 0  9994",
239                              "2 32260 053.3641 236.0940 0079746 026.4105 333.9774 02.00547771061402"));
240         prnToTLE.put(16,
241                      new TLE("1 27663U 03005A   16059.14440417 -.00000071  00000-0  00000+0 0  9996",
242                              "2 27663  56.7743   3.3691 0085346  17.6322 214.4333  2.00559487 95843"));
243         prnToTLE.put(17,
244                      new TLE("1 28874U 05038A   16059.21070933 -.00000024  00000-0  00000+0 0  9997",
245                              "2 28874  55.8916  61.9596 0112077 248.9647 205.7384  2.00567116 76379"));
246         prnToTLE.put(18,
247                      new TLE("1 26690U 01004A   16059.51332910  .00000008  00000-0  00000+0 0  9990",
248                              "2 26690  52.9999 177.6630 0169501 250.8579 153.6293  2.00563995110499"));
249         prnToTLE.put(19,
250                      new TLE("1 28190U 04009A   16058.12363503 -.00000030  00000-0  00000+0 0  9999",
251                              "2 28190  55.7230  64.8110 0105865  40.0254 321.4519  2.00572212 87500"));
252         prnToTLE.put(20,
253                      new TLE("1 26360U 00025A   16059.44770263  .00000005  00000-0  00000+0 0  9992",
254                              "2 26360  53.0712 174.6895 0046205  76.0615 334.4302  2.00559931115818"));
255         prnToTLE.put(21,
256                      new TLE("1 27704U 03010A   16059.50719524 -.00000019  00000-0  00000+0 0  9998",
257                              "2 27704  53.6134 117.9454 0234081 255.6874 199.2128  2.00564673 94659"));
258         prnToTLE.put(22,
259                      new TLE("1 28129U 03058A   16059.06680941  .00000008  00000-0  00000+0 0  9990",
260                              "2 28129  52.8771 177.7253 0079127 245.1376 114.0279  2.00398763 89357"));
261         prnToTLE.put(23,
262                      new TLE("1 28361U 04023A   16059.54310021  .00000046  00000-0  00000+0 0  9995",
263                              "2 28361  54.2347 239.3240 0106509 211.5355  11.7648  2.00557932 85613"));
264         prnToTLE.put(24,
265                      new TLE("1 38833U 12053A   16059.04618549 -.00000032  00000-0  00000+0 0  9999",
266                              "2 38833  54.4591 298.1383 0042253  18.7074 341.5041  2.00568407 24895"));
267         prnToTLE.put(25,
268                      new TLE("1 36585U 10022A   16059.29300735 -.00000074  00000-0  00000+0 0  9993",
269                              "2 36585  56.0738 359.4320 0050768  38.3425  49.1794  2.00578535 42134"));
270         prnToTLE.put(26,
271                      new TLE("1 40534U 15013A   16059.28299301 -.00000076  00000-0  00000+0 0  9994",
272                              "2 40534  55.0430 359.0082 0009349 342.4081  17.5685  2.00558853  6801"));
273         prnToTLE.put(27,
274                      new TLE("1 39166U 13023A   16059.40401153 -.00000025  00000-0  00000+0 0  9990",
275                              "2 39166  55.6020  59.1224 0032420   7.7969 352.2759  2.00568484 20414"));
276         prnToTLE.put(28,
277                      new TLE("1 26407U 00040A   16059.80383354 -.00000069 +00000-0 +00000-0 0  9994",
278                              "2 26407 056.6988 003.6328 0201499 267.0948 317.6209 02.00569902114508"));
279         prnToTLE.put(29,
280                      new TLE("1 32384U 07062A   16059.44770263 -.00000021  00000-0  00000+0 0  9992",
281                              "2 32384  55.9456  62.5022 0011922 319.9531 172.6730  2.00571577 60128"));
282         prnToTLE.put(30,
283                      new TLE("1 39533U 14008A   16059.40267873 -.00000038 +00000-0 +00000-0 0  9996",
284                              "2 39533 054.6126 303.3404 0017140 179.4267 180.6311 02.00568364014251"));
285         prnToTLE.put(31,
286                      new TLE("1 29486U 06042A   16059.50651990 -.00000032  00000-0  00000+0 0  9992",
287                              "2 29486  55.7041 301.2472 0084115 334.2804 254.9897  2.00560606 69098"));
288         prnToTLE.put(32,
289                      new TLE("1 41328U 16007A   16059.56873502  .00000049  00000-0  00000+0 0  9991",
290                              "2 41328  55.0137 239.0304 0002157 298.9074  61.0768  1.99172830   453"));
291 
292         for (final GNSSPropagator gpsPropagator : gpsPropagators) {
293             final int prn = gpsPropagator.getOrbitalElements().getPRN();
294             TLE tle = prnToTLE.get(prn);
295             TLEPropagator tlePropagator = TLEPropagator.selectExtrapolator(tle);
296             for (double dt = 0; dt < Constants.JULIAN_DAY; dt += 600) {
297                 final AbsoluteDate date = tlePropagator.getInitialState().getDate().shiftedBy(dt);
298                 final PVCoordinates gpsPV = gpsPropagator.getPVCoordinates(date, gpsPropagator.getECI());
299                 final PVCoordinates tlePV = tlePropagator.getPVCoordinates(date, gpsPropagator.getECI());
300                 Assertions.assertEquals(0.0,
301                                     Vector3D.distance(gpsPV.getPosition(), tlePV.getPosition()),
302                                     8400.0);
303             }
304         }
305     }
306 
307     @Test
308     void testDerivativesConsistency() {
309 
310         final Frame eme2000 = FramesFactory.getEME2000();
311         double errorP = 0;
312         double errorV = 0;
313         double errorA = 0;
314         for (final GPSAlmanac almanac : almanacs) {
315             final GNSSPropagator propagator = almanac.getPropagator();
316             GNSSOrbitalElements<?> elements = propagator.getOrbitalElements();
317             AbsoluteDate t0 = new GNSSDate(elements.getWeek(), elements.getTime(), SatelliteSystem.GPS).getDate();
318             for (double dt = 0; dt < Constants.JULIAN_DAY; dt += 600) {
319                 final AbsoluteDate central = t0.shiftedBy(dt);
320                 final PVCoordinates pv = propagator.getPVCoordinates(central, eme2000);
321                 final double h = 60.0;
322                 List<TimeStampedPVCoordinates> sample = new ArrayList<>();
323                 for (int i = -3; i <= 3; ++i) {
324                     sample.add(propagator.getPVCoordinates(central.shiftedBy(i * h), eme2000));
325                 }
326 
327                 // create interpolator
328                 final TimeInterpolator<TimeStampedPVCoordinates> interpolator =
329                         new TimeStampedPVCoordinatesHermiteInterpolator(sample.size(), CartesianDerivativesFilter.USE_P);
330 
331                 final PVCoordinates interpolated = interpolator.interpolate(central, sample);
332                 errorP = FastMath.max(errorP, Vector3D.distance(pv.getPosition(), interpolated.getPosition()));
333                 errorV = FastMath.max(errorV, Vector3D.distance(pv.getVelocity(), interpolated.getVelocity()));
334                 errorA = FastMath.max(errorA, Vector3D.distance(pv.getAcceleration(), interpolated.getAcceleration()));
335             }
336         }
337         Assertions.assertEquals(0.0, errorP, 3.8e-9);
338         Assertions.assertEquals(0.0, errorV, 3.5e-8);
339         Assertions.assertEquals(0.0, errorA, 1.1e-8);
340 
341     }
342 
343     @Test
344     void testPosition() {
345         // Initial GPS orbital elements (Ref: IGS)
346         final GPSLegacyNavigationMessage goe = new GPSLegacyNavigationMessage(DataContext.getDefault().getTimeScales(),
347                                                                               SatelliteSystem.GPS);
348         goe.setPRN(7);
349         goe.setWeek(0);
350         goe.setTime(288000);
351         goe.setSqrtA(5153.599830627441);
352         goe.setE(0.012442796607501805);
353         goe.setDeltaN0(4.419469802942352E-9);
354         goe.setI0(0.9558937988021613);
355         goe.setIDot(-2.4608167886110235E-10);
356         goe.setOmega0(1.0479401362158658);
357         goe.setOmegaDot(-7.967117576712062E-9);
358         goe.setPa(-2.4719019944000538);
359         goe.setM0(-1.0899023379614294);
360         goe.setCuc(4.3995678424835205E-6);
361         goe.setCus(1.002475619316101E-5);
362         goe.setCrc(183.40625);
363         goe.setCrs(87.03125);
364         goe.setCic(3.203749656677246E-7);
365         goe.setCis(4.0978193283081055E-8);
366 
367         // Date of the GPS orbital elements
368         final AbsoluteDate target = goe.getDate();
369         // Build the GPS propagator
370         final GNSSPropagator propagator = goe.getPropagator();
371         // Compute the PV coordinates at the date of the GPS orbital elements
372         final PVCoordinates pv = propagator.getPVCoordinates(target, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
373         // Computed position
374         final Vector3D computedPos = pv.getPosition();
375         // Expected position (reference from IGS file igu20484_00.sp3)
376         final Vector3D expectedPos = new Vector3D(-4920705.292, 24248099.200, 9236130.101);
377 
378         Assertions.assertEquals(0., Vector3D.distance(expectedPos, computedPos), 3.2);
379     }
380 
381     @Test
382     void testStmAndJacobian() {
383         // Initial GPS orbital elements (Ref: IGS)
384         final GPSLegacyNavigationMessage goe = new GPSLegacyNavigationMessage(DataContext.getDefault().getTimeScales(),
385                                                                               SatelliteSystem.GPS);
386         goe.setPRN(7);
387         goe.setWeek(0);
388         goe.setTime(288000);
389         goe.setSqrtA(5153.599830627441);
390         goe.setE(0.012442796607501805);
391         goe.setDeltaN0(4.419469802942352E-9);
392         goe.setI0(0.9558937988021613);
393         goe.setIDot(-2.4608167886110235E-10);
394         goe.setOmega0(1.0479401362158658);
395         goe.setOmegaDot(-7.967117576712062E-9);
396         goe.setPa(-2.4719019944000538);
397         goe.setM0(-1.0899023379614294);
398         goe.setCuc(4.3995678424835205E-6);
399         goe.setCus(1.002475619316101E-5);
400         goe.setCrc(183.40625);
401         goe.setCrs(87.03125);
402         goe.setCic(3.203749656677246E-7);
403         goe.setCis(4.0978193283081055E-8);
404         GNSSPropagator propagator = goe.getPropagator();
405 
406         // we want to compute the partial derivatives with respect to Crs and Crc parameters
407         Assertions.assertEquals(9, propagator.getOrbitalElements().getParameters().length);
408         propagator.getOrbitalElements().getParameterDriver(CommonGnssData.RADIUS_SINE).setSelected(true);
409         propagator.getOrbitalElements().getParameterDriver(CommonGnssData.RADIUS_COSINE).setSelected(true);
410         final DoubleArrayDictionary initialJacobianColumns = new DoubleArrayDictionary();
411         initialJacobianColumns.put(CommonGnssData.RADIUS_SINE,   new double[6]);
412         initialJacobianColumns.put(CommonGnssData.RADIUS_COSINE, new double[6]);
413         final MatricesHarvester harvester = propagator.setupMatricesComputation("stm", null, initialJacobianColumns);
414 
415         // harvester sorts the columns lexicographically, and wraps them as SpanXxx##
416         Assertions.assertEquals(2, harvester.getJacobiansColumnsNames().size());
417         Assertions.assertEquals("Span" + CommonGnssData.RADIUS_COSINE + "0", harvester.getJacobiansColumnsNames().get(0));
418         Assertions.assertEquals("Span" + CommonGnssData.RADIUS_SINE   + "0", harvester.getJacobiansColumnsNames().get(1));
419 
420         // propagate orbit
421         final SpacecraftState state = propagator.propagate(goe.getDate().shiftedBy(3600.0));
422 
423         // extract state transition matrix
424         final RealMatrix stm = harvester.getStateTransitionMatrix(state);
425         Assertions.assertEquals(OrbitType.CARTESIAN, harvester.getOrbitType());
426         Assertions.assertEquals(6, stm.getRowDimension());
427         Assertions.assertEquals(6, stm.getColumnDimension());
428         Assertions.assertEquals(1.202731937, stm.getEntry(0, 0), 1.0e-9);
429 
430         // extract Jacobian matrix
431         final RealMatrix jacobian = harvester.getParametersJacobian(state);
432         Assertions.assertEquals(6, jacobian.getRowDimension());
433         Assertions.assertEquals(2, jacobian.getColumnDimension());
434         Assertions.assertEquals(0.959547021, jacobian.getEntry(0, 0), 1.0e-9);
435 
436     }
437 
438     @Test
439     void testRebuildModel() {
440 
441         final Frame            eci              = FramesFactory.getEME2000();
442         final Frame            ecef             = FramesFactory.getITRF(IERSConventions.IERS_2010, false);
443         final double           mass             = Propagator.DEFAULT_MASS;
444         final AttitudeProvider attitudeProvider = FrameAlignedProvider.of(eci);
445 
446         // Initial GPS orbital elements (Ref: IGS)
447         final GPSLegacyNavigationMessage goe = new GPSLegacyNavigationMessage(DataContext.getDefault().getTimeScales(),
448                                                                               SatelliteSystem.GPS);
449         goe.setPRN(7);
450         goe.setWeek(0);
451         goe.setTime(288000);
452         goe.setSqrtA(5153.599830627441);
453         goe.setE(0.012442796607501805);
454         goe.setDeltaN0(4.419469802942352E-9);
455         goe.setI0(0.9558937988021613);
456         goe.setIDot(-2.4608167886110235E-10);
457         goe.setOmega0(1.0479401362158658);
458         goe.setOmegaDot(-7.967117576712062E-9);
459         goe.setPa(-2.4719019944000538);
460         goe.setM0(-1.0899023379614294);
461         goe.setCuc(4.3995678424835205E-6);
462         goe.setCus(1.002475619316101E-5);
463         goe.setCrc(183.40625);
464         goe.setCrs(87.03125);
465         goe.setCic(3.203749656677246E-7);
466         goe.setCis(4.0978193283081055E-8);
467         GNSSPropagator propagator = goe.getPropagator(DataContext.getDefault().getFrames(),
468                                                       attitudeProvider, eci, ecef, mass);
469 
470         final GNSSPropagator rebuilt = new GNSSPropagator(propagator.getInitialState(), goe,
471                                                           ecef, attitudeProvider, mass);
472         final GNSSOrbitalElements<?> oe2 = rebuilt.getOrbitalElements();
473         Assertions.assertEquals(0, goe.getDate().durationFrom(oe2),               1.0e-20);
474         Assertions.assertEquals(0,
475                                 Vector3D.distance(propagator.getInitialState().getPVCoordinates().getPosition(),
476                                                   rebuilt.getInitialState().getPVCoordinates().getPosition()),
477                                 3.8e-7);
478         Assertions.assertEquals(0,
479                                 Vector3D.distance(propagator.getInitialState().getPVCoordinates().getVelocity(),
480                                                   rebuilt.getInitialState().getPVCoordinates().getVelocity()),
481                                 4.0e-11);
482 
483         // general parameters
484         Assertions.assertEquals(goe.getMu(),            oe2.getMu(),            1.0e-20);
485         Assertions.assertEquals(goe.getCycleDuration(), oe2.getCycleDuration(), 1.0e-20);
486         Assertions.assertEquals(goe.getSystem(),        oe2.getSystem());
487         Assertions.assertEquals(goe.getPRN(),           oe2.getPRN());
488         Assertions.assertEquals(goe.getWeek(),          oe2.getWeek());
489 
490         // non-Keplerian parameters, which are just copied
491         Assertions.assertEquals(goe.getTime(),     oe2.getTime(),     1.0e-20);
492         Assertions.assertEquals(goe.getIDot(),     oe2.getIDot(),     1.0e-20);
493         Assertions.assertEquals(goe.getOmegaDot(), oe2.getOmegaDot(), 1.0e-20);
494         Assertions.assertEquals(goe.getCuc(),      oe2.getCuc(),      1.0e-20);
495         Assertions.assertEquals(goe.getCus(),      oe2.getCus(),      1.0e-20);
496         Assertions.assertEquals(goe.getCrc(),      oe2.getCrc(),      1.0e-20);
497         Assertions.assertEquals(goe.getCrs(),      oe2.getCrs(),      1.0e-20);
498         Assertions.assertEquals(goe.getCic(),      oe2.getCic(),      1.0e-20);
499         Assertions.assertEquals(goe.getCis(),      oe2.getCis(),      1.0e-20);
500 
501         // orbital parameters, those are rebuilt from the initial state
502         Assertions.assertEquals(goe.getSma(),    oe2.getSma(),                                               4.0e-8);
503         Assertions.assertEquals(goe.getE(),      oe2.getE(),                                                 1.e-15);
504         Assertions.assertEquals(goe.getI0(),     oe2.getI0(),                                                1.0e-20);
505         Assertions.assertEquals(goe.getPa(),     MathUtils.normalizeAngle(oe2.getPa(),     goe.getPa()),     1.e-13);
506         Assertions.assertEquals(goe.getOmega0(), MathUtils.normalizeAngle(oe2.getOmega0(), goe.getOmega0()), 1.7e-14);
507         Assertions.assertEquals(goe.getM0(),     MathUtils.normalizeAngle(oe2.getM0(),     goe.getM0()),     1.e-13);
508 
509     }
510 
511     @Test
512     void testIssue544() {
513         // Builds the GPSPropagator from the almanac
514         final GNSSPropagator propagator = new GNSSPropagatorBuilder(almanacs.get(0)).build();
515         // In order to test the issue, we voluntarily set a Double.NaN value in the date.
516         final AbsoluteDate date0 = new AbsoluteDate(2010, 5, 7, 7, 50, Double.NaN, TimeScalesFactory.getUTC());
517         final PVCoordinates pv0 = propagator.propagateInEcef(date0);
518         // Verify that an infinite loop did not occur
519         Assertions.assertEquals(Vector3D.NaN, pv0.getPosition());
520         Assertions.assertEquals(Vector3D.NaN, pv0.getVelocity());
521 
522     }
523 
524     @Test
525     void testFieldIssue544() {
526         // Builds the GPSPropagator from the almanac
527         final FieldGnssPropagator<Binary64> propagator =
528             new FieldGnssPropagatorBuilder<>(almanacs.get(0).toField(Binary64Field.getInstance())).build();
529         // In order to test the issue, we voluntarily set a Double.NaN value in the date.
530         final FieldAbsoluteDate<Binary64> date0 =
531             new FieldAbsoluteDate<>(Binary64Field.getInstance(),
532                                     2010, 5, 7, 7, 50, Double.NaN, TimeScalesFactory.getUTC());
533         final FieldPVCoordinates<Binary64> pv0 =
534             propagator.propagateInEcef(date0, propagator.getParameters(Binary64Field.getInstance()));
535         // Verify that an infinite loop did not occur
536         Assertions.assertTrue(pv0.getPosition().isNaN());
537         Assertions.assertTrue(pv0.getVelocity().isNaN());
538 
539     }
540 
541     /** Error with specific propagators & additional state provider throwing a NullPointerException when propagating */
542     @Test
543     void testIssue949() {
544         // GIVEN
545         // Setup propagator
546         final GNSSPropagator propagator = new GNSSPropagatorBuilder(almanacs.get(0)).build();
547 
548         // Setup additional data provider which use the initial state in its init method
549         final AdditionalDataProvider<double[]> additionalDataProvider = TestUtils.getAdditionalProviderWithInit();
550         propagator.addAdditionalDataProvider(additionalDataProvider);
551 
552         // WHEN & THEN
553         Assertions.assertDoesNotThrow(() -> propagator.propagate(new AbsoluteDate()), "No error should have been thrown");
554 
555     }
556 
557     @Test
558     void testConversion() {
559         GnssTestUtils.checkFieldConversion(almanacs.get(0));
560     }
561 
562 }