Added more unit tests for analytical orbit determination.

src/test/java/org/orekit/estimation/EcksteinHechlerContext.java

+/* Copyright 2002-2022 CS GROUP
+ * Licensed to CS GROUP (CS) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * CS licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.orekit.estimation;
+
+import java.util.List;
+import java.util.Map;
+
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.hipparchus.util.FastMath;
+import org.orekit.bodies.CelestialBody;
+import org.orekit.bodies.GeodeticPoint;
+import org.orekit.bodies.OneAxisEllipsoid;
+import org.orekit.estimation.measurements.GroundStation;
+import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
+import org.orekit.frames.TopocentricFrame;
+import org.orekit.models.earth.displacement.StationDisplacement;
+import org.orekit.orbits.CartesianOrbit;
+import org.orekit.orbits.Orbit;
+import org.orekit.orbits.PositionAngle;
+import org.orekit.propagation.conversion.EcksteinHechlerPropagatorBuilder;
+import org.orekit.time.TimeScale;
+import org.orekit.time.UT1Scale;
+import org.orekit.utils.IERSConventions;
+import org.orekit.utils.PVCoordinates;
+
+public class EcksteinHechlerContext implements StationDataProvider {
+
+    public IERSConventions                        conventions;
+    public OneAxisEllipsoid                       earth;
+    public CelestialBody                          sun;
+    public CelestialBody                          moon;
+    public UnnormalizedSphericalHarmonicsProvider gravity;
+    public TimeScale                              utc;
+    public UT1Scale                               ut1;
+    public Orbit                                  initialOrbit;
+    public StationDisplacement[]                  displacements;
+    public List<GroundStation>                    stations;
+    // Stations for turn-around range
+    // Map entry = primary station
+    // Map value = secondary station associated
+    public Map<GroundStation, GroundStation>     TARstations;
+
+    public EcksteinHechlerPropagatorBuilder createBuilder(final PositionAngle angleType, final boolean perfectStart, final double dP) {
+
+        final Orbit startOrbit;
+        if (perfectStart) {
+            // orbit estimation will start from a perfect orbit
+            startOrbit = initialOrbit;
+        } else {
+            // orbit estimation will start from a wrong point
+            final Vector3D initialPosition = initialOrbit.getPVCoordinates().getPosition();
+            final Vector3D initialVelocity = initialOrbit.getPVCoordinates().getVelocity();
+            final Vector3D wrongPosition   = initialPosition.add(new Vector3D(1000.0, 0, 0));
+            final Vector3D wrongVelocity   = initialVelocity.add(new Vector3D(0, 0, 0.01));
+            startOrbit                     = new CartesianOrbit(new PVCoordinates(wrongPosition, wrongVelocity),
+                                                                initialOrbit.getFrame(),
+                                                                initialOrbit.getDate(),
+                                                                initialOrbit.getMu());
+        }
+
+        // Initialize builder
+        final EcksteinHechlerPropagatorBuilder propagatorBuilder =
+                        new EcksteinHechlerPropagatorBuilder(startOrbit, gravity, angleType, dP);
+
+        // Return
+        return propagatorBuilder;
+
+    }
+
+    GroundStation createStation(double latitudeInDegrees, double longitudeInDegrees,
+                                double altitude, String name) {
+        final GeodeticPoint gp = new GeodeticPoint(FastMath.toRadians(latitudeInDegrees),
+                                                   FastMath.toRadians(longitudeInDegrees),
+                                                   altitude);
+        return new GroundStation(new TopocentricFrame(earth, gp, name),
+                                 ut1.getEOPHistory(), displacements);
+    }
+
+    @Override
+    public List<GroundStation> getStations() {
+        return stations;
+    }
+
+}

src/test/java/org/orekit/estimation/EcksteinHechlerEstimationTestUtils.java

+/* Copyright 2002-2022 CS GROUP
+ * Licensed to CS GROUP (CS) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * CS licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.orekit.estimation;
+
+import java.util.Arrays;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.hipparchus.linear.MatrixUtils;
+import org.hipparchus.linear.RealMatrix;
+import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresOptimizer.Optimum;
+import org.hipparchus.util.FastMath;
+import org.junit.Assert;
+import org.orekit.Utils;
+import org.orekit.bodies.CelestialBodyFactory;
+import org.orekit.bodies.OneAxisEllipsoid;
+import org.orekit.data.DataContext;
+import org.orekit.estimation.leastsquares.BatchLSEstimator;
+import org.orekit.estimation.measurements.EstimatedMeasurement;
+import org.orekit.estimation.measurements.GroundStation;
+import org.orekit.estimation.measurements.MeasurementCreator;
+import org.orekit.estimation.measurements.ObservedMeasurement;
+import org.orekit.estimation.sequential.KalmanEstimator;
+import org.orekit.forces.gravity.potential.AstronomicalAmplitudeReader;
+import org.orekit.forces.gravity.potential.FESCHatEpsilonReader;
+import org.orekit.forces.gravity.potential.GRGSFormatReader;
+import org.orekit.forces.gravity.potential.GravityFieldFactory;
+import org.orekit.forces.gravity.potential.OceanLoadDeformationCoefficients;
+import org.orekit.frames.EOPHistory;
+import org.orekit.frames.FramesFactory;
+import org.orekit.models.earth.displacement.StationDisplacement;
+import org.orekit.models.earth.displacement.TidalDisplacement;
+import org.orekit.orbits.KeplerianOrbit;
+import org.orekit.orbits.Orbit;
+import org.orekit.orbits.PositionAngle;
+import org.orekit.propagation.Propagator;
+import org.orekit.propagation.conversion.PropagatorBuilder;
+import org.orekit.time.AbsoluteDate;
+import org.orekit.time.TimeScalesFactory;
+import org.orekit.utils.Constants;
+import org.orekit.utils.IERSConventions;
+import org.orekit.utils.ParameterDriver;
+
+public class EcksteinHechlerEstimationTestUtils {
+
+    public static EcksteinHechlerContext smallEccentricContext(final String dataRoot) {
+
+        Utils.setDataRoot(dataRoot);
+        EcksteinHechlerContext context = new EcksteinHechlerContext();
+        context.conventions = IERSConventions.IERS_2010;
+        context.earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
+                                             Constants.WGS84_EARTH_FLATTENING,
+                                             FramesFactory.getITRF(context.conventions, true));
+        context.sun = CelestialBodyFactory.getSun();
+        context.moon = CelestialBodyFactory.getMoon();
+        final EOPHistory eopHistory = FramesFactory.getEOPHistory(context.conventions, true);
+        context.utc = TimeScalesFactory.getUTC();
+        context.ut1 = TimeScalesFactory.getUT1(eopHistory);
+        context.displacements = new StationDisplacement[] {
+            new TidalDisplacement(Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
+                                  Constants.JPL_SSD_SUN_EARTH_PLUS_MOON_MASS_RATIO,
+                                  Constants.JPL_SSD_EARTH_MOON_MASS_RATIO,
+                                  context.sun, context.moon,
+                                  context.conventions, false)
+        };
+        GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
+        AstronomicalAmplitudeReader aaReader =
+                        new AstronomicalAmplitudeReader("hf-fes2004.dat", 5, 2, 3, 1.0);
+        DataContext.getDefault().getDataProvidersManager().feed(aaReader.getSupportedNames(), aaReader);
+        Map<Integer, Double> map = aaReader.getAstronomicalAmplitudesMap();
+        GravityFieldFactory.addOceanTidesReader(new FESCHatEpsilonReader("fes2004-7x7.dat",
+                                                                         0.01, FastMath.toRadians(1.0),
+                                                                         OceanLoadDeformationCoefficients.IERS_2010,
+                                                                         map));
+        context.gravity = GravityFieldFactory.getUnnormalizedProvider(6, 0);
+        context.initialOrbit = new KeplerianOrbit(15000000.0, 0.001, 1.25,
+                                                  0.250, 1.375, 0.0625, PositionAngle.TRUE,
+                                                  FramesFactory.getEME2000(),
+                                                  new AbsoluteDate(2000, 2, 24, 11, 35, 47.0, context.utc),
+                                                  context.gravity.getMu());
+
+        context.stations = Arrays.asList(//context.createStation(-18.59146, -173.98363,   76.0, "Leimatu`a"),
+                                         context.createStation(-53.05388,  -75.01551, 1750.0, "Isla Desolación"),
+                                         context.createStation( 62.29639,   -7.01250,  880.0, "Slættaratindur")
+                                         //context.createStation( -4.01583,  103.12833, 3173.0, "Gunung Dempo")
+                        );
+
+        // Turn-around range stations
+        // Map entry = primary station
+        // Map value = secondary station associated
+        context.TARstations = new HashMap<GroundStation, GroundStation>();
+
+        context.TARstations.put(context.createStation(-53.05388,  -75.01551, 1750.0, "Isla Desolación"),
+                                context.createStation(-54.815833,  -68.317778, 6.0, "Ushuaïa"));
+
+        context.TARstations.put(context.createStation( 62.29639,   -7.01250,  880.0, "Slættaratindur"),
+                                context.createStation( 61.405833,   -6.705278,  470.0, "Sumba"));
+
+        return context;
+
+    }
+
+    public static Propagator createPropagator(final Orbit initialOrbit,
+                                              final PropagatorBuilder propagatorBuilder) {
+
+        // override orbital parameters
+        double[] orbitArray = new double[6];
+        propagatorBuilder.getOrbitType().mapOrbitToArray(initialOrbit,
+                                                         propagatorBuilder.getPositionAngle(),
+                                                         orbitArray, null);
+        for (int i = 0; i < orbitArray.length; ++i) {
+            propagatorBuilder.getOrbitalParametersDrivers().getDrivers().get(i).setValue(orbitArray[i]);
+        }
+
+        return propagatorBuilder.buildPropagator(propagatorBuilder.getSelectedNormalizedParameters());
+
+    }
+
+    public static List<ObservedMeasurement<?>> createMeasurements(final Propagator propagator,
+                                                                  final MeasurementCreator creator,
+                                                                  final double startPeriod, final double endPeriod,
+                                                                  final double step) {
+
+        propagator.setStepHandler(step, creator);
+        final double       period = propagator.getInitialState().getKeplerianPeriod();
+        final AbsoluteDate start  = propagator.getInitialState().getDate().shiftedBy(startPeriod * period);
+        final AbsoluteDate end    = propagator.getInitialState().getDate().shiftedBy(endPeriod   * period);
+        propagator.propagate(start, end);
+
+        final List<ObservedMeasurement<?>> measurements = creator.getMeasurements();
+
+        for (final ObservedMeasurement<?> measurement : measurements) {
+            for (final ParameterDriver driver : measurement.getParametersDrivers()) {
+                if (driver.getReferenceDate() == null) {
+                    driver.setReferenceDate(propagator.getInitialState().getDate());
+                }
+            }
+        }
+
+        return measurements;
+
+    }
+
+    /**
+     * Checker for batch LS estimator validation
+     * @param context Context used for the test
+     * @param estimator Batch LS estimator
+     * @param iterations Number of iterations expected
+     * @param evaluations Number of evaluations expected
+     * @param expectedRMS Expected RMS value
+     * @param rmsEps Tolerance on expected RMS
+     * @param expectedMax Expected weighted residual maximum
+     * @param maxEps Tolerance on weighted residual maximum
+     * @param expectedDeltaPos Expected position difference between estimated orbit and initial orbit
+     * @param posEps Tolerance on expected position difference
+     * @param expectedDeltaVel Expected velocity difference between estimated orbit and initial orbit
+     * @param velEps Tolerance on expected velocity difference
+     */
+    public static void checkFit(final EcksteinHechlerContext context, final BatchLSEstimator estimator,
+                                final int iterations, final int evaluations,
+                                final double expectedRMS,      final double rmsEps,
+                                final double expectedMax,      final double maxEps,
+                                final double expectedDeltaPos, final double posEps,
+                                final double expectedDeltaVel, final double velEps) {
+
+        final Orbit estimatedOrbit = estimator.estimate()[0].getInitialState().getOrbit();
+        final Vector3D estimatedPosition = estimatedOrbit.getPVCoordinates().getPosition();
+        final Vector3D estimatedVelocity = estimatedOrbit.getPVCoordinates().getVelocity();
+
+        Assert.assertEquals(iterations, estimator.getIterationsCount());
+        Assert.assertEquals(evaluations, estimator.getEvaluationsCount());
+        Optimum optimum = estimator.getOptimum();
+        Assert.assertEquals(iterations, optimum.getIterations());
+        Assert.assertEquals(evaluations, optimum.getEvaluations());
+
+        int    k   = 0;
+        double sum = 0;
+        double max = 0;
+        for (final Map.Entry<ObservedMeasurement<?>, EstimatedMeasurement<?>> entry :
+             estimator.getLastEstimations().entrySet()) {
+            final ObservedMeasurement<?>  m = entry.getKey();
+            final EstimatedMeasurement<?> e = entry.getValue();
+            final double[]    weight      = m.getBaseWeight();
+            final double[]    sigma       = m.getTheoreticalStandardDeviation();
+            final double[]    observed    = m.getObservedValue();
+            final double[]    theoretical = e.getEstimatedValue();
+            for (int i = 0; i < m.getDimension(); ++i) {
+                final double weightedResidual = weight[i] * (theoretical[i] - observed[i]) / sigma[i];
+                ++k;
+                sum += weightedResidual * weightedResidual;
+                max = FastMath.max(max, FastMath.abs(weightedResidual));
+            }
+        }
+
+        final double rms = FastMath.sqrt(sum / k);
+        final double deltaPos = Vector3D.distance(context.initialOrbit.getPVCoordinates().getPosition(), estimatedPosition);
+        final double deltaVel = Vector3D.distance(context.initialOrbit.getPVCoordinates().getVelocity(), estimatedVelocity);
+        Assert.assertEquals(expectedRMS,
+                            rms,
+                            rmsEps);
+        Assert.assertEquals(expectedMax,
+                            max,
+                            maxEps);
+        Assert.assertEquals(expectedDeltaPos,
+                            deltaPos,
+                            posEps);
+        Assert.assertEquals(expectedDeltaVel,
+                            deltaVel,
+                            velEps);
+
+    }
+    
+    /**
+     * Checker for Kalman estimator validation
+     * @param context context used for the test
+     * @param kalman Kalman filter
+     * @param measurements List of observed measurements to be processed by the Kalman
+     * @param refOrbit Reference orbits at last measurement date
+     * @param expectedDeltaPos Expected position difference between estimated orbit and reference orbit
+     * @param posEps Tolerance on expected position difference
+     * @param expectedDeltaVel Expected velocity difference between estimated orbit and reference orbit
+     * @param velEps Tolerance on expected velocity difference
+     * @param expectedSigmasPos Expected values for covariance matrix on position
+     * @param sigmaPosEps Tolerance on expected covariance matrix on position
+     * @param expectedSigmasVel Expected values for covariance matrix on velocity
+     * @param sigmaVelEps Tolerance on expected covariance matrix on velocity
+     */
+    public static void checkKalmanFit(final EcksteinHechlerContext context, final KalmanEstimator kalman,
+                                      final List<ObservedMeasurement<?>> measurements,
+                                      final Orbit refOrbit, final PositionAngle positionAngle,
+                                      final double expectedDeltaPos, final double posEps,
+                                      final double expectedDeltaVel, final double velEps,
+                                      final double[] expectedSigmasPos,final double sigmaPosEps,
+                                      final double[] expectedSigmasVel,final double sigmaVelEps)
+        {
+        checkKalmanFit(context, kalman, measurements,
+                       new Orbit[] { refOrbit },
+                       new PositionAngle[] { positionAngle },
+                       new double[] { expectedDeltaPos }, new double[] { posEps },
+                       new double[] { expectedDeltaVel }, new double[] { velEps },
+                       new double[][] { expectedSigmasPos }, new double[] { sigmaPosEps },
+                       new double[][] { expectedSigmasVel }, new double[] { sigmaVelEps });
+    }
+
+    /**
+     * Checker for Kalman estimator validation
+     * @param context context used for the test
+     * @param kalman Kalman filter
+     * @param measurements List of observed measurements to be processed by the Kalman
+     * @param refOrbit Reference orbits at last measurement date
+     * @param expectedDeltaPos Expected position difference between estimated orbit and reference orbits
+     * @param posEps Tolerance on expected position difference
+     * @param expectedDeltaVel Expected velocity difference between estimated orbit and reference orbits
+     * @param velEps Tolerance on expected velocity difference
+     * @param expectedSigmasPos Expected values for covariance matrix on position
+     * @param sigmaPosEps Tolerance on expected covariance matrix on position
+     * @param expectedSigmasVel Expected values for covariance matrix on velocity
+     * @param sigmaVelEps Tolerance on expected covariance matrix on velocity
+     */
+    public static void checkKalmanFit(final EcksteinHechlerContext context, final KalmanEstimator kalman,
+                                      final List<ObservedMeasurement<?>> measurements,
+                                      final Orbit[] refOrbit, final PositionAngle[] positionAngle,
+                                      final double[] expectedDeltaPos, final double[] posEps,
+                                      final double[] expectedDeltaVel, final double []velEps,
+                                      final double[][] expectedSigmasPos,final double[] sigmaPosEps,
+                                      final double[][] expectedSigmasVel,final double[] sigmaVelEps)
+                                                      {
+
+        // Add the measurements to the Kalman filter
+        Propagator[] estimated = kalman.processMeasurements(measurements);
+        
+        // Check the number of measurements processed by the filter
+        Assert.assertEquals(measurements.size(), kalman.getCurrentMeasurementNumber());
+
+        for (int k = 0; k < refOrbit.length; ++k) {
+            // Get the last estimation
+            final Orbit    estimatedOrbit    = estimated[k].getInitialState().getOrbit();
+            final Vector3D estimatedPosition = estimatedOrbit.getPVCoordinates().getPosition();
+            final Vector3D estimatedVelocity = estimatedOrbit.getPVCoordinates().getVelocity();        
+
+            // Get the last covariance matrix estimation
+            final RealMatrix estimatedP = kalman.getPhysicalEstimatedCovarianceMatrix();
+
+            // Convert the orbital part to Cartesian formalism
+            // Assuming all 6 orbital parameters are estimated by the filter
+            final double[][] dCdY = new double[6][6];
+            estimatedOrbit.getJacobianWrtParameters(positionAngle[k], dCdY);
+            final RealMatrix Jacobian = MatrixUtils.createRealMatrix(dCdY);
+            final RealMatrix estimatedCartesianP = 
+                            Jacobian.
+                            multiply(estimatedP.getSubMatrix(0, 5, 0, 5)).
+                            multiply(Jacobian.transpose());
+
+            // Get the final sigmas (ie.sqrt of the diagonal of the Cartesian orbital covariance matrix)
+            final double[] sigmas = new double[6];
+            for (int i = 0; i < 6; i++) {
+                sigmas[i] = FastMath.sqrt(estimatedCartesianP.getEntry(i, i));
+            }
+
+            // Check the final orbit estimation & PV sigmas
+            final double deltaPosK = Vector3D.distance(refOrbit[k].getPVCoordinates().getPosition(), estimatedPosition);
+            final double deltaVelK = Vector3D.distance(refOrbit[k].getPVCoordinates().getVelocity(), estimatedVelocity);
+            Assert.assertEquals(expectedDeltaPos[k], deltaPosK, posEps[k]);
+            Assert.assertEquals(expectedDeltaVel[k], deltaVelK, velEps[k]);
+
+            for (int i = 0; i < 3; i++) {
+                System.out.println(sigmas[i]);
+                System.out.println(sigmas[i+3]);
+                Assert.assertEquals(expectedSigmasPos[k][i], sigmas[i],   sigmaPosEps[k]);
+                Assert.assertEquals(expectedSigmasVel[k][i], sigmas[i+3], sigmaVelEps[k]);
+            }
+        }
+    }
+
+}

