1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package org.orekit.estimation.leastsquares;
18
19 import org.hipparchus.exception.LocalizedCoreFormats;
20 import org.hipparchus.geometry.euclidean.threed.Vector3D;
21 import org.hipparchus.linear.RealMatrix;
22 import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem.Evaluation;
23 import org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer;
24 import org.junit.jupiter.api.Assertions;
25 import org.junit.jupiter.api.Test;
26 import org.orekit.attitudes.LofOffset;
27 import org.orekit.errors.OrekitException;
28 import org.orekit.errors.OrekitMessages;
29 import org.orekit.estimation.Context;
30 import org.orekit.estimation.EstimationTestUtils;
31 import org.orekit.estimation.Force;
32 import org.orekit.estimation.measurements.AngularAzElMeasurementCreator;
33 import org.orekit.estimation.measurements.EstimationsProvider;
34 import org.orekit.estimation.measurements.GroundStation;
35 import org.orekit.estimation.measurements.InterSatellitesRangeMeasurementCreator;
36 import org.orekit.estimation.measurements.MultiplexedMeasurement;
37 import org.orekit.estimation.measurements.ObservedMeasurement;
38 import org.orekit.estimation.measurements.PVMeasurementCreator;
39 import org.orekit.estimation.measurements.Range;
40 import org.orekit.estimation.measurements.TwoWayRangeMeasurementCreator;
41 import org.orekit.estimation.measurements.RangeRateMeasurementCreator;
42 import org.orekit.estimation.measurements.modifiers.PhaseCentersRangeModifier;
43 import org.orekit.forces.gravity.NewtonianAttraction;
44 import org.orekit.forces.radiation.RadiationSensitive;
45 import org.orekit.frames.LOFType;
46 import org.orekit.gnss.antenna.FrequencyPattern;
47 import org.orekit.orbits.CartesianOrbit;
48 import org.orekit.orbits.KeplerianOrbit;
49 import org.orekit.orbits.Orbit;
50 import org.orekit.orbits.OrbitType;
51 import org.orekit.orbits.PositionAngleType;
52 import org.orekit.propagation.BoundedPropagator;
53 import org.orekit.propagation.EphemerisGenerator;
54 import org.orekit.propagation.Propagator;
55 import org.orekit.propagation.conversion.NumericalPropagatorBuilder;
56 import org.orekit.propagation.numerical.NumericalPropagator;
57 import org.orekit.time.AbsoluteDate;
58 import org.orekit.time.ChronologicalComparator;
59 import org.orekit.utils.ParameterDriver;
60 import org.orekit.utils.ParameterDriversList;
61 import org.orekit.utils.ParameterDriversList.DelegatingDriver;
62 import org.orekit.utils.TimeStampedPVCoordinates;
63
64 import java.util.ArrayList;
65 import java.util.List;
66 import java.util.stream.IntStream;
67
68 public class BatchLSEstimatorTest {
69
70
71
72
73 @Test
74 public void testKeplerPVMultipleDrag() {
75
76 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
77
78 final NumericalPropagatorBuilder propagatorBuilder =
79 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
80 1.0e-6, 60.0, 1.0, Force.DRAG);
81
82 for (ParameterDriver driver:propagatorBuilder.getPropagationParametersDrivers().getDrivers()) {
83 if (driver.getName().equals("drag coefficient")) {
84 driver.setSelected(true);
85 driver.addSpanAtDate(context.initialOrbit.getDate());
86 }
87 }
88
89
90 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
91 propagatorBuilder);
92 final List<ObservedMeasurement<?>> measurements =
93 EstimationTestUtils.createMeasurements(propagator,
94 new PVMeasurementCreator(),
95 -3.0, 3.0, 300.0);
96
97
98
99
100 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
101 propagatorBuilder);
102 for (final ObservedMeasurement<?> measurement : measurements) {
103 estimator.addMeasurement(measurement);
104 }
105 estimator.setParametersConvergenceThreshold(1.0e-2);
106 estimator.setMaxIterations(10);
107 estimator.setMaxEvaluations(20);
108
109 EstimationTestUtils.checkFit(context, estimator, 1, 2,
110 0.0, 7.8e-8,
111 0.0, 6.0e-7,
112 0.0, 3.2e-7,
113 0.0, 1.3e-10);
114
115
116 List<DelegatingDriver> Orbparameters = estimator.getOrbitalParametersDrivers(true).getDrivers();
117 Assertions.assertEquals(context.initialOrbit.getA(), Orbparameters.get(0).getValue() , 1.0e-8);
118 Assertions.assertEquals(context.initialOrbit.getE(), Orbparameters.get(1).getValue() , 1.0e-12);
119 Assertions.assertEquals(context.initialOrbit.getI(), Orbparameters.get(2).getValue() , 1.0e-12);
120
121 RealMatrix jacobian = estimator.getOptimum().getJacobian();
122 Assertions.assertEquals(8, jacobian.getColumnDimension());
123
124 }
125
126
127
128
129 @Test
130 public void testKeplerPV() {
131
132 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
133
134 final NumericalPropagatorBuilder propagatorBuilder =
135 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
136 1.0e-6, 60.0, 1.0);
137
138
139 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
140 propagatorBuilder);
141 final List<ObservedMeasurement<?>> measurements =
142 EstimationTestUtils.createMeasurements(propagator,
143 new PVMeasurementCreator(),
144 0.0, 1.0, 300.0);
145
146
147 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
148 propagatorBuilder);
149 for (final ObservedMeasurement<?> measurement : measurements) {
150 estimator.addMeasurement(measurement);
151 }
152 estimator.setParametersConvergenceThreshold(1.0e-2);
153 estimator.setMaxIterations(10);
154 estimator.setMaxEvaluations(20);
155
156 EstimationTestUtils.checkFit(context, estimator, 1, 4,
157 0.0, 2.2e-8,
158 0.0, 1.1e-7,
159 0.0, 1.4e-8,
160 0.0, 6.3e-12);
161
162 RealMatrix normalizedCovariances = estimator.getOptimum().getCovariances(1.0e-10);
163 RealMatrix physicalCovariances = estimator.getPhysicalCovariances(1.0e-10);
164 Assertions.assertEquals(6, normalizedCovariances.getRowDimension());
165 Assertions.assertEquals(6, normalizedCovariances.getColumnDimension());
166 Assertions.assertEquals(6, physicalCovariances.getRowDimension());
167 Assertions.assertEquals(6, physicalCovariances.getColumnDimension());
168 Assertions.assertEquals(0.00258, physicalCovariances.getEntry(0, 0), 1.0e-5);
169
170 }
171
172
173 @Test
174 public void testKeplerPVBackward() {
175
176 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
177
178 final NumericalPropagatorBuilder propagatorBuilder =
179 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
180 1.0e-6, 60.0, 1.0);
181
182
183 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
184 propagatorBuilder);
185 final List<ObservedMeasurement<?>> measurements =
186 EstimationTestUtils.createMeasurements(propagator,
187 new PVMeasurementCreator(),
188 0.0, -1.0, 300.0);
189
190
191 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
192 propagatorBuilder);
193 for (final ObservedMeasurement<?> measurement : measurements) {
194 estimator.addMeasurement(measurement);
195 }
196 estimator.setParametersConvergenceThreshold(1.0e-2);
197 estimator.setMaxIterations(10);
198 estimator.setMaxEvaluations(20);
199
200 EstimationTestUtils.checkFit(context, estimator, 1, 2,
201 0.0, 8.3e-9,
202 0.0, 5.3e-8,
203 0.0, 5.6e-9,
204 0.0, 1.6e-12);
205
206 RealMatrix normalizedCovariances = estimator.getOptimum().getCovariances(1.0e-10);
207 RealMatrix physicalCovariances = estimator.getPhysicalCovariances(1.0e-10);
208 Assertions.assertEquals(6, normalizedCovariances.getRowDimension());
209 Assertions.assertEquals(6, normalizedCovariances.getColumnDimension());
210 Assertions.assertEquals(6, physicalCovariances.getRowDimension());
211 Assertions.assertEquals(6, physicalCovariances.getColumnDimension());
212 Assertions.assertEquals(0.00258, physicalCovariances.getEntry(0, 0), 1.0e-5);
213
214 }
215
216
217
218
219 @Test
220 public void testKeplerRange() {
221
222 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
223
224 final NumericalPropagatorBuilder propagatorBuilder =
225 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
226 1.0e-6, 60.0, 1.0);
227
228
229 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
230 propagatorBuilder);
231 final List<ObservedMeasurement<?>> measurements =
232 EstimationTestUtils.createMeasurements(propagator,
233 new TwoWayRangeMeasurementCreator(context),
234 1.0, 3.0, 300.0);
235
236
237 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
238 propagatorBuilder);
239 for (final ObservedMeasurement<?> range : measurements) {
240 estimator.addMeasurement(range);
241 }
242 estimator.setParametersConvergenceThreshold(1.0e-2);
243 estimator.setMaxIterations(10);
244 estimator.setMaxEvaluations(20);
245 estimator.setObserver(new BatchLSObserver() {
246 int lastIter = 0;
247 int lastEval = 0;
248
249 @Override
250 public void evaluationPerformed(int iterationsCount, int evaluationscount,
251 Orbit[] orbits,
252 ParameterDriversList estimatedOrbitalParameters,
253 ParameterDriversList estimatedPropagatorParameters,
254 ParameterDriversList estimatedMeasurementsParameters,
255 EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) {
256 if (iterationsCount == lastIter) {
257 Assertions.assertEquals(lastEval + 1, evaluationscount);
258 } else {
259 Assertions.assertEquals(lastIter + 1, iterationsCount);
260 }
261 lastIter = iterationsCount;
262 lastEval = evaluationscount;
263 Assertions.assertEquals(measurements.size(), evaluationsProvider.getNumber());
264 try {
265 evaluationsProvider.getEstimatedMeasurement(-1);
266 Assertions.fail("an exception should have been thrown");
267 } catch (OrekitException oe) {
268 Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
269 }
270 try {
271 evaluationsProvider.getEstimatedMeasurement(measurements.size());
272 Assertions.fail("an exception should have been thrown");
273 } catch (OrekitException oe) {
274 Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
275 }
276 AbsoluteDate previous = AbsoluteDate.PAST_INFINITY;
277 for (int i = 0; i < evaluationsProvider.getNumber(); ++i) {
278 AbsoluteDate current = evaluationsProvider.getEstimatedMeasurement(i).getDate();
279 Assertions.assertTrue(current.compareTo(previous) >= 0);
280 previous = current;
281 }
282 }
283 });
284
285 ParameterDriver aDriver = estimator.getOrbitalParametersDrivers(true).getDrivers().get(0);
286 Assertions.assertEquals("a", aDriver.getName());
287 aDriver.setValue(aDriver.getValue() + 1.2);
288 aDriver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
289
290 EstimationTestUtils.checkFit(context, estimator, 2, 3,
291 0.0, 1.2e-6,
292 0.0, 2.8e-6,
293 0.0, 5.0e-7,
294 0.0, 2.3e-10);
295
296
297
298 for (final ParameterDriver driver : estimator.getOrbitalParametersDrivers(true).getDrivers()) {
299 if ("a".equals(driver.getName())) {
300
301 Assertions.assertEquals(0, driver.getReferenceDate().durationFrom(AbsoluteDate.GALILEO_EPOCH), 1.0e-15);
302 } else {
303
304 Assertions.assertEquals(0,
305 driver.getReferenceDate().durationFrom(propagatorBuilder.getInitialOrbitDate()), 1.0e-15);
306 }
307 }
308
309 }
310
311
312
313
314 @Test
315 public void testKeplerRangeWithOnBoardAntennaOffset() {
316
317 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
318
319 final NumericalPropagatorBuilder propagatorBuilder =
320 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
321 1.0e-6, 60.0, 1.0);
322 propagatorBuilder.setAttitudeProvider(new LofOffset(propagatorBuilder.getFrame(), LOFType.LVLH));
323 final Vector3D antennaPhaseCenter = new Vector3D(-1.2, 2.3, -0.7);
324
325
326 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
327 propagatorBuilder);
328 final List<ObservedMeasurement<?>> measurements =
329 EstimationTestUtils.createMeasurements(propagator,
330 new TwoWayRangeMeasurementCreator(context,
331 Vector3D.ZERO, null,
332 antennaPhaseCenter, null,
333 0),
334 1.0, 3.0, 300.0);
335
336
337 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
338 propagatorBuilder);
339 final PhaseCentersRangeModifier obaModifier = new PhaseCentersRangeModifier(FrequencyPattern.ZERO_CORRECTION,
340 new FrequencyPattern(antennaPhaseCenter,
341 null));
342 for (final ObservedMeasurement<?> range : measurements) {
343 ((Range) range).addModifier(obaModifier);
344 estimator.addMeasurement(range);
345 }
346 estimator.setParametersConvergenceThreshold(1.0e-2);
347 estimator.setMaxIterations(10);
348 estimator.setMaxEvaluations(20);
349 estimator.setObserver(new BatchLSObserver() {
350 int lastIter = 0;
351 int lastEval = 0;
352
353 @Override
354 public void evaluationPerformed(int iterationsCount, int evaluationscount,
355 Orbit[] orbits,
356 ParameterDriversList estimatedOrbitalParameters,
357 ParameterDriversList estimatedPropagatorParameters,
358 ParameterDriversList estimatedMeasurementsParameters,
359 EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) {
360 if (iterationsCount == lastIter) {
361 Assertions.assertEquals(lastEval + 1, evaluationscount);
362 } else {
363 Assertions.assertEquals(lastIter + 1, iterationsCount);
364 }
365 lastIter = iterationsCount;
366 lastEval = evaluationscount;
367 Assertions.assertEquals(measurements.size(), evaluationsProvider.getNumber());
368 try {
369 evaluationsProvider.getEstimatedMeasurement(-1);
370 Assertions.fail("an exception should have been thrown");
371 } catch (OrekitException oe) {
372 Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
373 }
374 try {
375 evaluationsProvider.getEstimatedMeasurement(measurements.size());
376 Assertions.fail("an exception should have been thrown");
377 } catch (OrekitException oe) {
378 Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
379 }
380 AbsoluteDate previous = AbsoluteDate.PAST_INFINITY;
381 for (int i = 0; i < evaluationsProvider.getNumber(); ++i) {
382 AbsoluteDate current = evaluationsProvider.getEstimatedMeasurement(i).getDate();
383 Assertions.assertTrue(current.compareTo(previous) >= 0);
384 previous = current;
385 }
386 }
387 });
388
389 ParameterDriver aDriver = estimator.getOrbitalParametersDrivers(true).getDrivers().get(0);
390 Assertions.assertEquals("a", aDriver.getName());
391 aDriver.setValue(aDriver.getValue() + 1.2);
392 aDriver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
393
394 EstimationTestUtils.checkFit(context, estimator, 2, 3,
395 0.0, 2.0e-5,
396 0.0, 5.2e-5,
397 0.0, 2.7e-5,
398 0.0, 1.1e-8);
399
400
401
402 for (final ParameterDriver driver : estimator.getOrbitalParametersDrivers(true).getDrivers()) {
403 if ("a".equals(driver.getName())) {
404
405 Assertions.assertEquals(0, driver.getReferenceDate().durationFrom(AbsoluteDate.GALILEO_EPOCH), 1.0e-15);
406 } else {
407
408 Assertions.assertEquals(0,
409 driver.getReferenceDate().durationFrom(propagatorBuilder.getInitialOrbitDate()), 1.0e-15);
410 }
411 }
412
413 }
414
415 @Test
416 public void testMultiSat() {
417
418 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
419
420 final NumericalPropagatorBuilder propagatorBuilder1 =
421 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
422 1.0e-6, 60.0, 1.0e-3);
423 final NumericalPropagatorBuilder propagatorBuilder2 =
424 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
425 1.0e-6, 60.0, 1.0e-3);
426
427
428 final TimeStampedPVCoordinates original = context.initialOrbit.getPVCoordinates();
429 final Orbit closeOrbit = new CartesianOrbit(new TimeStampedPVCoordinates(context.initialOrbit.getDate(),
430 original.getPosition().add(new Vector3D(1000, 2000, 3000)),
431 original.getVelocity().add(new Vector3D(-0.03, 0.01, 0.02))),
432 context.initialOrbit.getFrame(),
433 context.initialOrbit.getMu());
434 final Propagator closePropagator = EstimationTestUtils.createPropagator(closeOrbit,
435 propagatorBuilder2);
436 final EphemerisGenerator generator = closePropagator.getEphemerisGenerator();
437 closePropagator.propagate(context.initialOrbit.getDate().shiftedBy(3.5 * closeOrbit.getKeplerianPeriod()));
438 final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
439 Propagator propagator1 = EstimationTestUtils.createPropagator(context.initialOrbit,
440 propagatorBuilder1);
441
442 final double localClockOffset = 0.137e-6;
443 final double remoteClockOffset = 469.0e-6;
444 final List<ObservedMeasurement<?>> r12 =
445 EstimationTestUtils.createMeasurements(propagator1,
446 new InterSatellitesRangeMeasurementCreator(ephemeris,
447 localClockOffset,
448 remoteClockOffset),
449 1.0, 3.0, 300.0);
450
451
452 propagator1 = EstimationTestUtils.createPropagator(context.initialOrbit,
453 propagatorBuilder1);
454 final List<ObservedMeasurement<?>> r1 =
455 EstimationTestUtils.createMeasurements(propagator1,
456 new TwoWayRangeMeasurementCreator(context),
457 1.0, 3.0, 300.0);
458
459
460 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
461 propagatorBuilder1,
462 propagatorBuilder2);
463 for (final ObservedMeasurement<?> interSat : r12) {
464 estimator.addMeasurement(interSat);
465 }
466 for (final ObservedMeasurement<?> range : r1) {
467 estimator.addMeasurement(range);
468 }
469 estimator.setParametersConvergenceThreshold(1.0e-2);
470 estimator.setMaxIterations(10);
471 estimator.setMaxEvaluations(20);
472 estimator.setObserver(new BatchLSObserver() {
473 int lastIter = 0;
474 int lastEval = 0;
475
476 @Override
477 public void evaluationPerformed(int iterationsCount, int evaluationscount,
478 Orbit[] orbits,
479 ParameterDriversList estimatedOrbitalParameters,
480 ParameterDriversList estimatedPropagatorParameters,
481 ParameterDriversList estimatedMeasurementsParameters,
482 EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) {
483 if (iterationsCount == lastIter) {
484 Assertions.assertEquals(lastEval + 1, evaluationscount);
485 } else {
486 Assertions.assertEquals(lastIter + 1, iterationsCount);
487 }
488 lastIter = iterationsCount;
489 lastEval = evaluationscount;
490 Assertions.assertEquals(r12.size() + r1.size(), evaluationsProvider.getNumber());
491 try {
492 evaluationsProvider.getEstimatedMeasurement(-1);
493 Assertions.fail("an exception should have been thrown");
494 } catch (OrekitException oe) {
495 Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
496 }
497 try {
498 evaluationsProvider.getEstimatedMeasurement(r12.size() + r1.size());
499 Assertions.fail("an exception should have been thrown");
500 } catch (OrekitException oe) {
501 Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
502 }
503 AbsoluteDate previous = AbsoluteDate.PAST_INFINITY;
504 for (int i = 0; i < evaluationsProvider.getNumber(); ++i) {
505 AbsoluteDate current = evaluationsProvider.getEstimatedMeasurement(i).getDate();
506 Assertions.assertTrue(current.compareTo(previous) >= 0);
507 previous = current;
508 }
509 }
510 });
511
512 List<DelegatingDriver> parameters = estimator.getOrbitalParametersDrivers(true).getDrivers();
513 ParameterDriver a0Driver = parameters.get(0);
514 Assertions.assertEquals("a[0]", a0Driver.getName());
515 a0Driver.setValue(a0Driver.getValue() + 1.2);
516 a0Driver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
517
518 ParameterDriver a1Driver = parameters.get(6);
519 Assertions.assertEquals("a[1]", a1Driver.getName());
520 a1Driver.setValue(a1Driver.getValue() - 5.4);
521 a1Driver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
522
523 final Orbit before = new KeplerianOrbit(parameters.get( 6).getValue(),
524 parameters.get( 7).getValue(),
525 parameters.get( 8).getValue(),
526 parameters.get( 9).getValue(),
527 parameters.get(10).getValue(),
528 parameters.get(11).getValue(),
529 PositionAngleType.TRUE,
530 closeOrbit.getFrame(),
531 closeOrbit.getDate(),
532 closeOrbit.getMu());
533 Assertions.assertEquals(4.7246, Vector3D.distance(closeOrbit.getPosition(),
534 before.getPosition()), 1.0e-3);
535 Assertions.assertEquals(0.0010514, Vector3D.distance(closeOrbit.getPVCoordinates().getVelocity(),
536 before.getPVCoordinates().getVelocity()), 1.0e-6);
537 EstimationTestUtils.checkFit(context, estimator, 3, 4,
538 0.0, 2.9e-06,
539 0.0, 1.1e-05,
540 0.0, 8.3e-07,
541 0.0, 3.7e-10);
542
543 final Orbit determined = new KeplerianOrbit(parameters.get( 6).getValue(),
544 parameters.get( 7).getValue(),
545 parameters.get( 8).getValue(),
546 parameters.get( 9).getValue(),
547 parameters.get(10).getValue(),
548 parameters.get(11).getValue(),
549 PositionAngleType.TRUE,
550 closeOrbit.getFrame(),
551 closeOrbit.getDate(),
552 closeOrbit.getMu());
553 Assertions.assertEquals(0.0, Vector3D.distance(closeOrbit.getPosition(),
554 determined.getPosition()), 6.2e-6);
555 Assertions.assertEquals(0.0, Vector3D.distance(closeOrbit.getPVCoordinates().getVelocity(),
556 determined.getPVCoordinates().getVelocity()), 1.6e-9);
557
558
559
560 for (final ParameterDriver driver : estimator.getOrbitalParametersDrivers(true).getDrivers()) {
561 if (driver.getName().startsWith("a[")) {
562
563 Assertions.assertEquals(0, driver.getReferenceDate().durationFrom(AbsoluteDate.GALILEO_EPOCH), 1.0e-15);
564 } else {
565
566 Assertions.assertEquals(0,
567 driver.getReferenceDate().durationFrom(propagatorBuilder1.getInitialOrbitDate()), 1.0e-15);
568 }
569 }
570
571 }
572
573
574
575
576
577
578 @Test
579 public void testMultiSatWithParameters() {
580
581
582 final boolean muEstimated = true;
583 final boolean crEstimated1 = false;
584 final boolean caEstimated1 = true;
585 final boolean crEstimated2 = true;
586 final boolean caEstimated2 = false;
587
588
589
590 final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
591 final NumericalPropagatorBuilder propagatorBuilder1 =
592 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
593 1.0e-6, 60.0, 1.0e-3, Force.POTENTIAL, Force.SOLAR_RADIATION_PRESSURE);
594
595
596 String satName = "sat 1";
597 for (DelegatingDriver driver:propagatorBuilder1.getPropagationParametersDrivers().getDrivers()) {
598 if (driver.getName().equals("central attraction coefficient")) {
599 driver.setSelected(muEstimated);
600 }
601 if (driver.getName().equals(RadiationSensitive.REFLECTION_COEFFICIENT)) {
602 driver.setName(driver.getName() + " " + satName);
603 driver.setSelected(crEstimated1);
604 }
605 if (driver.getName().equals(RadiationSensitive.ABSORPTION_COEFFICIENT)) {
606 driver.setName(driver.getName() + " " + satName);
607 driver.setSelected(caEstimated1);
608 }
609 }
610
611
612 final Context context2 = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
613 final NumericalPropagatorBuilder propagatorBuilder2 =
614 context2.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
615 1.0e-6, 60.0, 1.0e-3, Force.POTENTIAL, Force.SOLAR_RADIATION_PRESSURE);
616
617
618 satName = "sat 2";
619 for (ParameterDriver driver:propagatorBuilder2.getPropagationParametersDrivers().getDrivers()) {
620 if (driver.getName().equals("central attraction coefficient")) {
621 driver.setSelected(muEstimated);
622 }
623 if (driver.getName().equals(RadiationSensitive.REFLECTION_COEFFICIENT)) {
624 driver.setName(driver.getName() + " " + satName);
625 driver.setSelected(crEstimated2);
626 }
627 if (driver.getName().equals(RadiationSensitive.ABSORPTION_COEFFICIENT)) {
628 driver.setName(driver.getName() + " " + satName);
629 driver.setSelected(caEstimated2);
630 }
631 }
632
633
634 final TimeStampedPVCoordinates original = context.initialOrbit.getPVCoordinates();
635 final Orbit closeOrbit = new CartesianOrbit(new TimeStampedPVCoordinates(context.initialOrbit.getDate(),
636 original.getPosition().add(new Vector3D(1000, 2000, 3000)),
637 original.getVelocity().add(new Vector3D(-0.03, 0.01, 0.02))),
638 context.initialOrbit.getFrame(),
639 context.initialOrbit.getMu());
640 final Propagator closePropagator = EstimationTestUtils.createPropagator(closeOrbit,
641 propagatorBuilder2);
642 final EphemerisGenerator generator = closePropagator.getEphemerisGenerator();
643 closePropagator.propagate(context.initialOrbit.getDate().shiftedBy(3.5 * closeOrbit.getKeplerianPeriod()));
644 final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
645 Propagator propagator1 = EstimationTestUtils.createPropagator(context.initialOrbit,
646 propagatorBuilder1);
647 final List<ObservedMeasurement<?>> r12 =
648 EstimationTestUtils.createMeasurements(propagator1,
649 new InterSatellitesRangeMeasurementCreator(ephemeris, 0., 0.),
650 1.0, 3.0, 120.0);
651
652
653 propagator1 = EstimationTestUtils.createPropagator(context.initialOrbit,
654 propagatorBuilder1);
655 final List<ObservedMeasurement<?>> r1 =
656 EstimationTestUtils.createMeasurements(propagator1,
657 new TwoWayRangeMeasurementCreator(context),
658 1.0, 3.0, 120.0);
659
660
661 propagator1 = EstimationTestUtils.createPropagator(context.initialOrbit,
662 propagatorBuilder1);
663 final List<ObservedMeasurement<?>> a1 =
664 EstimationTestUtils.createMeasurements(propagator1,
665 new AngularAzElMeasurementCreator(context),
666 1.0, 3.0, 120.0);
667
668
669 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
670 propagatorBuilder1,
671 propagatorBuilder2);
672 for (final ObservedMeasurement<?> interSat : r12) {
673 estimator.addMeasurement(interSat);
674 }
675 for (final ObservedMeasurement<?> range : r1) {
676 estimator.addMeasurement(range);
677 }
678 for (final ObservedMeasurement<?> angular : a1) {
679 estimator.addMeasurement(angular);
680 }
681 estimator.setParametersConvergenceThreshold(1.0e-3);
682 estimator.setMaxIterations(10);
683 estimator.setMaxEvaluations(20);
684 estimator.setObserver(new BatchLSObserver() {
685 int lastIter = 0;
686 int lastEval = 0;
687
688 @Override
689 public void evaluationPerformed(int iterationsCount, int evaluationscount,
690 Orbit[] orbits,
691 ParameterDriversList estimatedOrbitalParameters,
692 ParameterDriversList estimatedPropagatorParameters,
693 ParameterDriversList estimatedMeasurementsParameters,
694 EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) {
695 if (iterationsCount == lastIter) {
696 Assertions.assertEquals(lastEval + 1, evaluationscount);
697 } else {
698 Assertions.assertEquals(lastIter + 1, iterationsCount);
699 }
700 lastIter = iterationsCount;
701 lastEval = evaluationscount;
702
703 AbsoluteDate previous = AbsoluteDate.PAST_INFINITY;
704 for (int i = 0; i < evaluationsProvider.getNumber(); ++i) {
705 AbsoluteDate current = evaluationsProvider.getEstimatedMeasurement(i).getDate();
706 Assertions.assertTrue(current.compareTo(previous) >= 0);
707 previous = current;
708 }
709 }
710 });
711
712 List<DelegatingDriver> parameters = estimator.getOrbitalParametersDrivers(true).getDrivers();
713 ParameterDriver a0Driver = parameters.get(0);
714 Assertions.assertEquals("a[0]", a0Driver.getName());
715 a0Driver.setValue(a0Driver.getValue() + 1.2);
716 a0Driver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
717
718 ParameterDriver a1Driver = parameters.get(6);
719 Assertions.assertEquals("a[1]", a1Driver.getName());
720 a1Driver.setValue(a1Driver.getValue() - 5.4, null);
721 a1Driver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
722
723 final Orbit before = new KeplerianOrbit(parameters.get( 6).getValue(),
724 parameters.get( 7).getValue(),
725 parameters.get( 8).getValue(),
726 parameters.get( 9).getValue(),
727 parameters.get(10).getValue(),
728 parameters.get(11).getValue(),
729 PositionAngleType.TRUE,
730 closeOrbit.getFrame(),
731 closeOrbit.getDate(),
732 closeOrbit.getMu());
733 Assertions.assertEquals(4.7246, Vector3D.distance(closeOrbit.getPosition(),
734 before.getPosition()), 1.0e-3);
735 Assertions.assertEquals(0.0010514, Vector3D.distance(closeOrbit.getPVCoordinates().getVelocity(),
736 before.getPVCoordinates().getVelocity()), 1.0e-6);
737 EstimationTestUtils.checkFit(context, estimator, 4, 6,
738 0.0, 5.3e-06,
739 0.0, 1.4e-05,
740 0.0, 8.8e-07,
741 0.0, 3.6e-10);
742
743 final Orbit determined = new KeplerianOrbit(parameters.get( 6).getValue(),
744 parameters.get( 7).getValue(),
745 parameters.get( 8).getValue(),
746 parameters.get( 9).getValue(),
747 parameters.get(10).getValue(),
748 parameters.get(11).getValue(),
749 PositionAngleType.TRUE,
750 closeOrbit.getFrame(),
751 closeOrbit.getDate(),
752 closeOrbit.getMu());
753 Assertions.assertEquals(0.0, Vector3D.distance(closeOrbit.getPosition(),
754 determined.getPosition()), 5.3e-6);
755 Assertions.assertEquals(0.0, Vector3D.distance(closeOrbit.getPVCoordinates().getVelocity(),
756 determined.getPVCoordinates().getVelocity()), 2.9e-9);
757
758
759
760 for (final ParameterDriver driver : estimator.getOrbitalParametersDrivers(true).getDrivers()) {
761 if (driver.getName().startsWith("a[")) {
762
763 Assertions.assertEquals(0, driver.getReferenceDate().durationFrom(AbsoluteDate.GALILEO_EPOCH), 1.0e-15);
764 } else {
765
766 Assertions.assertEquals(0,
767 driver.getReferenceDate().durationFrom(propagatorBuilder1.getInitialOrbitDate()), 1.0e-15);
768 }
769 }
770
771 }
772
773
774
775
776
777 @Test
778 public void testIssue617() {
779
780 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
781
782 final NumericalPropagatorBuilder propagatorBuilder1 =
783 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
784 1.0e-6, 60.0, 1.0);
785 final NumericalPropagatorBuilder propagatorBuilder2 =
786 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
787 1.0e-6, 60.0, 1.0);
788
789
790 final TimeStampedPVCoordinates original = context.initialOrbit.getPVCoordinates();
791 final Orbit closeOrbit = new CartesianOrbit(new TimeStampedPVCoordinates(context.initialOrbit.getDate(),
792 original.getPosition().add(new Vector3D(1000, 2000, 3000)),
793 original.getVelocity().add(new Vector3D(-0.03, 0.01, 0.02))),
794 context.initialOrbit.getFrame(),
795 context.initialOrbit.getMu());
796 final Propagator closePropagator = EstimationTestUtils.createPropagator(closeOrbit,
797 propagatorBuilder2);
798 final EphemerisGenerator generator = closePropagator.getEphemerisGenerator();
799 closePropagator.propagate(context.initialOrbit.getDate().shiftedBy(3.5 * closeOrbit.getKeplerianPeriod()));
800 final BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
801 Propagator propagator1 = EstimationTestUtils.createPropagator(context.initialOrbit,
802 propagatorBuilder1);
803
804 final List<ObservedMeasurement<?>> r12 =
805 EstimationTestUtils.createMeasurements(propagator1,
806 new InterSatellitesRangeMeasurementCreator(ephemeris, 0., 0.),
807 1.0, 3.0, 300.0);
808
809
810 propagator1 = EstimationTestUtils.createPropagator(context.initialOrbit,
811 propagatorBuilder1);
812 final List<ObservedMeasurement<?>> r1 =
813 EstimationTestUtils.createMeasurements(propagator1,
814 new TwoWayRangeMeasurementCreator(context),
815 1.0, 3.0, 300.0);
816
817
818 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
819 propagatorBuilder1,
820 propagatorBuilder2);
821
822 final List<ObservedMeasurement<?>> independentMeasurements = new ArrayList<>();
823 independentMeasurements.addAll(r1);
824 independentMeasurements.addAll(r12);
825
826
827
828
829
830 final List<ObservedMeasurement<?>> multiplexed = multiplexMeasurements(independentMeasurements, 60.0);
831
832 for (final ObservedMeasurement<?> measurement : multiplexed) {
833 estimator.addMeasurement(measurement);
834 }
835
836 estimator.setParametersConvergenceThreshold(1.0e-2);
837 estimator.setMaxIterations(10);
838 estimator.setMaxEvaluations(20);
839 estimator.setObserver(new BatchLSObserver() {
840 int lastIter = 0;
841 int lastEval = 0;
842
843 @Override
844 public void evaluationPerformed(int iterationsCount, int evaluationscount,
845 Orbit[] orbits,
846 ParameterDriversList estimatedOrbitalParameters,
847 ParameterDriversList estimatedPropagatorParameters,
848 ParameterDriversList estimatedMeasurementsParameters,
849 EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) {
850 if (iterationsCount == lastIter) {
851 Assertions.assertEquals(lastEval + 1, evaluationscount);
852 } else {
853 Assertions.assertEquals(lastIter + 1, iterationsCount);
854 }
855 lastIter = iterationsCount;
856 lastEval = evaluationscount;
857 Assertions.assertEquals(multiplexed.size(), evaluationsProvider.getNumber());
858 try {
859 evaluationsProvider.getEstimatedMeasurement(-1);
860 Assertions.fail("an exception should have been thrown");
861 } catch (OrekitException oe) {
862 Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
863 }
864 try {
865 evaluationsProvider.getEstimatedMeasurement(r12.size() + r1.size());
866 Assertions.fail("an exception should have been thrown");
867 } catch (OrekitException oe) {
868 Assertions.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
869 }
870 AbsoluteDate previous = AbsoluteDate.PAST_INFINITY;
871 for (int i = 0; i < evaluationsProvider.getNumber(); ++i) {
872 AbsoluteDate current = evaluationsProvider.getEstimatedMeasurement(i).getDate();
873 Assertions.assertTrue(current.compareTo(previous) >= 0);
874 previous = current;
875 }
876 }
877 });
878
879 List<DelegatingDriver> parameters = estimator.getOrbitalParametersDrivers(true).getDrivers();
880 ParameterDriver a0Driver = parameters.get(0);
881 Assertions.assertEquals("a[0]", a0Driver.getName());
882 a0Driver.setValue(a0Driver.getValue() + 1.2);
883 a0Driver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
884
885 ParameterDriver a1Driver = parameters.get(6);
886 Assertions.assertEquals("a[1]", a1Driver.getName());
887 a1Driver.setValue(a1Driver.getValue() - 5.4);
888 a1Driver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
889
890 final Orbit before = new KeplerianOrbit(parameters.get( 6).getValue(),
891 parameters.get( 7).getValue(),
892 parameters.get( 8).getValue(),
893 parameters.get( 9).getValue(),
894 parameters.get(10).getValue(),
895 parameters.get(11).getValue(),
896 PositionAngleType.TRUE,
897 closeOrbit.getFrame(),
898 closeOrbit.getDate(),
899 closeOrbit.getMu());
900 Assertions.assertEquals(4.7246, Vector3D.distance(closeOrbit.getPosition(),
901 before.getPosition()), 1.0e-3);
902 Assertions.assertEquals(0.0010514, Vector3D.distance(closeOrbit.getPVCoordinates().getVelocity(),
903 before.getPVCoordinates().getVelocity()), 1.0e-6);
904 EstimationTestUtils.checkFit(context, estimator, 2, 3,
905 0.0, 2.9e-06,
906 0.0, 8.1e-06,
907 0.0, 7.1e-07,
908 0.0, 3.2e-10);
909
910 final Orbit determined = new KeplerianOrbit(parameters.get( 6).getValue(),
911 parameters.get( 7).getValue(),
912 parameters.get( 8).getValue(),
913 parameters.get( 9).getValue(),
914 parameters.get(10).getValue(),
915 parameters.get(11).getValue(),
916 PositionAngleType.TRUE,
917 closeOrbit.getFrame(),
918 closeOrbit.getDate(),
919 closeOrbit.getMu());
920 Assertions.assertEquals(0.0, Vector3D.distance(closeOrbit.getPosition(),
921 determined.getPosition()), 3.3e-6);
922 Assertions.assertEquals(0.0, Vector3D.distance(closeOrbit.getPVCoordinates().getVelocity(),
923 determined.getPVCoordinates().getVelocity()), 1.6e-9);
924
925
926
927 for (final ParameterDriver driver : estimator.getOrbitalParametersDrivers(true).getDrivers()) {
928 if (driver.getName().startsWith("a[")) {
929
930 Assertions.assertEquals(0, driver.getReferenceDate().durationFrom(AbsoluteDate.GALILEO_EPOCH), 1.0e-15);
931 } else {
932
933 Assertions.assertEquals(0,
934 driver.getReferenceDate().durationFrom(propagatorBuilder1.getInitialOrbitDate()), 1.0e-15);
935 }
936 }
937
938 }
939
940 @Test
941 public void testWrappedException() {
942
943 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
944
945 final NumericalPropagatorBuilder propagatorBuilder =
946 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
947 1.0e-6, 60.0, 1.0);
948
949
950 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
951 propagatorBuilder);
952 final List<ObservedMeasurement<?>> measurements =
953 EstimationTestUtils.createMeasurements(propagator,
954 new TwoWayRangeMeasurementCreator(context),
955 1.0, 3.0, 300.0);
956
957
958 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
959 propagatorBuilder);
960 for (final ObservedMeasurement<?> range : measurements) {
961 estimator.addMeasurement(range);
962 }
963 estimator.setParametersConvergenceThreshold(1.0e-2);
964 estimator.setMaxIterations(10);
965 estimator.setMaxEvaluations(20);
966 estimator.setObserver(new BatchLSObserver() {
967
968 @Override
969 public void evaluationPerformed(int iterationsCount, int evaluationscount,
970 Orbit[] orbits,
971 ParameterDriversList estimatedOrbitalParameters,
972 ParameterDriversList estimatedPropagatorParameters,
973 ParameterDriversList estimatedMeasurementsParameters,
974 EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) {
975 throw new DummyException();
976 }
977 });
978
979 try {
980 EstimationTestUtils.checkFit(context, estimator, 3, 4,
981 0.0, 1.5e-6,
982 0.0, 3.2e-6,
983 0.0, 3.8e-7,
984 0.0, 1.5e-10);
985 Assertions.fail("an exception should have been thrown");
986 } catch (DummyException de) {
987
988 }
989
990 }
991
992 private static class DummyException extends OrekitException {
993 private static final long serialVersionUID = 1L;
994 public DummyException() {
995 super(OrekitMessages.INTERNAL_ERROR);
996 }
997 }
998
999
1000
1001
1002 @Test
1003 public void testKeplerRangeRate() {
1004
1005 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
1006
1007 final NumericalPropagatorBuilder propagatorBuilder =
1008 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
1009 1.0e-6, 60.0, 1.0);
1010
1011
1012 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
1013 propagatorBuilder);
1014 final double groundClockDrift = 4.8e-9;
1015 for (final GroundStation station : context.stations) {
1016 station.getClockDriftDriver().setValue(groundClockDrift);
1017 }
1018 final double satClkDrift = 3.2e-10;
1019 final List<ObservedMeasurement<?>> measurements1 =
1020 EstimationTestUtils.createMeasurements(propagator,
1021 new RangeRateMeasurementCreator(context, false, satClkDrift),
1022 1.0, 3.0, 300.0);
1023
1024 final List<ObservedMeasurement<?>> measurements = new ArrayList<ObservedMeasurement<?>>();
1025 measurements.addAll(measurements1);
1026
1027
1028 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
1029 propagatorBuilder);
1030 for (final ObservedMeasurement<?> rangerate : measurements) {
1031 estimator.addMeasurement(rangerate);
1032 }
1033 estimator.setParametersConvergenceThreshold(1.0e-3);
1034 estimator.setMaxIterations(10);
1035 estimator.setMaxEvaluations(20);
1036
1037 EstimationTestUtils.checkFit(context, estimator, 1, 2,
1038 0.0, 5.3e-7,
1039 0.0, 1.3e-6,
1040 0.0, 8.4e-4,
1041 0.0, 5.1e-7);
1042 }
1043
1044
1045
1046
1047 @Test
1048 public void testKeplerRangeAndRangeRate() {
1049
1050 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
1051
1052 final NumericalPropagatorBuilder propagatorBuilder =
1053 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
1054 1.0e-6, 60.0, 1.0);
1055
1056
1057 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
1058 propagatorBuilder);
1059
1060 final List<ObservedMeasurement<?>> measurementsRange =
1061 EstimationTestUtils.createMeasurements(propagator,
1062 new TwoWayRangeMeasurementCreator(context),
1063 1.0, 3.0, 300.0);
1064 final double groundClockDrift = 4.8e-9;
1065 for (final GroundStation station : context.stations) {
1066 station.getClockDriftDriver().setValue(groundClockDrift);
1067 }
1068 final double satClkDrift = 3.2e-10;
1069 final List<ObservedMeasurement<?>> measurementsRangeRate =
1070 EstimationTestUtils.createMeasurements(propagator,
1071 new RangeRateMeasurementCreator(context, false, satClkDrift),
1072 1.0, 3.0, 300.0);
1073
1074
1075 final List<ObservedMeasurement<?>> measurements = new ArrayList<ObservedMeasurement<?>>();
1076 measurements.addAll(measurementsRange);
1077 measurements.addAll(measurementsRangeRate);
1078
1079
1080 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
1081 propagatorBuilder);
1082 for (final ObservedMeasurement<?> meas : measurements) {
1083 estimator.addMeasurement(meas);
1084 }
1085 estimator.setParametersConvergenceThreshold(1.0e-3);
1086 estimator.setMaxIterations(10);
1087 estimator.setMaxEvaluations(20);
1088
1089
1090 EstimationTestUtils.checkFit(context, estimator, 1, 2,
1091 0.0, 4.6e7,
1092 0.0, 1.8e-6,
1093 0.0, 5.8e-7,
1094 0.0, 2.7e-10);
1095 }
1096
1097
1098
1099
1100 @Test
1101 public void testIssue359() {
1102
1103 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
1104
1105 final NumericalPropagatorBuilder propagatorBuilder =
1106 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
1107 1.0e-6, 60.0, 1.0);
1108
1109
1110
1111 propagatorBuilder.getPropagationParametersDrivers().getDrivers().get(0).setSelected(true);
1112
1113 final NumericalPropagator propagator = (NumericalPropagator) EstimationTestUtils.createPropagator(context.initialOrbit,
1114 propagatorBuilder);
1115 final List<ObservedMeasurement<?>> measurements =
1116 EstimationTestUtils.createMeasurements(propagator,
1117 new PVMeasurementCreator(),
1118 0.0, 1.0, 300.0);
1119
1120
1121 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
1122 propagatorBuilder);
1123
1124 for (final ObservedMeasurement<?> measurement : measurements) {
1125 estimator.addMeasurement(measurement);
1126 }
1127
1128 ParameterDriversList estimatedParameters = estimator.getPropagatorParametersDrivers(true);
1129
1130 final String driverName = NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT;
1131 Assertions.assertTrue(propagator.getAllForceModels().get(0) instanceof NewtonianAttraction);
1132 Assertions.assertTrue(propagatorBuilder.getAllForceModels().get(0) instanceof NewtonianAttraction);
1133 Assertions.assertNotNull(estimatedParameters.findByName(driverName));
1134 Assertions.assertTrue(propagator.getAllForceModels().get(0).getParameterDriver(driverName).isSelected());
1135 Assertions.assertTrue(propagatorBuilder.getAllForceModels().get(0).getParameterDriver(driverName).isSelected());
1136 }
1137
1138 @Test
1139 public void testEstimateOnlyOneOrbitalParameter() {
1140 doTestEstimateOnlySomeOrbitalParameters(new boolean[]{ true, false, false, false, false, false });
1141 doTestEstimateOnlySomeOrbitalParameters(new boolean[]{ false, true, false, false, false, false });
1142 doTestEstimateOnlySomeOrbitalParameters(new boolean[]{ false, false, true, false, false, false });
1143 doTestEstimateOnlySomeOrbitalParameters(new boolean[]{ false, false, false, true, false, false });
1144 doTestEstimateOnlySomeOrbitalParameters(new boolean[]{ false, false, false, false, true, false });
1145 doTestEstimateOnlySomeOrbitalParameters(new boolean[]{ false, false, false, false, false, true });
1146 }
1147
1148 @Test
1149 public void testEstimateOnlyFewOrbitalParameters() {
1150 doTestEstimateOnlySomeOrbitalParameters(new boolean[]{ false, true, false, true, false, false });
1151 }
1152
1153 private void doTestEstimateOnlySomeOrbitalParameters(boolean[] orbitalParametersEstimated) {
1154
1155 Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
1156
1157 final NumericalPropagatorBuilder propagatorBuilder =
1158 context.createBuilder(OrbitType.KEPLERIAN, PositionAngleType.TRUE, true,
1159 1.0e-6, 60.0, 1.0e-3);
1160
1161
1162 final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit,
1163 propagatorBuilder);
1164 final List<ObservedMeasurement<?>> measurements =
1165 EstimationTestUtils.createMeasurements(propagator,
1166 new PVMeasurementCreator(),
1167 0.0, 1.0, 300.0);
1168
1169 List<ParameterDriversList.DelegatingDriver> orbitalElementsDrivers = propagatorBuilder.getOrbitalParametersDrivers().getDrivers();
1170 IntStream.range(0, orbitalParametersEstimated.length).forEach(i -> {
1171 final ParameterDriver driver = orbitalElementsDrivers.get(i);
1172 if (orbitalParametersEstimated[i]) {
1173 driver.setSelected(true);
1174 driver.setValue(1.0001 * driver.getReferenceValue());
1175 } else {
1176 driver.setSelected(false);
1177 }
1178 });
1179
1180
1181 final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
1182 propagatorBuilder);
1183 for (final ObservedMeasurement<?> measurement : measurements) {
1184 estimator.addMeasurement(measurement);
1185 }
1186 estimator.setParametersConvergenceThreshold(1.0e-2);
1187 estimator.setMaxIterations(10);
1188 estimator.setMaxEvaluations(20);
1189
1190
1191 estimator.estimate();
1192
1193
1194 IntStream.range(0, orbitalParametersEstimated.length).forEach(i -> {
1195 if (orbitalParametersEstimated[i]) {
1196 final ParameterDriver driver = orbitalElementsDrivers.get(i);
1197 Assertions.assertEquals(driver.getReferenceValue(), driver.getValue(),
1198 driver.getReferenceValue() * 1.0e-3);
1199 }
1200 });
1201
1202 }
1203
1204
1205
1206
1207
1208
1209 private List<ObservedMeasurement<?>> multiplexMeasurements(final List<ObservedMeasurement<?>> independentMeasurements,
1210 final double tol) {
1211 final List<ObservedMeasurement<?>> multiplexed = new ArrayList<>();
1212 independentMeasurements.sort(new ChronologicalComparator());
1213 List<ObservedMeasurement<?>> clump = new ArrayList<>();
1214 for (final ObservedMeasurement<?> measurement : independentMeasurements) {
1215 if (!clump.isEmpty() && measurement.getDate().durationFrom(clump.get(0).getDate()) > tol) {
1216
1217
1218 if (clump.size() == 1) {
1219 multiplexed.add(clump.get(0));
1220 } else {
1221 multiplexed.add(new MultiplexedMeasurement(clump));
1222 }
1223
1224
1225 clump = new ArrayList<>();
1226
1227 }
1228 clump.add(measurement);
1229 }
1230
1231 if (clump.size() == 1) {
1232 multiplexed.add(clump.get(0));
1233 } else {
1234 multiplexed.add(new MultiplexedMeasurement(clump));
1235 }
1236 return multiplexed;
1237 }
1238
1239 }
1240
1241