1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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.BistaticRangeTroposphericDelayModifier;
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.Constants;
37 import org.orekit.utils.Differentiation;
38 import org.orekit.utils.ParameterDriver;
39 import org.orekit.utils.ParameterFunction;
40
41 public class BistaticRangeTest {
42
43
44
45
46
47 @Test
48 public void testValues() {
49
50 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
51
52
53 final NumericalPropagatorBuilder propagatorBuilder =
54 context.createBuilder(OrbitType.EQUINOCTIAL, PositionAngleType.TRUE, false,
55 1.0e-6, 60.0, 0.001);
56 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
57 propagatorBuilder);
58 final List<ObservedMeasurement<?>> measurements =
59 EstimationTestUtils.createMeasurements(propagator,
60 new BistaticRangeMeasurementCreator(context),
61 1.0, 3.0, 300.0);
62 propagator.clearStepHandlers();
63
64
65 final StreamingStatistics diffStat = new StreamingStatistics();
66
67 for (final ObservedMeasurement<?> measurement : measurements) {
68
69
70 final AbsoluteDate datemeas = measurement.getDate();
71 SpacecraftState state = propagator.propagate(datemeas);
72
73
74 final EstimatedMeasurementBase<?> estimated = measurement.estimateWithoutDerivatives(new SpacecraftState[] { state });
75
76
77 diffStat.addValue(FastMath.abs(estimated.getEstimatedValue()[0] - measurement.getObservedValue()[0]));
78 }
79
80
81 Assertions.assertEquals(0.0, diffStat.getMean(), 1.59e-7);
82 Assertions.assertEquals(0.0, diffStat.getStandardDeviation(), 1.3e-7);
83
84
85 Assertions.assertEquals(BistaticRange.MEASUREMENT_TYPE, measurements.get(0).getMeasurementType());
86 }
87
88
89
90
91
92 @Test
93 public void testStateDerivatives() {
94
95 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
96
97
98 final NumericalPropagatorBuilder propagatorBuilder =
99 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
100 1.0e-6, 60.0, 0.001);
101 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
102 propagatorBuilder);
103 final List<ObservedMeasurement<?>> measurements =
104 EstimationTestUtils.createMeasurements(propagator,
105 new BistaticRangeMeasurementCreator(context),
106 1.0, 3.0, 300.0);
107 propagator.clearStepHandlers();
108
109 double maxRelativeError = 0;
110 for (final ObservedMeasurement<?> measurement : measurements) {
111
112 final AbsoluteDate date = measurement.getDate().shiftedBy(1);
113 final SpacecraftState state = propagator.propagate(date);
114
115 final EstimatedMeasurement<?> estimated = measurement.estimate(0, 0, new SpacecraftState[] { state });
116 Assertions.assertEquals(3, estimated.getParticipants().length);
117 final double[][] jacobian = estimated.getStateDerivatives(0);
118
119 final double[][] finiteDifferencesJacobian =
120 Differentiation.differentiate(state1 -> measurement.
121 estimateWithoutDerivatives(new SpacecraftState[] { state1 }).
122 getEstimatedValue(), 1, propagator.getAttitudeProvider(),
123 OrbitType.CARTESIAN, PositionAngleType.TRUE, 15.0, 3).value(state);
124
125 Assertions.assertEquals(finiteDifferencesJacobian.length, jacobian.length);
126 Assertions.assertEquals(finiteDifferencesJacobian[0].length, jacobian[0].length);
127
128 for (int i = 0; i < jacobian.length; ++i) {
129 for (int j = 0; j < jacobian[i].length; ++j) {
130
131 maxRelativeError = FastMath.max(maxRelativeError,
132 FastMath.abs((finiteDifferencesJacobian[i][j] - jacobian[i][j]) /
133 finiteDifferencesJacobian[i][j]));
134 }
135 }
136 }
137
138 Assertions.assertEquals(0, maxRelativeError, 2.7e-5);
139
140 }
141
142
143
144
145
146 @Test
147 public void testStateDerivativesWithModifier() {
148
149 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
150
151
152 final NumericalPropagatorBuilder propagatorBuilder =
153 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
154 1.0e-6, 60.0, 0.001);
155 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
156 propagatorBuilder);
157 final double clockOffset = 4.8e-9;
158 for (final GroundStation station : Arrays.asList(context.BRRstations.getKey(),
159 context.BRRstations.getValue())) {
160 station.getClockOffsetDriver().setValue(clockOffset);
161 }
162 final List<ObservedMeasurement<?>> measurements =
163 EstimationTestUtils.createMeasurements(propagator,
164 new BistaticRangeMeasurementCreator(context),
165 1.0, 3.0, 300.0);
166 propagator.clearStepHandlers();
167
168 final BistaticRangeTroposphericDelayModifier modifier =
169 new BistaticRangeTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
170
171 double maxRelativeError = 0;
172 for (final ObservedMeasurement<?> measurement : measurements) {
173
174 ((BistaticRange) measurement).addModifier(modifier);
175
176 final AbsoluteDate date = measurement.getDate().shiftedBy(1);
177 final SpacecraftState state = propagator.propagate(date);
178
179 final double[][] jacobian = measurement.estimate(0, 0, new SpacecraftState[] { state }).getStateDerivatives(0);
180
181 final double[][] finiteDifferencesJacobian =
182 Differentiation.differentiate(state1 -> measurement.
183 estimateWithoutDerivatives(new SpacecraftState[] { state1 }).
184 getEstimatedValue(), 1, propagator.getAttitudeProvider(),
185 OrbitType.CARTESIAN, PositionAngleType.TRUE, 15.0, 3).value(state);
186
187 Assertions.assertEquals(finiteDifferencesJacobian.length, jacobian.length);
188 Assertions.assertEquals(finiteDifferencesJacobian[0].length, jacobian[0].length);
189
190 for (int i = 0; i < jacobian.length; ++i) {
191 for (int j = 0; j < jacobian[i].length; ++j) {
192
193 maxRelativeError = FastMath.max(maxRelativeError,
194 FastMath.abs((finiteDifferencesJacobian[i][j] - jacobian[i][j]) /
195 finiteDifferencesJacobian[i][j]));
196 }
197 }
198 }
199
200 Assertions.assertEquals(0, maxRelativeError, 2.7e-5);
201
202 }
203
204
205
206
207
208 @Test
209 public void testParameterDerivatives() {
210
211 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
212
213
214 final NumericalPropagatorBuilder propagatorBuilder =
215 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
216 1.0e-6, 60.0, 0.001);
217
218 final BistaticRangeMeasurementCreator creator = new BistaticRangeMeasurementCreator(context);
219 final double clockOffset = 4.8e-9;
220 for (final GroundStation station : Arrays.asList(context.BRRstations.getKey(),
221 context.BRRstations.getValue())) {
222 station.getClockOffsetDriver().setValue(clockOffset);
223 station.getClockOffsetDriver().setSelected(true);
224 station.getEastOffsetDriver().setSelected(true);
225 station.getNorthOffsetDriver().setSelected(true);
226 station.getZenithOffsetDriver().setSelected(true);
227 }
228
229 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
230 propagatorBuilder);
231 final List<ObservedMeasurement<?>> measurements =
232 EstimationTestUtils.createMeasurements(propagator,
233 creator,
234 1.0, 3.0, 300.0);
235 propagator.clearStepHandlers();
236
237 double maxRelativeError = 0;
238 for (final ObservedMeasurement<?> measurement : measurements) {
239
240
241 final GroundStation emitterParameter = ((BistaticRange) measurement).getEmitterStation();
242 final GroundStation receiverParameter = ((BistaticRange) measurement).getReceiverStation();
243
244
245
246
247
248
249
250
251 final AbsoluteDate date = measurement.getDate().shiftedBy(0.05);
252 final SpacecraftState state = propagator.propagate(date);
253 final ParameterDriver[] drivers = new ParameterDriver[] {
254 emitterParameter.getEastOffsetDriver(),
255 emitterParameter.getNorthOffsetDriver(),
256 emitterParameter.getZenithOffsetDriver(),
257 receiverParameter.getEastOffsetDriver(),
258 receiverParameter.getNorthOffsetDriver(),
259 receiverParameter.getZenithOffsetDriver(),
260 };
261 for (int i = 0; i < drivers.length; ++i) {
262 final double[] gradient = measurement.estimate(0, 0, new SpacecraftState[] { state }).getParameterDerivatives(drivers[i]);
263 Assertions.assertEquals(1, measurement.getDimension());
264 Assertions.assertEquals(1, gradient.length);
265
266 final ParameterFunction dMkdP =
267 Differentiation.differentiate(new ParameterFunction() {
268
269 @Override
270 public double value(final ParameterDriver parameterDriver, final AbsoluteDate date) {
271 return measurement.
272 estimateWithoutDerivatives(new SpacecraftState[] { state }).
273 getEstimatedValue()[0];
274 }
275 }, 3, 20.0 * drivers[i].getScale());
276 final double ref = dMkdP.value(drivers[i], date);
277 maxRelativeError = FastMath.max(maxRelativeError, FastMath.abs((ref - gradient[0]) / ref));
278 }
279 }
280
281 Assertions.assertEquals(0, maxRelativeError, 1.9e-7);
282
283 }
284
285
286
287
288
289 @Test
290 public void testParameterDerivativesWithModifier() {
291
292 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
293
294
295 final NumericalPropagatorBuilder propagatorBuilder =
296 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
297 1.0e-6, 60.0, 0.001);
298 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
299 propagatorBuilder);
300
301 final GroundStation emitter = context.BRRstations.getKey();
302 emitter.getEastOffsetDriver().setSelected(true);
303 emitter.getNorthOffsetDriver().setSelected(true);
304 emitter.getZenithOffsetDriver().setSelected(true);
305
306 final double clockOffset = 4.8e-9;
307 final GroundStation receiver = context.BRRstations.getValue();
308 receiver.getClockOffsetDriver().setValue(clockOffset);
309 receiver.getClockOffsetDriver().setSelected(true);
310 receiver.getEastOffsetDriver().setSelected(true);
311 receiver.getNorthOffsetDriver().setSelected(true);
312 receiver.getZenithOffsetDriver().setSelected(true);
313
314 final BistaticRangeMeasurementCreator creator = new BistaticRangeMeasurementCreator(context);
315 final List<ObservedMeasurement<?>> measurements =
316 EstimationTestUtils.createMeasurements(propagator,
317 creator,
318 1.0, 3.0, 300.0);
319 propagator.clearStepHandlers();
320
321 final BistaticRangeTroposphericDelayModifier modifier =
322 new BistaticRangeTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
323
324 double maxRelativeError = 0;
325 for (final ObservedMeasurement<?> measurement : measurements) {
326
327 ((BistaticRange) measurement).addModifier(modifier);
328
329
330 final GroundStation emitterParameter = ((BistaticRange) measurement).getEmitterStation();
331 final GroundStation receiverParameter = ((BistaticRange) measurement).getReceiverStation();
332
333
334
335
336
337
338
339
340 final AbsoluteDate date = measurement.getDate().shiftedBy(0.05);
341 final SpacecraftState state = propagator.propagate(date);
342 final ParameterDriver[] drivers = new ParameterDriver[] {
343 emitterParameter.getEastOffsetDriver(),
344 emitterParameter.getNorthOffsetDriver(),
345 emitterParameter.getZenithOffsetDriver(),
346 receiverParameter.getClockOffsetDriver(),
347 receiverParameter.getEastOffsetDriver(),
348 receiverParameter.getNorthOffsetDriver(),
349 receiverParameter.getZenithOffsetDriver(),
350 };
351 for (int i = 0; i < drivers.length; ++i) {
352 final double[] gradient = measurement.estimate(0, 0, new SpacecraftState[] { state }).getParameterDerivatives(drivers[i]);
353 Assertions.assertEquals(1, measurement.getDimension());
354 Assertions.assertEquals(1, gradient.length);
355
356 final ParameterFunction dMkdP =
357 Differentiation.differentiate(new ParameterFunction() {
358
359 @Override
360 public double value(final ParameterDriver parameterDriver, final AbsoluteDate date) {
361 return measurement.estimate(0, 0, new SpacecraftState[] { state }).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.9e-6);
370
371 }
372
373
374
375
376
377 @Test
378 public void testIssue1418() {
379 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
380
381
382 final double emitterClockOffset = 5e-9;
383 final double receiverClockOffset = 10e-9;
384 final GroundStation emitter = context.BRRstations.getKey();
385 final GroundStation receiver = context.BRRstations.getValue();
386 emitter.getClockOffsetDriver().setValue(emitterClockOffset);
387 receiver.getClockOffsetDriver().setValue(receiverClockOffset);
388
389
390 final NumericalPropagatorBuilder propagatorBuilder =
391 context.createBuilder(OrbitType.EQUINOCTIAL, PositionAngleType.TRUE, false,
392 1.0e-6, 60.0, 0.001);
393 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
394 propagatorBuilder);
395 final List<ObservedMeasurement<?>> measurements =
396 EstimationTestUtils.createMeasurements(propagator,
397 new BistaticRangeMeasurementCreator(context),
398 1.0, 3.0, 300.0);
399 propagator.clearStepHandlers();
400
401
402 final StreamingStatistics diffStat = new StreamingStatistics();
403
404 for (final ObservedMeasurement<?> measurement : measurements) {
405
406
407 final AbsoluteDate datemeas = measurement.getDate();
408 SpacecraftState state = propagator.propagate(datemeas);
409
410
411 final EstimatedMeasurementBase<?> estimated = measurement.estimateWithoutDerivatives(new SpacecraftState[] { state });
412
413
414 diffStat.addValue(FastMath.abs(estimated.getEstimatedValue()[0] - measurement.getObservedValue()[0]));
415 }
416
417
418 final double clockOffsetShift = (receiverClockOffset - emitterClockOffset) * Constants.SPEED_OF_LIGHT;
419 Assertions.assertEquals(clockOffsetShift, diffStat.getMean(), 1e-6);
420
421
422 }
423 }