src/test/java/org/orekit/estimation/KeplerianContext.java

+package org.orekit.estimation;
+
+import java.util.List;
+
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.hipparchus.util.FastMath;
+import org.orekit.bodies.CelestialBody;
+import org.orekit.bodies.GeodeticPoint;
+import org.orekit.bodies.OneAxisEllipsoid;
+import org.orekit.estimation.measurements.GroundStation;
+import org.orekit.frames.TopocentricFrame;
+import org.orekit.models.earth.displacement.StationDisplacement;
+import org.orekit.orbits.CartesianOrbit;
+import org.orekit.orbits.Orbit;
+import org.orekit.orbits.PositionAngle;
+import org.orekit.propagation.conversion.KeplerianPropagatorBuilder;
+import org.orekit.time.TimeScale;
+import org.orekit.time.UT1Scale;
+import org.orekit.utils.IERSConventions;
+import org.orekit.utils.PVCoordinates;
+
+public class KeplerianContext implements StationDataProvider {
+
+    public IERSConventions                        conventions;
+    public OneAxisEllipsoid                       earth;
+    public CelestialBody                          sun;
+    public CelestialBody                          moon;
+    public TimeScale                              utc;
+    public UT1Scale                               ut1;
+    public Orbit                                  initialOrbit;
+    public StationDisplacement[]                  displacements;
+    public List<GroundStation>                    stations;
+
+    public KeplerianPropagatorBuilder createBuilder(final PositionAngle angleType, final boolean perfectStart, final double dP) {
+
+        final Orbit startOrbit;
+        if (perfectStart) {
+            // orbit estimation will start from a perfect orbit
+            startOrbit = initialOrbit;
+        } else {
+            // orbit estimation will start from a wrong point
+            final Vector3D initialPosition = initialOrbit.getPVCoordinates().getPosition();
+            final Vector3D initialVelocity = initialOrbit.getPVCoordinates().getVelocity();
+            final Vector3D wrongPosition   = initialPosition.add(new Vector3D(1000.0, 0, 0));
+            final Vector3D wrongVelocity   = initialVelocity.add(new Vector3D(0, 0, 0.01));
+            startOrbit                     = new CartesianOrbit(new PVCoordinates(wrongPosition, wrongVelocity),
+                                                                initialOrbit.getFrame(),
+                                                                initialOrbit.getDate(),
+                                                                initialOrbit.getMu());
+        }
+
+        // Initialize builder
+        final KeplerianPropagatorBuilder propagatorBuilder =
+                        new KeplerianPropagatorBuilder(startOrbit, angleType, dP);
+
+        // Return
+        return propagatorBuilder;
+
+    }
+
+    GroundStation createStation(double latitudeInDegrees, double longitudeInDegrees,
+                                double altitude, String name) {
+        final GeodeticPoint gp = new GeodeticPoint(FastMath.toRadians(latitudeInDegrees),
+                                                   FastMath.toRadians(longitudeInDegrees),
+                                                   altitude);
+        return new GroundStation(new TopocentricFrame(earth, gp, name),
+                                 ut1.getEOPHistory(), displacements);
+    }
+
+    @Override
+    public List<GroundStation> getStations() {
+        return stations;
+    }
+
+}

