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