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.estimation.measurements.gnss;
18  
19  import java.util.ArrayList;
20  import java.util.Comparator;
21  import java.util.List;
22  import java.util.Locale;
23  
24  import org.hipparchus.geometry.euclidean.threed.Vector3D;
25  import org.hipparchus.stat.descriptive.moment.Mean;
26  import org.hipparchus.stat.descriptive.rank.Max;
27  import org.hipparchus.stat.descriptive.rank.Median;
28  import org.hipparchus.stat.descriptive.rank.Min;
29  import org.hipparchus.util.FastMath;
30  import org.junit.Assert;
31  import org.junit.Test;
32  import org.orekit.Utils;
33  import org.orekit.estimation.Context;
34  import org.orekit.estimation.EstimationTestUtils;
35  import org.orekit.estimation.measurements.EstimatedMeasurement;
36  import org.orekit.estimation.measurements.ObservableSatellite;
37  import org.orekit.estimation.measurements.ObservedMeasurement;
38  import org.orekit.gnss.Frequency;
39  import org.orekit.orbits.CartesianOrbit;
40  import org.orekit.orbits.Orbit;
41  import org.orekit.orbits.OrbitType;
42  import org.orekit.orbits.PositionAngle;
43  import org.orekit.propagation.BoundedPropagator;
44  import org.orekit.propagation.EphemerisGenerator;
45  import org.orekit.propagation.Propagator;
46  import org.orekit.propagation.SpacecraftState;
47  import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
48  import org.orekit.time.AbsoluteDate;
49  import org.orekit.utils.Constants;
50  import org.orekit.utils.Differentiation;
51  import org.orekit.utils.ParameterDriver;
52  import org.orekit.utils.ParameterFunction;
53  import org.orekit.utils.StateFunction;
54  import org.orekit.utils.TimeStampedPVCoordinates;
55  
56  public class OneWayGNSSPhaseTest {
57  
58      private static final Frequency FREQUENCY = Frequency.G01;
59  
60      /**
61       * Test the values of the phase comparing the observed values and the estimated values
62       * Both are calculated with a different algorithm
63       */
64      @Test
65      public void testValues() {
66          boolean printResults = false;
67          if (printResults) {
68              System.out.println("\nTest One-way GNSS phase Values\n");
69          }
70          // Run test
71          this.genericTestValues(printResults);
72      }
73  
74      /**
75       * Test the values of the state derivatives using a numerical
76       * finite differences calculation as a reference
77       */
78      @Test
79      public void testStateDerivatives() {
80  
81          boolean printResults = false;
82          if (printResults) {
83              System.out.println("\nTest One-way GNSS phase State Derivatives - Finite Differences Comparison\n");
84          }
85          // Run test
86          double refErrorsPMedian = 1.9e-10;
87          double refErrorsPMean   = 5.3e-10;
88          double refErrorsPMax    = 3.7e-08;
89          double refErrorsVMedian = 4.7e-04;
90          double refErrorsVMean   = 1.6e-03;
91          double refErrorsVMax    = 1.1e-01;
92          this.genericTestStateDerivatives(printResults, 0,
93                                           refErrorsPMedian, refErrorsPMean, refErrorsPMax,
94                                           refErrorsVMedian, refErrorsVMean, refErrorsVMax);
95      }
96  
97      /**
98       * Test the values of the parameters' derivatives using a numerical
99       * finite differences calculation as a reference
100      */
101     @Test
102     public void testParameterDerivatives() {
103 
104         // Print the results ?
105         boolean printResults = false;
106 
107         if (printResults) {
108             System.out.println("\nTest One-way GNSS phase Derivatives - Finite Differences Comparison\n");
109         }
110         // Run test
111         double refErrorsMedian = 1.0e-15;
112         double refErrorsMean   = 1.0e-15;
113         double refErrorsMax    = 1.0e-15;
114         this.genericTestParameterDerivatives(printResults,
115                                              refErrorsMedian, refErrorsMean, refErrorsMax);
116 
117     }
118 
119     /**
120      * Generic test function for values of the one-way GNSS phase
121      * @param printResults Print the results ?
122      */
123     void genericTestValues(final boolean printResults) {
124 
125         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
126 
127         final NumericalPropagatorBuilder propagatorBuilder =
128                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true,
129                                               1.0e-6, 60.0, 0.001);
130 
131         // Create perfect inter-satellites range measurements
132         final TimeStampedPVCoordinates original = context.initialOrbit.getPVCoordinates();
133         final Orbit closeOrbit = new CartesianOrbit(new TimeStampedPVCoordinates(context.initialOrbit.getDate(),
134                                                                                  original.getPosition().add(new Vector3D(1000, 2000, 3000)),
135                                                                                  original.getVelocity().add(new Vector3D(-0.03, 0.01, 0.02))),
136                                                     context.initialOrbit.getFrame(),
137                                                     context.initialOrbit.getMu());
138         final Propagator closePropagator = EstimationTestUtils.createPropagator(closeOrbit,
139                                                                                 propagatorBuilder);
140         final EphemerisGenerator generator = closePropagator.getEphemerisGenerator();
141         closePropagator.propagate(context.initialOrbit.getDate().shiftedBy(3.5 * closeOrbit.getKeplerianPeriod()));
142         final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
143         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
144                                                                            propagatorBuilder);
145 
146         final int    ambiguity         = 1234;
147         final double localClockOffset  = 0.137e-6;
148         final double remoteClockOffset = 469.0e-6;
149         final List<ObservedMeasurement<?>> measurements =
150                         EstimationTestUtils.createMeasurements(propagator,
151                                                                new OneWayGNSSPhaseCreator(ephemeris, FREQUENCY, ambiguity, localClockOffset, remoteClockOffset),
152                                                                1.0, 3.0, 300.0);
153 
154         // Lists for results' storage - Used only for derivatives with respect to state
155         // "final" value to be seen by "handleStep" function of the propagator
156         final List<Double> absoluteErrors = new ArrayList<Double>();
157         final List<Double> relativeErrors = new ArrayList<Double>();
158 
159         // Use a lambda function to implement "handleStep" function
160         propagator.setStepHandler(interpolator -> {
161 
162             for (final ObservedMeasurement<?> measurement : measurements) {
163 
164                 //  Play test if the measurement date is between interpolator previous and current date
165                 if ((measurement.getDate().durationFrom(interpolator.getPreviousState().getDate()) > 0.) &&
166                     (measurement.getDate().durationFrom(interpolator.getCurrentState().getDate())  <=  0.)
167                    ) {
168                     // We intentionally propagate to a date which is close to the
169                     // real spacecraft state but is *not* the accurate date, by
170                     // compensating only part of the downlink delay. This is done
171                     // in order to validate the partial derivatives with respect
172                     // to velocity.
173                     final double          meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
174                     final AbsoluteDate    date      = measurement.getDate().shiftedBy(-0.75 * meanDelay);
175                     final SpacecraftState state     = interpolator.getInterpolatedState(date);
176 
177                     // Values of the phase & errors
178                     final double phaseObserved  = measurement.getObservedValue()[0];
179                     final EstimatedMeasurement<?> estimated = measurement.estimate(0, 0,
180                                                                                    new SpacecraftState[] {
181                                                                                        state,
182                                                                                        ephemeris.propagate(state.getDate())
183                                                                                    });
184                     Assert.assertEquals(FREQUENCY.getWavelength(), ((OneWayGNSSPhase) measurement).getWavelength(), 1.0e-15);
185                     final double phaseEstimated = estimated.getEstimatedValue()[0];
186                     final double absoluteError  = phaseEstimated-phaseObserved;
187                     absoluteErrors.add(absoluteError);
188                     relativeErrors.add(FastMath.abs(absoluteError)/FastMath.abs(phaseObserved));
189 
190                     // Print results on console ?
191                     if (printResults) {
192                         final AbsoluteDate measurementDate = measurement.getDate();
193 
194                         System.out.format(Locale.US, "%-23s  %-23s  %19.6f  %19.6f  %13.6e  %13.6e%n",
195                                          measurementDate, date,
196                                          phaseObserved, phaseEstimated,
197                                          FastMath.abs(phaseEstimated-phaseObserved),
198                                          FastMath.abs((phaseEstimated-phaseObserved)/phaseObserved));
199                     }
200 
201                 } // End if measurement date between previous and current interpolator step
202             } // End for loop on the measurements
203         }); // End lambda function handlestep
204 
205         // Print results on console ? Header
206         if (printResults) {
207             System.out.format(Locale.US, "%-23s  %-23s  %19s  %19s  %19s  %19s%n",
208                               "Measurement Date", "State Date",
209                               "Phase observed [cycles]", "Phase estimated [cycles]",
210                               "ΔPhase [cycles]", "rel ΔPhase");
211         }
212 
213         // Rewind the propagator to initial date
214         propagator.propagate(context.initialOrbit.getDate());
215 
216         // Sort measurements chronologically
217         measurements.sort(Comparator.naturalOrder());
218 
219         // Propagate to final measurement's date
220         propagator.propagate(measurements.get(measurements.size()-1).getDate());
221 
222         // Convert lists to double array
223         final double[] absErrors = absoluteErrors.stream().mapToDouble(Double::doubleValue).toArray();
224         final double[] relErrors = relativeErrors.stream().mapToDouble(Double::doubleValue).toArray();
225 
226         // Statistics' assertion
227         final double absErrorsMedian = new Median().evaluate(absErrors);
228         final double absErrorsMin    = new Min().evaluate(absErrors);
229         final double absErrorsMax    = new Max().evaluate(absErrors);
230         final double relErrorsMedian = new Median().evaluate(relErrors);
231         final double relErrorsMax    = new Max().evaluate(relErrors);
232 
233         // Print the results on console ? Final results
234         if (printResults) {
235             System.out.println();
236             System.out.println("Absolute errors median: " +  absErrorsMedian);
237             System.out.println("Absolute errors min   : " +  absErrorsMin);
238             System.out.println("Absolute errors max   : " +  absErrorsMax);
239             System.out.println("Relative errors median: " +  relErrorsMedian);
240             System.out.println("Relative errors max   : " +  relErrorsMax);
241         }
242 
243         Assert.assertEquals(0.0, absErrorsMedian, 6.5e-7);
244         Assert.assertEquals(0.0, absErrorsMin,    3.1e-6);
245         Assert.assertEquals(0.0, absErrorsMax,    9.0e-7);
246         Assert.assertEquals(0.0, relErrorsMedian, 5.9e-12);
247         Assert.assertEquals(0.0, relErrorsMax,    1.4e-10);
248     }
249 
250     void genericTestStateDerivatives(final boolean printResults, final int index,
251                                      final double refErrorsPMedian, final double refErrorsPMean, final double refErrorsPMax,
252                                      final double refErrorsVMedian, final double refErrorsVMean, final double refErrorsVMax) {
253 
254         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
255 
256         final NumericalPropagatorBuilder propagatorBuilder =
257                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true,
258                                               1.0e-6, 60.0, 0.001);
259 
260         // Create perfect one-way GNSS range measurements
261         final TimeStampedPVCoordinates original = context.initialOrbit.getPVCoordinates();
262         final Orbit closeOrbit = new CartesianOrbit(new TimeStampedPVCoordinates(context.initialOrbit.getDate(),
263                                                                                  original.getPosition().add(new Vector3D(1000, 2000, 3000)),
264                                                                                  original.getVelocity().add(new Vector3D(-0.03, 0.01, 0.02))),
265                                                     context.initialOrbit.getFrame(),
266                                                     context.initialOrbit.getMu());
267         final Propagator closePropagator = EstimationTestUtils.createPropagator(closeOrbit,
268                                                                                 propagatorBuilder);
269         final EphemerisGenerator generator = closePropagator.getEphemerisGenerator();
270         closePropagator.propagate(context.initialOrbit.getDate().shiftedBy(3.5 * closeOrbit.getKeplerianPeriod()));
271         final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
272         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
273                                                                            propagatorBuilder);
274 
275         final int    ambiguity         = 1234;
276         final double localClockOffset  = 0.137e-6;
277         final double remoteClockOffset = 469.0e-6;
278         final List<ObservedMeasurement<?>> measurements =
279                         EstimationTestUtils.createMeasurements(propagator,
280                                                                new OneWayGNSSPhaseCreator(ephemeris, FREQUENCY, ambiguity, localClockOffset, remoteClockOffset),
281                                                                1.0, 3.0, 300.0);
282 
283         // Lists for results' storage - Used only for derivatives with respect to state
284         // "final" value to be seen by "handleStep" function of the propagator
285         final List<Double> errorsP = new ArrayList<Double>();
286         final List<Double> errorsV = new ArrayList<Double>();
287 
288         // Use a lambda function to implement "handleStep" function
289         propagator.setStepHandler(interpolator -> {
290 
291             for (final ObservedMeasurement<?> measurement : measurements) {
292 
293                 //  Play test if the measurement date is between interpolator previous and current date
294                 if ((measurement.getDate().durationFrom(interpolator.getPreviousState().getDate()) > 0.) &&
295                     (measurement.getDate().durationFrom(interpolator.getCurrentState().getDate())  <=  0.)
296                    ) {
297 
298                     // We intentionally propagate to a date which is close to the
299                     // real spacecraft state but is *not* the accurate date, by
300                     // compensating only part of the downlink delay. This is done
301                     // in order to validate the partial derivatives with respect
302                     // to velocity.
303                     final double            meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
304                     final AbsoluteDate      date      = measurement.getDate().shiftedBy(-0.75 * meanDelay);
305                     final SpacecraftState[] states    = {
306                         interpolator.getInterpolatedState(date),
307                         ephemeris.propagate(date)
308                     };
309                     final double[][]      jacobian  = measurement.estimate(0, 0, states).getStateDerivatives(index);
310 
311                     // Jacobian reference value
312                     final double[][] jacobianRef;
313 
314                     // Compute a reference value using finite differences
315                     jacobianRef = Differentiation.differentiate(new StateFunction() {
316                         public double[] value(final SpacecraftState state) {
317                             final SpacecraftState[] s = states.clone();
318                             s[index] = state;
319                             return measurement.estimate(0, 0, s).getEstimatedValue();
320                         }
321                     }, measurement.getDimension(), propagator.getAttitudeProvider(),
322                        OrbitType.CARTESIAN, PositionAngle.TRUE, 2.0, 3).value(states[index]);
323 
324                     Assert.assertEquals(jacobianRef.length, jacobian.length);
325                     Assert.assertEquals(jacobianRef[0].length, jacobian[0].length);
326 
327                     // Errors & relative errors on the Jacobian
328                     double [][] dJacobian         = new double[jacobian.length][jacobian[0].length];
329                     double [][] dJacobianRelative = new double[jacobian.length][jacobian[0].length];
330                     for (int i = 0; i < jacobian.length; ++i) {
331                         for (int j = 0; j < jacobian[i].length; ++j) {
332                             dJacobian[i][j] = jacobian[i][j] - jacobianRef[i][j];
333                             dJacobianRelative[i][j] = FastMath.abs(dJacobian[i][j]/jacobianRef[i][j]);
334 
335                             if (j < 3) { errorsP.add(dJacobianRelative[i][j]);
336                             } else { errorsV.add(dJacobianRelative[i][j]); }
337                         }
338                     }
339                     // Print values in console ?
340                     if (printResults) {
341                         System.out.format(Locale.US, "%-23s  %-23s  " +
342                                         "%10.3e  %10.3e  %10.3e  " +
343                                         "%10.3e  %10.3e  %10.3e  " +
344                                         "%10.3e  %10.3e  %10.3e  " +
345                                         "%10.3e  %10.3e  %10.3e%n",
346                                         measurement.getDate(), date,
347                                         dJacobian[0][0], dJacobian[0][1], dJacobian[0][2],
348                                         dJacobian[0][3], dJacobian[0][4], dJacobian[0][5],
349                                         dJacobianRelative[0][0], dJacobianRelative[0][1], dJacobianRelative[0][2],
350                                         dJacobianRelative[0][3], dJacobianRelative[0][4], dJacobianRelative[0][5]);
351                     }
352                 } // End if measurement date between previous and current interpolator step
353             } // End for loop on the measurements
354         });
355 
356         // Print results on console ?
357         if (printResults) {
358             System.out.format(Locale.US, "%-23s  %-23s  " +
359                             "%10s  %10s  %10s  " +
360                             "%10s  %10s  %10s  " +
361                             "%10s  %10s  %10s  " +
362                             "%10s  %10s  %10s%n",
363                             "Measurement Date", "State Date",
364                             "ΔdPx", "ΔdPy", "ΔdPz", "ΔdVx", "ΔdVy", "ΔdVz",
365                             "rel ΔdPx", "rel ΔdPy", "rel ΔdPz",
366                             "rel ΔdVx", "rel ΔdVy", "rel ΔdVz");
367         }
368 
369         // Rewind the propagator to initial date
370         propagator.propagate(context.initialOrbit.getDate());
371 
372         // Sort measurements, primarily chronologically
373         measurements.sort(Comparator.naturalOrder());
374 
375         // Propagate to final measurement's date
376         propagator.propagate(measurements.get(measurements.size()-1).getDate());
377 
378         // Convert lists to double[] and evaluate some statistics
379         final double relErrorsP[] = errorsP.stream().mapToDouble(Double::doubleValue).toArray();
380         final double relErrorsV[] = errorsV.stream().mapToDouble(Double::doubleValue).toArray();
381 
382         final double errorsPMedian = new Median().evaluate(relErrorsP);
383         final double errorsPMean   = new Mean().evaluate(relErrorsP);
384         final double errorsPMax    = new Max().evaluate(relErrorsP);
385         final double errorsVMedian = new Median().evaluate(relErrorsV);
386         final double errorsVMean   = new Mean().evaluate(relErrorsV);
387         final double errorsVMax    = new Max().evaluate(relErrorsV);
388 
389         // Print the results on console ?
390         if (printResults) {
391             System.out.println();
392             System.out.format(Locale.US, "Relative errors dR/dP -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n",
393                               errorsPMedian, errorsPMean, errorsPMax);
394             System.out.format(Locale.US, "Relative errors dR/dV -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n",
395                               errorsVMedian, errorsVMean, errorsVMax);
396         }
397 
398         Assert.assertEquals(0.0, errorsPMedian, refErrorsPMedian);
399         Assert.assertEquals(0.0, errorsPMean, refErrorsPMean);
400         Assert.assertEquals(0.0, errorsPMax, refErrorsPMax);
401         Assert.assertEquals(0.0, errorsVMedian, refErrorsVMedian);
402         Assert.assertEquals(0.0, errorsVMean, refErrorsVMean);
403         Assert.assertEquals(0.0, errorsVMax, refErrorsVMax);
404     }
405 
406     void genericTestParameterDerivatives(final boolean printResults,
407                                          final double refErrorsMedian, final double refErrorsMean, final double refErrorsMax) {
408 
409         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
410 
411         final NumericalPropagatorBuilder propagatorBuilder =
412                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true,
413                                               1.0e-6, 60.0, 0.001);
414 
415         // Create perfect one-way GNSS range measurements
416         final TimeStampedPVCoordinates original = context.initialOrbit.getPVCoordinates();
417         final Orbit closeOrbit = new CartesianOrbit(new TimeStampedPVCoordinates(context.initialOrbit.getDate(),
418                                                                                  original.getPosition().add(new Vector3D(1000, 2000, 3000)),
419                                                                                  original.getVelocity().add(new Vector3D(-0.03, 0.01, 0.02))),
420                                                     context.initialOrbit.getFrame(),
421                                                     context.initialOrbit.getMu());
422         final Propagator closePropagator = EstimationTestUtils.createPropagator(closeOrbit, propagatorBuilder);
423         final EphemerisGenerator generator = closePropagator.getEphemerisGenerator();
424         closePropagator.propagate(context.initialOrbit.getDate().shiftedBy(3.5 * closeOrbit.getKeplerianPeriod()));
425         final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
426 
427         // Create perfect phase measurements
428         final int    ambiguity         = 1234;
429         final double localClockOffset  = 0.137e-6;
430         final double remoteClockOffset = 469.0e-6;
431         final OneWayGNSSPhaseCreator creator = new OneWayGNSSPhaseCreator(ephemeris, FREQUENCY, ambiguity, localClockOffset, remoteClockOffset);
432         creator.getLocalSatellite().getClockOffsetDriver().setSelected(true);
433 
434         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
435                                                                            propagatorBuilder);
436         final List<ObservedMeasurement<?>> measurements =
437                         EstimationTestUtils.createMeasurements(propagator, creator, 1.0, 3.0, 300.0);
438 
439         // List to store the results
440         final List<Double> relErrorList = new ArrayList<Double>();
441 
442         // Use a lambda function to implement "handleStep" function
443         propagator.setStepHandler(interpolator -> {
444 
445             for (final ObservedMeasurement<?> measurement : measurements) {
446 
447                 //  Play test if the measurement date is between interpolator previous and current date
448                 if ((measurement.getDate().durationFrom(interpolator.getPreviousState().getDate()) > 0.) &&
449                     (measurement.getDate().durationFrom(interpolator.getCurrentState().getDate())  <=  0.)) {
450 
451                     // We intentionally propagate to a date which is close to the
452                     // real spacecraft state but is *not* the accurate date, by
453                     // compensating only part of the downlink delay. This is done
454                     // in order to validate the partial derivatives with respect
455                     // to velocity. If we had chosen the proper state date, the
456                     // phase would have depended only on the current position but
457                     // not on the current velocity.
458                     final double          meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
459                     final AbsoluteDate    date      = measurement.getDate().shiftedBy(-0.75 * meanDelay);
460                     final SpacecraftState[] states    = {
461                         interpolator.getInterpolatedState(date),
462                         ephemeris.propagate(date)
463                     };
464                     final ParameterDriver[] drivers = new ParameterDriver[] {
465                         measurement.getSatellites().get(0).getClockOffsetDriver(),
466                     };
467 
468                     for (int i = 0; i < drivers.length; ++i) {
469                         final double[] gradient  = measurement.estimate(0, 0, states).getParameterDerivatives(drivers[i]);
470                         Assert.assertEquals(1, measurement.getDimension());
471                         Assert.assertEquals(1, gradient.length);
472 
473                         // Compute a reference value using finite differences
474                         final ParameterFunction dMkdP =
475                                         Differentiation.differentiate(new ParameterFunction() {
476                                             /** {@inheritDoc} */
477                                             @Override
478                                             public double value(final ParameterDriver parameterDriver) {
479                                                 return measurement.estimate(0, 0, states).getEstimatedValue()[0];
480                                             }
481                                         }, 3, 20.0 * drivers[i].getScale());
482                         final double ref = dMkdP.value(drivers[i]);
483 
484                         if (printResults) {
485                             System.out.format(Locale.US, "%10.3e  %10.3e  ", gradient[0]-ref, FastMath.abs((gradient[0]-ref)/ref));
486                         }
487 
488                         final double relError = FastMath.abs((ref-gradient[0])/ref);
489                         relErrorList.add(relError);
490                     }
491                     if (printResults) {
492                         System.out.format(Locale.US, "%n");
493                     }
494 
495                 } // End if measurement date between previous and current interpolator step
496             } // End for loop on the measurements
497         });
498 
499         // Rewind the propagator to initial date
500         propagator.propagate(context.initialOrbit.getDate());
501 
502         // Sort measurements chronologically
503         measurements.sort(Comparator.naturalOrder());
504 
505         // Print results ? Header
506         if (printResults) {
507             System.out.format(Locale.US, "%-15s  %-23s  %-23s  " +
508                               "%10s  %10s  %10s  %10s  %10s  %10s  %10s  %10s  %10s  %10s%n",
509                               "Station", "Measurement Date", "State Date",
510                               "Δt",    "rel Δt",
511                               "ΔdQx",  "rel ΔdQx",
512                               "ΔdQy",  "rel ΔdQy",
513                               "ΔdQz",  "rel ΔdQz",
514                               "Δtsat", "rel Δtsat");
515          }
516 
517         // Propagate to final measurement's date
518         propagator.propagate(measurements.get(measurements.size()-1).getDate());
519 
520         // Convert error list to double[]
521         final double relErrors[] = relErrorList.stream().mapToDouble(Double::doubleValue).toArray();
522 
523         // Compute statistics
524         final double relErrorsMedian = new Median().evaluate(relErrors);
525         final double relErrorsMean   = new Mean().evaluate(relErrors);
526         final double relErrorsMax    = new Max().evaluate(relErrors);
527 
528         // Print the results on console ?
529         if (printResults) {
530             System.out.println();
531             System.out.format(Locale.US, "Relative errors dR/dQ -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n",
532                               relErrorsMedian, relErrorsMean, relErrorsMax);
533         }
534 
535         Assert.assertEquals(0.0, relErrorsMedian, refErrorsMedian);
536         Assert.assertEquals(0.0, relErrorsMean, refErrorsMean);
537         Assert.assertEquals(0.0, relErrorsMax, refErrorsMax);
538 
539     }
540 
541     @Test
542     public void testIssue734() {
543 
544         Utils.setDataRoot("regular-data");
545 
546         // Create a phase measurement. Remote is set to null since it not used by the test
547         final OneWayGNSSPhase phase = new OneWayGNSSPhase(null, 635.0e-6, AbsoluteDate.J2000_EPOCH, 467614.701,
548                                                           Frequency.G01.getWavelength(), 0.02, 1.0, new ObservableSatellite(0));
549 
550         // First check
551         Assert.assertEquals(0.0, phase.getAmbiguityDriver().getValue(), Double.MIN_VALUE);
552         Assert.assertFalse(phase.getAmbiguityDriver().isSelected());
553 
554         // Perform some changes in ambiguity driver
555         phase.getAmbiguityDriver().setValue(1234.0);
556         phase.getAmbiguityDriver().setSelected(true);
557 
558         // Second check
559         Assert.assertEquals(1234.0, phase.getAmbiguityDriver().getValue(), Double.MIN_VALUE);
560         Assert.assertTrue(phase.getAmbiguityDriver().isSelected());
561         for (ParameterDriver driver : phase.getParametersDrivers()) {
562             // Verify if the current driver corresponds to the phase ambiguity
563             if (driver.getName() == Phase.AMBIGUITY_NAME) {
564                 Assert.assertEquals(1234.0, phase.getAmbiguityDriver().getValue(), Double.MIN_VALUE);
565                 Assert.assertTrue(phase.getAmbiguityDriver().isSelected());
566             }
567         }
568 
569     }
570 
571 }