src/test/java/org/orekit/estimation/KeplerianEstimationTestUtils.java

+/* Copyright 2002-2022 CS GROUP
+ * Licensed to CS GROUP (CS) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * CS licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.orekit.estimation;
+
+import java.util.Arrays;
+import java.util.List;
+import java.util.Map;
+
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.hipparchus.linear.MatrixUtils;
+import org.hipparchus.linear.RealMatrix;
+import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresOptimizer.Optimum;
+import org.hipparchus.util.FastMath;
+import org.junit.Assert;
+import org.orekit.Utils;
+import org.orekit.bodies.CelestialBodyFactory;
+import org.orekit.bodies.OneAxisEllipsoid;
+import org.orekit.data.DataContext;
+import org.orekit.estimation.leastsquares.BatchLSEstimator;
+import org.orekit.estimation.measurements.EstimatedMeasurement;
+import org.orekit.estimation.measurements.MeasurementCreator;
+import org.orekit.estimation.measurements.ObservedMeasurement;
+import org.orekit.estimation.sequential.KalmanEstimator;
+import org.orekit.forces.gravity.potential.AstronomicalAmplitudeReader;
+import org.orekit.forces.gravity.potential.FESCHatEpsilonReader;
+import org.orekit.forces.gravity.potential.GRGSFormatReader;
+import org.orekit.forces.gravity.potential.GravityFieldFactory;
+import org.orekit.forces.gravity.potential.OceanLoadDeformationCoefficients;
+import org.orekit.frames.EOPHistory;
+import org.orekit.frames.FramesFactory;
+import org.orekit.models.earth.displacement.StationDisplacement;
+import org.orekit.models.earth.displacement.TidalDisplacement;
+import org.orekit.orbits.KeplerianOrbit;
+import org.orekit.orbits.Orbit;
+import org.orekit.orbits.PositionAngle;
+import org.orekit.propagation.Propagator;
+import org.orekit.propagation.conversion.PropagatorBuilder;
+import org.orekit.time.AbsoluteDate;
+import org.orekit.time.TimeScalesFactory;
+import org.orekit.utils.Constants;
+import org.orekit.utils.IERSConventions;
+import org.orekit.utils.ParameterDriver;
+
+/** Utility class for orbit determination tests. */
+public class KeplerianEstimationTestUtils {
+
+    public static KeplerianContext eccentricContext(final String dataRoot) {
+
+        Utils.setDataRoot(dataRoot);
+        KeplerianContext context = new KeplerianContext();
+        context.conventions = IERSConventions.IERS_2010;
+        context.earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
+                                             Constants.WGS84_EARTH_FLATTENING,
+                                             FramesFactory.getITRF(context.conventions, true));
+        context.sun = CelestialBodyFactory.getSun();
+        context.moon = CelestialBodyFactory.getMoon();
+        final EOPHistory eopHistory = FramesFactory.getEOPHistory(context.conventions, true);
+        context.utc = TimeScalesFactory.getUTC();
+        context.ut1 = TimeScalesFactory.getUT1(eopHistory);
+        context.displacements = new StationDisplacement[] {
+            new TidalDisplacement(Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS,
+                                  Constants.JPL_SSD_SUN_EARTH_PLUS_MOON_MASS_RATIO,
+                                  Constants.JPL_SSD_EARTH_MOON_MASS_RATIO,
+                                  context.sun, context.moon,
+                                  context.conventions, false)
+        };
+        GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
+        AstronomicalAmplitudeReader aaReader =
+                        new AstronomicalAmplitudeReader("hf-fes2004.dat", 5, 2, 3, 1.0);
+        DataContext.getDefault().getDataProvidersManager().feed(aaReader.getSupportedNames(), aaReader);
+        Map<Integer, Double> map = aaReader.getAstronomicalAmplitudesMap();
+        GravityFieldFactory.addOceanTidesReader(new FESCHatEpsilonReader("fes2004-7x7.dat",
+                                                                         0.01, FastMath.toRadians(1.0),
+                                                                         OceanLoadDeformationCoefficients.IERS_2010,
+                                                                         map));
+        context.initialOrbit = new KeplerianOrbit(15000000.0, 0.125, 1.25,
+                                                  0.250, 1.375, 0.0625, PositionAngle.TRUE,
+                                                  FramesFactory.getEME2000(),
+                                                  new AbsoluteDate(2000, 2, 24, 11, 35, 47.0, context.utc),
+                                                  Constants.WGS84_EARTH_MU);
+
+        context.stations = Arrays.asList(//context.createStation(-18.59146, -173.98363,   76.0, "Leimatu`a"),
+                                         context.createStation(-53.05388,  -75.01551, 1750.0, "Isla Desolación"),
+                                         context.createStation( 62.29639,   -7.01250,  880.0, "Slættaratindur")
+                                         //context.createStation( -4.01583,  103.12833, 3173.0, "Gunung Dempo")
+                        );
+
+        return context;
+
+    }
+
+    public static Propagator createPropagator(final Orbit initialOrbit,
+                                              final PropagatorBuilder propagatorBuilder) {
+
+        // override orbital parameters
+        double[] orbitArray = new double[6];
+        propagatorBuilder.getOrbitType().mapOrbitToArray(initialOrbit,
+                                                         propagatorBuilder.getPositionAngle(),
+                                                         orbitArray, null);
+        for (int i = 0; i < orbitArray.length; ++i) {
+            propagatorBuilder.getOrbitalParametersDrivers().getDrivers().get(i).setValue(orbitArray[i]);
+        }
+
+        return propagatorBuilder.buildPropagator(propagatorBuilder.getSelectedNormalizedParameters());
+
+    }
+
+    public static List<ObservedMeasurement<?>> createMeasurements(final Propagator propagator,
+                                                                  final MeasurementCreator creator,
+                                                                  final double startPeriod, final double endPeriod,
+                                                                  final double step) {
+
+        propagator.setStepHandler(step, creator);
+        final double       period = propagator.getInitialState().getKeplerianPeriod();
+        final AbsoluteDate start  = propagator.getInitialState().getDate().shiftedBy(startPeriod * period);
+        final AbsoluteDate end    = propagator.getInitialState().getDate().shiftedBy(endPeriod   * period);
+        propagator.propagate(start, end);
+
+        final List<ObservedMeasurement<?>> measurements = creator.getMeasurements();
+
+        for (final ObservedMeasurement<?> measurement : measurements) {
+            for (final ParameterDriver driver : measurement.getParametersDrivers()) {
+                if (driver.getReferenceDate() == null) {
+                    driver.setReferenceDate(propagator.getInitialState().getDate());
+                }
+            }
+        }
+
+        return measurements;
+
+    }
+
+    /**
+     * Checker for batch LS estimator validation
+     * @param context Context used for the test
+     * @param estimator Batch LS estimator
+     * @param iterations Number of iterations expected
+     * @param evaluations Number of evaluations expected
+     * @param expectedRMS Expected RMS value
+     * @param rmsEps Tolerance on expected RMS
+     * @param expectedMax Expected weighted residual maximum
+     * @param maxEps Tolerance on weighted residual maximum
+     * @param expectedDeltaPos Expected position difference between estimated orbit and initial orbit
+     * @param posEps Tolerance on expected position difference
+     * @param expectedDeltaVel Expected velocity difference between estimated orbit and initial orbit
+     * @param velEps Tolerance on expected velocity difference
+     */
+    public static void checkFit(final KeplerianContext context, final BatchLSEstimator estimator,
+                                final int iterations, final int evaluations,
+                                final double expectedRMS,      final double rmsEps,
+                                final double expectedMax,      final double maxEps,
+                                final double expectedDeltaPos, final double posEps,
+                                final double expectedDeltaVel, final double velEps) {
+
+        final Orbit estimatedOrbit = estimator.estimate()[0].getInitialState().getOrbit();
+        final Vector3D estimatedPosition = estimatedOrbit.getPVCoordinates().getPosition();
+        final Vector3D estimatedVelocity = estimatedOrbit.getPVCoordinates().getVelocity();
+
+        Assert.assertEquals(iterations, estimator.getIterationsCount());
+        Assert.assertEquals(evaluations, estimator.getEvaluationsCount());
+        Optimum optimum = estimator.getOptimum();
+        Assert.assertEquals(iterations, optimum.getIterations());
+        Assert.assertEquals(evaluations, optimum.getEvaluations());
+
+        int    k   = 0;
+        double sum = 0;
+        double max = 0;
+        for (final Map.Entry<ObservedMeasurement<?>, EstimatedMeasurement<?>> entry :
+             estimator.getLastEstimations().entrySet()) {
+            final ObservedMeasurement<?>  m = entry.getKey();
+            final EstimatedMeasurement<?> e = entry.getValue();
+            final double[]    weight      = m.getBaseWeight();
+            final double[]    sigma       = m.getTheoreticalStandardDeviation();
+            final double[]    observed    = m.getObservedValue();
+            final double[]    theoretical = e.getEstimatedValue();
+            for (int i = 0; i < m.getDimension(); ++i) {
+                final double weightedResidual = weight[i] * (theoretical[i] - observed[i]) / sigma[i];
+                ++k;
+                sum += weightedResidual * weightedResidual;
+                max = FastMath.max(max, FastMath.abs(weightedResidual));
+            }
+        }
+
+        final double rms = FastMath.sqrt(sum / k);
+        final double deltaPos = Vector3D.distance(context.initialOrbit.getPVCoordinates().getPosition(), estimatedPosition);
+        final double deltaVel = Vector3D.distance(context.initialOrbit.getPVCoordinates().getVelocity(), estimatedVelocity);
+        Assert.assertEquals(expectedRMS,
+                            rms,
+                            rmsEps);
+        Assert.assertEquals(expectedMax,
+                            max,
+                            maxEps);
+        Assert.assertEquals(expectedDeltaPos,
+                            deltaPos,
+                            posEps);
+        Assert.assertEquals(expectedDeltaVel,
+                            deltaVel,
+                            velEps);
+
+    }
+    
+    /**
+     * Checker for Kalman estimator validation
+     * @param context context used for the test
+     * @param kalman Kalman filter
+     * @param measurements List of observed measurements to be processed by the Kalman
+     * @param refOrbit Reference orbits at last measurement date
+     * @param expectedDeltaPos Expected position difference between estimated orbit and reference orbit
+     * @param posEps Tolerance on expected position difference
+     * @param expectedDeltaVel Expected velocity difference between estimated orbit and reference orbit
+     * @param velEps Tolerance on expected velocity difference
+     * @param expectedSigmasPos Expected values for covariance matrix on position
+     * @param sigmaPosEps Tolerance on expected covariance matrix on position
+     * @param expectedSigmasVel Expected values for covariance matrix on velocity
+     * @param sigmaVelEps Tolerance on expected covariance matrix on velocity
+     */
+    public static void checkKalmanFit(final KeplerianContext context, final KalmanEstimator kalman,
+                                      final List<ObservedMeasurement<?>> measurements,
+                                      final Orbit refOrbit, final PositionAngle positionAngle,
+                                      final double expectedDeltaPos, final double posEps,
+                                      final double expectedDeltaVel, final double velEps,
+                                      final double[] expectedSigmasPos,final double sigmaPosEps,
+                                      final double[] expectedSigmasVel,final double sigmaVelEps)
+        {
+        checkKalmanFit(context, kalman, measurements,
+                       new Orbit[] { refOrbit },
+                       new PositionAngle[] { positionAngle },
+                       new double[] { expectedDeltaPos }, new double[] { posEps },
+                       new double[] { expectedDeltaVel }, new double[] { velEps },
+                       new double[][] { expectedSigmasPos }, new double[] { sigmaPosEps },
+                       new double[][] { expectedSigmasVel }, new double[] { sigmaVelEps });
+    }
+
+    /**
+     * Checker for Kalman estimator validation
+     * @param context context used for the test
+     * @param kalman Kalman filter
+     * @param measurements List of observed measurements to be processed by the Kalman
+     * @param refOrbit Reference orbits at last measurement date
+     * @param expectedDeltaPos Expected position difference between estimated orbit and reference orbits
+     * @param posEps Tolerance on expected position difference
+     * @param expectedDeltaVel Expected velocity difference between estimated orbit and reference orbits
+     * @param velEps Tolerance on expected velocity difference
+     * @param expectedSigmasPos Expected values for covariance matrix on position
+     * @param sigmaPosEps Tolerance on expected covariance matrix on position
+     * @param expectedSigmasVel Expected values for covariance matrix on velocity
+     * @param sigmaVelEps Tolerance on expected covariance matrix on velocity
+     */
+    public static void checkKalmanFit(final KeplerianContext context, final KalmanEstimator kalman,
+                                      final List<ObservedMeasurement<?>> measurements,
+                                      final Orbit[] refOrbit, final PositionAngle[] positionAngle,
+                                      final double[] expectedDeltaPos, final double[] posEps,
+                                      final double[] expectedDeltaVel, final double []velEps,
+                                      final double[][] expectedSigmasPos,final double[] sigmaPosEps,
+                                      final double[][] expectedSigmasVel,final double[] sigmaVelEps)
+                                                      {
+
+        // Add the measurements to the Kalman filter
+        Propagator[] estimated = kalman.processMeasurements(measurements);
+        
+        // Check the number of measurements processed by the filter
+        Assert.assertEquals(measurements.size(), kalman.getCurrentMeasurementNumber());
+
+        for (int k = 0; k < refOrbit.length; ++k) {
+            // Get the last estimation
+            final Orbit    estimatedOrbit    = estimated[k].getInitialState().getOrbit();
+            final Vector3D estimatedPosition = estimatedOrbit.getPVCoordinates().getPosition();
+            final Vector3D estimatedVelocity = estimatedOrbit.getPVCoordinates().getVelocity();        
+
+            // Get the last covariance matrix estimation
+            final RealMatrix estimatedP = kalman.getPhysicalEstimatedCovarianceMatrix();
+
+            // Convert the orbital part to Cartesian formalism
+            // Assuming all 6 orbital parameters are estimated by the filter
+            final double[][] dCdY = new double[6][6];
+            estimatedOrbit.getJacobianWrtParameters(positionAngle[k], dCdY);
+            final RealMatrix Jacobian = MatrixUtils.createRealMatrix(dCdY);
+            final RealMatrix estimatedCartesianP = 
+                            Jacobian.
+                            multiply(estimatedP.getSubMatrix(0, 5, 0, 5)).
+                            multiply(Jacobian.transpose());
+
+            // Get the final sigmas (ie.sqrt of the diagonal of the Cartesian orbital covariance matrix)
+            final double[] sigmas = new double[6];
+            for (int i = 0; i < 6; i++) {
+                sigmas[i] = FastMath.sqrt(estimatedCartesianP.getEntry(i, i));
+            }
+
+            // Check the final orbit estimation & PV sigmas
+            final double deltaPosK = Vector3D.distance(refOrbit[k].getPVCoordinates().getPosition(), estimatedPosition);
+            final double deltaVelK = Vector3D.distance(refOrbit[k].getPVCoordinates().getVelocity(), estimatedVelocity);
+            Assert.assertEquals(expectedDeltaPos[k], deltaPosK, posEps[k]);
+            Assert.assertEquals(expectedDeltaVel[k], deltaVelK, velEps[k]);
+
+            for (int i = 0; i < 3; i++) {
+                System.out.println(sigmas[i]);
+                System.out.println(sigmas[i+3]);
+                Assert.assertEquals(expectedSigmasPos[k][i], sigmas[i],   sigmaPosEps[k]);
+                Assert.assertEquals(expectedSigmasVel[k][i], sigmas[i+3], sigmaVelEps[k]);
+            }
+        }
+    }
+
+}
+

