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;
18  
19  import java.util.ArrayList;
20  import java.util.List;
21  import java.util.Locale;
22  import java.util.Map;
23  
24  import org.hipparchus.stat.descriptive.moment.Mean;
25  import org.hipparchus.stat.descriptive.rank.Max;
26  import org.hipparchus.stat.descriptive.rank.Median;
27  import org.hipparchus.stat.descriptive.rank.Min;
28  import org.hipparchus.util.FastMath;
29  import org.junit.jupiter.api.Assertions;
30  import org.junit.jupiter.api.Test;
31  import org.orekit.estimation.Context;
32  import org.orekit.estimation.EstimationTestUtils;
33  import org.orekit.estimation.measurements.modifiers.TurnAroundRangeTroposphericDelayModifier;
34  import org.orekit.models.earth.troposphere.ModifiedSaastamoinenModel;
35  import org.orekit.orbits.OrbitType;
36  import org.orekit.orbits.PositionAngleType;
37  import org.orekit.propagation.Propagator;
38  import org.orekit.propagation.SpacecraftState;
39  import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
40  import org.orekit.time.AbsoluteDate;
41  import org.orekit.utils.Constants;
42  import org.orekit.utils.Differentiation;
43  import org.orekit.utils.ParameterDriver;
44  import org.orekit.utils.ParameterFunction;
45  
46  public class TurnAroundRangeAnalyticTest {
47  
48  
49      /**
50       * Test the values of the TAR (Turn-Around Range) comparing the observed values and the estimated values
51       * Both are calculated with a different algorithm
52       */
53      @Test
54      public void testValues() {
55          boolean printResults = false;
56          if (printResults) {
57              System.out.println("\nTest TAR Analytical Values\n");
58          }
59          // Run test
60          genericTestValues(printResults);
61      }
62  
63      /**
64       * Test the values of the analytic state derivatives computation
65       * using TurnAroundRange function as a comparison
66       */
67      @Test
68      public void testStateDerivatives() {
69  
70          boolean printResults = false;
71          if (printResults) {
72              System.out.println("\nTest TAR Analytical State Derivatives - Derivative Structure comparison\n");
73          }
74          // Run test
75          boolean isModifier = false;
76          boolean isFiniteDifferences  = false;
77          genericTestStateDerivatives(isModifier, isFiniteDifferences, printResults,
78                                      1.4e-6, 1.4e-6, 2.6e-6, 9.7e-7, 1.0e-6, 2.1e-6);
79      }
80  
81      /**
82       * Test the values of the analytic state derivatives
83       * using a numerical finite differences calculation as a reference
84       */
85      @Test
86      public void testStateDerivativesFiniteDifferences() {
87  
88          boolean printResults = false;
89          if (printResults) {
90              System.out.println("\nTest TAR Analytical State Derivatives - Finite Differences comparison\n");
91          }
92          // Run test
93          boolean isModifier = false;
94          boolean isFiniteDifferences  = true;
95          genericTestStateDerivatives(isModifier, isFiniteDifferences, printResults,
96                                      7.0e-9, 2.5e-8, 3.9e-7, 8.2e-5, 3.1e-4, 7.3e-3);
97      }
98  
99      /**
100      * Test the values of the state derivatives with modifier
101      * using TurnAroundRange function as a comparison
102      */
103     @Test
104     public void testStateDerivativesWithModifier() {
105 
106         boolean printResults = false;
107         if (printResults) {
108             System.out.println("\nTest TAR Analytical State Derivatives with Modifier - Derivative Structure comparison\n");
109         }
110         // Run test
111         boolean isModifier = true;
112         boolean isFiniteDifferences  = false;
113         genericTestStateDerivatives(isModifier, isFiniteDifferences, printResults,
114                                     1.5e-6, 2.0e-6, 2.3e-5, 1.0e-6, 1.0e-6, 2.1e-6);
115     }
116 
117     /**
118      * Test the values of the state derivatives with modifier
119      * using a numerical finite differences calculation as a reference
120      */
121     @Test
122     public void testStateDerivativesWithModifierFiniteDifferences() {
123 
124         boolean printResults = false;
125         if (printResults) {
126             System.out.println("\nTest TAR Analytical State Derivatives with Modifier - Finite Differences comparison\n");
127         }
128         // Run test
129         boolean isModifier = true;
130         boolean isFiniteDifferences  = true;
131         genericTestStateDerivatives(isModifier, isFiniteDifferences, printResults,
132                                     3.1e-8, 9.9e-8, 1.8e-6, 8.5e-5, 3.2e-4, 1.2e-2);
133     }
134 
135     /**
136      * Test the values of the analytic parameters derivatives computation
137      * using TurnAroundRange function as a comparison
138      */
139     @Test
140     public void testParameterDerivatives() {
141 
142         // Print the results ?
143         boolean printResults = false;
144 
145         if (printResults) {
146             System.out.println("\nTest TAR Analytical Parameter Derivatives - Derivative Structure comparison\n");
147         }
148         // Run test
149         boolean isModifier = false;
150         boolean isFiniteDifferences  = false;
151         genericTestParameterDerivatives(isModifier, isFiniteDifferences, printResults,
152                                         4.4e-06, 6.9e-06, 1.3e-04, 3.4e-6, 4.9e-6, 3.6e-5);
153 
154     }
155 
156     /**
157      * Test the values of the analytic parameters derivatives
158      * using a numerical finite differences calculation as a reference
159      */
160     @Test
161     public void testParameterDerivativesFiniteDifferences() {
162 
163         // Print the results ?
164         boolean printResults = false;
165 
166         if (printResults) {
167             System.out.println("\nTest TAR Analytical Parameter Derivatives - Finite Differences comparison\n");
168         }
169         // Run test
170         boolean isModifier = false;
171         boolean isFiniteDifferences  = true;
172         genericTestParameterDerivatives(isModifier, isFiniteDifferences, printResults,
173                                         3.0e-06, 5.9e-06, 1.3e-04, 2.9e-6, 5.0e-6, 3.9e-5);
174 
175     }
176 
177     /**
178      * Test the values of the analytic parameters derivatives with modifiers computation
179      * using TurnAroundRange function as a comparison
180      */
181     @Test
182     public void testParameterDerivativesWithModifier() {
183 
184         // Print the results ?
185         boolean printResults = false;
186 
187         if (printResults) {
188             System.out.println("\nTest TAR Analytical Parameter Derivatives with Modifier - Derivative Structure comparison\n");
189         }
190         // Run test
191         boolean isModifier = true;
192         boolean isFiniteDifferences  = false;
193         genericTestParameterDerivatives(isModifier, isFiniteDifferences, printResults,
194                                         4.4e-06, 6.9e-06, 1.3e-04, 3.4e-6, 4.9e-6, 3.6e-5);
195 
196     }
197 
198     /**
199      * Test the values of the analytic parameters derivatives with modifiers computation
200      * using a numerical finite differences calculation as a reference
201      */
202     @Test
203     public void testParameterDerivativesWithModifierFiniteDifferences() {
204 
205         // Print the results ?
206         boolean printResults = false;
207 
208         if (printResults) {
209             System.out.println("\nTest TAR Analytical Parameter Derivatives with Modifier - Finite Differences comparison\n");
210         }
211         // Run test
212         boolean isModifier = true;
213         boolean isFiniteDifferences  = true;
214         genericTestParameterDerivatives(isModifier, isFiniteDifferences, printResults,
215                                         2.8e-06, 5.9e-06, 1.2e-04, 4.3e-6, 2.3e-5, 2.2e-4);
216 
217     }
218 
219     /**
220      * Generic test function for values of the TAR
221      * @param printResults Print the results ?
222      */
223     void genericTestValues(final boolean printResults) {
224 
225         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
226         //Context context = EstimationTestUtils.geoStationnaryContext();
227 
228         final NumericalPropagatorBuilder propagatorBuilder =
229                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
230                                               1.0e-6, 60.0, 0.001);
231 
232         // create perfect range measurements
233         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
234                                                                            propagatorBuilder);
235         final List<ObservedMeasurement<?>> measurements =
236                         EstimationTestUtils.createMeasurements(propagator,
237                                                                new TurnAroundRangeMeasurementCreator(context),
238                                                                1.0, 3.0, 300.0);
239         propagator.clearStepHandlers();
240 
241         double[] absoluteErrors = new double[measurements.size()];
242         double[] relativeErrors = new double[measurements.size()];
243         int index = 0;
244 
245         // Print the results ? Header
246         if (printResults) {
247            System.out.format(Locale.US, "%-15s  %-15s  %-23s  %-23s  %17s  %17s  %13s %13s%n",
248                               "Primary Station", "secondary Station",
249                               "Measurement Date", "State Date",
250                               "TAR observed [m]", "TAR estimated [m]",
251                               "|ΔTAR| [m]", "rel |ΔTAR|");
252         }
253 
254         // Loop on the measurements
255         for (final ObservedMeasurement<?> measurement : measurements) {
256 
257             final double          meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
258             final AbsoluteDate    date      = measurement.getDate().shiftedBy(meanDelay);
259             final SpacecraftState state     = propagator.propagate(date);
260 
261             // Values of the TAR & errors
262             final double TARobserved  = measurement.getObservedValue()[0];
263             final double TARestimated = new TurnAroundRangeAnalytic((TurnAroundRange) measurement).
264                                                                     theoreticalEvaluationAnalytic(0, 0, propagator.getInitialState(), state).getEstimatedValue()[0];
265 
266             absoluteErrors[index] = TARestimated-TARobserved;
267             relativeErrors[index] = FastMath.abs(absoluteErrors[index])/FastMath.abs(TARobserved);
268             index++;
269 
270             // Print results ? Values
271             if (printResults) {
272                 final AbsoluteDate measurementDate = measurement.getDate();
273 
274                 String primaryStationName = ((TurnAroundRange) measurement).getPrimaryStation().getBaseFrame().getName();
275                 String secondaryStationName = ((TurnAroundRange) measurement).getSecondaryStation().getBaseFrame().getName();
276                 System.out.format(Locale.US, "%-15s  %-15s  %-23s  %-23s  %17.6f  %17.6f  %13.6e %13.6e%n",
277                                   primaryStationName, secondaryStationName, measurementDate, date,
278                                  TARobserved, TARestimated,
279                                  FastMath.abs(TARestimated-TARobserved),
280                                  FastMath.abs((TARestimated-TARobserved)/TARobserved));
281             }
282 
283         }
284 
285         // Compute some statistics
286         final double absErrorsMedian = new Median().evaluate(absoluteErrors);
287         final double absErrorsMin    = new Min().evaluate(absoluteErrors);
288         final double absErrorsMax    = new Max().evaluate(absoluteErrors);
289         final double relErrorsMedian = new Median().evaluate(relativeErrors);
290         final double relErrorsMax    = new Max().evaluate(relativeErrors);
291 
292         // Print the results on console ? Final results
293         if (printResults) {
294             System.out.println();
295             System.out.println("Absolute errors median: " +  absErrorsMedian);
296             System.out.println("Absolute errors min   : " +  absErrorsMin);
297             System.out.println("Absolute errors max   : " +  absErrorsMax);
298             System.out.println("Relative errors median: " +  relErrorsMedian);
299             System.out.println("Relative errors max   : " +  relErrorsMax);
300         }
301 
302         // Assert statistical errors
303         Assertions.assertEquals(0.0, absErrorsMedian, 8.5e-08);
304         Assertions.assertEquals(0.0, absErrorsMin,    9.7e-08);
305         Assertions.assertEquals(0.0, absErrorsMax,    2.0e-07);
306         Assertions.assertEquals(0.0, relErrorsMedian, 5.1e-15);
307         Assertions.assertEquals(0.0, relErrorsMax,    1.2e-14);
308  
309         // Test measurement type
310         final TurnAroundRangeAnalytic taRangeAnalytic =
311                         new TurnAroundRangeAnalytic((TurnAroundRange) measurements.get(0));
312         Assertions.assertEquals(TurnAroundRangeAnalytic.MEASUREMENT_TYPE, taRangeAnalytic.getMeasurementType());
313     }
314 
315     /**
316      * Generic test function for derivatives with respect to state
317      * @param isModifier Use of atmospheric modifiers
318      * @param isFiniteDifferences Finite differences reference calculation if true, TurnAroundRange class otherwise
319      * @param printResults Print the results ?
320       */
321     void genericTestStateDerivatives(final boolean isModifier, final boolean isFiniteDifferences, final boolean printResults,
322                                      final double refErrorsPMedian, final double refErrorsPMean,
323                                      final double refErrorsPMax, final double refErrorsVMedian,
324                                      final double refErrorsVMean, final double refErrorsVMax)
325                     {
326 
327         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
328         //Context context = EstimationTestUtils.geoStationnaryContext();
329 
330         final NumericalPropagatorBuilder propagatorBuilder =
331                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
332                                               1.0e-6, 60.0, 0.001);
333 
334         // create perfect range2 measurements
335         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
336                                                                            propagatorBuilder);
337         final List<ObservedMeasurement<?>> measurements =
338                         EstimationTestUtils.createMeasurements(propagator,
339                                                                new TurnAroundRangeMeasurementCreator(context),
340                                                                1.0, 3.0, 300.0);
341         propagator.clearStepHandlers();
342 
343         double[] errorsP = new double[3 * measurements.size()];
344         double[] errorsV = new double[3 * measurements.size()];
345         int indexP = 0;
346         int indexV = 0;
347 
348         // Print the results ? Header
349         if (printResults) {
350             System.out.format(Locale.US, "%-15s  %-15s  %-23s  %-23s  " +
351                             "%10s  %10s  %10s  " +
352                             "%10s  %10s  %10s  " +
353                             "%10s  %10s  %10s  " +
354                             "%10s  %10s  %10s%n",
355                             "Primary Station", "secondary Station",
356                             "Measurement Date", "State Date",
357                             "ΔdPx", "ΔdPy", "ΔdPz", "ΔdVx", "ΔdVy", "ΔdVz",
358                             "rel ΔdPx", "rel ΔdPy", "rel ΔdPz",
359                             "rel ΔdVx", "rel ΔdVy", "rel ΔdVz");
360         }
361 
362         // Loop on the measurements
363         for (final ObservedMeasurement<?> measurement : measurements) {
364 
365             // Add modifiers if test implies it
366             final TurnAroundRangeTroposphericDelayModifier modifier =
367                             new TurnAroundRangeTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
368             if (isModifier) {
369                 ((TurnAroundRange) measurement).addModifier(modifier);
370             }
371 
372             // We intentionally propagate to a date which is close to the
373             // real spacecraft state but is *not* the accurate date, by
374             // compensating only part of the downlink delay. This is done
375             // in order to validate the partial derivatives with respect
376             // to velocity. If we had chosen the proper state date, the
377             // range would have depended only on the current position but
378             // not on the current velocity.
379             final double          meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
380             final AbsoluteDate    date      = measurement.getDate().shiftedBy(-0.75 * meanDelay);
381             final SpacecraftState state     = propagator.propagate(date);
382             final EstimatedMeasurement<TurnAroundRange> TAR = new TurnAroundRangeAnalytic((TurnAroundRange)measurement).
383                             theoreticalEvaluationAnalytic(0, 0, propagator.getInitialState(), state);
384             if (isModifier) {
385                 modifier.modify(TAR);
386             }
387             final double[][] jacobian  = TAR.getStateDerivatives(0);
388 
389             // Jacobian reference value
390             final double[][] jacobianRef;
391 
392             if (isFiniteDifferences) {
393                 // Compute a reference value using finite differences
394                 jacobianRef = Differentiation.differentiate(
395                     state1 -> measurement.
396                            estimateWithoutDerivatives(new SpacecraftState[] { state1 }).
397                            getEstimatedValue(), measurement.getDimension(), propagator.getAttitudeProvider(),
398                     OrbitType.CARTESIAN, PositionAngleType.TRUE, 2.0, 3).value(state);
399             } else {
400                 // Compute a reference value using TurnAroundRange class function
401                 jacobianRef = ((TurnAroundRange) measurement).theoreticalEvaluation(0, 0, new SpacecraftState[] { state }).getStateDerivatives(0);
402             }
403 
404 //            //Test: Test point by point with the debugger
405 //            if (!isFiniteDifferences && !isModifier) {
406 //                final EstimatedMeasurement<TurnAroundRange> test =
407 //                                new TurnAroundRangeAnalytic((TurnAroundRange)measurement).theoreticalEvaluationValidation(0, 0, state);
408 //            }
409 //            //Test
410 
411             Assertions.assertEquals(jacobianRef.length, jacobian.length);
412             Assertions.assertEquals(jacobianRef[0].length, jacobian[0].length);
413 
414             double [][] dJacobian         = new double[jacobian.length][jacobian[0].length];
415             double [][] dJacobianRelative = new double[jacobian.length][jacobian[0].length];
416             for (int i = 0; i < jacobian.length; ++i) {
417                 for (int j = 0; j < jacobian[i].length; ++j) {
418                     dJacobian[i][j] = jacobian[i][j] - jacobianRef[i][j];
419                     dJacobianRelative[i][j] = FastMath.abs(dJacobian[i][j]/jacobianRef[i][j]);
420 
421                     if (j < 3) {
422                         errorsP[indexP++] = dJacobianRelative[i][j];
423                     } else {
424                         errorsV[indexV++] = dJacobianRelative[i][j];
425                     }
426                 }
427             }
428             // Print results on the console ? Print the Jacobian
429             if (printResults) {
430                 String primaryStationName = ((TurnAroundRange) measurement).getPrimaryStation().getBaseFrame().getName();
431                 String secondaryStationName  = ((TurnAroundRange) measurement).getSecondaryStation().getBaseFrame().getName();
432                 System.out.format(Locale.US, "%-15s  %-15s  %-23s  %-23s  " +
433                                 "%10.3e  %10.3e  %10.3e  " +
434                                 "%10.3e  %10.3e  %10.3e  " +
435                                 "%10.3e  %10.3e  %10.3e  " +
436                                 "%10.3e  %10.3e  %10.3e%n",
437                                 primaryStationName, secondaryStationName, measurement.getDate(), date,
438                                 dJacobian[0][0], dJacobian[0][1], dJacobian[0][2],
439                                 dJacobian[0][3], dJacobian[0][4], dJacobian[0][5],
440                                 dJacobianRelative[0][0], dJacobianRelative[0][1], dJacobianRelative[0][2],
441                                 dJacobianRelative[0][3], dJacobianRelative[0][4], dJacobianRelative[0][5]);
442             }
443         } // End loop on the measurements
444 
445         // Compute some statistics
446         final double errorsPMedian = new Median().evaluate(errorsP);
447         final double errorsPMean   = new Mean().evaluate(errorsP);
448         final double errorsPMax    = new Max().evaluate(errorsP);
449         final double errorsVMedian = new Median().evaluate(errorsV);
450         final double errorsVMean   = new Mean().evaluate(errorsV);
451         final double errorsVMax    = new Max().evaluate(errorsV);
452 
453 
454         // Print the results on console ? Final results
455         if (printResults) {
456             System.out.println();
457             System.out.format(Locale.US, "Relative errors dR/dP -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n",
458                               errorsPMedian, errorsPMean, errorsPMax);
459             System.out.format(Locale.US, "Relative errors dR/dV -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n",
460                               errorsVMedian, errorsVMean, errorsVMax);
461         }
462 
463         // Assert the results / max values depend on the test
464         Assertions.assertEquals(0.0, errorsPMedian, refErrorsPMedian);
465         Assertions.assertEquals(0.0, errorsPMean, refErrorsPMean);
466         Assertions.assertEquals(0.0, errorsPMax, refErrorsPMax);
467         Assertions.assertEquals(0.0, errorsVMedian, refErrorsVMedian);
468         Assertions.assertEquals(0.0, errorsVMean, refErrorsVMean);
469         Assertions.assertEquals(0.0, errorsVMax, refErrorsVMax);
470     }
471 
472 
473     /**
474      * Generic test function for derivatives with respect to parameters (station's position in station's topocentric frame)
475      * @param isModifier Use of atmospheric modifiers
476      * @param isFiniteDifferences Finite differences reference calculation if true, TurnAroundRange class otherwise
477      * @param printResults Print the results ?
478      */
479     void genericTestParameterDerivatives(final boolean isModifier, final boolean isFiniteDifferences, final boolean printResults,
480                                          final double refErrorQMMedian, final double refErrorQMMean, final double refErrorQMMax,
481                                          final double refErrorQSMedian, final double refErrorQSMean, final double refErrorQSMax)
482                     {
483 
484         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
485 
486         final NumericalPropagatorBuilder propagatorBuilder =
487                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
488                                               1.0e-6, 60.0, 0.001);
489 
490         // Create perfect TAR measurements
491         for (Map.Entry<GroundStation, GroundStation> entry : context.TARstations.entrySet()) {
492             final GroundStation    primaryStation = entry.getKey();
493             final GroundStation    secondaryStation  = entry.getValue();
494             primaryStation.getClockOffsetDriver().setSelected(false);
495             primaryStation.getEastOffsetDriver().setSelected(true);
496             primaryStation.getNorthOffsetDriver().setSelected(true);
497             primaryStation.getZenithOffsetDriver().setSelected(true);
498             secondaryStation.getClockOffsetDriver().setSelected(false);
499             secondaryStation.getEastOffsetDriver().setSelected(true);
500             secondaryStation.getNorthOffsetDriver().setSelected(true);
501             secondaryStation.getZenithOffsetDriver().setSelected(true);
502         }
503         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
504                                                                            propagatorBuilder);
505         final List<ObservedMeasurement<?>> measurements =
506                         EstimationTestUtils.createMeasurements(propagator,
507                                                                new TurnAroundRangeMeasurementCreator(context),
508                                                                1.0, 3.0, 300.0);
509         propagator.clearStepHandlers();
510 
511         // Print results on console ? Header
512         if (printResults) {
513             System.out.format(Locale.US, "%-15s %-15s %-23s  %-23s  " +
514                             "%10s  %10s  %10s  " +
515                             "%10s  %10s  %10s  " +
516                             "%10s  %10s  %10s  " +
517                             "%10s  %10s  %10s%n",
518                             "Primary Station", "secondary Station",
519                             "Measurement Date", "State Date",
520                             "ΔdQMx", "rel ΔdQMx",
521                             "ΔdQMy", "rel ΔdQMy",
522                             "ΔdQMz", "rel ΔdQMz",
523                             "ΔdQSx", "rel ΔdQSx",
524                             "ΔdQSy", "rel ΔdQSy",
525                             "ΔdQSz", "rel ΔdQSz");
526          }
527 
528         // List to store the results for primary and secondary station
529         final List<Double> relErrorQMList = new ArrayList<>();
530         final List<Double> relErrorQSList = new ArrayList<>();
531 
532         // Loop on the measurements
533         for (final ObservedMeasurement<?> measurement : measurements) {
534 
535             // Add modifiers if test implies it
536             final TurnAroundRangeTroposphericDelayModifier modifier =
537                 new TurnAroundRangeTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
538             if (isModifier) {
539                 ((TurnAroundRange) measurement).addModifier(modifier);
540             }
541 
542             // parameter corresponding to station position offset
543             final GroundStation primaryStationParameter = ((TurnAroundRange) measurement).getPrimaryStation();
544             final GroundStation secondaryStationParameter = ((TurnAroundRange) measurement).getSecondaryStation();
545 
546             // We intentionally propagate to a date which is close to the
547             // real spacecraft state but is *not* the accurate date, by
548             // compensating only part of the downlink delay. This is done
549             // in order to validate the partial derivatives with respect
550             // to velocity. If we had chosen the proper state date, the
551             // range would have depended only on the current position but
552             // not on the current velocity.
553             final double          meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
554             final AbsoluteDate    date      = measurement.getDate().shiftedBy(-0.75 * meanDelay);
555             final SpacecraftState state     = propagator.propagate(date);
556             final ParameterDriver[] drivers = new ParameterDriver[] {
557                                                                      primaryStationParameter.getEastOffsetDriver(),
558                                                                      primaryStationParameter.getNorthOffsetDriver(),
559                                                                      primaryStationParameter.getZenithOffsetDriver(),
560                                                                      secondaryStationParameter.getEastOffsetDriver(),
561                                                                      secondaryStationParameter.getNorthOffsetDriver(),
562                                                                      secondaryStationParameter.getZenithOffsetDriver(),
563             };
564 
565             // Print results on console ? Stations' names
566             if (printResults) {
567                 String primaryStationName  = primaryStationParameter.getBaseFrame().getName();
568                 String secondaryStationName  = secondaryStationParameter.getBaseFrame().getName();
569                 System.out.format(Locale.US, "%-15s %-15s %-23s  %-23s  ",
570                                   primaryStationName, secondaryStationName, measurement.getDate(), date);
571             }
572 
573             // Loop on the parameters
574             for (int i = 0; i < 6; ++i) {
575                 // Analytical computation of the parameters derivatives
576                 final EstimatedMeasurement<TurnAroundRange> TAR =
577                                 new TurnAroundRangeAnalytic((TurnAroundRange) measurement).
578                                 theoreticalEvaluationAnalytic(0, 0, propagator.getInitialState(), state);
579                 // Optional modifier addition
580                 if (isModifier) {
581                   modifier.modify(TAR);
582                 }
583                 final double[] gradient  = TAR.getParameterDerivatives(drivers[i], new AbsoluteDate());
584 
585                 Assertions.assertEquals(1, measurement.getDimension());
586                 Assertions.assertEquals(1, gradient.length);
587 
588                 // Reference value
589                 double ref;
590                 if (isFiniteDifferences) {
591                     // Compute a reference value using finite differences
592                     final ParameterFunction dMkdP =
593                                     Differentiation.differentiate(new ParameterFunction() {
594                                         /** {@inheritDoc} */
595                                         @Override
596                                         public double value(final ParameterDriver parameterDriver, AbsoluteDate date) {
597                                             return measurement.
598                                                    estimateWithoutDerivatives(new SpacecraftState[] { state }).
599                                                    getEstimatedValue()[0];
600                                         }
601                                     }, 3, 20.0 * drivers[i].getScale());
602                     ref = dMkdP.value(drivers[i], date);
603                 } else {
604                     // Compute a reference value using TurnAroundRange function
605                     ref = measurement.estimate(0, 0, new SpacecraftState[] { state }).getParameterDerivatives(drivers[i], new AbsoluteDate())[0];
606                 }
607 
608                 // Deltas
609                 double dGradient         = gradient[0] - ref;
610                 double dGradientRelative = FastMath.abs(dGradient/ref);
611 
612                 // Print results on console ? Gradient difference
613                 if (printResults) {
614                     System.out.format(Locale.US, "%10.3e  %10.3e  ", dGradient, dGradientRelative);
615                 }
616 
617                 // Add relative error to the list
618                 if (i<3) { relErrorQMList.add(dGradientRelative);}
619                 else     { relErrorQSList.add(dGradientRelative);}
620 
621             } // End for loop on the parameters
622             if (printResults) {
623                 System.out.format(Locale.US, "%n");
624             }
625         } // End for loop on the measurements
626 
627         // Convert error list to double[]
628         final double[] relErrorQM = relErrorQMList.stream().mapToDouble(Double::doubleValue).toArray();
629         final double[] relErrorQS = relErrorQSList.stream().mapToDouble(Double::doubleValue).toArray();
630 
631         // Compute statistics
632         final double relErrorsQMMedian = new Median().evaluate(relErrorQM);
633         final double relErrorsQMMean   = new Mean().evaluate(relErrorQM);
634         final double relErrorsQMMax    = new Max().evaluate(relErrorQM);
635 
636         final double relErrorsQSMedian = new Median().evaluate(relErrorQS);
637         final double relErrorsQSMean   = new Mean().evaluate(relErrorQS);
638         final double relErrorsQSMax    = new Max().evaluate(relErrorQS);
639 
640 
641         // Print the results on console ?
642         if (printResults) {
643             System.out.println();
644             System.out.format(Locale.US, "Relative errors dR/dQ primary station -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n",
645                               relErrorsQMMedian, relErrorsQMMean, relErrorsQMMax);
646             System.out.format(Locale.US, "Relative errors dR/dQ secondary station  -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n",
647                               relErrorsQSMedian, relErrorsQSMean, relErrorsQSMax);
648         }
649 
650         // Check values
651         Assertions.assertEquals(0.0, relErrorsQMMedian, refErrorQMMedian);
652         Assertions.assertEquals(0.0, relErrorsQMMean, refErrorQMMean);
653         Assertions.assertEquals(0.0, relErrorsQMMax, refErrorQMMax);
654         Assertions.assertEquals(0.0, relErrorsQSMedian, refErrorQSMedian);
655         Assertions.assertEquals(0.0, relErrorsQSMean, refErrorQSMean);
656         Assertions.assertEquals(0.0, relErrorsQSMax, refErrorQSMax);
657 
658     }
659 
660 }