1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.estimation.measurements.gnss;
18
19 import java.util.ArrayList;
20 import java.util.Comparator;
21 import java.util.List;
22 import java.util.Locale;
23
24 import org.hipparchus.geometry.euclidean.threed.Vector3D;
25 import org.hipparchus.stat.descriptive.moment.Mean;
26 import org.hipparchus.stat.descriptive.rank.Max;
27 import org.hipparchus.stat.descriptive.rank.Median;
28 import org.hipparchus.stat.descriptive.rank.Min;
29 import org.hipparchus.util.FastMath;
30 import org.junit.Assert;
31 import org.junit.Test;
32 import org.orekit.Utils;
33 import org.orekit.estimation.Context;
34 import org.orekit.estimation.EstimationTestUtils;
35 import org.orekit.estimation.measurements.EstimatedMeasurement;
36 import org.orekit.estimation.measurements.ObservableSatellite;
37 import org.orekit.estimation.measurements.ObservedMeasurement;
38 import org.orekit.gnss.Frequency;
39 import org.orekit.orbits.CartesianOrbit;
40 import org.orekit.orbits.Orbit;
41 import org.orekit.orbits.OrbitType;
42 import org.orekit.orbits.PositionAngle;
43 import org.orekit.propagation.BoundedPropagator;
44 import org.orekit.propagation.EphemerisGenerator;
45 import org.orekit.propagation.Propagator;
46 import org.orekit.propagation.SpacecraftState;
47 import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
48 import org.orekit.time.AbsoluteDate;
49 import org.orekit.utils.Constants;
50 import org.orekit.utils.Differentiation;
51 import org.orekit.utils.ParameterDriver;
52 import org.orekit.utils.ParameterFunction;
53 import org.orekit.utils.StateFunction;
54 import org.orekit.utils.TimeStampedPVCoordinates;
55
56 public class OneWayGNSSPhaseTest {
57
58 private static final Frequency FREQUENCY = Frequency.G01;
59
60
61
62
63
64 @Test
65 public void testValues() {
66 boolean printResults = false;
67 if (printResults) {
68 System.out.println("\nTest One-way GNSS phase Values\n");
69 }
70
71 this.genericTestValues(printResults);
72 }
73
74
75
76
77
78 @Test
79 public void testStateDerivatives() {
80
81 boolean printResults = false;
82 if (printResults) {
83 System.out.println("\nTest One-way GNSS phase State Derivatives - Finite Differences Comparison\n");
84 }
85
86 double refErrorsPMedian = 1.9e-10;
87 double refErrorsPMean = 5.3e-10;
88 double refErrorsPMax = 3.7e-08;
89 double refErrorsVMedian = 4.7e-04;
90 double refErrorsVMean = 1.6e-03;
91 double refErrorsVMax = 1.1e-01;
92 this.genericTestStateDerivatives(printResults, 0,
93 refErrorsPMedian, refErrorsPMean, refErrorsPMax,
94 refErrorsVMedian, refErrorsVMean, refErrorsVMax);
95 }
96
97
98
99
100
101 @Test
102 public void testParameterDerivatives() {
103
104
105 boolean printResults = false;
106
107 if (printResults) {
108 System.out.println("\nTest One-way GNSS phase Derivatives - Finite Differences Comparison\n");
109 }
110
111 double refErrorsMedian = 1.0e-15;
112 double refErrorsMean = 1.0e-15;
113 double refErrorsMax = 1.0e-15;
114 this.genericTestParameterDerivatives(printResults,
115 refErrorsMedian, refErrorsMean, refErrorsMax);
116
117 }
118
119
120
121
122
123 void genericTestValues(final boolean printResults) {
124
125 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
126
127 final NumericalPropagatorBuilder propagatorBuilder =
128 context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true,
129 1.0e-6, 60.0, 0.001);
130
131
132 final TimeStampedPVCoordinates original = context.initialOrbit.getPVCoordinates();
133 final Orbit closeOrbit = new CartesianOrbit(new TimeStampedPVCoordinates(context.initialOrbit.getDate(),
134 original.getPosition().add(new Vector3D(1000, 2000, 3000)),
135 original.getVelocity().add(new Vector3D(-0.03, 0.01, 0.02))),
136 context.initialOrbit.getFrame(),
137 context.initialOrbit.getMu());
138 final Propagator closePropagator = EstimationTestUtils.createPropagator(closeOrbit,
139 propagatorBuilder);
140 final EphemerisGenerator generator = closePropagator.getEphemerisGenerator();
141 closePropagator.propagate(context.initialOrbit.getDate().shiftedBy(3.5 * closeOrbit.getKeplerianPeriod()));
142 final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
143 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
144 propagatorBuilder);
145
146 final int ambiguity = 1234;
147 final double localClockOffset = 0.137e-6;
148 final double remoteClockOffset = 469.0e-6;
149 final List<ObservedMeasurement<?>> measurements =
150 EstimationTestUtils.createMeasurements(propagator,
151 new OneWayGNSSPhaseCreator(ephemeris, FREQUENCY, ambiguity, localClockOffset, remoteClockOffset),
152 1.0, 3.0, 300.0);
153
154
155
156 final List<Double> absoluteErrors = new ArrayList<Double>();
157 final List<Double> relativeErrors = new ArrayList<Double>();
158
159
160 propagator.setStepHandler(interpolator -> {
161
162 for (final ObservedMeasurement<?> measurement : measurements) {
163
164
165 if ((measurement.getDate().durationFrom(interpolator.getPreviousState().getDate()) > 0.) &&
166 (measurement.getDate().durationFrom(interpolator.getCurrentState().getDate()) <= 0.)
167 ) {
168
169
170
171
172
173 final double meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
174 final AbsoluteDate date = measurement.getDate().shiftedBy(-0.75 * meanDelay);
175 final SpacecraftState state = interpolator.getInterpolatedState(date);
176
177
178 final double phaseObserved = measurement.getObservedValue()[0];
179 final EstimatedMeasurement<?> estimated = measurement.estimate(0, 0,
180 new SpacecraftState[] {
181 state,
182 ephemeris.propagate(state.getDate())
183 });
184 Assert.assertEquals(FREQUENCY.getWavelength(), ((OneWayGNSSPhase) measurement).getWavelength(), 1.0e-15);
185 final double phaseEstimated = estimated.getEstimatedValue()[0];
186 final double absoluteError = phaseEstimated-phaseObserved;
187 absoluteErrors.add(absoluteError);
188 relativeErrors.add(FastMath.abs(absoluteError)/FastMath.abs(phaseObserved));
189
190
191 if (printResults) {
192 final AbsoluteDate measurementDate = measurement.getDate();
193
194 System.out.format(Locale.US, "%-23s %-23s %19.6f %19.6f %13.6e %13.6e%n",
195 measurementDate, date,
196 phaseObserved, phaseEstimated,
197 FastMath.abs(phaseEstimated-phaseObserved),
198 FastMath.abs((phaseEstimated-phaseObserved)/phaseObserved));
199 }
200
201 }
202 }
203 });
204
205
206 if (printResults) {
207 System.out.format(Locale.US, "%-23s %-23s %19s %19s %19s %19s%n",
208 "Measurement Date", "State Date",
209 "Phase observed [cycles]", "Phase estimated [cycles]",
210 "ΔPhase [cycles]", "rel ΔPhase");
211 }
212
213
214 propagator.propagate(context.initialOrbit.getDate());
215
216
217 measurements.sort(Comparator.naturalOrder());
218
219
220 propagator.propagate(measurements.get(measurements.size()-1).getDate());
221
222
223 final double[] absErrors = absoluteErrors.stream().mapToDouble(Double::doubleValue).toArray();
224 final double[] relErrors = relativeErrors.stream().mapToDouble(Double::doubleValue).toArray();
225
226
227 final double absErrorsMedian = new Median().evaluate(absErrors);
228 final double absErrorsMin = new Min().evaluate(absErrors);
229 final double absErrorsMax = new Max().evaluate(absErrors);
230 final double relErrorsMedian = new Median().evaluate(relErrors);
231 final double relErrorsMax = new Max().evaluate(relErrors);
232
233
234 if (printResults) {
235 System.out.println();
236 System.out.println("Absolute errors median: " + absErrorsMedian);
237 System.out.println("Absolute errors min : " + absErrorsMin);
238 System.out.println("Absolute errors max : " + absErrorsMax);
239 System.out.println("Relative errors median: " + relErrorsMedian);
240 System.out.println("Relative errors max : " + relErrorsMax);
241 }
242
243 Assert.assertEquals(0.0, absErrorsMedian, 6.5e-7);
244 Assert.assertEquals(0.0, absErrorsMin, 3.1e-6);
245 Assert.assertEquals(0.0, absErrorsMax, 9.0e-7);
246 Assert.assertEquals(0.0, relErrorsMedian, 5.9e-12);
247 Assert.assertEquals(0.0, relErrorsMax, 1.4e-10);
248 }
249
250 void genericTestStateDerivatives(final boolean printResults, final int index,
251 final double refErrorsPMedian, final double refErrorsPMean, final double refErrorsPMax,
252 final double refErrorsVMedian, final double refErrorsVMean, final double refErrorsVMax) {
253
254 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
255
256 final NumericalPropagatorBuilder propagatorBuilder =
257 context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true,
258 1.0e-6, 60.0, 0.001);
259
260
261 final TimeStampedPVCoordinates original = context.initialOrbit.getPVCoordinates();
262 final Orbit closeOrbit = new CartesianOrbit(new TimeStampedPVCoordinates(context.initialOrbit.getDate(),
263 original.getPosition().add(new Vector3D(1000, 2000, 3000)),
264 original.getVelocity().add(new Vector3D(-0.03, 0.01, 0.02))),
265 context.initialOrbit.getFrame(),
266 context.initialOrbit.getMu());
267 final Propagator closePropagator = EstimationTestUtils.createPropagator(closeOrbit,
268 propagatorBuilder);
269 final EphemerisGenerator generator = closePropagator.getEphemerisGenerator();
270 closePropagator.propagate(context.initialOrbit.getDate().shiftedBy(3.5 * closeOrbit.getKeplerianPeriod()));
271 final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
272 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
273 propagatorBuilder);
274
275 final int ambiguity = 1234;
276 final double localClockOffset = 0.137e-6;
277 final double remoteClockOffset = 469.0e-6;
278 final List<ObservedMeasurement<?>> measurements =
279 EstimationTestUtils.createMeasurements(propagator,
280 new OneWayGNSSPhaseCreator(ephemeris, FREQUENCY, ambiguity, localClockOffset, remoteClockOffset),
281 1.0, 3.0, 300.0);
282
283
284
285 final List<Double> errorsP = new ArrayList<Double>();
286 final List<Double> errorsV = new ArrayList<Double>();
287
288
289 propagator.setStepHandler(interpolator -> {
290
291 for (final ObservedMeasurement<?> measurement : measurements) {
292
293
294 if ((measurement.getDate().durationFrom(interpolator.getPreviousState().getDate()) > 0.) &&
295 (measurement.getDate().durationFrom(interpolator.getCurrentState().getDate()) <= 0.)
296 ) {
297
298
299
300
301
302
303 final double meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
304 final AbsoluteDate date = measurement.getDate().shiftedBy(-0.75 * meanDelay);
305 final SpacecraftState[] states = {
306 interpolator.getInterpolatedState(date),
307 ephemeris.propagate(date)
308 };
309 final double[][] jacobian = measurement.estimate(0, 0, states).getStateDerivatives(index);
310
311
312 final double[][] jacobianRef;
313
314
315 jacobianRef = Differentiation.differentiate(new StateFunction() {
316 public double[] value(final SpacecraftState state) {
317 final SpacecraftState[] s = states.clone();
318 s[index] = state;
319 return measurement.estimate(0, 0, s).getEstimatedValue();
320 }
321 }, measurement.getDimension(), propagator.getAttitudeProvider(),
322 OrbitType.CARTESIAN, PositionAngle.TRUE, 2.0, 3).value(states[index]);
323
324 Assert.assertEquals(jacobianRef.length, jacobian.length);
325 Assert.assertEquals(jacobianRef[0].length, jacobian[0].length);
326
327
328 double [][] dJacobian = new double[jacobian.length][jacobian[0].length];
329 double [][] dJacobianRelative = new double[jacobian.length][jacobian[0].length];
330 for (int i = 0; i < jacobian.length; ++i) {
331 for (int j = 0; j < jacobian[i].length; ++j) {
332 dJacobian[i][j] = jacobian[i][j] - jacobianRef[i][j];
333 dJacobianRelative[i][j] = FastMath.abs(dJacobian[i][j]/jacobianRef[i][j]);
334
335 if (j < 3) { errorsP.add(dJacobianRelative[i][j]);
336 } else { errorsV.add(dJacobianRelative[i][j]); }
337 }
338 }
339
340 if (printResults) {
341 System.out.format(Locale.US, "%-23s %-23s " +
342 "%10.3e %10.3e %10.3e " +
343 "%10.3e %10.3e %10.3e " +
344 "%10.3e %10.3e %10.3e " +
345 "%10.3e %10.3e %10.3e%n",
346 measurement.getDate(), date,
347 dJacobian[0][0], dJacobian[0][1], dJacobian[0][2],
348 dJacobian[0][3], dJacobian[0][4], dJacobian[0][5],
349 dJacobianRelative[0][0], dJacobianRelative[0][1], dJacobianRelative[0][2],
350 dJacobianRelative[0][3], dJacobianRelative[0][4], dJacobianRelative[0][5]);
351 }
352 }
353 }
354 });
355
356
357 if (printResults) {
358 System.out.format(Locale.US, "%-23s %-23s " +
359 "%10s %10s %10s " +
360 "%10s %10s %10s " +
361 "%10s %10s %10s " +
362 "%10s %10s %10s%n",
363 "Measurement Date", "State Date",
364 "ΔdPx", "ΔdPy", "ΔdPz", "ΔdVx", "ΔdVy", "ΔdVz",
365 "rel ΔdPx", "rel ΔdPy", "rel ΔdPz",
366 "rel ΔdVx", "rel ΔdVy", "rel ΔdVz");
367 }
368
369
370 propagator.propagate(context.initialOrbit.getDate());
371
372
373 measurements.sort(Comparator.naturalOrder());
374
375
376 propagator.propagate(measurements.get(measurements.size()-1).getDate());
377
378
379 final double relErrorsP[] = errorsP.stream().mapToDouble(Double::doubleValue).toArray();
380 final double relErrorsV[] = errorsV.stream().mapToDouble(Double::doubleValue).toArray();
381
382 final double errorsPMedian = new Median().evaluate(relErrorsP);
383 final double errorsPMean = new Mean().evaluate(relErrorsP);
384 final double errorsPMax = new Max().evaluate(relErrorsP);
385 final double errorsVMedian = new Median().evaluate(relErrorsV);
386 final double errorsVMean = new Mean().evaluate(relErrorsV);
387 final double errorsVMax = new Max().evaluate(relErrorsV);
388
389
390 if (printResults) {
391 System.out.println();
392 System.out.format(Locale.US, "Relative errors dR/dP -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n",
393 errorsPMedian, errorsPMean, errorsPMax);
394 System.out.format(Locale.US, "Relative errors dR/dV -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n",
395 errorsVMedian, errorsVMean, errorsVMax);
396 }
397
398 Assert.assertEquals(0.0, errorsPMedian, refErrorsPMedian);
399 Assert.assertEquals(0.0, errorsPMean, refErrorsPMean);
400 Assert.assertEquals(0.0, errorsPMax, refErrorsPMax);
401 Assert.assertEquals(0.0, errorsVMedian, refErrorsVMedian);
402 Assert.assertEquals(0.0, errorsVMean, refErrorsVMean);
403 Assert.assertEquals(0.0, errorsVMax, refErrorsVMax);
404 }
405
406 void genericTestParameterDerivatives(final boolean printResults,
407 final double refErrorsMedian, final double refErrorsMean, final double refErrorsMax) {
408
409 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
410
411 final NumericalPropagatorBuilder propagatorBuilder =
412 context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true,
413 1.0e-6, 60.0, 0.001);
414
415
416 final TimeStampedPVCoordinates original = context.initialOrbit.getPVCoordinates();
417 final Orbit closeOrbit = new CartesianOrbit(new TimeStampedPVCoordinates(context.initialOrbit.getDate(),
418 original.getPosition().add(new Vector3D(1000, 2000, 3000)),
419 original.getVelocity().add(new Vector3D(-0.03, 0.01, 0.02))),
420 context.initialOrbit.getFrame(),
421 context.initialOrbit.getMu());
422 final Propagator closePropagator = EstimationTestUtils.createPropagator(closeOrbit, propagatorBuilder);
423 final EphemerisGenerator generator = closePropagator.getEphemerisGenerator();
424 closePropagator.propagate(context.initialOrbit.getDate().shiftedBy(3.5 * closeOrbit.getKeplerianPeriod()));
425 final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
426
427
428 final int ambiguity = 1234;
429 final double localClockOffset = 0.137e-6;
430 final double remoteClockOffset = 469.0e-6;
431 final OneWayGNSSPhaseCreator creator = new OneWayGNSSPhaseCreator(ephemeris, FREQUENCY, ambiguity, localClockOffset, remoteClockOffset);
432 creator.getLocalSatellite().getClockOffsetDriver().setSelected(true);
433
434 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
435 propagatorBuilder);
436 final List<ObservedMeasurement<?>> measurements =
437 EstimationTestUtils.createMeasurements(propagator, creator, 1.0, 3.0, 300.0);
438
439
440 final List<Double> relErrorList = new ArrayList<Double>();
441
442
443 propagator.setStepHandler(interpolator -> {
444
445 for (final ObservedMeasurement<?> measurement : measurements) {
446
447
448 if ((measurement.getDate().durationFrom(interpolator.getPreviousState().getDate()) > 0.) &&
449 (measurement.getDate().durationFrom(interpolator.getCurrentState().getDate()) <= 0.)) {
450
451
452
453
454
455
456
457
458 final double meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
459 final AbsoluteDate date = measurement.getDate().shiftedBy(-0.75 * meanDelay);
460 final SpacecraftState[] states = {
461 interpolator.getInterpolatedState(date),
462 ephemeris.propagate(date)
463 };
464 final ParameterDriver[] drivers = new ParameterDriver[] {
465 measurement.getSatellites().get(0).getClockOffsetDriver(),
466 };
467
468 for (int i = 0; i < drivers.length; ++i) {
469 final double[] gradient = measurement.estimate(0, 0, states).getParameterDerivatives(drivers[i]);
470 Assert.assertEquals(1, measurement.getDimension());
471 Assert.assertEquals(1, gradient.length);
472
473
474 final ParameterFunction dMkdP =
475 Differentiation.differentiate(new ParameterFunction() {
476
477 @Override
478 public double value(final ParameterDriver parameterDriver) {
479 return measurement.estimate(0, 0, states).getEstimatedValue()[0];
480 }
481 }, 3, 20.0 * drivers[i].getScale());
482 final double ref = dMkdP.value(drivers[i]);
483
484 if (printResults) {
485 System.out.format(Locale.US, "%10.3e %10.3e ", gradient[0]-ref, FastMath.abs((gradient[0]-ref)/ref));
486 }
487
488 final double relError = FastMath.abs((ref-gradient[0])/ref);
489 relErrorList.add(relError);
490 }
491 if (printResults) {
492 System.out.format(Locale.US, "%n");
493 }
494
495 }
496 }
497 });
498
499
500 propagator.propagate(context.initialOrbit.getDate());
501
502
503 measurements.sort(Comparator.naturalOrder());
504
505
506 if (printResults) {
507 System.out.format(Locale.US, "%-15s %-23s %-23s " +
508 "%10s %10s %10s %10s %10s %10s %10s %10s %10s %10s%n",
509 "Station", "Measurement Date", "State Date",
510 "Δt", "rel Δt",
511 "ΔdQx", "rel ΔdQx",
512 "ΔdQy", "rel ΔdQy",
513 "ΔdQz", "rel ΔdQz",
514 "Δtsat", "rel Δtsat");
515 }
516
517
518 propagator.propagate(measurements.get(measurements.size()-1).getDate());
519
520
521 final double relErrors[] = relErrorList.stream().mapToDouble(Double::doubleValue).toArray();
522
523
524 final double relErrorsMedian = new Median().evaluate(relErrors);
525 final double relErrorsMean = new Mean().evaluate(relErrors);
526 final double relErrorsMax = new Max().evaluate(relErrors);
527
528
529 if (printResults) {
530 System.out.println();
531 System.out.format(Locale.US, "Relative errors dR/dQ -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n",
532 relErrorsMedian, relErrorsMean, relErrorsMax);
533 }
534
535 Assert.assertEquals(0.0, relErrorsMedian, refErrorsMedian);
536 Assert.assertEquals(0.0, relErrorsMean, refErrorsMean);
537 Assert.assertEquals(0.0, relErrorsMax, refErrorsMax);
538
539 }
540
541 @Test
542 public void testIssue734() {
543
544 Utils.setDataRoot("regular-data");
545
546
547 final OneWayGNSSPhase phase = new OneWayGNSSPhase(null, 635.0e-6, AbsoluteDate.J2000_EPOCH, 467614.701,
548 Frequency.G01.getWavelength(), 0.02, 1.0, new ObservableSatellite(0));
549
550
551 Assert.assertEquals(0.0, phase.getAmbiguityDriver().getValue(), Double.MIN_VALUE);
552 Assert.assertFalse(phase.getAmbiguityDriver().isSelected());
553
554
555 phase.getAmbiguityDriver().setValue(1234.0);
556 phase.getAmbiguityDriver().setSelected(true);
557
558
559 Assert.assertEquals(1234.0, phase.getAmbiguityDriver().getValue(), Double.MIN_VALUE);
560 Assert.assertTrue(phase.getAmbiguityDriver().isSelected());
561 for (ParameterDriver driver : phase.getParametersDrivers()) {
562
563 if (driver.getName() == Phase.AMBIGUITY_NAME) {
564 Assert.assertEquals(1234.0, phase.getAmbiguityDriver().getValue(), Double.MIN_VALUE);
565 Assert.assertTrue(phase.getAmbiguityDriver().isSelected());
566 }
567 }
568
569 }
570
571 }