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.Arrays;
20  import java.util.List;
21  
22  import org.hipparchus.stat.descriptive.StreamingStatistics;
23  import org.hipparchus.util.FastMath;
24  import org.junit.jupiter.api.Assertions;
25  import org.junit.jupiter.api.Test;
26  import org.orekit.estimation.Context;
27  import org.orekit.estimation.EstimationTestUtils;
28  import org.orekit.estimation.measurements.modifiers.BistaticRangeRateTroposphericDelayModifier;
29  import org.orekit.models.earth.troposphere.ModifiedSaastamoinenModel;
30  import org.orekit.orbits.OrbitType;
31  import org.orekit.orbits.PositionAngleType;
32  import org.orekit.propagation.Propagator;
33  import org.orekit.propagation.SpacecraftState;
34  import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
35  import org.orekit.time.AbsoluteDate;
36  import org.orekit.utils.Differentiation;
37  import org.orekit.utils.ParameterDriver;
38  import org.orekit.utils.ParameterFunction;
39  
40  public class BistaticRangeRateTest {
41  
42      /**
43       * Compare observed values and estimated values.
44       * Both are calculated with a different algorithm.
45       */
46      @Test
47      public void testValues() {
48  
49          Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
50  
51          // Create perfect measurements
52          final NumericalPropagatorBuilder propagatorBuilder =
53                          context.createBuilder(OrbitType.EQUINOCTIAL, PositionAngleType.TRUE, false,
54                                                1.0e-6, 60.0, 0.001);
55          final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
56                                                                             propagatorBuilder);
57          final List<ObservedMeasurement<?>> measurements =
58                          EstimationTestUtils.createMeasurements(propagator,
59                                                                 new BistaticRangeRateMeasurementCreator(context),
60                                                                 1.0, 3.0, 300.0);
61          propagator.clearStepHandlers();
62  
63          // Prepare statistics for values difference
64          final StreamingStatistics diffStat = new StreamingStatistics();
65  
66          for (final ObservedMeasurement<?> measurement : measurements) {
67  
68              // Propagate to measurement date
69              final AbsoluteDate datemeas  = measurement.getDate();
70              SpacecraftState    state     = propagator.propagate(datemeas);
71  
72              // Estimate the measurement value
73              final EstimatedMeasurementBase<?> estimated = measurement.estimateWithoutDerivatives(new SpacecraftState[] { state });
74  
75              // Store the difference between estimated and observed values in the stats
76              diffStat.addValue(FastMath.abs(estimated.getEstimatedValue()[0] - measurement.getObservedValue()[0]));
77          }
78  
79          // Mean and std errors check
80          Assertions.assertEquals(0.0, diffStat.getMean(), 1.3e-7);
81          Assertions.assertEquals(0.0, diffStat.getStandardDeviation(), 1.2e-7);
82          // Test measurement type
83          Assertions.assertEquals(BistaticRangeRate.MEASUREMENT_TYPE, measurements.get(0).getMeasurementType());
84      }
85  
86      /**
87       * Test the values of the state derivatives using a numerical
88       * finite differences calculation as a reference.
89       */
90      @Test
91      public void testStateDerivatives() {
92  
93          Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
94  
95          // create perfect measurements
96          final NumericalPropagatorBuilder propagatorBuilder =
97                          context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
98                                                1.0e-6, 60.0, 0.001);
99          final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
100                                                                            propagatorBuilder);
101         final List<ObservedMeasurement<?>> measurements =
102                         EstimationTestUtils.createMeasurements(propagator,
103                                                                new BistaticRangeRateMeasurementCreator(context),
104                                                                1.0, 3.0, 300.0);
105         propagator.clearStepHandlers();
106 
107         double maxRelativeError = 0;
108         for (final ObservedMeasurement<?> measurement : measurements) {
109 
110             final AbsoluteDate    date  = measurement.getDate().shiftedBy(1);
111             final SpacecraftState state = propagator.propagate(date);
112 
113             final EstimatedMeasurement<?> estimated = measurement.estimate(0, 0, new SpacecraftState[] { state });
114             Assertions.assertEquals(3, estimated.getParticipants().length);
115             final double[][] jacobian = estimated.getStateDerivatives(0);
116 
117             final double[][] finiteDifferencesJacobian =
118                     Differentiation.differentiate(state1 -> measurement.
119                            estimateWithoutDerivatives(new SpacecraftState[] { state1 }).
120                            getEstimatedValue(), 1, propagator.getAttitudeProvider(),
121                                                   OrbitType.CARTESIAN, PositionAngleType.TRUE, 15.0, 3).value(state);
122 
123             Assertions.assertEquals(finiteDifferencesJacobian.length, jacobian.length);
124             Assertions.assertEquals(finiteDifferencesJacobian[0].length, jacobian[0].length);
125 
126             for (int i = 0; i < jacobian.length; ++i) {
127                 for (int j = 0; j < jacobian[i].length; ++j) {
128                     // check the values returned by getStateDerivatives() are correct
129                     maxRelativeError = FastMath.max(maxRelativeError,
130                                                     FastMath.abs((finiteDifferencesJacobian[i][j] - jacobian[i][j]) /
131                                                                  finiteDifferencesJacobian[i][j]));
132                 }
133             }
134         }
135 
136         Assertions.assertEquals(0, maxRelativeError, 3.7e-6);
137 
138     }
139 
140     /**
141      * Test the values of the state derivatives, using a numerical
142      * finite differences calculation as a reference, with modifiers (tropospheric corrections).
143      */
144     @Test
145     public void testStateDerivativesWithModifier() {
146 
147         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
148 
149         // create perfect measurements
150         final NumericalPropagatorBuilder propagatorBuilder =
151                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
152                                               1.0e-6, 60.0, 0.001);
153         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
154                                                                            propagatorBuilder);
155         final double clockOffset = 4.8e-9;
156         for (final GroundStation station : Arrays.asList(context.BRRstations.getKey(),
157                                                          context.BRRstations.getValue())) {
158             station.getClockOffsetDriver().setValue(clockOffset);
159         }
160         final List<ObservedMeasurement<?>> measurements =
161                         EstimationTestUtils.createMeasurements(propagator,
162                                                                new BistaticRangeRateMeasurementCreator(context),
163                                                                1.0, 3.0, 300.0);
164         propagator.clearStepHandlers();
165 
166         final BistaticRangeRateTroposphericDelayModifier modifier =
167             new BistaticRangeRateTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
168 
169         double maxRelativeError = 0;
170         for (final ObservedMeasurement<?> measurement : measurements) {
171 
172             ((BistaticRangeRate) measurement).addModifier(modifier);
173 
174             final AbsoluteDate    date  = measurement.getDate().shiftedBy(1);
175             final SpacecraftState state = propagator.propagate(date);
176 
177             final double[][] jacobian = measurement.estimate(0, 0, new SpacecraftState[] { state }).getStateDerivatives(0);
178 
179             final double[][] finiteDifferencesJacobian =
180                     Differentiation.differentiate(state1 -> measurement.
181                            estimateWithoutDerivatives(new SpacecraftState[] { state1 }).
182                            getEstimatedValue(), 1, propagator.getAttitudeProvider(),
183                                                   OrbitType.CARTESIAN, PositionAngleType.TRUE, 15.0, 3).value(state);
184 
185             Assertions.assertEquals(finiteDifferencesJacobian.length, jacobian.length);
186             Assertions.assertEquals(finiteDifferencesJacobian[0].length, jacobian[0].length);
187 
188             for (int i = 0; i < jacobian.length; ++i) {
189                 for (int j = 0; j < jacobian[i].length; ++j) {
190                     // check the values returned by getStateDerivatives() are correct
191                     maxRelativeError = FastMath.max(maxRelativeError,
192                                                     FastMath.abs((finiteDifferencesJacobian[i][j] - jacobian[i][j]) /
193                                                                  finiteDifferencesJacobian[i][j]));
194                 }
195             }
196         }
197 
198         Assertions.assertEquals(0, maxRelativeError, 3.7e-6);
199 
200     }
201 
202     /**
203      * Test the values of the parameters' derivatives using a numerical
204      * finite differences calculation as a reference.
205      */
206     @Test
207     public void testParameterDerivatives() {
208 
209         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
210 
211         // create perfect measurements
212         final NumericalPropagatorBuilder propagatorBuilder =
213                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
214                                               1.0e-6, 60.0, 0.001);
215 
216         final BistaticRangeRateMeasurementCreator creator = new BistaticRangeRateMeasurementCreator(context);
217         final double clockOffset = 4.8e-9;
218         for (final GroundStation station : Arrays.asList(context.BRRstations.getKey(),
219                                                          context.BRRstations.getValue())) {
220             station.getClockOffsetDriver().setValue(clockOffset);
221             station.getClockOffsetDriver().setSelected(true);
222             station.getEastOffsetDriver().setSelected(true);
223             station.getNorthOffsetDriver().setSelected(true);
224             station.getZenithOffsetDriver().setSelected(true);
225         }
226 
227         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
228                                                                            propagatorBuilder);
229         final List<ObservedMeasurement<?>> measurements =
230                         EstimationTestUtils.createMeasurements(propagator,
231                                                                creator,
232                                                                1.0, 3.0, 300.0);
233         propagator.clearStepHandlers();
234 
235         double maxRelativeError = 0;
236         for (final ObservedMeasurement<?> measurement : measurements) {
237 
238             // parameter corresponding to station position offset
239             final GroundStation emitterParameter = ((BistaticRangeRate) measurement).getEmitterStation();
240             final GroundStation receiverParameter = ((BistaticRangeRate) measurement).getReceiverStation();
241 
242             // We intentionally propagate to a date which is close to the
243             // real spacecraft state but is *not* the accurate date, by
244             // compensating only part of the downlink delay. This is done
245             // in order to validate the partial derivatives with respect
246             // to velocity. If we had chosen the proper state date, the
247             // range would have depended only on the current position but
248             // not on the current velocity.
249             final AbsoluteDate    date  = measurement.getDate().shiftedBy(0.05);
250             final SpacecraftState state = propagator.propagate(date);
251             final ParameterDriver[] drivers = new ParameterDriver[] {
252                 emitterParameter.getEastOffsetDriver(),
253                 emitterParameter.getNorthOffsetDriver(),
254                 emitterParameter.getZenithOffsetDriver(),
255                 receiverParameter.getEastOffsetDriver(),
256                 receiverParameter.getNorthOffsetDriver(),
257                 receiverParameter.getZenithOffsetDriver(),
258             };
259             for (int i = 0; i < drivers.length; ++i) {
260                 final double[] gradient  = measurement.estimate(0, 0, new SpacecraftState[] { state }).getParameterDerivatives(drivers[i]);
261                 Assertions.assertEquals(1, measurement.getDimension());
262                 Assertions.assertEquals(1, gradient.length);
263 
264                 final ParameterFunction dMkdP =
265                                 Differentiation.differentiate(new ParameterFunction() {
266                                     /** {@inheritDoc} */
267                                     @Override
268                                     public double value(final ParameterDriver parameterDriver, AbsoluteDate date) {
269                                         return measurement.
270                                                estimateWithoutDerivatives(new SpacecraftState[] { state }).
271                                                getEstimatedValue()[0];
272                                     }
273                                 }, 3, 20.0 * drivers[i].getScale());
274                 final double ref = dMkdP.value(drivers[i], date);
275                 maxRelativeError = FastMath.max(maxRelativeError, FastMath.abs((ref - gradient[0]) / ref));
276             }
277         }
278 
279         Assertions.assertEquals(0, maxRelativeError, 9.3e-8);
280 
281     }
282 
283     /**
284      * Test the values of the parameters' derivatives, using a numerical
285      * finite differences calculation as a reference, with modifiers (tropospheric corrections).
286      */
287     @Test
288     public void testParameterDerivativesWithModifier() {
289 
290         Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
291 
292         // create perfect measurements
293         final NumericalPropagatorBuilder propagatorBuilder =
294                         context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
295                                               1.0e-6, 60.0, 0.001);
296         final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
297                                                                            propagatorBuilder);
298 
299         final GroundStation emitter = context.BRRstations.getKey();
300         emitter.getEastOffsetDriver().setSelected(true);
301         emitter.getNorthOffsetDriver().setSelected(true);
302         emitter.getZenithOffsetDriver().setSelected(true);
303 
304         final double clockOffset = 4.8e-9;
305         final GroundStation receiver = context.BRRstations.getValue();
306         receiver.getClockOffsetDriver().setValue(clockOffset);
307         receiver.getClockOffsetDriver().setSelected(true);
308         receiver.getEastOffsetDriver().setSelected(true);
309         receiver.getNorthOffsetDriver().setSelected(true);
310         receiver.getZenithOffsetDriver().setSelected(true);
311 
312         final BistaticRangeRateMeasurementCreator creator = new BistaticRangeRateMeasurementCreator(context);
313         final List<ObservedMeasurement<?>> measurements =
314                         EstimationTestUtils.createMeasurements(propagator,
315                                                                creator,
316                                                                1.0, 3.0, 300.0);
317         propagator.clearStepHandlers();
318 
319         final BistaticRangeRateTroposphericDelayModifier modifier =
320             new BistaticRangeRateTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
321 
322         double maxRelativeError = 0;
323         for (final ObservedMeasurement<?> measurement : measurements) {
324 
325             ((BistaticRangeRate) measurement).addModifier(modifier);
326 
327             // parameter corresponding to station position offset
328             final GroundStation emitterParameter  = ((BistaticRangeRate) measurement).getEmitterStation();
329             final GroundStation receiverParameter = ((BistaticRangeRate) measurement).getReceiverStation();
330 
331             // We intentionally propagate to a date which is close to the
332             // real spacecraft state but is *not* the accurate date, by
333             // compensating only part of the downlink delay. This is done
334             // in order to validate the partial derivatives with respect
335             // to velocity. If we had chosen the proper state date, the
336             // range would have depended only on the current position but
337             // not on the current velocity.
338             final AbsoluteDate    date  = measurement.getDate().shiftedBy(0.05);
339             final SpacecraftState state = propagator.propagate(date);
340             final ParameterDriver[] drivers = new ParameterDriver[] {
341                 emitterParameter.getEastOffsetDriver(),
342                 emitterParameter.getNorthOffsetDriver(),
343                 emitterParameter.getZenithOffsetDriver(),
344                 receiverParameter.getClockOffsetDriver(),
345                 receiverParameter.getEastOffsetDriver(),
346                 receiverParameter.getNorthOffsetDriver(),
347                 receiverParameter.getZenithOffsetDriver(),
348             };
349             for (int i = 0; i < drivers.length; ++i) {
350                 final double[] gradient = measurement.estimate(0, 0, new SpacecraftState[] { state }).getParameterDerivatives(drivers[i]);
351                 Assertions.assertEquals(1, measurement.getDimension());
352                 Assertions.assertEquals(1, gradient.length);
353 
354                 final ParameterFunction dMkdP =
355                                 Differentiation.differentiate(new ParameterFunction() {
356                                     /** {@inheritDoc} */
357                                     @Override
358                                     public double value(final ParameterDriver parameterDriver, AbsoluteDate date) {
359                                         return measurement.
360                                                estimateWithoutDerivatives(new SpacecraftState[] { state }).
361                                                getEstimatedValue()[0];
362                                     }
363                                 }, 3, 20.0 * drivers[i].getScale());
364                 final double ref = dMkdP.value(drivers[i], date);
365                 maxRelativeError = FastMath.max(maxRelativeError, FastMath.abs((ref - gradient[0]) / ref));
366             }
367         }
368 
369         Assertions.assertEquals(0, maxRelativeError, 2.0e-5);
370 
371     }
372 
373 }