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.List;
20  
21  import org.hipparchus.stat.descriptive.StreamingStatistics;
22  import org.hipparchus.util.FastMath;
23  import org.junit.jupiter.api.Assertions;
24  import org.junit.jupiter.api.Test;
25  import org.orekit.estimation.Context;
26  import org.orekit.estimation.EstimationTestUtils;
27  import org.orekit.estimation.measurements.modifiers.RangeRateTroposphericDelayModifier;
28  import org.orekit.models.earth.troposphere.EstimatedModel;
29  import org.orekit.models.earth.troposphere.GlobalMappingFunctionModel;
30  import org.orekit.models.earth.troposphere.ModifiedSaastamoinenModel;
31  import org.orekit.orbits.OrbitType;
32  import org.orekit.orbits.PositionAngleType;
33  import org.orekit.propagation.Propagator;
34  import org.orekit.propagation.SpacecraftState;
35  import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
36  import org.orekit.time.AbsoluteDate;
37  import org.orekit.utils.Constants;
38  import org.orekit.utils.Differentiation;
39  import org.orekit.utils.ParameterDriver;
40  import org.orekit.utils.ParameterFunction;
41  
42  public class RangeRateTest {
43  
44      /** Compare observed values and estimated values.
45       *  Both are calculated with a different algorithm.
46       *  One-way measurements.
47       */
48      @Test
49      public void testValuesOneWay() {
50  
51          Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
52  
53          final NumericalPropagatorBuilder propagatorBuilder =
54                          context.createBuilder(OrbitType.EQUINOCTIAL, PositionAngleType.TRUE, false,
55                                                1.0e-6, 60.0, 0.001);
56  
57          // Create perfect right-ascension/declination measurements
58          final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
59                                                                             propagatorBuilder);
60          final double satClkDrift = 3.2e-10;
61          final List<ObservedMeasurement<?>> measurements =
62                          EstimationTestUtils.createMeasurements(propagator,
63                                                                 new RangeRateMeasurementCreator(context, false, satClkDrift),
64                                                                 1.0, 3.0, 300.0);
65  
66          propagator.clearStepHandlers();
67  
68          // Prepare statistics for values difference
69          final StreamingStatistics diffStat = new StreamingStatistics();
70  
71          for (final ObservedMeasurement<?> measurement : measurements) {
72  
73              // Propagate to measurement date
74              final AbsoluteDate datemeas  = measurement.getDate();
75              SpacecraftState    state     = propagator.propagate(datemeas);
76  
77              // Estimate the AZEL value
78              final EstimatedMeasurementBase<?> estimated = measurement.estimateWithoutDerivatives(new SpacecraftState[] { state });
79  
80              // Store the difference between estimated and observed values in the stats
81              diffStat.addValue(FastMath.abs(estimated.getEstimatedValue()[0] - measurement.getObservedValue()[0]));
82          }
83  
84          // Mean and std errors check
85          Assertions.assertEquals(0.0, diffStat.getMean(), 6.5e-8);
86          Assertions.assertEquals(0.0, diffStat.getStandardDeviation(), 5.5e-8);
87  
88          // Test measurement type
89          Assertions.assertEquals(RangeRate.MEASUREMENT_TYPE, measurements.get(0).getMeasurementType());
90      }
91  
92      /** Compare observed values and estimated values.
93       *  Both are calculated with a different algorithm.
94       *  Two-ways measurements.
95       */
96      @Test
97      public void testValuesTwoWays() {
98  
99          Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
100 
101         final NumericalPropagatorBuilder propagatorBuilder =
102                         context.createBuilder(OrbitType.EQUINOCTIAL, PositionAngleType.TRUE, false,
103                                               1.0e-6, 60.0, 0.001);
104 
105         // Create perfect right-ascension/declination measurements
106         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
107                                                                            propagatorBuilder);
108         final double satClkDrift = 3.2e-10;
109         final List<ObservedMeasurement<?>> measurements =
110                         EstimationTestUtils.createMeasurements(propagator,
111                                                                new RangeRateMeasurementCreator(context, true, satClkDrift),
112                                                                1.0, 3.0, 300.0);
113 
114         propagator.clearStepHandlers();
115 
116         // Prepare statistics for values difference
117         final StreamingStatistics diffStat = new StreamingStatistics();
118 
119         for (final ObservedMeasurement<?> measurement : measurements) {
120 
121             // Propagate to measurement date
122             final AbsoluteDate datemeas  = measurement.getDate();
123             SpacecraftState    state     = propagator.propagate(datemeas);
124 
125             // Estimate the AZEL value
126             final EstimatedMeasurementBase<?> estimated = measurement.estimateWithoutDerivatives(new SpacecraftState[] { state });
127 
128             // Store the difference between estimated and observed values in the stats
129             diffStat.addValue(FastMath.abs(estimated.getEstimatedValue()[0] - measurement.getObservedValue()[0]));
130         }
131 
132         // Mean and std errors check
133         Assertions.assertEquals(0.0, diffStat.getMean(), 6.5e-8);
134         Assertions.assertEquals(0.0, diffStat.getStandardDeviation(), 5.5e-8);
135     }
136 
137     /** Test the values of the state derivatives using a numerical
138      * finite differences calculation as a reference.
139      * One way measurements.
140      */
141     @Test
142     public void testStateDerivativesOneWay() {
143 
144         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
145 
146         final NumericalPropagatorBuilder propagatorBuilder =
147                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
148                                               1.0e-6, 60.0, 0.001);
149 
150         // create perfect range rate measurements
151         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
152                                                                            propagatorBuilder);
153         final double satClkDrift = 3.2e-10;
154         final List<ObservedMeasurement<?>> measurements =
155                         EstimationTestUtils.createMeasurements(propagator,
156                                                                new RangeRateMeasurementCreator(context, false, satClkDrift),
157                                                                1.0, 3.0, 300.0);
158         for (final ObservedMeasurement<?> m : measurements) {
159             Assertions.assertFalse(((RangeRate) m).isTwoWay());
160         }
161         propagator.clearStepHandlers();
162 
163         double maxRelativeError = 0;
164         for (final ObservedMeasurement<?> measurement : measurements) {
165 
166             final double          meanDelay = 1; // measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
167             final AbsoluteDate    date      = measurement.getDate().shiftedBy(-0.75 * meanDelay);
168             final SpacecraftState state = propagator.propagate(date);
169 
170             final EstimatedMeasurement<?> estimated = measurement.estimate(0, 0, new SpacecraftState[] { state });
171             Assertions.assertEquals(2, estimated.getParticipants().length);
172             final double[][] jacobian = estimated.getStateDerivatives(0);
173 
174             final double[][] finiteDifferencesJacobian =
175                     Differentiation.differentiate(state1 ->
176                                                       measurement.
177                                                           estimateWithoutDerivatives(new SpacecraftState[] { state1 }).
178                                                           getEstimatedValue(),
179                                                   1, propagator.getAttitudeProvider(),
180                OrbitType.CARTESIAN, PositionAngleType.TRUE, 15.0, 3).value(state);
181 
182             Assertions.assertEquals(finiteDifferencesJacobian.length, jacobian.length);
183             Assertions.assertEquals(finiteDifferencesJacobian[0].length, jacobian[0].length);
184 
185             for (int i = 0; i < jacobian.length; ++i) {
186                 for (int j = 0; j < jacobian[i].length; ++j) {
187                     // check the values returned by getStateDerivatives() are correct
188                     maxRelativeError = FastMath.max(maxRelativeError,
189                                                     FastMath.abs((finiteDifferencesJacobian[i][j] - jacobian[i][j]) /
190                                                                  finiteDifferencesJacobian[i][j]));
191                 }
192             }
193 
194         }
195         Assertions.assertEquals(0, maxRelativeError, 8.1e-6);
196 
197     }
198 
199     /** Test the values of the state derivatives using a numerical
200      * finite differences calculation as a reference.
201      * Two-ways measurements.
202      */
203     @Test
204     public void testStateDerivativesTwoWays() {
205 
206         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
207 
208         final NumericalPropagatorBuilder propagatorBuilder =
209                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
210                                               1.0e-6, 60.0, 0.001);
211 
212         // create perfect range rate measurements
213         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
214                                                                            propagatorBuilder);
215         final double satClkDrift = 3.2e-10;
216         final List<ObservedMeasurement<?>> measurements =
217                         EstimationTestUtils.createMeasurements(propagator,
218                                                                new RangeRateMeasurementCreator(context, true, satClkDrift),
219                                                                1.0, 3.0, 300.0);
220         for (final ObservedMeasurement<?> m : measurements) {
221             Assertions.assertTrue(((RangeRate) m).isTwoWay());
222         }
223         propagator.clearStepHandlers();
224 
225         double maxRelativeError = 0;
226         for (final ObservedMeasurement<?> measurement : measurements) {
227 
228             //
229             //final AbsoluteDate date = measurement.getDate();
230             final double          meanDelay = 1; // measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
231             final AbsoluteDate    date      = measurement.getDate().shiftedBy(-0.75 * meanDelay);
232             final SpacecraftState state     = propagator.propagate(date);
233 
234             final EstimatedMeasurement<?> estimated = measurement.estimate(0, 0, new SpacecraftState[] { state });
235             Assertions.assertEquals(3, estimated.getParticipants().length);
236             final double[][] jacobian = estimated.getStateDerivatives(0);
237 
238             final double[][] finiteDifferencesJacobian =
239                     Differentiation.differentiate(state1 ->
240                                                       measurement.
241                                                           estimateWithoutDerivatives(new SpacecraftState[] { state1 }).
242                                                           getEstimatedValue(),
243                                                   1, propagator.getAttitudeProvider(),
244                OrbitType.CARTESIAN, PositionAngleType.TRUE, 15.0, 3).value(state);
245 
246             Assertions.assertEquals(finiteDifferencesJacobian.length, jacobian.length);
247             Assertions.assertEquals(finiteDifferencesJacobian[0].length, jacobian[0].length);
248 
249             for (int i = 0; i < jacobian.length; ++i) {
250                 for (int j = 0; j < jacobian[i].length; ++j) {
251                     // check the values returned by getStateDerivatives() are correct
252                     maxRelativeError = FastMath.max(maxRelativeError,
253                                                     FastMath.abs((finiteDifferencesJacobian[i][j] - jacobian[i][j]) /
254                                                                  finiteDifferencesJacobian[i][j]));
255                 }
256             }
257 
258         }
259         Assertions.assertEquals(0, maxRelativeError, 8.1e-6);
260 
261     }
262 
263     /** Test the values of the parameters' derivatives using a numerical
264      * finite differences calculation as a reference.
265      * One-way measurements.
266      */
267     @Test
268     public void testParameterDerivativesOneWay() {
269 
270         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
271 
272         final NumericalPropagatorBuilder propagatorBuilder =
273                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
274                                               1.0e-6, 60.0, 0.001);
275 
276         // create perfect range rate measurements
277         final double groundClockDrift =  4.8e-9;
278         for (final GroundStation station : context.stations) {
279             station.getClockDriftDriver().setValue(groundClockDrift);
280         }
281         final double satClkDrift = 3.2e-10;
282         final RangeRateMeasurementCreator creator = new RangeRateMeasurementCreator(context, false, satClkDrift);
283         creator.getSatellite().getClockDriftDriver().setSelected(true);
284         for (final GroundStation station : context.stations) {
285             station.getClockOffsetDriver().setSelected(true);
286             station.getClockDriftDriver().setSelected(true);
287             station.getEastOffsetDriver().setSelected(true);
288             station.getNorthOffsetDriver().setSelected(true);
289             station.getZenithOffsetDriver().setSelected(true);
290         }
291         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
292                                                                            propagatorBuilder);
293 
294 
295 
296         final List<ObservedMeasurement<?>> measurements =
297                         EstimationTestUtils.createMeasurements(propagator,
298                                                                creator,
299                                                                1.0, 3.0, 300.0);
300         propagator.clearStepHandlers();
301 
302         double maxRelativeError = 0;
303         for (final ObservedMeasurement<?> measurement : measurements) {
304 
305             // parameter corresponding to station position offset
306             final GroundStation stationParameter = ((RangeRate) measurement).getStation();
307 
308             // We intentionally propagate to a date which is close to the
309             // real spacecraft state but is *not* the accurate date, by
310             // compensating only part of the downlink delay. This is done
311             // in order to validate the partial derivatives with respect
312             // to velocity. If we had chosen the proper state date, the
313             // range would have depended only on the current position but
314             // not on the current velocity.
315             final double          meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
316             final AbsoluteDate    date      = measurement.getDate().shiftedBy(-0.75 * meanDelay);
317             final SpacecraftState state     = propagator.propagate(date);
318             final ParameterDriver[] drivers = new ParameterDriver[] {
319                 stationParameter.getClockDriftDriver(),
320                 stationParameter.getEastOffsetDriver(),
321                 stationParameter.getNorthOffsetDriver(),
322                 stationParameter.getZenithOffsetDriver(),
323                 measurement.getSatellites().get(0).getClockDriftDriver()
324             };
325             for (int i = 0; i < drivers.length; ++i) {
326                 final double[] gradient  = measurement.estimate(0, 0, new SpacecraftState[] { state }).getParameterDerivatives(drivers[i]);
327                 Assertions.assertEquals(1, measurement.getDimension());
328                 Assertions.assertEquals(1, gradient.length);
329 
330                 final ParameterFunction dMkdP =
331                                 Differentiation.differentiate(new ParameterFunction() {
332                                     /** {@inheritDoc} */
333                                     @Override
334                                     public double value(final ParameterDriver parameterDriver, AbsoluteDate date) {
335                                         return measurement.
336                                                estimateWithoutDerivatives(new SpacecraftState[] { state }).
337                                                getEstimatedValue()[0];
338                                     }
339                                 }, 3, 20.0 * drivers[i].getScale());
340                 final double ref = dMkdP.value(drivers[i], date);
341                 maxRelativeError = FastMath.max(maxRelativeError, FastMath.abs((ref - gradient[0]) / ref));
342             }
343 
344         }
345         Assertions.assertEquals(0, maxRelativeError, 1.2e-6);
346 
347     }
348 
349     /** Test the values of the parameters' derivatives using a numerical
350      * finite differences calculation as a reference.
351      * Two-ways measurements.
352      */
353     @Test
354     public void testParameterDerivativesTwoWays() {
355 
356         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
357 
358         final NumericalPropagatorBuilder propagatorBuilder =
359                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
360                                               1.0e-6, 60.0, 0.001);
361 
362         // create perfect range rate measurements
363         final double groundClockDrift =  4.8e-9;
364         for (final GroundStation station : context.stations) {
365             station.getClockDriftDriver().setValue(groundClockDrift);
366         }
367         final double satClkDrift = 3.2e-10;
368         final RangeRateMeasurementCreator creator = new RangeRateMeasurementCreator(context, false, satClkDrift);
369         for (final GroundStation station : context.stations) {
370             station.getClockOffsetDriver().setSelected(true);
371             station.getEastOffsetDriver().setSelected(true);
372             station.getNorthOffsetDriver().setSelected(true);
373             station.getZenithOffsetDriver().setSelected(true);
374         }
375 
376         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
377                                                                            propagatorBuilder);
378 
379 
380         final List<ObservedMeasurement<?>> measurements =
381                         EstimationTestUtils.createMeasurements(propagator,
382                                                                creator,
383                                                                1.0, 3.0, 300.0);
384         propagator.clearStepHandlers();
385 
386         double maxRelativeError = 0;
387         for (final ObservedMeasurement<?> measurement : measurements) {
388 
389             // parameter corresponding to station position offset
390             final GroundStation stationParameter = ((RangeRate) measurement).getStation();
391 
392             // We intentionally propagate to a date which is close to the
393             // real spacecraft state but is *not* the accurate date, by
394             // compensating only part of the downlink delay. This is done
395             // in order to validate the partial derivatives with respect
396             // to velocity. If we had chosen the proper state date, the
397             // range would have depended only on the current position but
398             // not on the current velocity.
399             final double          meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
400             final AbsoluteDate    date      = measurement.getDate().shiftedBy(-0.75 * meanDelay);
401             final SpacecraftState state     = propagator.propagate(date);
402             final ParameterDriver[] drivers = new ParameterDriver[] {
403                 stationParameter.getEastOffsetDriver(),
404                 stationParameter.getNorthOffsetDriver(),
405                 stationParameter.getZenithOffsetDriver(),
406             };
407             for (int i = 0; i < drivers.length; ++i) {
408                 final double[] gradient  = measurement.estimate(0, 0, new SpacecraftState[] { state }).getParameterDerivatives(drivers[i]);
409                 Assertions.assertEquals(1, measurement.getDimension());
410                 Assertions.assertEquals(1, gradient.length);
411 
412                 final ParameterFunction dMkdP =
413                                 Differentiation.differentiate(new ParameterFunction() {
414                                     /** {@inheritDoc} */
415                                     @Override
416                                     public double value(final ParameterDriver parameterDriver, AbsoluteDate date) {
417                                         return measurement.
418                                                estimateWithoutDerivatives(new SpacecraftState[] { state }).
419                                                getEstimatedValue()[0];
420                                     }
421                                 }, 3, 20.0 * drivers[i].getScale());
422                 final double ref = dMkdP.value(drivers[i], date);
423                 maxRelativeError = FastMath.max(maxRelativeError, FastMath.abs((ref - gradient[0]) / ref));
424             }
425 
426         }
427         Assertions.assertEquals(0, maxRelativeError, 1.2e-6);
428 
429     }
430 
431     /** Test the values of the state derivatives using a numerical
432      * finite differences calculation as a reference.
433      * One-way measurements with modifiers (tropospheric corrections).
434      */
435     @Test
436     public void testStateDerivativesWithModifier() {
437 
438         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
439 
440         final NumericalPropagatorBuilder propagatorBuilder =
441                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
442                                               1.0e-6, 60.0, 0.001);
443 
444         // create perfect range rate measurements
445         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
446                                                                            propagatorBuilder);
447         final double groundClockDrift =  4.8e-9;
448         for (final GroundStation station : context.stations) {
449             station.getClockDriftDriver().setValue(groundClockDrift);
450         }
451         final double satClkDrift = 3.2e-10;
452         final List<ObservedMeasurement<?>> measurements =
453                         EstimationTestUtils.createMeasurements(propagator,
454                                                                new RangeRateMeasurementCreator(context, false, satClkDrift),
455                                                                1.0, 3.0, 300.0);
456         propagator.clearStepHandlers();
457 
458         double maxRelativeError = 0;
459         for (final ObservedMeasurement<?> measurement : measurements) {
460 
461             final RangeRateTroposphericDelayModifier modifier =
462                 new RangeRateTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel(), true);
463             ((RangeRate) measurement).addModifier(modifier);
464 
465             //
466             //final AbsoluteDate date = measurement.getDate();
467             final double          meanDelay = 1; // measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
468             final AbsoluteDate    date      = measurement.getDate().shiftedBy(-0.75 * meanDelay);
469             final SpacecraftState state = propagator.propagate(date);
470 
471             final double[][] jacobian = measurement.estimate(0, 0, new SpacecraftState[] { state }).getStateDerivatives(0);
472 
473             final double[][] finiteDifferencesJacobian =
474                     Differentiation.differentiate( state1 -> measurement.
475                             estimate(0, 0, new SpacecraftState[] { state1 }).
476                                                        getEstimatedValue(), 1, propagator.getAttitudeProvider(),
477                                                    OrbitType.CARTESIAN, PositionAngleType.TRUE, 15.0, 3).value(state);
478 
479             Assertions.assertEquals(finiteDifferencesJacobian.length, jacobian.length);
480             Assertions.assertEquals(finiteDifferencesJacobian[0].length, jacobian[0].length);
481 
482             for (int i = 0; i < jacobian.length; ++i) {
483                 for (int j = 0; j < jacobian[i].length; ++j) {
484                     // check the values returned by getStateDerivatives() are correct
485                     maxRelativeError = FastMath.max(maxRelativeError,
486                                                     FastMath.abs((finiteDifferencesJacobian[i][j] - jacobian[i][j]) /
487                                                                  finiteDifferencesJacobian[i][j]));
488                 }
489             }
490 
491         }
492         Assertions.assertEquals(0, maxRelativeError, 1.5e-7);
493 
494     }
495 
496     /** Test the values of the state derivatives using a numerical
497      * finite differences calculation as a reference.
498      * One-way measurements with estimated modifiers (tropospheric corrections).
499      */
500     @Test
501     public void testStateDerivativesWithEstimatedModifier() {
502 
503         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
504 
505         final NumericalPropagatorBuilder propagatorBuilder =
506                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
507                                               1.0e-6, 60.0, 0.001);
508 
509         // create perfect range rate measurements
510         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
511                                                                            propagatorBuilder);
512 
513         final double groundClockDrift =  4.8e-9;
514         for (final GroundStation station : context.stations) {
515             station.getClockDriftDriver().setValue(groundClockDrift);
516         }
517         final double satClkDrift = 3.2e-10;
518         final List<ObservedMeasurement<?>> measurements =
519                         EstimationTestUtils.createMeasurements(propagator,
520                                                                new RangeRateMeasurementCreator(context, false, satClkDrift),
521                                                                1.0, 3.0, 300.0);
522         propagator.clearStepHandlers();
523 
524         double maxRelativeError = 0;
525         for (final ObservedMeasurement<?> measurement : measurements) {
526 
527             // Add modifiers if test implies it
528             final GlobalMappingFunctionModel mappingFunction = new GlobalMappingFunctionModel();
529             final EstimatedModel             tropoModel      = new EstimatedModel(mappingFunction, 5.0);
530 
531             final RangeRateTroposphericDelayModifier modifier = new RangeRateTroposphericDelayModifier(tropoModel, true);
532             ((RangeRate) measurement).addModifier(modifier);
533 
534             //
535             //final AbsoluteDate date = measurement.getDate();
536             final double          meanDelay = 1; // measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
537             final AbsoluteDate    date      = measurement.getDate().shiftedBy(-0.75 * meanDelay);
538             final SpacecraftState state = propagator.propagate(date);
539 
540             final double[][] jacobian = measurement.estimate(0, 0, new SpacecraftState[] { state }).getStateDerivatives(0);
541 
542             final double[][] finiteDifferencesJacobian =
543                     Differentiation.differentiate(state1 -> measurement.
544                                                       estimate(0, 0, new SpacecraftState[] { state1 }).
545                                                       getEstimatedValue(), 1, propagator.getAttitudeProvider(),
546                         OrbitType.CARTESIAN, PositionAngleType.TRUE, 15.0, 3).value(state);
547 
548             Assertions.assertEquals(finiteDifferencesJacobian.length, jacobian.length);
549             Assertions.assertEquals(finiteDifferencesJacobian[0].length, jacobian[0].length);
550 
551             for (int i = 0; i < jacobian.length; ++i) {
552                 for (int j = 0; j < jacobian[i].length; ++j) {
553                     // check the values returned by getStateDerivatives() are correct
554                     maxRelativeError = FastMath.max(maxRelativeError,
555                                                     FastMath.abs((finiteDifferencesJacobian[i][j] - jacobian[i][j]) /
556                                                                  finiteDifferencesJacobian[i][j]));
557                 }
558             }
559 
560         }
561         Assertions.assertEquals(0, maxRelativeError, 3.4e-7);
562 
563     }
564 
565     /** Test the values of the parameters' derivatives using a numerical
566      * finite differences calculation as a reference.
567      * One-way measurements with modifiers (tropospheric corrections).
568      */
569     @Test
570     public void testParameterDerivativesWithModifier() {
571 
572         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
573 
574         final NumericalPropagatorBuilder propagatorBuilder =
575                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
576                                               1.0e-6, 60.0, 0.001);
577 
578         // create perfect range rate measurements
579         final double groundClockDrift =  4.8e-9;
580         for (final GroundStation station : context.stations) {
581             station.getClockDriftDriver().setValue(groundClockDrift);
582         }
583         final double satClkDrift = 3.2e-10;
584         final RangeRateMeasurementCreator creator = new RangeRateMeasurementCreator(context, false, satClkDrift);
585         creator.getSatellite().getClockDriftDriver().setSelected(true);
586         for (final GroundStation station : context.stations) {
587             station.getClockOffsetDriver().setSelected(true);
588             station.getClockDriftDriver().setSelected(true);
589             station.getEastOffsetDriver().setSelected(true);
590             station.getNorthOffsetDriver().setSelected(true);
591             station.getZenithOffsetDriver().setSelected(true);
592         }
593         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
594                                                                            propagatorBuilder);
595 
596         final List<ObservedMeasurement<?>> measurements =
597                         EstimationTestUtils.createMeasurements(propagator,
598                                                                creator,
599                                                                1.0, 3.0, 300.0);
600         propagator.clearStepHandlers();
601 
602         double maxRelativeError = 0;
603         for (final ObservedMeasurement<?> measurement : measurements) {
604 
605             final RangeRateTroposphericDelayModifier modifier =
606                 new RangeRateTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel(), true);
607             ((RangeRate) measurement).addModifier(modifier);
608 
609             // parameter corresponding to station position offset
610             final GroundStation stationParameter = ((RangeRate) measurement).getStation();
611 
612             // We intentionally propagate to a date which is close to the
613             // real spacecraft state but is *not* the accurate date, by
614             // compensating only part of the downlink delay. This is done
615             // in order to validate the partial derivatives with respect
616             // to velocity. If we had chosen the proper state date, the
617             // range would have depended only on the current position but
618             // not on the current velocity.
619             final double          meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
620             final AbsoluteDate    date      = measurement.getDate().shiftedBy(-0.75 * meanDelay);
621             final SpacecraftState state     = propagator.propagate(date);
622             final ParameterDriver[] drivers = new ParameterDriver[] {
623                 stationParameter.getClockDriftDriver(),
624                 stationParameter.getEastOffsetDriver(),
625                 stationParameter.getNorthOffsetDriver(),
626                 stationParameter.getZenithOffsetDriver(),
627                 measurement.getSatellites().get(0).getClockDriftDriver()
628             };
629             for (int i = 0; i < drivers.length; ++i) {
630                 final double[] gradient  = measurement.estimate(0, 0, new SpacecraftState[] { state }).getParameterDerivatives(drivers[i]);
631                 Assertions.assertEquals(1, measurement.getDimension());
632                 Assertions.assertEquals(1, gradient.length);
633 
634                 final ParameterFunction dMkdP =
635                                 Differentiation.differentiate(new ParameterFunction() {
636                                     /** {@inheritDoc} */
637                                     @Override
638                                     public double value(final ParameterDriver parameterDriver, AbsoluteDate date) {
639                                         return measurement.
640                                                estimateWithoutDerivatives(new SpacecraftState[] { state }).
641                                                getEstimatedValue()[0];
642                                     }
643                                 }, 3, 20.0 * drivers[i].getScale());
644                 final double ref = dMkdP.value(drivers[i], date);
645                 maxRelativeError = FastMath.max(maxRelativeError, FastMath.abs((ref - gradient[0]) / ref));
646             }
647 
648         }
649         Assertions.assertEquals(0, maxRelativeError, 1.2e-6);
650 
651     }
652 
653     /** Test the values of the parameters' derivatives using a numerical
654      * finite differences calculation as a reference.
655      * One-way measurements with estimated modifiers (tropospheric corrections).
656      */
657     @Test
658     public void testParameterDerivativesWithEstimatedModifier() {
659 
660         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
661 
662         final NumericalPropagatorBuilder propagatorBuilder =
663                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
664                                               1.0e-6, 60.0, 0.001);
665 
666         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
667                                                                            propagatorBuilder);
668 
669         final double groundClockDrift =  4.8e-9;
670         for (final GroundStation station : context.stations) {
671             station.getClockDriftDriver().setValue(groundClockDrift);
672         }
673         final double satClkDrift = 3.2e-10;
674         final List<ObservedMeasurement<?>> measurements =
675                         EstimationTestUtils.createMeasurements(propagator,
676                                                                new RangeRateMeasurementCreator(context, false, satClkDrift),
677                                                                1.0, 3.0, 300.0);
678         propagator.clearStepHandlers();
679 
680         double maxRelativeError = 0;
681         for (final ObservedMeasurement<?> measurement : measurements) {
682 
683             // Add modifiers if test implies it
684             final GlobalMappingFunctionModel mappingFunction = new GlobalMappingFunctionModel();
685             final EstimatedModel             tropoModel     = new EstimatedModel(mappingFunction, 10.0);
686 
687             final List<ParameterDriver> parameters = tropoModel.getParametersDrivers();
688             for (ParameterDriver driver : parameters) {
689                 driver.setSelected(true);
690             }
691 
692             final RangeRateTroposphericDelayModifier modifier = new RangeRateTroposphericDelayModifier(tropoModel, true);
693             ((RangeRate) measurement).addModifier(modifier);
694 
695             // We intentionally propagate to a date which is close to the
696             // real spacecraft state but is *not* the accurate date, by
697             // compensating only part of the downlink delay. This is done
698             // in order to validate the partial derivatives with respect
699             // to velocity. If we had chosen the proper state date, the
700             // range would have depended only on the current position but
701             // not on the current velocity.
702             final double          meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
703             final AbsoluteDate    date      = measurement.getDate().shiftedBy(-0.75 * meanDelay);
704             final SpacecraftState state     = propagator.propagate(date);
705 
706             final ParameterDriver[] drivers = new ParameterDriver[] {
707                 parameters.get(0)
708             };
709             for (int i = 0; i < 1; ++i) {
710                 final double[] gradient  = measurement.estimate(0, 0, new SpacecraftState[] { state }).getParameterDerivatives(drivers[i]);
711                 Assertions.assertEquals(1, measurement.getDimension());
712                 Assertions.assertEquals(1, gradient.length);
713 
714                 final ParameterFunction dMkdP =
715                                 Differentiation.differentiate(new ParameterFunction() {
716                                     /** {@inheritDoc} */
717                                     @Override
718                                     public double value(final ParameterDriver parameterDriver, AbsoluteDate date) {
719                                         return measurement.
720                                                estimateWithoutDerivatives(new SpacecraftState[] { state }).
721                                                getEstimatedValue()[0];
722                                     }
723                                 }, 3, 0.1 * drivers[i].getScale());
724                 final double ref = dMkdP.value(drivers[i], date);
725                 maxRelativeError = FastMath.max(maxRelativeError, FastMath.abs((ref - gradient[0]) / ref));
726             }
727 
728         }
729         Assertions.assertEquals(0, maxRelativeError, 2.2e-7);
730 
731     }
732 
733 }
734 
735