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.Comparator;
21  import java.util.List;
22  import java.util.Locale;
23  
24  import org.hipparchus.geometry.euclidean.threed.Vector3D;
25  import org.hipparchus.stat.descriptive.moment.Mean;
26  import org.hipparchus.stat.descriptive.rank.Max;
27  import org.hipparchus.stat.descriptive.rank.Median;
28  import org.hipparchus.stat.descriptive.rank.Min;
29  import org.hipparchus.util.FastMath;
30  import org.junit.Assert;
31  import org.junit.Test;
32  import org.orekit.estimation.Context;
33  import org.orekit.estimation.EstimationTestUtils;
34  import org.orekit.estimation.measurements.modifiers.RangeTroposphericDelayModifier;
35  import org.orekit.models.earth.troposphere.EstimatedTroposphericModel;
36  import org.orekit.models.earth.troposphere.NiellMappingFunctionModel;
37  import org.orekit.models.earth.troposphere.SaastamoinenModel;
38  import org.orekit.orbits.OrbitType;
39  import org.orekit.orbits.PositionAngle;
40  import org.orekit.propagation.Propagator;
41  import org.orekit.propagation.SpacecraftState;
42  import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
43  import org.orekit.time.AbsoluteDate;
44  import org.orekit.utils.Constants;
45  import org.orekit.utils.Differentiation;
46  import org.orekit.utils.ParameterDriver;
47  import org.orekit.utils.ParameterFunction;
48  import org.orekit.utils.StateFunction;
49  import org.orekit.utils.TimeStampedPVCoordinates;
50  
51  public class RangeTest {
52  
53      /**
54       * Test the values of the range comparing the observed values and the estimated values
55       * Both are calculated with a different algorithm
56       */
57      @Test
58      public void testValues() {
59          boolean printResults = false;
60          if (printResults) {
61              System.out.println("\nTest Range Values\n");
62          }
63          // Run test
64          this.genericTestValues(printResults);
65      }
66  
67      /**
68       * Test the values of the state derivatives using a numerical
69       * finite differences calculation as a reference
70       */
71      @Test
72      public void testStateDerivatives() {
73  
74          boolean printResults = false;
75          if (printResults) {
76              System.out.println("\nTest Range State Derivatives - Finite Differences Comparison\n");
77          }
78          // Run test
79          boolean isModifier = false;
80          double refErrorsPMedian = 6.5e-10;
81          double refErrorsPMean   = 4.1e-09;
82          double refErrorsPMax    = 2.1e-07;
83          double refErrorsVMedian = 2.2e-04;
84          double refErrorsVMean   = 6.2e-04;
85          double refErrorsVMax    = 1.3e-02;
86          this.genericTestStateDerivatives(isModifier, printResults,
87                                           refErrorsPMedian, refErrorsPMean, refErrorsPMax,
88                                           refErrorsVMedian, refErrorsVMean, refErrorsVMax);
89      }
90  
91      /**
92       * Test the values of the state derivatives with modifier using a numerical
93       * finite differences calculation as a reference
94       */
95      @Test
96      public void testStateDerivativesWithModifier() {
97  
98          boolean printResults = false;
99          if (printResults) {
100             System.out.println("\nTest Range State Derivatives with Modifier - Finite Differences Comparison\n");
101         }
102         // Run test
103         boolean isModifier = true;
104         double refErrorsPMedian = 6.2e-10;
105         double refErrorsPMean   = 3.5e-09;
106         double refErrorsPMax    = 1.6e-07;
107         double refErrorsVMedian = 2.2e-04;
108         double refErrorsVMean   = 6.2e-04;
109         double refErrorsVMax    = 1.3e-02;
110         this.genericTestStateDerivatives(isModifier, printResults,
111                                          refErrorsPMedian, refErrorsPMean, refErrorsPMax,
112                                          refErrorsVMedian, refErrorsVMean, refErrorsVMax);
113     }
114 
115     /**
116      * Test the values of the parameters' derivatives using a numerical
117      * finite differences calculation as a reference
118      */
119     @Test
120     public void testParameterDerivatives() {
121 
122         // Print the results ?
123         boolean printResults = false;
124 
125         if (printResults) {
126             System.out.println("\nTest Range Parameter Derivatives - Finite Differences Comparison\n");
127         }
128         // Run test
129         boolean isModifier = false;
130         double refErrorsMedian = 5.7e-9;
131         double refErrorsMean   = 6.4e-8;
132         double refErrorsMax    = 5.1e-6;
133         this.genericTestParameterDerivatives(isModifier, printResults,
134                                              refErrorsMedian, refErrorsMean, refErrorsMax);
135 
136     }
137 
138     /**
139      * Test the values of the parameters' derivatives with modifier, using a numerical
140      * finite differences calculation as a reference
141      */
142     @Test
143     public void testParameterDerivativesWithModifier() {
144 
145         // Print the results ?
146         boolean printResults = false;
147 
148         if (printResults) {
149             System.out.println("\nTest Range Parameter Derivatives with Modifier - Finite Differences Comparison\n");
150         }
151         // Run test
152         boolean isModifier = true;
153         double refErrorsMedian = 2.9e-8;
154         double refErrorsMean   = 4.9e-7;
155         double refErrorsMax    = 1.5e-5;
156         this.genericTestParameterDerivatives(isModifier, printResults,
157                                              refErrorsMedian, refErrorsMean, refErrorsMax);
158 
159     }
160 
161     /**
162      * Test the values of the parameters' derivatives with estimated modifier, using a numerical
163      * finite differences calculation as a reference
164      */
165     @Test
166     public void testParameterDerivativesWithEstimatedModifier() {
167 
168         // Print the results ?
169         boolean printResults = false;
170 
171         if (printResults) {
172             System.out.println("\nTest Range Parameter Derivatives with Estimated Modifier - Finite Differences Comparison\n");
173         }
174         // Run test
175         boolean isModifier = true;
176         double refErrorsMedian = 1.2e-9;
177         double refErrorsMean   = 1.9e-9;
178         double refErrorsMax    = 6.6e-9;
179         this.genericTestEstimatedParameterDerivatives(isModifier, printResults,
180                                                       refErrorsMedian, refErrorsMean, refErrorsMax);
181 
182     }
183 
184     /**
185      * Generic test function for values of the range
186      * @param printResults Print the results ?
187      */
188     void genericTestValues(final boolean printResults) {
189 
190         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
191 
192         final NumericalPropagatorBuilder propagatorBuilder =
193                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true,
194                                               1.0e-6, 60.0, 0.001);
195 
196         // Create perfect range measurements
197         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
198                                                                            propagatorBuilder);
199         final double groundClockOffset = 1.234e-3;
200         for (final GroundStation station : context.stations) {
201             station.getClockOffsetDriver().setValue(groundClockOffset);
202         }
203         final List<ObservedMeasurement<?>> measurements =
204                         EstimationTestUtils.createMeasurements(propagator,
205                                                                new RangeMeasurementCreator(context, Vector3D.ZERO),
206                                                                1.0, 3.0, 300.0);
207 
208         // Lists for results' storage - Used only for derivatives with respect to state
209         // "final" value to be seen by "handleStep" function of the propagator
210         final List<Double> absoluteErrors = new ArrayList<Double>();
211         final List<Double> relativeErrors = new ArrayList<Double>();
212 
213         // Set step handler
214         // Use a lambda function to implement "handleStep" function
215         propagator.setStepHandler(interpolator -> {
216 
217             for (final ObservedMeasurement<?> measurement : measurements) {
218 
219                 //  Play test if the measurement date is between interpolator previous and current date
220                 if ((measurement.getDate().durationFrom(interpolator.getPreviousState().getDate()) > 0.) &&
221                     (measurement.getDate().durationFrom(interpolator.getCurrentState().getDate())  <=  0.)
222                    ) {
223                     // We intentionally propagate to a date which is close to the
224                     // real spacecraft state but is *not* the accurate date, by
225                     // compensating only part of the downlink delay. This is done
226                     // in order to validate the partial derivatives with respect
227                     // to velocity. If we had chosen the proper state date, the
228                     // range would have depended only on the current position but
229                     // not on the current velocity.
230                     final double          meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
231                     final AbsoluteDate    date      = measurement.getDate().shiftedBy(-0.75 * meanDelay);
232                     final SpacecraftState state     = interpolator.getInterpolatedState(date);
233 
234                     // Values of the Range & errors
235                     final double RangeObserved  = measurement.getObservedValue()[0];
236                     final EstimatedMeasurement<?> estimated = measurement.estimate(0, 0, new SpacecraftState[] { state });
237 
238                     final TimeStampedPVCoordinates[] participants = estimated.getParticipants();
239                     Assert.assertEquals(3, participants.length);
240                     Assert.assertEquals(0.5 * Constants.SPEED_OF_LIGHT * participants[2].getDate().durationFrom(participants[0].getDate()),
241                                         estimated.getEstimatedValue()[0],
242                                         2.0e-8);
243 
244                     // the real state used for estimation is adjusted according to downlink delay
245                     double adjustment = state.getDate().durationFrom(estimated.getStates()[0].getDate());
246                     Assert.assertTrue(adjustment > 0.006);
247                     Assert.assertTrue(adjustment < 0.011);
248 
249                     final double RangeEstimated = estimated.getEstimatedValue()[0];
250                     final double absoluteError = RangeEstimated-RangeObserved;
251                     absoluteErrors.add(absoluteError);
252                     relativeErrors.add(FastMath.abs(absoluteError)/FastMath.abs(RangeObserved));
253 
254                     // Print results on console ?
255                     if (printResults) {
256                         final AbsoluteDate measurementDate = measurement.getDate();
257                         String stationName = ((Range) measurement).getStation().getBaseFrame().getName();
258 
259                         System.out.format(Locale.US, "%-15s  %-23s  %-23s  %19.6f  %19.6f  %13.6e  %13.6e%n",
260                                          stationName, measurementDate, date,
261                                          RangeObserved, RangeEstimated,
262                                          FastMath.abs(RangeEstimated-RangeObserved),
263                                          FastMath.abs((RangeEstimated-RangeObserved)/RangeObserved));
264                     }
265 
266                 } // End if measurement date between previous and current interpolator step
267             } // End for loop on the measurements
268         }); // End lambda function handlestep
269 
270         // Print results on console ? Header
271         if (printResults) {
272             System.out.format(Locale.US, "%-15s  %-23s  %-23s  %19s  %19s  %13s  %13s%n",
273                               "Station", "Measurement Date", "State Date",
274                               "Range observed [m]", "Range estimated [m]",
275                               "ΔRange [m]", "rel ΔRange");
276         }
277 
278         // Rewind the propagator to initial date
279         propagator.propagate(context.initialOrbit.getDate());
280 
281         // Sort measurements chronologically
282         measurements.sort(Comparator.naturalOrder());
283 
284         // Propagate to final measurement's date
285         propagator.propagate(measurements.get(measurements.size()-1).getDate());
286 
287         // Convert lists to double array
288         final double[] absErrors = absoluteErrors.stream().mapToDouble(Double::doubleValue).toArray();
289         final double[] relErrors = relativeErrors.stream().mapToDouble(Double::doubleValue).toArray();
290 
291         // Statistics' assertion
292         final double absErrorsMedian = new Median().evaluate(absErrors);
293         final double absErrorsMin    = new Min().evaluate(absErrors);
294         final double absErrorsMax    = new Max().evaluate(absErrors);
295         final double relErrorsMedian = new Median().evaluate(relErrors);
296         final double relErrorsMax    = new Max().evaluate(relErrors);
297 
298         // Print the results on console ? Final results
299         if (printResults) {
300             System.out.println();
301             System.out.println("Absolute errors median: " +  absErrorsMedian);
302             System.out.println("Absolute errors min   : " +  absErrorsMin);
303             System.out.println("Absolute errors max   : " +  absErrorsMax);
304             System.out.println("Relative errors median: " +  relErrorsMedian);
305             System.out.println("Relative errors max   : " +  relErrorsMax);
306         }
307 
308         Assert.assertEquals(0.0, absErrorsMedian, 4.9e-8);
309         Assert.assertEquals(0.0, absErrorsMin,    2.2e-7);
310         Assert.assertEquals(0.0, absErrorsMax,    2.1e-7);
311         Assert.assertEquals(0.0, relErrorsMedian, 1.0e-14);
312         Assert.assertEquals(0.0, relErrorsMax,    2.6e-14);
313 
314 
315     }
316 
317     void genericTestStateDerivatives(final boolean isModifier, final boolean printResults,
318                                      final double refErrorsPMedian, final double refErrorsPMean, final double refErrorsPMax,
319                                      final double refErrorsVMedian, final double refErrorsVMean, final double refErrorsVMax) {
320 
321         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
322 
323         final NumericalPropagatorBuilder propagatorBuilder =
324                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true,
325                                               1.0e-6, 60.0, 0.001);
326 
327         // Create perfect range measurements
328         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
329                                                                            propagatorBuilder);
330         final List<ObservedMeasurement<?>> measurements =
331                         EstimationTestUtils.createMeasurements(propagator,
332                                                                new RangeMeasurementCreator(context),
333                                                                1.0, 3.0, 300.0);
334 
335         // Lists for results' storage - Used only for derivatives with respect to state
336         // "final" value to be seen by "handleStep" function of the propagator
337         final List<Double> errorsP = new ArrayList<Double>();
338         final List<Double> errorsV = new ArrayList<Double>();
339 
340         // Set step handler
341         // Use a lambda function to implement "handleStep" function
342         propagator.setStepHandler(interpolator -> {
343 
344             for (final ObservedMeasurement<?> measurement : measurements) {
345 
346                 //  Play test if the measurement date is between interpolator previous and current date
347                 if ((measurement.getDate().durationFrom(interpolator.getPreviousState().getDate()) > 0.) &&
348                     (measurement.getDate().durationFrom(interpolator.getCurrentState().getDate())  <=  0.)
349                    ) {
350 
351                     // Add modifiers if test implies it
352                     final RangeTroposphericDelayModifier modifier = new RangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
353                     if (isModifier) {
354                         ((Range) measurement).addModifier(modifier);
355                     }
356 
357                     // We intentionally propagate to a date which is close to the
358                     // real spacecraft state but is *not* the accurate date, by
359                     // compensating only part of the downlink delay. This is done
360                     // in order to validate the partial derivatives with respect
361                     // to velocity. If we had chosen the proper state date, the
362                     // range would have depended only on the current position but
363                     // not on the current velocity.
364                     final double          meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
365                     final AbsoluteDate    date      = measurement.getDate().shiftedBy(-0.75 * meanDelay);
366                     final SpacecraftState state     = interpolator.getInterpolatedState(date);
367                     final double[][]      jacobian  = measurement.estimate(0, 0, new SpacecraftState[] { state }).getStateDerivatives(0);
368 
369                     // Jacobian reference value
370                     final double[][] jacobianRef;
371 
372                     // Compute a reference value using finite differences
373                     jacobianRef = Differentiation.differentiate(new StateFunction() {
374                         public double[] value(final SpacecraftState state) {
375                             return measurement.estimate(0, 0, new SpacecraftState[] { state }).getEstimatedValue();
376                         }
377                     }, measurement.getDimension(), propagator.getAttitudeProvider(),
378                        OrbitType.CARTESIAN, PositionAngle.TRUE, 2.0, 3).value(state);
379 
380                     Assert.assertEquals(jacobianRef.length, jacobian.length);
381                     Assert.assertEquals(jacobianRef[0].length, jacobian[0].length);
382 
383                     // Errors & relative errors on the Jacobian
384                     double [][] dJacobian         = new double[jacobian.length][jacobian[0].length];
385                     double [][] dJacobianRelative = new double[jacobian.length][jacobian[0].length];
386                     for (int i = 0; i < jacobian.length; ++i) {
387                         for (int j = 0; j < jacobian[i].length; ++j) {
388                             dJacobian[i][j] = jacobian[i][j] - jacobianRef[i][j];
389                             dJacobianRelative[i][j] = FastMath.abs(dJacobian[i][j]/jacobianRef[i][j]);
390 
391                             if (j < 3) { errorsP.add(dJacobianRelative[i][j]);
392                             } else { errorsV.add(dJacobianRelative[i][j]); }
393                         }
394                     }
395                     // Print values in console ?
396                     if (printResults) {
397                         String stationName  = ((Range) measurement).getStation().getBaseFrame().getName();
398                         System.out.format(Locale.US, "%-15s  %-23s  %-23s  " +
399                                         "%10.3e  %10.3e  %10.3e  " +
400                                         "%10.3e  %10.3e  %10.3e  " +
401                                         "%10.3e  %10.3e  %10.3e  " +
402                                         "%10.3e  %10.3e  %10.3e%n",
403                                         stationName, measurement.getDate(), date,
404                                         dJacobian[0][0], dJacobian[0][1], dJacobian[0][2],
405                                         dJacobian[0][3], dJacobian[0][4], dJacobian[0][5],
406                                         dJacobianRelative[0][0], dJacobianRelative[0][1], dJacobianRelative[0][2],
407                                         dJacobianRelative[0][3], dJacobianRelative[0][4], dJacobianRelative[0][5]);
408                     }
409                 } // End if measurement date between previous and current interpolator step
410             } // End for loop on the measurements
411         });
412 
413         // Print results on console ?
414         if (printResults) {
415             System.out.format(Locale.US, "%-15s  %-23s  %-23s  " +
416                             "%10s  %10s  %10s  " +
417                             "%10s  %10s  %10s  " +
418                             "%10s  %10s  %10s  " +
419                             "%10s  %10s  %10s%n",
420                             "Station", "Measurement Date", "State Date",
421                             "ΔdPx", "ΔdPy", "ΔdPz", "ΔdVx", "ΔdVy", "ΔdVz",
422                             "rel ΔdPx", "rel ΔdPy", "rel ΔdPz",
423                             "rel ΔdVx", "rel ΔdVy", "rel ΔdVz");
424         }
425 
426         // Rewind the propagator to initial date
427         propagator.propagate(context.initialOrbit.getDate());
428 
429         // Sort measurements chronologically
430         measurements.sort(Comparator.naturalOrder());
431 
432         // Propagate to final measurement's date
433         propagator.propagate(measurements.get(measurements.size()-1).getDate());
434 
435         // Convert lists to double[] and evaluate some statistics
436         final double relErrorsP[] = errorsP.stream().mapToDouble(Double::doubleValue).toArray();
437         final double relErrorsV[] = errorsV.stream().mapToDouble(Double::doubleValue).toArray();
438 
439         final double errorsPMedian = new Median().evaluate(relErrorsP);
440         final double errorsPMean   = new Mean().evaluate(relErrorsP);
441         final double errorsPMax    = new Max().evaluate(relErrorsP);
442         final double errorsVMedian = new Median().evaluate(relErrorsV);
443         final double errorsVMean   = new Mean().evaluate(relErrorsV);
444         final double errorsVMax    = new Max().evaluate(relErrorsV);
445 
446         // Print the results on console ?
447         if (printResults) {
448             System.out.println();
449             System.out.format(Locale.US, "Relative errors dR/dP -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n",
450                               errorsPMedian, errorsPMean, errorsPMax);
451             System.out.format(Locale.US, "Relative errors dR/dV -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n",
452                               errorsVMedian, errorsVMean, errorsVMax);
453         }
454 
455         Assert.assertEquals(0.0, errorsPMedian, refErrorsPMedian);
456         Assert.assertEquals(0.0, errorsPMean, refErrorsPMean);
457         Assert.assertEquals(0.0, errorsPMax, refErrorsPMax);
458         Assert.assertEquals(0.0, errorsVMedian, refErrorsVMedian);
459         Assert.assertEquals(0.0, errorsVMean, refErrorsVMean);
460         Assert.assertEquals(0.0, errorsVMax, refErrorsVMax);
461     }
462 
463     void genericTestParameterDerivatives(final boolean isModifier, final boolean printResults,
464                                          final double refErrorsMedian, final double refErrorsMean, final double refErrorsMax) {
465 
466         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
467 
468         final NumericalPropagatorBuilder propagatorBuilder =
469                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true,
470                                               1.0e-6, 60.0, 0.001);
471 
472         // Create perfect range measurements
473         for (final GroundStation station : context.stations) {
474             station.getClockOffsetDriver().setSelected(true);
475             station.getEastOffsetDriver().setSelected(true);
476             station.getNorthOffsetDriver().setSelected(true);
477             station.getZenithOffsetDriver().setSelected(true);
478         }
479         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
480                                                                            propagatorBuilder);
481         final List<ObservedMeasurement<?>> measurements =
482                         EstimationTestUtils.createMeasurements(propagator,
483                                                                new RangeMeasurementCreator(context),
484                                                                1.0, 3.0, 300.0);
485 
486         // List to store the results
487         final List<Double> relErrorList = new ArrayList<Double>();
488 
489         // Set step handler
490         // Use a lambda function to implement "handleStep" function
491         propagator.setStepHandler(interpolator -> {
492 
493             for (final ObservedMeasurement<?> measurement : measurements) {
494 
495                 //  Play test if the measurement date is between interpolator previous and current date
496                 if ((measurement.getDate().durationFrom(interpolator.getPreviousState().getDate()) > 0.) &&
497                     (measurement.getDate().durationFrom(interpolator.getCurrentState().getDate())  <=  0.)
498                    ) {
499 
500                     // Add modifiers if test implies it
501                     final RangeTroposphericDelayModifier modifier = new RangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
502                     if (isModifier) {
503                         ((Range) measurement).addModifier(modifier);
504                     }
505 
506                     // Parameter corresponding to station position offset
507                     final GroundStation station = ((Range) measurement).getStation();
508 
509                     // We intentionally propagate to a date which is close to the
510                     // real spacecraft state but is *not* the accurate date, by
511                     // compensating only part of the downlink delay. This is done
512                     // in order to validate the partial derivatives with respect
513                     // to velocity. If we had chosen the proper state date, the
514                     // range would have depended only on the current position but
515                     // not on the current velocity.
516                     final double          meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
517                     final AbsoluteDate    date      = measurement.getDate().shiftedBy(-0.75 * meanDelay);
518                     final SpacecraftState state     = interpolator.getInterpolatedState(date);
519                     final ParameterDriver[] drivers = new ParameterDriver[] {
520                         station.getClockOffsetDriver(),
521                         station.getEastOffsetDriver(),
522                         station.getNorthOffsetDriver(),
523                         station.getZenithOffsetDriver()
524                     };
525 
526                     if (printResults) {
527                         String stationName  = ((Range) measurement).getStation().getBaseFrame().getName();
528                         System.out.format(Locale.US, "%-15s  %-23s  %-23s  ",
529                                           stationName, measurement.getDate(), date);
530                     }
531 
532                     for (int i = 0; i < drivers.length; ++i) {
533                         final double[] gradient  = measurement.estimate(0, 0, new SpacecraftState[] { state }).getParameterDerivatives(drivers[i]);
534                         Assert.assertEquals(1, measurement.getDimension());
535                         Assert.assertEquals(1, gradient.length);
536 
537                         // Compute a reference value using finite differences
538                         final ParameterFunction dMkdP =
539                                         Differentiation.differentiate(new ParameterFunction() {
540                                             /** {@inheritDoc} */
541                                             @Override
542                                             public double value(final ParameterDriver parameterDriver) {
543                                                 return measurement.estimate(0, 0, new SpacecraftState[] { state }).getEstimatedValue()[0];
544                                             }
545                                         }, 3, 20.0 * drivers[i].getScale());
546                         final double ref = dMkdP.value(drivers[i]);
547 
548                         if (printResults) {
549                             System.out.format(Locale.US, "%10.3e  %10.3e  ", gradient[0]-ref, FastMath.abs((gradient[0]-ref)/ref));
550                         }
551 
552                         final double relError = FastMath.abs((ref-gradient[0])/ref);
553                         relErrorList.add(relError);
554                         Assert.assertEquals(ref, gradient[0], 6.1e-5 * FastMath.abs(ref));
555                     }
556                     if (printResults) {
557                         System.out.format(Locale.US, "%n");
558                     }
559 
560                 } // End if measurement date between previous and current interpolator step
561             } // End for loop on the measurements
562         });
563 
564         // Rewind the propagator to initial date
565         propagator.propagate(context.initialOrbit.getDate());
566 
567         // Sort measurements chronologically
568         measurements.sort(Comparator.naturalOrder());
569 
570         // Print results ? Header
571         if (printResults) {
572             System.out.format(Locale.US, "%-15s  %-23s  %-23s  " +
573                               "%10s  %10s  %10s  %10s  %10s  %10s  %10s  %10s%n",
574                               "Station", "Measurement Date", "State Date",
575                               "Δt",   "rel Δt",
576                               "ΔdQx", "rel ΔdQx",
577                               "ΔdQy", "rel ΔdQy",
578                               "ΔdQz", "rel ΔdQz");
579          }
580 
581         // Propagate to final measurement's date
582         propagator.propagate(measurements.get(measurements.size()-1).getDate());
583 
584         // Convert error list to double[]
585         final double relErrors[] = relErrorList.stream().mapToDouble(Double::doubleValue).toArray();
586 
587         // Compute statistics
588         final double relErrorsMedian = new Median().evaluate(relErrors);
589         final double relErrorsMean   = new Mean().evaluate(relErrors);
590         final double relErrorsMax    = new Max().evaluate(relErrors);
591 
592         // Print the results on console ?
593         if (printResults) {
594             System.out.println();
595             System.out.format(Locale.US, "Relative errors dR/dQ -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n",
596                               relErrorsMedian, relErrorsMean, relErrorsMax);
597         }
598 
599         Assert.assertEquals(0.0, relErrorsMedian, refErrorsMedian);
600         Assert.assertEquals(0.0, relErrorsMean, refErrorsMean);
601         Assert.assertEquals(0.0, relErrorsMax, refErrorsMax);
602 
603     }
604 
605     void genericTestEstimatedParameterDerivatives(final boolean isModifier, final boolean printResults,
606                                                   final double refErrorsMedian, final double refErrorsMean, final double refErrorsMax)
607                     {
608 
609         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
610 
611         final NumericalPropagatorBuilder propagatorBuilder =
612                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true,
613                                               1.0e-6, 60.0, 0.001);
614 
615         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
616                                                                            propagatorBuilder);
617         final List<ObservedMeasurement<?>> measurements =
618                         EstimationTestUtils.createMeasurements(propagator,
619                                                                new RangeMeasurementCreator(context),
620                                                                1.0, 3.0, 300.0);
621 
622         // List to store the results
623         final List<Double> relErrorList = new ArrayList<Double>();
624 
625         // Set step handler
626         // Use a lambda function to implement "handleStep" function
627         propagator.setStepHandler(interpolator -> {
628 
629             for (final ObservedMeasurement<?> measurement : measurements) {
630 
631                 //  Play test if the measurement date is between interpolator previous and current date
632                 if ((measurement.getDate().durationFrom(interpolator.getPreviousState().getDate()) > 0.) &&
633                     (measurement.getDate().durationFrom(interpolator.getCurrentState().getDate())  <=  0.)
634                    ) {
635 
636                     String stationName  = ((Range) measurement).getStation().getBaseFrame().getName();
637 
638                     // Add modifiers if test implies it
639                     final NiellMappingFunctionModel mappingFunction = new NiellMappingFunctionModel();
640                     final EstimatedTroposphericModel tropoModel     = new EstimatedTroposphericModel(mappingFunction, 5.0);
641                     
642                     final List<ParameterDriver> parameters = tropoModel.getParametersDrivers();
643                     for (ParameterDriver driver : parameters) {
644                         driver.setSelected(true);
645                     }
646 
647                     parameters.get(0).setName(stationName + "/" + EstimatedTroposphericModel.TOTAL_ZENITH_DELAY);
648                     final RangeTroposphericDelayModifier modifier = new RangeTroposphericDelayModifier(tropoModel);
649                     if (isModifier) {
650                         ((Range) measurement).addModifier(modifier);
651                     }
652 
653                     // We intentionally propagate to a date which is close to the
654                     // real spacecraft state but is *not* the accurate date, by
655                     // compensating only part of the downlink delay. This is done
656                     // in order to validate the partial derivatives with respect
657                     // to velocity. If we had chosen the proper state date, the
658                     // range would have depended only on the current position but
659                     // not on the current velocity.
660                     final double          meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
661                     final AbsoluteDate    date      = measurement.getDate().shiftedBy(-0.75 * meanDelay);
662                     final SpacecraftState state     = interpolator.getInterpolatedState(date);
663                     final ParameterDriver[] drivers = new ParameterDriver[] {
664                         parameters.get(0)
665                     };
666 
667                     if (printResults) {
668                         System.out.format(Locale.US, "%-15s  %-23s  %-23s  ",
669                                           stationName, measurement.getDate(), date);
670                     }
671 
672                     for (int i = 0; i < 1; ++i) {
673                         final double[] gradient  = measurement.estimate(0, 0, new SpacecraftState[] { state }).getParameterDerivatives(drivers[i]);
674                         Assert.assertEquals(1, measurement.getDimension());
675                         Assert.assertEquals(1, gradient.length);
676 
677                         // Compute a reference value using finite differences
678                         final ParameterFunction dMkdP =
679                                         Differentiation.differentiate(new ParameterFunction() {
680                                             /** {@inheritDoc} */
681                                             @Override
682                                             public double value(final ParameterDriver parameterDriver) {
683                                                 return measurement.estimate(0, 0, new SpacecraftState[] { state }).getEstimatedValue()[0];
684                                             }
685                                         }, 3, 0.1 * drivers[i].getScale());
686                         final double ref = dMkdP.value(drivers[i]);
687 
688                         if (printResults) {
689                             System.out.format(Locale.US, "%10.3e  %10.3e  ", gradient[0]-ref, FastMath.abs((gradient[0]-ref)/ref));
690                         }
691 
692                         final double relError = FastMath.abs((ref-gradient[0])/ref);
693                         relErrorList.add(relError);
694 //                        Assert.assertEquals(ref, gradient[0], 6.1e-5 * FastMath.abs(ref));
695                     }
696                     if (printResults) {
697                         System.out.format(Locale.US, "%n");
698                     }
699 
700                 } // End if measurement date between previous and current interpolator step
701             } // End for loop on the measurements
702         });
703 
704         // Rewind the propagator to initial date
705         propagator.propagate(context.initialOrbit.getDate());
706 
707         // Sort measurements chronologically
708         measurements.sort(Comparator.naturalOrder());
709 
710         // Print results ? Header
711         if (printResults) {
712             System.out.format(Locale.US, "%-15s  %-23s  %-23s  " +
713                             "%10s  %10s  %10s  " +
714                             "%10s  %10s  %10s%n ",
715                             "Station", "Measurement Date", "State Date",
716                             "ΔdDelay", "rel ΔdDelay",
717                             "ΔdNorthG", "rel ΔdNorthG",
718                             "ΔdEastG", "rel ΔdEastG");
719          }
720 
721         // Propagate to final measurement's date
722         propagator.propagate(measurements.get(measurements.size()-1).getDate());
723 
724         // Convert error list to double[]
725         final double relErrors[] = relErrorList.stream().mapToDouble(Double::doubleValue).toArray();
726 
727         // Compute statistics
728         final double relErrorsMedian = new Median().evaluate(relErrors);
729         final double relErrorsMean   = new Mean().evaluate(relErrors);
730         final double relErrorsMax    = new Max().evaluate(relErrors);
731 
732         // Print the results on console ?
733         if (printResults) {
734             System.out.println();
735             System.out.format(Locale.US, "Relative errors dR/dQ -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n",
736                               relErrorsMedian, relErrorsMean, relErrorsMax);
737         }
738 
739         Assert.assertEquals(0.0, relErrorsMedian, refErrorsMedian);
740         Assert.assertEquals(0.0, relErrorsMean, refErrorsMean);
741         Assert.assertEquals(0.0, relErrorsMax, refErrorsMax);
742 
743     }
744 
745 }