src/test/java/org/orekit/estimation/leastsquares/EcksteinHechlerBatchLSEstimatorTest.java

+/* Copyright 2002-2022 CS GROUP
+ * Licensed to CS GROUP (CS) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * CS licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.orekit.estimation.leastsquares;
+
+import java.util.ArrayList;
+import java.util.List;
+
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.hipparchus.linear.RealMatrix;
+import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem.Evaluation;
+import org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer;
+import org.junit.Assert;
+import org.junit.Test;
+import org.orekit.attitudes.LofOffset;
+import org.orekit.errors.OrekitException;
+import org.orekit.errors.OrekitMessages;
+import org.orekit.estimation.EcksteinHechlerContext;
+import org.orekit.estimation.EcksteinHechlerEstimationTestUtils;
+import org.orekit.estimation.measurements.EstimationsProvider;
+import org.orekit.estimation.measurements.GroundStation;
+import org.orekit.estimation.measurements.ObservedMeasurement;
+import org.orekit.estimation.measurements.PVMeasurementCreator;
+import org.orekit.estimation.measurements.Range;
+import org.orekit.estimation.measurements.RangeMeasurementCreator;
+import org.orekit.estimation.measurements.RangeRateMeasurementCreator;
+import org.orekit.estimation.measurements.modifiers.OnBoardAntennaRangeModifier;
+import org.orekit.frames.LOFType;
+import org.orekit.orbits.Orbit;
+import org.orekit.orbits.PositionAngle;
+import org.orekit.propagation.Propagator;
+import org.orekit.propagation.conversion.EcksteinHechlerPropagatorBuilder;
+import org.orekit.utils.ParameterDriversList;
+
+public class EcksteinHechlerBatchLSEstimatorTest {
+
+    /**
+     * Perfect PV measurements with a perfect start
+     */
+    @Test
+    public void testPV() {
+
+        EcksteinHechlerContext context = EcksteinHechlerEstimationTestUtils.smallEccentricContext("regular-data:potential:tides");
+
+        final EcksteinHechlerPropagatorBuilder propagatorBuilder = context.createBuilder(PositionAngle.MEAN, true, 1.0);
+
+        // create perfect PV measurements
+        final Propagator propagator = EcksteinHechlerEstimationTestUtils.createPropagator(context.initialOrbit,
+                                                                                          propagatorBuilder);
+        final List<ObservedMeasurement<?>> measurements =
+                        EcksteinHechlerEstimationTestUtils.createMeasurements(propagator,
+                                                                             new PVMeasurementCreator(),
+                                                                             0.0, 1.0, 300.0);
+
+        // create orbit estimator
+        final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
+                                                                propagatorBuilder);
+        for (final ObservedMeasurement<?> measurement : measurements) {
+            estimator.addMeasurement(measurement);
+        }
+        estimator.setParametersConvergenceThreshold(1.0e-2);
+        estimator.setMaxIterations(10);
+        estimator.setMaxEvaluations(20);
+
+        EcksteinHechlerEstimationTestUtils.checkFit(context, estimator, 1, 1,
+                                                    0.0, 1.0e-15,
+                                                    0.0, 1.0e-15,
+                                                    0.0, 1.0e-15,
+                                                    0.0, 1.0e-15);
+
+        RealMatrix normalizedCovariances = estimator.getOptimum().getCovariances(1.0e-10);
+        RealMatrix physicalCovariances   = estimator.getPhysicalCovariances(1.0e-10);
+        Assert.assertEquals(6,       normalizedCovariances.getRowDimension());
+        Assert.assertEquals(6,       normalizedCovariances.getColumnDimension());
+        Assert.assertEquals(6,       physicalCovariances.getRowDimension());
+        Assert.assertEquals(6,       physicalCovariances.getColumnDimension());
+
+    }
+
+    /** Test PV measurements generation and backward propagation in least-square orbit determination. */
+    @Test
+    public void testKeplerPVBackward() {
+
+        EcksteinHechlerContext context = EcksteinHechlerEstimationTestUtils.smallEccentricContext("regular-data:potential:tides");
+
+        final EcksteinHechlerPropagatorBuilder propagatorBuilder =
+                        context.createBuilder(PositionAngle.MEAN, true, 1.0);
+
+        // create perfect PV measurements
+        final Propagator propagator = EcksteinHechlerEstimationTestUtils.createPropagator(context.initialOrbit,
+                                                                                         propagatorBuilder);
+        final List<ObservedMeasurement<?>> measurements =
+                        EcksteinHechlerEstimationTestUtils.createMeasurements(propagator,
+                                                               new PVMeasurementCreator(),
+                                                               0.0, -1.0, 300.0);
+
+        // create orbit estimator
+        final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
+                                                                propagatorBuilder);
+        for (final ObservedMeasurement<?> measurement : measurements) {
+            estimator.addMeasurement(measurement);
+        }
+        estimator.setParametersConvergenceThreshold(1.0e-2);
+        estimator.setMaxIterations(10);
+        estimator.setMaxEvaluations(20);
+
+        EcksteinHechlerEstimationTestUtils.checkFit(context, estimator, 1, 1,
+                                                    0.0, 1.0e-15,
+                                                    0.0, 1.0e-15,
+                                                    0.0, 1.0e-15,
+                                                    0.0, 1.0e-15);
+
+        RealMatrix normalizedCovariances = estimator.getOptimum().getCovariances(1.0e-10);
+        RealMatrix physicalCovariances   = estimator.getPhysicalCovariances(1.0e-10);
+        Assert.assertEquals(6,       normalizedCovariances.getRowDimension());
+        Assert.assertEquals(6,       normalizedCovariances.getColumnDimension());
+        Assert.assertEquals(6,       physicalCovariances.getRowDimension());
+        Assert.assertEquals(6,       physicalCovariances.getColumnDimension());
+
+    }
+
+    /**
+     * Perfect range measurements with a perfect start
+     */
+    @Test
+    public void testKeplerRange() {
+
+        EcksteinHechlerContext context = EcksteinHechlerEstimationTestUtils.smallEccentricContext("regular-data:potential:tides");
+
+        final EcksteinHechlerPropagatorBuilder propagatorBuilder =
+                        context.createBuilder(PositionAngle.MEAN, true, 1.0);
+
+        // create perfect range measurements
+        final Propagator propagator = EcksteinHechlerEstimationTestUtils.createPropagator(context.initialOrbit,
+                                                                           propagatorBuilder);
+        final List<ObservedMeasurement<?>> measurements =
+                        EcksteinHechlerEstimationTestUtils.createMeasurements(propagator,
+                                                               new RangeMeasurementCreator(context),
+                                                               1.0, 3.0, 300.0);
+
+        // create orbit estimator
+        final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
+                                                                propagatorBuilder);
+        for (final ObservedMeasurement<?> range : measurements) {
+            estimator.addMeasurement(range);
+        }
+        estimator.setParametersConvergenceThreshold(1.0e-2);
+        estimator.setMaxIterations(10);
+        estimator.setMaxEvaluations(20);
+
+        EcksteinHechlerEstimationTestUtils.checkFit(context, estimator, 1, 10,
+                                                    0.0, 1.1e-6,
+                                                    0.0, 3.0e-6,
+                                                    0.0, 5.6e-8,
+                                                    0.0, 5.4e-11);
+
+    }
+
+    /**
+     * Perfect range measurements with a perfect start and an on-board antenna range offset 
+     */
+    @Test
+    public void testKeplerRangeWithOnBoardAntennaOffset() {
+
+        EcksteinHechlerContext context = EcksteinHechlerEstimationTestUtils.smallEccentricContext("regular-data:potential:tides");
+
+        final EcksteinHechlerPropagatorBuilder propagatorBuilder =
+                        context.createBuilder(PositionAngle.MEAN, true, 1.0);
+        propagatorBuilder.setAttitudeProvider(new LofOffset(propagatorBuilder.getFrame(), LOFType.LVLH));
+        final Vector3D antennaPhaseCenter = new Vector3D(-1.2, 2.3, -0.7);
+
+        // create perfect range measurements with antenna offset
+        final Propagator propagator = EcksteinHechlerEstimationTestUtils.createPropagator(context.initialOrbit,
+                                                                           propagatorBuilder);
+        final List<ObservedMeasurement<?>> measurements =
+                        EcksteinHechlerEstimationTestUtils.createMeasurements(propagator,
+                                                               new RangeMeasurementCreator(context, antennaPhaseCenter),
+                                                               1.0, 3.0, 300.0);
+
+        // create orbit estimator
+        final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
+                                                                propagatorBuilder);
+        final OnBoardAntennaRangeModifier obaModifier = new OnBoardAntennaRangeModifier(antennaPhaseCenter);
+        for (final ObservedMeasurement<?> range : measurements) {
+            ((Range) range).addModifier(obaModifier);
+            estimator.addMeasurement(range);
+        }
+        estimator.setParametersConvergenceThreshold(1.0e-2);
+        estimator.setMaxIterations(10);
+        estimator.setMaxEvaluations(20);
+
+        EcksteinHechlerEstimationTestUtils.checkFit(context, estimator, 1, 11,
+                                                   0.0, 4.3e-5,
+                                                   0.0, 1.2e-4,
+                                                   0.0, 3.0e-8,
+                                                   0.0, 2.5e-11);
+
+    }
+
+    /**
+     * Perfect range rate measurements with a perfect start
+     */
+    @Test
+    public void testKeplerRangeRate() {
+
+        EcksteinHechlerContext context = EcksteinHechlerEstimationTestUtils.smallEccentricContext("regular-data:potential:tides");
+
+        final EcksteinHechlerPropagatorBuilder propagatorBuilder =
+                        context.createBuilder(PositionAngle.MEAN, true, 1.0);
+
+        // create perfect range rate measurements
+        final Propagator propagator = EcksteinHechlerEstimationTestUtils.createPropagator(context.initialOrbit,
+                                                                                         propagatorBuilder);
+        final double groundClockDrift =  4.8e-9;
+        for (final GroundStation station : context.stations) {
+            station.getClockDriftDriver().setValue(groundClockDrift);
+        }
+        final double satClkDrift = 3.2e-10;
+        final List<ObservedMeasurement<?>> measurements1 =
+                        EcksteinHechlerEstimationTestUtils.createMeasurements(propagator,
+                                                               new RangeRateMeasurementCreator(context, false, satClkDrift),
+                                                               1.0, 3.0, 300.0);
+
+        final List<ObservedMeasurement<?>> measurements = new ArrayList<ObservedMeasurement<?>>();
+        measurements.addAll(measurements1);
+
+        // create orbit estimator
+        final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
+                                                                propagatorBuilder);
+        for (final ObservedMeasurement<?> rangerate : measurements) {
+            estimator.addMeasurement(rangerate);
+        }
+        estimator.setParametersConvergenceThreshold(1.0e-3);
+        estimator.setMaxIterations(10);
+        estimator.setMaxEvaluations(20);
+
+        EcksteinHechlerEstimationTestUtils.checkFit(context, estimator, 1, 10,
+                                                   0.0, 7.4e-8,
+                                                   0.0, 1.3e-7,
+                                                   0.0, 3.2e-5,
+                                                   0.0, 3.2e-8);
+    }
+
+    @Test
+    public void testWrappedException() {
+
+        EcksteinHechlerContext context = EcksteinHechlerEstimationTestUtils.smallEccentricContext("regular-data:potential:tides");
+
+        final EcksteinHechlerPropagatorBuilder propagatorBuilder =
+                        context.createBuilder(PositionAngle.MEAN, true, 1.0);
+
+        // create perfect range measurements
+        final Propagator propagator = EcksteinHechlerEstimationTestUtils.createPropagator(context.initialOrbit,
+                                                                                         propagatorBuilder);
+        final List<ObservedMeasurement<?>> measurements =
+                        EcksteinHechlerEstimationTestUtils.createMeasurements(propagator,
+                                                                             new RangeMeasurementCreator(context),
+                                                                             1.0, 3.0, 300.0);
+
+        // create orbit estimator
+        final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
+                                                                propagatorBuilder);
+        for (final ObservedMeasurement<?> range : measurements) {
+            estimator.addMeasurement(range);
+        }
+        estimator.setParametersConvergenceThreshold(1.0e-2);
+        estimator.setMaxIterations(10);
+        estimator.setMaxEvaluations(20);
+        estimator.setObserver(new BatchLSObserver() {
+            /** {@inheritDoc} */
+            @Override
+            public void evaluationPerformed(int iterationsCount, int evaluationscount,
+                                           Orbit[] orbits,
+                                           ParameterDriversList estimatedOrbitalParameters,
+                                           ParameterDriversList estimatedPropagatorParameters,
+                                           ParameterDriversList estimatedMeasurementsParameters,
+                                           EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) {
+                throw new DummyException();
+            }
+        });
+
+        try {
+            EcksteinHechlerEstimationTestUtils.checkFit(context, estimator, 3, 4,
+                                                       0.0, 1.5e-6,
+                                                       0.0, 3.2e-6,
+                                                       0.0, 3.8e-7,
+                                                       0.0, 1.5e-10);
+            Assert.fail("an exception should have been thrown");
+        } catch (DummyException de) {
+            // expected
+        }
+
+    }
+
+    private static class DummyException extends OrekitException {
+        private static final long serialVersionUID = 1L;
+        public DummyException() {
+            super(OrekitMessages.INTERNAL_ERROR);
+        }
+    }
+
+}

