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.ArrayList;
20 import java.util.List;
21 import java.util.Locale;
22 import java.util.Map;
23
24 import org.hipparchus.stat.descriptive.moment.Mean;
25 import org.hipparchus.stat.descriptive.rank.Max;
26 import org.hipparchus.stat.descriptive.rank.Median;
27 import org.hipparchus.stat.descriptive.rank.Min;
28 import org.hipparchus.util.FastMath;
29 import org.junit.jupiter.api.Assertions;
30 import org.junit.jupiter.api.Test;
31 import org.orekit.estimation.Context;
32 import org.orekit.estimation.EstimationTestUtils;
33 import org.orekit.estimation.measurements.modifiers.TurnAroundRangeTroposphericDelayModifier;
34 import org.orekit.models.earth.troposphere.ModifiedSaastamoinenModel;
35 import org.orekit.orbits.OrbitType;
36 import org.orekit.orbits.PositionAngleType;
37 import org.orekit.propagation.Propagator;
38 import org.orekit.propagation.SpacecraftState;
39 import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
40 import org.orekit.time.AbsoluteDate;
41 import org.orekit.utils.Constants;
42 import org.orekit.utils.Differentiation;
43 import org.orekit.utils.ParameterDriver;
44 import org.orekit.utils.ParameterFunction;
45 import org.orekit.utils.TimeStampedPVCoordinates;
46
47 public class TurnAroundRangeTest {
48
49
50
51
52
53
54 @Test
55 public void testValues() {
56 boolean printResults = false;
57 if (printResults) {
58 System.out.println("\nTest TAR Values\n");
59 }
60
61 this.genericTestValues(printResults);
62 }
63
64
65
66
67
68 @Test
69 public void testStateDerivatives() {
70
71 boolean printResults = false;
72 if (printResults) {
73 System.out.println("\nTest TAR State Derivatives - Finite Differences comparison\n");
74 }
75
76 boolean isModifier = false;
77 double refErrorsPMedian = 1.4e-6;
78 double refErrorsPMean = 1.4e-06;
79 double refErrorsPMax = 2.6e-06;
80 double refErrorsVMedian = 8.2e-05;
81 double refErrorsVMean = 3.6e-04;
82 double refErrorsVMax = 1.4e-02;
83
84 this.genericTestStateDerivatives(isModifier, printResults,
85 refErrorsPMedian, refErrorsPMean, refErrorsPMax,
86 refErrorsVMedian, refErrorsVMean, refErrorsVMax);
87 }
88
89
90
91
92
93 @Test
94 public void testStateDerivativesWithModifier() {
95
96 boolean printResults = false;
97 if (printResults) {
98 System.out.println("\nTest TAR State Derivatives with Modifier - Finite Differences comparison\n");
99 }
100
101 boolean isModifier = true;
102 double refErrorsPMedian = 1.4e-6;
103 double refErrorsPMean = 1.4e-06;
104 double refErrorsPMax = 2.6e-06;
105 double refErrorsVMedian = 8.2e-05;
106 double refErrorsVMean = 3.6e-04;
107 double refErrorsVMax = 1.4e-2;
108
109 this.genericTestStateDerivatives(isModifier, printResults,
110 refErrorsPMedian, refErrorsPMean, refErrorsPMax,
111 refErrorsVMedian, refErrorsVMean, refErrorsVMax);
112 }
113
114
115
116
117
118 @Test
119 public void testParameterDerivatives() {
120
121
122 boolean printResults = false;
123
124 if (printResults) {
125 System.out.println("\nTest TAR Parameter Derivatives - Finite Differences comparison\n");
126 }
127
128 boolean isModifier = false;
129 double refErrorQMMedian = 2.6e-6;
130 double refErrorQMMean = 2.5e-6;
131 double refErrorQMMax = 5.6e-6;
132 double refErrorQSMedian = 3.7e-7;
133 double refErrorQSMean = 3.6e-7;
134 double refErrorQSMax = 7.7e-7;
135 this.genericTestParameterDerivatives(isModifier, printResults,
136 refErrorQMMedian, refErrorQMMean, refErrorQMMax,
137 refErrorQSMedian, refErrorQSMean, refErrorQSMax);
138
139 }
140
141
142
143
144
145 @Test
146 public void testParameterDerivativesWithModifier() {
147
148
149 boolean printResults = false;
150
151 if (printResults) {
152 System.out.println("\nTest TAR Parameter Derivatives with Modifier - Finite Differences comparison\n");
153 }
154
155 boolean isModifier = true;
156 double refErrorQMMedian = 1.3e-8;
157 double refErrorQMMean = 8.9e-8;
158 double refErrorQMMax = 5.1e-6;
159 double refErrorQSMedian = 3.4e-8;
160 double refErrorQSMean = 1.9e-5;
161 double refErrorQSMax = 2.2e-4;
162 this.genericTestParameterDerivatives(isModifier, printResults,
163 refErrorQMMedian, refErrorQMMean, refErrorQMMax,
164 refErrorQSMedian, refErrorQSMean, refErrorQSMax);
165 }
166
167
168
169
170
171 void genericTestValues(final boolean printResults) {
172
173 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
174
175
176 final NumericalPropagatorBuilder propagatorBuilder =
177 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
178 1.0e-6, 60.0, 0.001);
179
180
181 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
182 propagatorBuilder);
183 final List<ObservedMeasurement<?>> measurements =
184 EstimationTestUtils.createMeasurements(propagator,
185 new TurnAroundRangeMeasurementCreator(context),
186 1.0, 3.0, 300.0);
187 propagator.clearStepHandlers();
188
189 double[] absoluteErrors = new double[measurements.size()];
190 double[] relativeErrors = new double[measurements.size()];
191 int index = 0;
192
193
194 if (printResults) {
195 System.out.format(Locale.US, "%-15s %-15s %-23s %-23s %17s %17s %13s %13s%n",
196 "Primary Station", "secondary Station",
197 "Measurement Date", "State Date",
198 "TAR observed [m]", "TAR estimated [m]",
199 "|ΔTAR| [m]", "rel |ΔTAR|");
200 }
201
202
203 for (final ObservedMeasurement<?> measurement : measurements) {
204
205 final double meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
206 final AbsoluteDate date = measurement.getDate().shiftedBy(meanDelay);
207 final SpacecraftState state = propagator.propagate(date);
208
209
210 final double TARobserved = measurement.getObservedValue()[0];
211 final EstimatedMeasurementBase<?> estimated = measurement.estimateWithoutDerivatives(new SpacecraftState[] { state });
212 final double TARestimated = estimated.getEstimatedValue()[0];
213
214 final TimeStampedPVCoordinates[] participants = estimated.getParticipants();
215 Assertions.assertEquals(5, participants.length);
216 Assertions.assertEquals(0.5 * Constants.SPEED_OF_LIGHT * participants[4].getDate().durationFrom(participants[0].getDate()),
217 estimated.getEstimatedValue()[0],
218 2.3e-8);
219
220 absoluteErrors[index] = TARestimated-TARobserved;
221 relativeErrors[index] = FastMath.abs(absoluteErrors[index])/FastMath.abs(TARobserved);
222 index++;
223
224
225 if (printResults) {
226 final AbsoluteDate measurementDate = measurement.getDate();
227
228 String primaryStationName = ((TurnAroundRange) measurement).getPrimaryStation().getBaseFrame().getName();
229 String secondaryStationName = ((TurnAroundRange) measurement).getSecondaryStation().getBaseFrame().getName();
230 System.out.format(Locale.US, "%-15s %-15s %-23s %-23s %17.6f %17.6f %13.6e %13.6e%n",
231 primaryStationName, secondaryStationName, measurementDate, date,
232 TARobserved, TARestimated,
233 FastMath.abs(TARestimated-TARobserved),
234 FastMath.abs((TARestimated-TARobserved)/TARobserved));
235 }
236
237 }
238
239
240 final double absErrorsMedian = new Median().evaluate(absoluteErrors);
241 final double absErrorsMin = new Min().evaluate(absoluteErrors);
242 final double absErrorsMax = new Max().evaluate(absoluteErrors);
243 final double relErrorsMedian = new Median().evaluate(relativeErrors);
244 final double relErrorsMax = new Max().evaluate(relativeErrors);
245
246
247 if (printResults) {
248 System.out.println();
249 System.out.println("Absolute errors median: " + absErrorsMedian);
250 System.out.println("Absolute errors min : " + absErrorsMin);
251 System.out.println("Absolute errors max : " + absErrorsMax);
252 System.out.println("Relative errors median: " + relErrorsMedian);
253 System.out.println("Relative errors max : " + relErrorsMax);
254 }
255
256
257 Assertions.assertEquals(0.0, absErrorsMedian, 1.5e-7);
258 Assertions.assertEquals(0.0, absErrorsMin, 5.0e-7);
259 Assertions.assertEquals(0.0, absErrorsMax, 4.9e-7);
260 Assertions.assertEquals(0.0, relErrorsMedian, 9.2e-15);
261 Assertions.assertEquals(0.0, relErrorsMax , 2.9e-14);
262
263
264 Assertions.assertEquals(TurnAroundRange.MEASUREMENT_TYPE, measurements.get(0).getMeasurementType());
265 }
266
267 void genericTestStateDerivatives(final boolean isModifier, final boolean printResults,
268 final double refErrorsPMedian, final double refErrorsPMean, final double refErrorsPMax,
269 final double refErrorsVMedian, final double refErrorsVMean, final double refErrorsVMax) {
270
271 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
272
273
274 final NumericalPropagatorBuilder propagatorBuilder =
275 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
276 1.0e-6, 60.0, 0.001);
277
278
279 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
280 propagatorBuilder);
281 final List<ObservedMeasurement<?>> measurements =
282 EstimationTestUtils.createMeasurements(propagator,
283 new TurnAroundRangeMeasurementCreator(context),
284 1.0, 3.0, 300.0);
285 propagator.clearStepHandlers();
286
287 double[] errorsP = new double[3 * measurements.size()];
288 double[] errorsV = new double[3 * measurements.size()];
289 int indexP = 0;
290 int indexV = 0;
291
292
293 if (printResults) {
294 System.out.format(Locale.US, "%-15s %-15s %-23s %-23s " +
295 "%10s %10s %10s " +
296 "%10s %10s %10s " +
297 "%10s %10s %10s " +
298 "%10s %10s %10s%n",
299 "Primary Station", "secondary Station",
300 "Measurement Date", "State Date",
301 "ΔdPx", "ΔdPy", "ΔdPz", "ΔdVx", "ΔdVy", "ΔdVz",
302 "rel ΔdPx", "rel ΔdPy", "rel ΔdPz",
303 "rel ΔdVx", "rel ΔdVy", "rel ΔdVz");
304 }
305
306
307 for (final ObservedMeasurement<?> measurement : measurements) {
308
309
310 final TurnAroundRangeTroposphericDelayModifier modifier =
311 new TurnAroundRangeTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
312 if (isModifier) {
313 ((TurnAroundRange) measurement).addModifier(modifier);
314 }
315
316
317
318
319
320
321
322
323 final double meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
324 final AbsoluteDate date = measurement.getDate().shiftedBy(-0.75 * meanDelay);
325 final SpacecraftState state = propagator.propagate(date);
326 final double[][] jacobian = measurement.estimate(0, 0, new SpacecraftState[] { state }).getStateDerivatives(0);
327
328
329 final double[][] jacobianRef;
330
331
332 jacobianRef = Differentiation.differentiate(state1 ->
333 measurement.
334 estimateWithoutDerivatives(new SpacecraftState[] { state1 }).
335 getEstimatedValue(),
336 measurement.getDimension(), propagator.getAttitudeProvider(),
337 OrbitType.CARTESIAN, PositionAngleType.TRUE, 2.0, 3).value(state);
338
339 Assertions.assertEquals(jacobianRef.length, jacobian.length);
340 Assertions.assertEquals(jacobianRef[0].length, jacobian[0].length);
341
342 double [][] dJacobian = new double[jacobian.length][jacobian[0].length];
343 double [][] dJacobianRelative = new double[jacobian.length][jacobian[0].length];
344 for (int i = 0; i < jacobian.length; ++i) {
345 for (int j = 0; j < jacobian[i].length; ++j) {
346 dJacobian[i][j] = jacobian[i][j] - jacobianRef[i][j];
347 dJacobianRelative[i][j] = FastMath.abs(dJacobian[i][j]/jacobianRef[i][j]);
348
349 if (j < 3) {
350 errorsP[indexP++] = dJacobianRelative[i][j];
351 } else {
352 errorsV[indexV++] = dJacobianRelative[i][j];
353 }
354 }
355 }
356
357 if (printResults) {
358 String primaryStationName = ((TurnAroundRange) measurement).getPrimaryStation().getBaseFrame().getName();
359 String secondaryStationName = ((TurnAroundRange) measurement).getSecondaryStation().getBaseFrame().getName();
360 System.out.format(Locale.US, "%-15s %-15s %-23s %-23s " +
361 "%10.3e %10.3e %10.3e " +
362 "%10.3e %10.3e %10.3e " +
363 "%10.3e %10.3e %10.3e " +
364 "%10.3e %10.3e %10.3e%n",
365 primaryStationName, secondaryStationName, measurement.getDate(), date,
366 dJacobian[0][0], dJacobian[0][1], dJacobian[0][2],
367 dJacobian[0][3], dJacobian[0][4], dJacobian[0][5],
368 dJacobianRelative[0][0], dJacobianRelative[0][1], dJacobianRelative[0][2],
369 dJacobianRelative[0][3], dJacobianRelative[0][4], dJacobianRelative[0][5]);
370 }
371 }
372
373
374 final double errorsPMedian = new Median().evaluate(errorsP);
375 final double errorsPMean = new Mean().evaluate(errorsP);
376 final double errorsPMax = new Max().evaluate(errorsP);
377 final double errorsVMedian = new Median().evaluate(errorsV);
378 final double errorsVMean = new Mean().evaluate(errorsV);
379 final double errorsVMax = new Max().evaluate(errorsV);
380
381
382
383 if (printResults) {
384 System.out.println();
385 System.out.format(Locale.US, "Relative errors dR/dP -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n",
386 errorsPMedian, errorsPMean, errorsPMax);
387 System.out.format(Locale.US, "Relative errors dR/dV -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n",
388 errorsVMedian, errorsVMean, errorsVMax);
389 }
390
391 Assertions.assertEquals(0.0, errorsPMedian, refErrorsPMedian);
392 Assertions.assertEquals(0.0, errorsPMean, refErrorsPMean);
393 Assertions.assertEquals(0.0, errorsPMax, refErrorsPMax);
394 Assertions.assertEquals(0.0, errorsVMedian, refErrorsVMedian);
395 Assertions.assertEquals(0.0, errorsVMean, refErrorsVMean);
396 Assertions.assertEquals(0.0, errorsVMax, refErrorsVMax);
397 }
398
399
400 void genericTestParameterDerivatives(final boolean isModifier, final boolean printResults,
401 final double refErrorQMMedian, final double refErrorQMMean, final double refErrorQMMax,
402 final double refErrorQSMedian, final double refErrorQSMean, final double refErrorQSMax) {
403
404 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
405
406 final NumericalPropagatorBuilder propagatorBuilder =
407 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
408 1.0e-6, 60.0, 0.001);
409
410
411 for (Map.Entry<GroundStation, GroundStation> entry : context.TARstations.entrySet()) {
412 final GroundStation primaryStation = entry.getKey();
413 final GroundStation secondaryStation = entry.getValue();
414 primaryStation.getClockOffsetDriver().setSelected(true);
415 primaryStation.getEastOffsetDriver().setSelected(true);
416 primaryStation.getNorthOffsetDriver().setSelected(true);
417 primaryStation.getZenithOffsetDriver().setSelected(true);
418 secondaryStation.getClockOffsetDriver().setSelected(false);
419 secondaryStation.getEastOffsetDriver().setSelected(true);
420 secondaryStation.getNorthOffsetDriver().setSelected(true);
421 secondaryStation.getZenithOffsetDriver().setSelected(true);
422 }
423 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
424 propagatorBuilder);
425 final List<ObservedMeasurement<?>> measurements =
426 EstimationTestUtils.createMeasurements(propagator,
427 new TurnAroundRangeMeasurementCreator(context),
428 1.0, 3.0, 300.0);
429 propagator.clearStepHandlers();
430
431
432 if (printResults) {
433 System.out.format(Locale.US, "%-15s %-15s %-23s %-23s " +
434 "%10s %10s %10s " +
435 "%10s %10s %10s " +
436 "%10s %10s %10s " +
437 "%10s %10s %10s%n",
438 "Primary Station", "secondary Station",
439 "Measurement Date", "State Date",
440 "ΔdQMx", "rel ΔdQMx",
441 "ΔdQMy", "rel ΔdQMy",
442 "ΔdQMz", "rel ΔdQMz",
443 "ΔdQSx", "rel ΔdQSx",
444 "ΔdQSy", "rel ΔdQSy",
445 "ΔdQSz", "rel ΔdQSz");
446 }
447
448
449 final List<Double> relErrorQMList = new ArrayList<>();
450 final List<Double> relErrorQSList = new ArrayList<>();
451
452
453 for (final ObservedMeasurement<?> measurement : measurements) {
454
455
456 final TurnAroundRangeTroposphericDelayModifier modifier =
457 new TurnAroundRangeTroposphericDelayModifier(ModifiedSaastamoinenModel.getStandardModel());
458 if (isModifier) {
459 ((TurnAroundRange) measurement).addModifier(modifier);
460 }
461
462
463 final GroundStation primaryStationParameter = ((TurnAroundRange) measurement).getPrimaryStation();
464 final GroundStation secondaryStationParameter = ((TurnAroundRange) measurement).getSecondaryStation();
465
466
467
468
469
470
471
472
473 final double meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
474 final AbsoluteDate date = measurement.getDate().shiftedBy(-0.75 * meanDelay);
475 final SpacecraftState state = propagator.propagate(date);
476 final ParameterDriver[] drivers = new ParameterDriver[] {
477 primaryStationParameter.getEastOffsetDriver(),
478 primaryStationParameter.getNorthOffsetDriver(),
479 primaryStationParameter.getZenithOffsetDriver(),
480 secondaryStationParameter.getEastOffsetDriver(),
481 secondaryStationParameter.getNorthOffsetDriver(),
482 secondaryStationParameter.getZenithOffsetDriver(),
483 };
484
485
486 if (printResults) {
487 String primaryStationName = primaryStationParameter.getBaseFrame().getName();
488 String secondaryStationName = secondaryStationParameter.getBaseFrame().getName();
489 System.out.format(Locale.US, "%-15s %-15s %-23s %-23s ",
490 primaryStationName, secondaryStationName, measurement.getDate(), date);
491 }
492
493
494 for (int i = 0; i < 6; ++i) {
495 final double[] gradient = measurement.estimate(0, 0, new SpacecraftState[] { state }).getParameterDerivatives(drivers[i], new AbsoluteDate());
496 Assertions.assertEquals(1, measurement.getDimension());
497 Assertions.assertEquals(1, gradient.length);
498
499
500 final ParameterFunction dMkdP =
501 Differentiation.differentiate(new ParameterFunction() {
502
503 @Override
504 public double value(final ParameterDriver parameterDriver, AbsoluteDate date) {
505 return measurement.
506 estimateWithoutDerivatives(new SpacecraftState[] { state }).
507 getEstimatedValue()[0];
508 }
509 }, 3, 20.0 * drivers[i].getScale());
510 final double ref = dMkdP.value(drivers[i], date);
511
512
513 double dGradient = gradient[0] - ref;
514 double dGradientRelative = FastMath.abs(dGradient/ref);
515
516
517 if (printResults) {
518 System.out.format(Locale.US, "%10.3e %10.3e ", dGradient, dGradientRelative);
519 }
520
521
522 if (i<3) { relErrorQMList.add(dGradientRelative);}
523 else { relErrorQSList.add(dGradientRelative);}
524
525 }
526 if (printResults) {
527 System.out.format(Locale.US, "%n");
528 }
529 }
530
531
532 final double[] relErrorQM = relErrorQMList.stream().mapToDouble(Double::doubleValue).toArray();
533 final double[] relErrorQS = relErrorQSList.stream().mapToDouble(Double::doubleValue).toArray();
534
535
536 final double relErrorsQMMedian = new Median().evaluate(relErrorQM);
537 final double relErrorsQMMean = new Mean().evaluate(relErrorQM);
538 final double relErrorsQMMax = new Max().evaluate(relErrorQM);
539
540 final double relErrorsQSMedian = new Median().evaluate(relErrorQS);
541 final double relErrorsQSMean = new Mean().evaluate(relErrorQS);
542 final double relErrorsQSMax = new Max().evaluate(relErrorQS);
543
544
545
546 if (printResults) {
547 System.out.println();
548 System.out.format(Locale.US, "Relative errors dR/dQ primary station -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n",
549 relErrorsQMMedian, relErrorsQMMean, relErrorsQMMax);
550 System.out.format(Locale.US, "Relative errors dR/dQ secondary station -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n",
551 relErrorsQSMedian, relErrorsQSMean, relErrorsQSMax);
552 }
553
554
555 Assertions.assertEquals(0.0, relErrorsQMMedian, refErrorQMMedian);
556 Assertions.assertEquals(0.0, relErrorsQMMean, refErrorQMMean);
557 Assertions.assertEquals(0.0, relErrorsQMMax, refErrorQMMax);
558 Assertions.assertEquals(0.0, relErrorsQSMedian, refErrorQSMedian);
559 Assertions.assertEquals(0.0, relErrorsQSMean, refErrorQSMean);
560 Assertions.assertEquals(0.0, relErrorsQSMax, refErrorQSMax);
561
562 }
563 }