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