src/test/java/org/orekit/estimation/leastsquares/KeplerianBatchLSEstimatorTest.java

+/* Copyright 2002-2022 CS GROUP
+ * Licensed to CS GROUP (CS) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * CS licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.orekit.estimation.leastsquares;
+
+import java.util.ArrayList;
+import java.util.List;
+
+import org.hipparchus.geometry.euclidean.threed.Vector3D;
+import org.hipparchus.linear.RealMatrix;
+import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem.Evaluation;
+import org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer;
+import org.junit.Assert;
+import org.junit.Test;
+import org.orekit.attitudes.LofOffset;
+import org.orekit.errors.OrekitException;
+import org.orekit.errors.OrekitMessages;
+import org.orekit.estimation.KeplerianContext;
+import org.orekit.estimation.KeplerianEstimationTestUtils;
+import org.orekit.estimation.measurements.EstimationsProvider;
+import org.orekit.estimation.measurements.GroundStation;
+import org.orekit.estimation.measurements.ObservedMeasurement;
+import org.orekit.estimation.measurements.PVMeasurementCreator;
+import org.orekit.estimation.measurements.Range;
+import org.orekit.estimation.measurements.RangeMeasurementCreator;
+import org.orekit.estimation.measurements.RangeRateMeasurementCreator;
+import org.orekit.estimation.measurements.modifiers.OnBoardAntennaRangeModifier;
+import org.orekit.frames.LOFType;
+import org.orekit.orbits.Orbit;
+import org.orekit.orbits.PositionAngle;
+import org.orekit.propagation.Propagator;
+import org.orekit.propagation.conversion.KeplerianPropagatorBuilder;
+import org.orekit.utils.ParameterDriversList;
+
+public class KeplerianBatchLSEstimatorTest {
+
+    /**
+     * Perfect PV measurements with a perfect start
+     */
+    @Test
+    public void testPV() {
+
+        KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
+
+        final KeplerianPropagatorBuilder propagatorBuilder = context.createBuilder(PositionAngle.MEAN, true, 1.0);
+
+        // create perfect PV measurements
+        final Propagator propagator = KeplerianEstimationTestUtils.createPropagator(context.initialOrbit,
+                                                                                         propagatorBuilder);
+        final List<ObservedMeasurement<?>> measurements =
+                        KeplerianEstimationTestUtils.createMeasurements(propagator,
+                                                                             new PVMeasurementCreator(),
+                                                                             0.0, 1.0, 300.0);
+
+        // create orbit estimator
+        final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
+                                                                propagatorBuilder);
+        for (final ObservedMeasurement<?> measurement : measurements) {
+            estimator.addMeasurement(measurement);
+        }
+        estimator.setParametersConvergenceThreshold(1.0e-2);
+        estimator.setMaxIterations(10);
+        estimator.setMaxEvaluations(20);
+
+        KeplerianEstimationTestUtils.checkFit(context, estimator, 1, 1,
+                                                   0.0, 1.0e-15,
+                                                   0.0, 1.0e-15,
+                                                   0.0, 1.0e-15,
+                                                   0.0, 1.0e-15);
+
+        RealMatrix normalizedCovariances = estimator.getOptimum().getCovariances(1.0e-10);
+        RealMatrix physicalCovariances   = estimator.getPhysicalCovariances(1.0e-10);
+        Assert.assertEquals(6,       normalizedCovariances.getRowDimension());
+        Assert.assertEquals(6,       normalizedCovariances.getColumnDimension());
+        Assert.assertEquals(6,       physicalCovariances.getRowDimension());
+        Assert.assertEquals(6,       physicalCovariances.getColumnDimension());
+
+    }
+
+    /** Test PV measurements generation and backward propagation in least-square orbit determination. */
+    @Test
+    public void testKeplerPVBackward() {
+
+        KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
+
+        final KeplerianPropagatorBuilder propagatorBuilder =
+                        context.createBuilder(PositionAngle.MEAN, true, 1.0);
+
+        // create perfect PV measurements
+        final Propagator propagator = KeplerianEstimationTestUtils.createPropagator(context.initialOrbit,
+                                                                                         propagatorBuilder);
+        final List<ObservedMeasurement<?>> measurements =
+                        KeplerianEstimationTestUtils.createMeasurements(propagator,
+                                                               new PVMeasurementCreator(),
+                                                               0.0, -1.0, 300.0);
+
+        // create orbit estimator
+        final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
+                                                                propagatorBuilder);
+        for (final ObservedMeasurement<?> measurement : measurements) {
+            estimator.addMeasurement(measurement);
+        }
+        estimator.setParametersConvergenceThreshold(1.0e-2);
+        estimator.setMaxIterations(10);
+        estimator.setMaxEvaluations(20);
+
+        KeplerianEstimationTestUtils.checkFit(context, estimator, 1, 1,
+                                                   0.0, 1.0e-15,
+                                                   0.0, 1.0e-15,
+                                                   0.0, 1.0e-15,
+                                                   0.0, 1.0e-15);
+
+        RealMatrix normalizedCovariances = estimator.getOptimum().getCovariances(1.0e-10);
+        RealMatrix physicalCovariances   = estimator.getPhysicalCovariances(1.0e-10);
+        Assert.assertEquals(6,       normalizedCovariances.getRowDimension());
+        Assert.assertEquals(6,       normalizedCovariances.getColumnDimension());
+        Assert.assertEquals(6,       physicalCovariances.getRowDimension());
+        Assert.assertEquals(6,       physicalCovariances.getColumnDimension());
+        Assert.assertEquals(0.0,     physicalCovariances.getEntry(0, 0), 1.7e-15);
+
+    }
+
+    /**
+     * Perfect range measurements with a perfect start
+     */
+    @Test
+    public void testKeplerRange() {
+
+        KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
+
+        final KeplerianPropagatorBuilder propagatorBuilder =
+                        context.createBuilder(PositionAngle.MEAN, true, 1.0);
+
+        // create perfect range measurements
+        final Propagator propagator = KeplerianEstimationTestUtils.createPropagator(context.initialOrbit,
+                                                                           propagatorBuilder);
+        final List<ObservedMeasurement<?>> measurements =
+                        KeplerianEstimationTestUtils.createMeasurements(propagator,
+                                                               new RangeMeasurementCreator(context),
+                                                               1.0, 3.0, 300.0);
+
+        // create orbit estimator
+        final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
+                                                                propagatorBuilder);
+        for (final ObservedMeasurement<?> range : measurements) {
+            estimator.addMeasurement(range);
+        }
+        estimator.setParametersConvergenceThreshold(1.0e-2);
+        estimator.setMaxIterations(10);
+        estimator.setMaxEvaluations(20);
+
+        KeplerianEstimationTestUtils.checkFit(context, estimator, 1, 5,
+                                                   0.0, 8.4e-7,
+                                                   0.0, 2.0e-6,
+                                                   0.0, 7.7e-9,
+                                                   0.0, 4.9e-12);
+
+    }
+
+    /**
+     * Perfect range measurements with a perfect start and an on-board antenna range offset 
+     */
+    @Test
+    public void testKeplerRangeWithOnBoardAntennaOffset() {
+
+        KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
+
+        final KeplerianPropagatorBuilder propagatorBuilder =
+                        context.createBuilder(PositionAngle.MEAN, true, 1.0);
+        propagatorBuilder.setAttitudeProvider(new LofOffset(propagatorBuilder.getFrame(), LOFType.LVLH));
+        final Vector3D antennaPhaseCenter = new Vector3D(-1.2, 2.3, -0.7);
+
+        // create perfect range measurements with antenna offset
+        final Propagator propagator = KeplerianEstimationTestUtils.createPropagator(context.initialOrbit,
+                                                                           propagatorBuilder);
+        final List<ObservedMeasurement<?>> measurements =
+                        KeplerianEstimationTestUtils.createMeasurements(propagator,
+                                                               new RangeMeasurementCreator(context, antennaPhaseCenter),
+                                                               1.0, 3.0, 300.0);
+
+        // create orbit estimator
+        final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
+                                                                propagatorBuilder);
+        final OnBoardAntennaRangeModifier obaModifier = new OnBoardAntennaRangeModifier(antennaPhaseCenter);
+        for (final ObservedMeasurement<?> range : measurements) {
+            ((Range) range).addModifier(obaModifier);
+            estimator.addMeasurement(range);
+        }
+        estimator.setParametersConvergenceThreshold(1.0e-2);
+        estimator.setMaxIterations(10);
+        estimator.setMaxEvaluations(20);
+
+        KeplerianEstimationTestUtils.checkFit(context, estimator, 1, 11,
+                                                   0.0, 5.9e-5,
+                                                   0.0, 1.5e-4,
+                                                   0.0, 7.2e-9,
+                                                   0.0, 7.8e-12);
+
+    }
+
+    /**
+     * Perfect range rate measurements with a perfect start
+     */
+    @Test
+    public void testKeplerRangeRate() {
+
+        KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
+
+        final KeplerianPropagatorBuilder propagatorBuilder =
+                        context.createBuilder(PositionAngle.MEAN, true, 1.0);
+
+        // create perfect range rate measurements
+        final Propagator propagator = KeplerianEstimationTestUtils.createPropagator(context.initialOrbit,
+                                                                                         propagatorBuilder);
+        final double groundClockDrift =  4.8e-9;
+        for (final GroundStation station : context.stations) {
+            station.getClockDriftDriver().setValue(groundClockDrift);
+        }
+        final double satClkDrift = 3.2e-10;
+        final List<ObservedMeasurement<?>> measurements1 =
+                        KeplerianEstimationTestUtils.createMeasurements(propagator,
+                                                               new RangeRateMeasurementCreator(context, false, satClkDrift),
+                                                               1.0, 3.0, 300.0);
+
+        final List<ObservedMeasurement<?>> measurements = new ArrayList<ObservedMeasurement<?>>();
+        measurements.addAll(measurements1);
+
+        // create orbit estimator
+        final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
+                                                                propagatorBuilder);
+        for (final ObservedMeasurement<?> rangerate : measurements) {
+            estimator.addMeasurement(rangerate);
+        }
+        estimator.setParametersConvergenceThreshold(1.0e-3);
+        estimator.setMaxIterations(10);
+        estimator.setMaxEvaluations(20);
+
+        KeplerianEstimationTestUtils.checkFit(context, estimator, 1, 3,
+                                                   0.0, 5.6e-7,
+                                                   0.0, 7.7e-7,
+                                                   0.0, 8.7e-5,
+                                                   0.0, 3.5e-8);
+    }
+
+    @Test
+    public void testWrappedException() {
+
+        KeplerianContext context = KeplerianEstimationTestUtils.eccentricContext("regular-data:potential:tides");
+
+        final KeplerianPropagatorBuilder propagatorBuilder =
+                        context.createBuilder(PositionAngle.MEAN, true, 1.0);
+
+        // create perfect range measurements
+        final Propagator propagator = KeplerianEstimationTestUtils.createPropagator(context.initialOrbit,
+                                                                                         propagatorBuilder);
+        final List<ObservedMeasurement<?>> measurements =
+                        KeplerianEstimationTestUtils.createMeasurements(propagator,
+                                                                             new RangeMeasurementCreator(context),
+                                                                             1.0, 3.0, 300.0);
+
+        // create orbit estimator
+        final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(),
+                                                                propagatorBuilder);
+        for (final ObservedMeasurement<?> range : measurements) {
+            estimator.addMeasurement(range);
+        }
+        estimator.setParametersConvergenceThreshold(1.0e-2);
+        estimator.setMaxIterations(10);
+        estimator.setMaxEvaluations(20);
+        estimator.setObserver(new BatchLSObserver() {
+            /** {@inheritDoc} */
+            @Override
+            public void evaluationPerformed(int iterationsCount, int evaluationscount,
+                                           Orbit[] orbits,
+                                           ParameterDriversList estimatedOrbitalParameters,
+                                           ParameterDriversList estimatedPropagatorParameters,
+                                           ParameterDriversList estimatedMeasurementsParameters,
+                                           EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) {
+                throw new DummyException();
+            }
+        });
+
+        try {
+            KeplerianEstimationTestUtils.checkFit(context, estimator, 3, 4,
+                                                       0.0, 1.5e-6,
+                                                       0.0, 3.2e-6,
+                                                       0.0, 3.8e-7,
+                                                       0.0, 1.5e-10);
+            Assert.fail("an exception should have been thrown");
+        } catch (DummyException de) {
+            // expected
+        }
+
+    }
+
+    private static class DummyException extends OrekitException {
+        private static final long serialVersionUID = 1L;
+        public DummyException() {
+            super(OrekitMessages.INTERNAL_ERROR);
+        }
+    }
+
+}

src/test/java/org/orekit/estimation/sequential/BrouwerLyddaneKalmanEstimatorTest.java

+/* Copyright 2002-2022 CS GROUP
+ * Licensed to CS GROUP (CS) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * CS licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *   http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
 package org.orekit.estimation.sequential;
 
 import java.util.List;