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