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