Search in sources :

Example 1 with EstimationsProvider

use of org.orekit.estimation.measurements.EstimationsProvider in project Orekit by CS-SI.

the class OrbitDeterminationTest method run.

private ResultOD run(final File input, final boolean print) throws IOException, IllegalArgumentException, OrekitException, ParseException {
    // read input parameters
    KeyValueFileParser<ParameterKey> parser = new KeyValueFileParser<ParameterKey>(ParameterKey.class);
    parser.parseInput(input.getAbsolutePath(), new FileInputStream(input));
    // log file
    final RangeLog rangeLog = new RangeLog();
    final RangeRateLog rangeRateLog = new RangeRateLog();
    final AzimuthLog azimuthLog = new AzimuthLog();
    final ElevationLog elevationLog = new ElevationLog();
    final PositionLog positionLog = new PositionLog();
    final VelocityLog velocityLog = new VelocityLog();
    // gravity field
    GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-5c.gfc", true));
    final NormalizedSphericalHarmonicsProvider gravityField = createGravityField(parser);
    // Orbit initial guess
    final Orbit initialGuess = createOrbit(parser, gravityField.getMu());
    // IERS conventions
    final IERSConventions conventions;
    if (!parser.containsKey(ParameterKey.IERS_CONVENTIONS)) {
        conventions = IERSConventions.IERS_2010;
    } else {
        conventions = IERSConventions.valueOf("IERS_" + parser.getInt(ParameterKey.IERS_CONVENTIONS));
    }
    // central body
    final OneAxisEllipsoid body = createBody(parser);
    // propagator builder
    final NumericalPropagatorBuilder propagatorBuilder = createPropagatorBuilder(parser, conventions, gravityField, body, initialGuess);
    // estimator
    final BatchLSEstimator estimator = createEstimator(parser, propagatorBuilder);
    // measurements
    final List<ObservedMeasurement<?>> measurements = new ArrayList<ObservedMeasurement<?>>();
    for (final String fileName : parser.getStringsList(ParameterKey.MEASUREMENTS_FILES, ',')) {
        measurements.addAll(readMeasurements(new File(input.getParentFile(), fileName), createStationsData(parser, body), createPVData(parser), createSatRangeBias(parser), createWeights(parser), createRangeOutliersManager(parser), createRangeRateOutliersManager(parser), createAzElOutliersManager(parser), createPVOutliersManager(parser)));
    }
    for (ObservedMeasurement<?> measurement : measurements) {
        estimator.addMeasurement(measurement);
    }
    if (print) {
        estimator.setObserver(new BatchLSObserver() {

            private PVCoordinates previousPV;

            {
                previousPV = initialGuess.getPVCoordinates();
                final String header = "iteration evaluations      ΔP(m)        ΔV(m/s)           RMS          nb Range    nb Range-rate  nb Angular     nb PV%n";
                System.out.format(Locale.US, header);
            }

            /**
             * {@inheritDoc}
             */
            @Override
            public void evaluationPerformed(final int iterationsCount, final int evaluationsCount, final Orbit[] orbits, final ParameterDriversList estimatedOrbitalParameters, final ParameterDriversList estimatedPropagatorParameters, final ParameterDriversList estimatedMeasurementsParameters, final EstimationsProvider evaluationsProvider, final LeastSquaresProblem.Evaluation lspEvaluation) {
                PVCoordinates currentPV = orbits[0].getPVCoordinates();
                final String format0 = "    %2d         %2d                                 %16.12f     %s       %s     %s     %s%n";
                final String format = "    %2d         %2d      %13.6f %12.9f %16.12f     %s       %s     %s     %s%n";
                final EvaluationCounter<Range> rangeCounter = new EvaluationCounter<Range>();
                final EvaluationCounter<RangeRate> rangeRateCounter = new EvaluationCounter<RangeRate>();
                final EvaluationCounter<AngularAzEl> angularCounter = new EvaluationCounter<AngularAzEl>();
                final EvaluationCounter<PV> pvCounter = new EvaluationCounter<PV>();
                for (final Map.Entry<ObservedMeasurement<?>, EstimatedMeasurement<?>> entry : estimator.getLastEstimations().entrySet()) {
                    if (entry.getKey() instanceof Range) {
                        @SuppressWarnings("unchecked") EstimatedMeasurement<Range> evaluation = (EstimatedMeasurement<Range>) entry.getValue();
                        rangeCounter.add(evaluation);
                    } else if (entry.getKey() instanceof RangeRate) {
                        @SuppressWarnings("unchecked") EstimatedMeasurement<RangeRate> evaluation = (EstimatedMeasurement<RangeRate>) entry.getValue();
                        rangeRateCounter.add(evaluation);
                    } else if (entry.getKey() instanceof AngularAzEl) {
                        @SuppressWarnings("unchecked") EstimatedMeasurement<AngularAzEl> evaluation = (EstimatedMeasurement<AngularAzEl>) entry.getValue();
                        angularCounter.add(evaluation);
                    } else if (entry.getKey() instanceof PV) {
                        @SuppressWarnings("unchecked") EstimatedMeasurement<PV> evaluation = (EstimatedMeasurement<PV>) entry.getValue();
                        pvCounter.add(evaluation);
                    }
                }
                if (evaluationsCount == 1) {
                    System.out.format(Locale.US, format0, iterationsCount, evaluationsCount, lspEvaluation.getRMS(), rangeCounter.format(8), rangeRateCounter.format(8), angularCounter.format(8), pvCounter.format(8));
                } else {
                    System.out.format(Locale.US, format, iterationsCount, evaluationsCount, Vector3D.distance(previousPV.getPosition(), currentPV.getPosition()), Vector3D.distance(previousPV.getVelocity(), currentPV.getVelocity()), lspEvaluation.getRMS(), rangeCounter.format(8), rangeRateCounter.format(8), angularCounter.format(8), pvCounter.format(8));
                }
                previousPV = currentPV;
            }
        });
    }
    Orbit estimated = estimator.estimate()[0].getInitialState().getOrbit();
    // compute some statistics
    for (final Map.Entry<ObservedMeasurement<?>, EstimatedMeasurement<?>> entry : estimator.getLastEstimations().entrySet()) {
        if (entry.getKey() instanceof Range) {
            @SuppressWarnings("unchecked") EstimatedMeasurement<Range> evaluation = (EstimatedMeasurement<Range>) entry.getValue();
            rangeLog.add(evaluation);
        } else if (entry.getKey() instanceof RangeRate) {
            @SuppressWarnings("unchecked") EstimatedMeasurement<RangeRate> evaluation = (EstimatedMeasurement<RangeRate>) entry.getValue();
            rangeRateLog.add(evaluation);
        } else if (entry.getKey() instanceof AngularAzEl) {
            @SuppressWarnings("unchecked") EstimatedMeasurement<AngularAzEl> evaluation = (EstimatedMeasurement<AngularAzEl>) entry.getValue();
            azimuthLog.add(evaluation);
            elevationLog.add(evaluation);
        } else if (entry.getKey() instanceof PV) {
            @SuppressWarnings("unchecked") EstimatedMeasurement<PV> evaluation = (EstimatedMeasurement<PV>) entry.getValue();
            positionLog.add(evaluation);
            velocityLog.add(evaluation);
        }
    }
    // parmaters and measurements.
    final ParameterDriversList propagatorParameters = estimator.getPropagatorParametersDrivers(true);
    final ParameterDriversList measurementsParameters = estimator.getMeasurementsParametersDrivers(true);
    // instation of results
    return new ResultOD(propagatorParameters, measurementsParameters, estimator.getIterationsCount(), estimator.getEvaluationsCount(), estimated.getPVCoordinates(), rangeLog.createStatisticsSummary(), rangeRateLog.createStatisticsSummary(), azimuthLog.createStatisticsSummary(), elevationLog.createStatisticsSummary(), positionLog.createStatisticsSummary(), velocityLog.createStatisticsSummary(), estimator.getPhysicalCovariances(1.0e-10));
}
Also used : OneAxisEllipsoid(org.orekit.bodies.OneAxisEllipsoid) ICGEMFormatReader(org.orekit.forces.gravity.potential.ICGEMFormatReader) PV(org.orekit.estimation.measurements.PV) ArrayList(java.util.ArrayList) PVCoordinates(org.orekit.utils.PVCoordinates) TimeStampedPVCoordinates(org.orekit.utils.TimeStampedPVCoordinates) ObservedMeasurement(org.orekit.estimation.measurements.ObservedMeasurement) EstimatedMeasurement(org.orekit.estimation.measurements.EstimatedMeasurement) CartesianOrbit(org.orekit.orbits.CartesianOrbit) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) Orbit(org.orekit.orbits.Orbit) CircularOrbit(org.orekit.orbits.CircularOrbit) EquinoctialOrbit(org.orekit.orbits.EquinoctialOrbit) Range(org.orekit.estimation.measurements.Range) FileInputStream(java.io.FileInputStream) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) RangeRate(org.orekit.estimation.measurements.RangeRate) File(java.io.File) AngularAzEl(org.orekit.estimation.measurements.AngularAzEl) Map(java.util.Map) HashMap(java.util.HashMap) ParameterDriversList(org.orekit.utils.ParameterDriversList) LeastSquaresProblem(org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem) NormalizedSphericalHarmonicsProvider(org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider) EstimationsProvider(org.orekit.estimation.measurements.EstimationsProvider) KeyValueFileParser(org.orekit.KeyValueFileParser) IERSConventions(org.orekit.utils.IERSConventions) GeodeticPoint(org.orekit.bodies.GeodeticPoint)

Example 2 with EstimationsProvider

use of org.orekit.estimation.measurements.EstimationsProvider in project Orekit by CS-SI.

the class BatchLSEstimatorTest method testMultiSatWithParameters.

/**
 * A modified version of the previous test with a selection of propagation drivers to estimate
 *  One common (µ)
 *  Some specifics for each satellite (Cr and Ca)
 *
 * @throws OrekitException
 */
@Test
public void testMultiSatWithParameters() throws OrekitException {
    // Test: Set the propagator drivers to estimate for each satellite
    final boolean muEstimated = true;
    final boolean crEstimated1 = true;
    final boolean caEstimated1 = true;
    final boolean crEstimated2 = true;
    final boolean caEstimated2 = false;
    // Builder sat 1
    final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    final NumericalPropagatorBuilder propagatorBuilder1 = context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 1.0, Force.POTENTIAL, Force.SOLAR_RADIATION_PRESSURE);
    // Adding selection of parameters
    String satName = "sat 1";
    for (DelegatingDriver driver : propagatorBuilder1.getPropagationParametersDrivers().getDrivers()) {
        if (driver.getName().equals("central attraction coefficient")) {
            driver.setSelected(muEstimated);
        }
        if (driver.getName().equals(RadiationSensitive.REFLECTION_COEFFICIENT)) {
            driver.setName(driver.getName() + " " + satName);
            driver.setSelected(crEstimated1);
        }
        if (driver.getName().equals(RadiationSensitive.ABSORPTION_COEFFICIENT)) {
            driver.setName(driver.getName() + " " + satName);
            driver.setSelected(caEstimated1);
        }
    }
    // Builder for sat 2
    final Context context2 = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    final NumericalPropagatorBuilder propagatorBuilder2 = context2.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 1.0, Force.POTENTIAL, Force.SOLAR_RADIATION_PRESSURE);
    // Adding selection of parameters
    satName = "sat 2";
    for (ParameterDriver driver : propagatorBuilder2.getPropagationParametersDrivers().getDrivers()) {
        if (driver.getName().equals("central attraction coefficient")) {
            driver.setSelected(muEstimated);
        }
        if (driver.getName().equals(RadiationSensitive.REFLECTION_COEFFICIENT)) {
            driver.setName(driver.getName() + " " + satName);
            driver.setSelected(crEstimated2);
        }
        if (driver.getName().equals(RadiationSensitive.ABSORPTION_COEFFICIENT)) {
            driver.setName(driver.getName() + " " + satName);
            driver.setSelected(caEstimated2);
        }
    }
    // Create perfect inter-satellites range measurements
    final TimeStampedPVCoordinates original = context.initialOrbit.getPVCoordinates();
    final Orbit closeOrbit = new CartesianOrbit(new TimeStampedPVCoordinates(context.initialOrbit.getDate(), original.getPosition().add(new Vector3D(1000, 2000, 3000)), original.getVelocity().add(new Vector3D(-0.03, 0.01, 0.02))), context.initialOrbit.getFrame(), context.initialOrbit.getMu());
    final Propagator closePropagator = EstimationTestUtils.createPropagator(closeOrbit, propagatorBuilder2);
    closePropagator.setEphemerisMode();
    closePropagator.propagate(context.initialOrbit.getDate().shiftedBy(3.5 * closeOrbit.getKeplerianPeriod()));
    final BoundedPropagator ephemeris = closePropagator.getGeneratedEphemeris();
    Propagator propagator1 = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder1);
    final List<ObservedMeasurement<?>> r12 = EstimationTestUtils.createMeasurements(propagator1, new InterSatellitesRangeMeasurementCreator(ephemeris), 1.0, 3.0, 300.0);
    // create perfect range measurements for first satellite
    propagator1 = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder1);
    final List<ObservedMeasurement<?>> r1 = EstimationTestUtils.createMeasurements(propagator1, new RangeMeasurementCreator(context), 1.0, 3.0, 300.0);
    // create orbit estimator
    final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(), propagatorBuilder1, propagatorBuilder2);
    for (final ObservedMeasurement<?> interSat : r12) {
        estimator.addMeasurement(interSat);
    }
    for (final ObservedMeasurement<?> range : r1) {
        estimator.addMeasurement(range);
    }
    estimator.setParametersConvergenceThreshold(1.0e-2);
    estimator.setMaxIterations(10);
    estimator.setMaxEvaluations(20);
    estimator.setObserver(new BatchLSObserver() {

        int lastIter = 0;

        int lastEval = 0;

        /**
         * {@inheritDoc}
         */
        @Override
        public void evaluationPerformed(int iterationsCount, int evaluationscount, Orbit[] orbits, ParameterDriversList estimatedOrbitalParameters, ParameterDriversList estimatedPropagatorParameters, ParameterDriversList estimatedMeasurementsParameters, EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) throws OrekitException {
            if (iterationsCount == lastIter) {
                Assert.assertEquals(lastEval + 1, evaluationscount);
            } else {
                Assert.assertEquals(lastIter + 1, iterationsCount);
            }
            lastIter = iterationsCount;
            lastEval = evaluationscount;
            AbsoluteDate previous = AbsoluteDate.PAST_INFINITY;
            for (int i = 0; i < evaluationsProvider.getNumber(); ++i) {
                AbsoluteDate current = evaluationsProvider.getEstimatedMeasurement(i).getDate();
                Assert.assertTrue(current.compareTo(previous) >= 0);
                previous = current;
            }
        }
    });
    List<DelegatingDriver> parameters = estimator.getOrbitalParametersDrivers(true).getDrivers();
    ParameterDriver a0Driver = parameters.get(0);
    Assert.assertEquals("a[0]", a0Driver.getName());
    a0Driver.setValue(a0Driver.getValue() + 1.2);
    a0Driver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
    ParameterDriver a1Driver = parameters.get(6);
    Assert.assertEquals("a[1]", a1Driver.getName());
    a1Driver.setValue(a1Driver.getValue() - 5.4);
    a1Driver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
    final Orbit before = new KeplerianOrbit(parameters.get(6).getValue(), parameters.get(7).getValue(), parameters.get(8).getValue(), parameters.get(9).getValue(), parameters.get(10).getValue(), parameters.get(11).getValue(), PositionAngle.TRUE, closeOrbit.getFrame(), closeOrbit.getDate(), closeOrbit.getMu());
    Assert.assertEquals(4.7246, Vector3D.distance(closeOrbit.getPVCoordinates().getPosition(), before.getPVCoordinates().getPosition()), 1.0e-3);
    Assert.assertEquals(0.0010514, Vector3D.distance(closeOrbit.getPVCoordinates().getVelocity(), before.getPVCoordinates().getVelocity()), 1.0e-6);
    EstimationTestUtils.checkFit(context, estimator, 4, 5, 0.0, 6.0e-06, 0.0, 1.7e-05, 0.0, 4.4e-07, 0.0, 1.7e-10);
    final Orbit determined = new KeplerianOrbit(parameters.get(6).getValue(), parameters.get(7).getValue(), parameters.get(8).getValue(), parameters.get(9).getValue(), parameters.get(10).getValue(), parameters.get(11).getValue(), PositionAngle.TRUE, closeOrbit.getFrame(), closeOrbit.getDate(), closeOrbit.getMu());
    Assert.assertEquals(0.0, Vector3D.distance(closeOrbit.getPVCoordinates().getPosition(), determined.getPVCoordinates().getPosition()), 5.8e-6);
    Assert.assertEquals(0.0, Vector3D.distance(closeOrbit.getPVCoordinates().getVelocity(), determined.getPVCoordinates().getVelocity()), 3.5e-9);
    // got a default one
    for (final ParameterDriver driver : estimator.getOrbitalParametersDrivers(true).getDrivers()) {
        if (driver.getName().startsWith("a[")) {
            // user-specified reference date
            Assert.assertEquals(0, driver.getReferenceDate().durationFrom(AbsoluteDate.GALILEO_EPOCH), 1.0e-15);
        } else {
            // default reference date
            Assert.assertEquals(0, driver.getReferenceDate().durationFrom(propagatorBuilder1.getInitialOrbitDate()), 1.0e-15);
        }
    }
}
Also used : CartesianOrbit(org.orekit.orbits.CartesianOrbit) TimeStampedPVCoordinates(org.orekit.utils.TimeStampedPVCoordinates) AbsoluteDate(org.orekit.time.AbsoluteDate) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) ParameterDriversList(org.orekit.utils.ParameterDriversList) BoundedPropagator(org.orekit.propagation.BoundedPropagator) Propagator(org.orekit.propagation.Propagator) OrekitException(org.orekit.errors.OrekitException) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) DelegatingDriver(org.orekit.utils.ParameterDriversList.DelegatingDriver) BoundedPropagator(org.orekit.propagation.BoundedPropagator) ObservedMeasurement(org.orekit.estimation.measurements.ObservedMeasurement) EstimationsProvider(org.orekit.estimation.measurements.EstimationsProvider) Context(org.orekit.estimation.Context) Evaluation(org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem.Evaluation) Orbit(org.orekit.orbits.Orbit) CartesianOrbit(org.orekit.orbits.CartesianOrbit) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) ParameterDriver(org.orekit.utils.ParameterDriver) LevenbergMarquardtOptimizer(org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) InterSatellitesRangeMeasurementCreator(org.orekit.estimation.measurements.InterSatellitesRangeMeasurementCreator) RangeMeasurementCreator(org.orekit.estimation.measurements.RangeMeasurementCreator) InterSatellitesRangeMeasurementCreator(org.orekit.estimation.measurements.InterSatellitesRangeMeasurementCreator) Test(org.junit.Test)

Example 3 with EstimationsProvider

use of org.orekit.estimation.measurements.EstimationsProvider in project Orekit by CS-SI.

the class BatchLSEstimatorTest method testWrappedException.

@Test
public void testWrappedException() throws OrekitException {
    Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 1.0);
    // create perfect range measurements
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.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) throws DummyException {
            throw new DummyException();
        }
    });
    try {
        EstimationTestUtils.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
    }
}
Also used : Context(org.orekit.estimation.Context) Evaluation(org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem.Evaluation) Orbit(org.orekit.orbits.Orbit) CartesianOrbit(org.orekit.orbits.CartesianOrbit) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) LevenbergMarquardtOptimizer(org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer) ParameterDriversList(org.orekit.utils.ParameterDriversList) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) BoundedPropagator(org.orekit.propagation.BoundedPropagator) Propagator(org.orekit.propagation.Propagator) RangeMeasurementCreator(org.orekit.estimation.measurements.RangeMeasurementCreator) InterSatellitesRangeMeasurementCreator(org.orekit.estimation.measurements.InterSatellitesRangeMeasurementCreator) ObservedMeasurement(org.orekit.estimation.measurements.ObservedMeasurement) EstimationsProvider(org.orekit.estimation.measurements.EstimationsProvider) Test(org.junit.Test)

Example 4 with EstimationsProvider

use of org.orekit.estimation.measurements.EstimationsProvider in project Orekit by CS-SI.

the class OrbitDetermination method run.

private void run(final File input, final File home) throws IOException, IllegalArgumentException, OrekitException, ParseException {
    // read input parameters
    KeyValueFileParser<ParameterKey> parser = new KeyValueFileParser<ParameterKey>(ParameterKey.class);
    try (final FileInputStream fis = new FileInputStream(input)) {
        parser.parseInput(input.getAbsolutePath(), fis);
    }
    // log file
    final String baseName;
    final PrintStream logStream;
    if (parser.containsKey(ParameterKey.OUTPUT_BASE_NAME) && parser.getString(ParameterKey.OUTPUT_BASE_NAME).length() > 0) {
        baseName = parser.getString(ParameterKey.OUTPUT_BASE_NAME);
        logStream = new PrintStream(new File(home, baseName + "-log.out"), "UTF-8");
    } else {
        baseName = null;
        logStream = null;
    }
    final RangeLog rangeLog = new RangeLog(home, baseName);
    final RangeRateLog rangeRateLog = new RangeRateLog(home, baseName);
    final AzimuthLog azimuthLog = new AzimuthLog(home, baseName);
    final ElevationLog elevationLog = new ElevationLog(home, baseName);
    final PositionLog positionLog = new PositionLog(home, baseName);
    final VelocityLog velocityLog = new VelocityLog(home, baseName);
    try {
        // gravity field
        final NormalizedSphericalHarmonicsProvider gravityField = createGravityField(parser);
        // Orbit initial guess
        final Orbit initialGuess = createOrbit(parser, gravityField.getMu());
        // IERS conventions
        final IERSConventions conventions;
        if (!parser.containsKey(ParameterKey.IERS_CONVENTIONS)) {
            conventions = IERSConventions.IERS_2010;
        } else {
            conventions = IERSConventions.valueOf("IERS_" + parser.getInt(ParameterKey.IERS_CONVENTIONS));
        }
        // central body
        final OneAxisEllipsoid body = createBody(parser);
        // propagator builder
        final NumericalPropagatorBuilder propagatorBuilder = createPropagatorBuilder(parser, conventions, gravityField, body, initialGuess);
        // estimator
        final BatchLSEstimator estimator = createEstimator(parser, propagatorBuilder);
        // measurements
        final List<ObservedMeasurement<?>> measurements = new ArrayList<ObservedMeasurement<?>>();
        for (final String fileName : parser.getStringsList(ParameterKey.MEASUREMENTS_FILES, ',')) {
            measurements.addAll(readMeasurements(new File(input.getParentFile(), fileName), createStationsData(parser, body), createPVData(parser), createSatRangeBias(parser), createWeights(parser), createRangeOutliersManager(parser), createRangeRateOutliersManager(parser), createAzElOutliersManager(parser), createPVOutliersManager(parser)));
        }
        for (ObservedMeasurement<?> measurement : measurements) {
            estimator.addMeasurement(measurement);
        }
        // estimate orbit
        estimator.setObserver(new BatchLSObserver() {

            private PVCoordinates previousPV;

            {
                previousPV = initialGuess.getPVCoordinates();
                final String header = "iteration evaluations      ΔP(m)        ΔV(m/s)           RMS          nb Range    nb Range-rate  nb Angular     nb PV%n";
                System.out.format(Locale.US, header);
                if (logStream != null) {
                    logStream.format(Locale.US, header);
                }
            }

            /**
             * {@inheritDoc}
             */
            @Override
            public void evaluationPerformed(final int iterationsCount, final int evaluationsCount, final Orbit[] orbits, final ParameterDriversList estimatedOrbitalParameters, final ParameterDriversList estimatedPropagatorParameters, final ParameterDriversList estimatedMeasurementsParameters, final EstimationsProvider evaluationsProvider, final LeastSquaresProblem.Evaluation lspEvaluation) {
                PVCoordinates currentPV = orbits[0].getPVCoordinates();
                final String format0 = "    %2d         %2d                                 %16.12f     %s       %s     %s     %s%n";
                final String format = "    %2d         %2d      %13.6f %12.9f %16.12f     %s       %s     %s     %s%n";
                final EvaluationCounter<Range> rangeCounter = new EvaluationCounter<Range>();
                final EvaluationCounter<RangeRate> rangeRateCounter = new EvaluationCounter<RangeRate>();
                final EvaluationCounter<AngularAzEl> angularCounter = new EvaluationCounter<AngularAzEl>();
                final EvaluationCounter<PV> pvCounter = new EvaluationCounter<PV>();
                for (final Map.Entry<ObservedMeasurement<?>, EstimatedMeasurement<?>> entry : estimator.getLastEstimations().entrySet()) {
                    if (entry.getKey() instanceof Range) {
                        @SuppressWarnings("unchecked") EstimatedMeasurement<Range> evaluation = (EstimatedMeasurement<Range>) entry.getValue();
                        rangeCounter.add(evaluation);
                    } else if (entry.getKey() instanceof RangeRate) {
                        @SuppressWarnings("unchecked") EstimatedMeasurement<RangeRate> evaluation = (EstimatedMeasurement<RangeRate>) entry.getValue();
                        rangeRateCounter.add(evaluation);
                    } else if (entry.getKey() instanceof AngularAzEl) {
                        @SuppressWarnings("unchecked") EstimatedMeasurement<AngularAzEl> evaluation = (EstimatedMeasurement<AngularAzEl>) entry.getValue();
                        angularCounter.add(evaluation);
                    } else if (entry.getKey() instanceof PV) {
                        @SuppressWarnings("unchecked") EstimatedMeasurement<PV> evaluation = (EstimatedMeasurement<PV>) entry.getValue();
                        pvCounter.add(evaluation);
                    }
                }
                if (evaluationsCount == 1) {
                    System.out.format(Locale.US, format0, iterationsCount, evaluationsCount, lspEvaluation.getRMS(), rangeCounter.format(8), rangeRateCounter.format(8), angularCounter.format(8), pvCounter.format(8));
                    if (logStream != null) {
                        logStream.format(Locale.US, format0, iterationsCount, evaluationsCount, lspEvaluation.getRMS(), rangeCounter.format(8), rangeRateCounter.format(8), angularCounter.format(8), pvCounter.format(8));
                    }
                } else {
                    System.out.format(Locale.US, format, iterationsCount, evaluationsCount, Vector3D.distance(previousPV.getPosition(), currentPV.getPosition()), Vector3D.distance(previousPV.getVelocity(), currentPV.getVelocity()), lspEvaluation.getRMS(), rangeCounter.format(8), rangeRateCounter.format(8), angularCounter.format(8), pvCounter.format(8));
                    if (logStream != null) {
                        logStream.format(Locale.US, format, iterationsCount, evaluationsCount, Vector3D.distance(previousPV.getPosition(), currentPV.getPosition()), Vector3D.distance(previousPV.getVelocity(), currentPV.getVelocity()), lspEvaluation.getRMS(), rangeCounter.format(8), rangeRateCounter.format(8), angularCounter.format(8), pvCounter.format(8));
                    }
                }
                previousPV = currentPV;
            }
        });
        Orbit estimated = estimator.estimate()[0].getInitialState().getOrbit();
        // compute some statistics
        for (final Map.Entry<ObservedMeasurement<?>, EstimatedMeasurement<?>> entry : estimator.getLastEstimations().entrySet()) {
            if (entry.getKey() instanceof Range) {
                @SuppressWarnings("unchecked") EstimatedMeasurement<Range> evaluation = (EstimatedMeasurement<Range>) entry.getValue();
                rangeLog.add(evaluation);
            } else if (entry.getKey() instanceof RangeRate) {
                @SuppressWarnings("unchecked") EstimatedMeasurement<RangeRate> evaluation = (EstimatedMeasurement<RangeRate>) entry.getValue();
                rangeRateLog.add(evaluation);
            } else if (entry.getKey() instanceof AngularAzEl) {
                @SuppressWarnings("unchecked") EstimatedMeasurement<AngularAzEl> evaluation = (EstimatedMeasurement<AngularAzEl>) entry.getValue();
                azimuthLog.add(evaluation);
                elevationLog.add(evaluation);
            } else if (entry.getKey() instanceof PV) {
                @SuppressWarnings("unchecked") EstimatedMeasurement<PV> evaluation = (EstimatedMeasurement<PV>) entry.getValue();
                positionLog.add(evaluation);
                velocityLog.add(evaluation);
            }
        }
        System.out.println("Estimated orbit: " + estimated);
        if (logStream != null) {
            logStream.println("Estimated orbit: " + estimated);
        }
        final ParameterDriversList orbitalParameters = estimator.getOrbitalParametersDrivers(true);
        final ParameterDriversList propagatorParameters = estimator.getPropagatorParametersDrivers(true);
        final ParameterDriversList measurementsParameters = estimator.getMeasurementsParametersDrivers(true);
        int length = 0;
        for (final ParameterDriver parameterDriver : orbitalParameters.getDrivers()) {
            length = FastMath.max(length, parameterDriver.getName().length());
        }
        for (final ParameterDriver parameterDriver : propagatorParameters.getDrivers()) {
            length = FastMath.max(length, parameterDriver.getName().length());
        }
        for (final ParameterDriver parameterDriver : measurementsParameters.getDrivers()) {
            length = FastMath.max(length, parameterDriver.getName().length());
        }
        displayParametersChanges(System.out, "Estimated orbital parameters changes: ", false, length, orbitalParameters);
        if (logStream != null) {
            displayParametersChanges(logStream, "Estimated orbital parameters changes: ", false, length, orbitalParameters);
        }
        displayParametersChanges(System.out, "Estimated propagator parameters changes: ", true, length, propagatorParameters);
        if (logStream != null) {
            displayParametersChanges(logStream, "Estimated propagator parameters changes: ", true, length, propagatorParameters);
        }
        displayParametersChanges(System.out, "Estimated measurements parameters changes: ", true, length, measurementsParameters);
        if (logStream != null) {
            displayParametersChanges(logStream, "Estimated measurements parameters changes: ", true, length, measurementsParameters);
        }
        System.out.println("Number of iterations: " + estimator.getIterationsCount());
        System.out.println("Number of evaluations: " + estimator.getEvaluationsCount());
        rangeLog.displaySummary(System.out);
        rangeRateLog.displaySummary(System.out);
        azimuthLog.displaySummary(System.out);
        elevationLog.displaySummary(System.out);
        positionLog.displaySummary(System.out);
        velocityLog.displaySummary(System.out);
        if (logStream != null) {
            logStream.println("Number of iterations: " + estimator.getIterationsCount());
            logStream.println("Number of evaluations: " + estimator.getEvaluationsCount());
            rangeLog.displaySummary(logStream);
            rangeRateLog.displaySummary(logStream);
            azimuthLog.displaySummary(logStream);
            elevationLog.displaySummary(logStream);
            positionLog.displaySummary(logStream);
            velocityLog.displaySummary(logStream);
        }
        rangeLog.displayResiduals();
        rangeRateLog.displayResiduals();
        azimuthLog.displayResiduals();
        elevationLog.displayResiduals();
        positionLog.displayResiduals();
        velocityLog.displayResiduals();
    } finally {
        if (logStream != null) {
            logStream.close();
        }
        rangeLog.close();
        rangeRateLog.close();
        azimuthLog.close();
        elevationLog.close();
        positionLog.close();
        velocityLog.close();
    }
}
Also used : OneAxisEllipsoid(org.orekit.bodies.OneAxisEllipsoid) PV(org.orekit.estimation.measurements.PV) ArrayList(java.util.ArrayList) PVCoordinates(org.orekit.utils.PVCoordinates) ObservedMeasurement(org.orekit.estimation.measurements.ObservedMeasurement) EstimatedMeasurement(org.orekit.estimation.measurements.EstimatedMeasurement) BatchLSObserver(org.orekit.estimation.leastsquares.BatchLSObserver) EquinoctialOrbit(org.orekit.orbits.EquinoctialOrbit) CartesianOrbit(org.orekit.orbits.CartesianOrbit) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) Orbit(org.orekit.orbits.Orbit) CircularOrbit(org.orekit.orbits.CircularOrbit) Range(org.orekit.estimation.measurements.Range) FileInputStream(java.io.FileInputStream) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) RangeRate(org.orekit.estimation.measurements.RangeRate) File(java.io.File) AngularAzEl(org.orekit.estimation.measurements.AngularAzEl) Map(java.util.Map) HashMap(java.util.HashMap) BatchLSEstimator(org.orekit.estimation.leastsquares.BatchLSEstimator) ParameterDriversList(org.orekit.utils.ParameterDriversList) LeastSquaresProblem(org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem) NormalizedSphericalHarmonicsProvider(org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider) EstimationsProvider(org.orekit.estimation.measurements.EstimationsProvider) PrintStream(java.io.PrintStream) KeyValueFileParser(fr.cs.examples.KeyValueFileParser) IERSConventions(org.orekit.utils.IERSConventions) ParameterDriver(org.orekit.utils.ParameterDriver) GeodeticPoint(org.orekit.bodies.GeodeticPoint)

Example 5 with EstimationsProvider

use of org.orekit.estimation.measurements.EstimationsProvider in project Orekit by CS-SI.

the class BatchLSEstimatorTest method testKeplerRange.

/**
 * Perfect range measurements with a biased start
 * @throws OrekitException
 */
@Test
public void testKeplerRange() throws OrekitException {
    Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 1.0);
    // create perfect range measurements
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.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() {

        int lastIter = 0;

        int lastEval = 0;

        /**
         * {@inheritDoc}
         */
        @Override
        public void evaluationPerformed(int iterationsCount, int evaluationscount, Orbit[] orbits, ParameterDriversList estimatedOrbitalParameters, ParameterDriversList estimatedPropagatorParameters, ParameterDriversList estimatedMeasurementsParameters, EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) throws OrekitException {
            if (iterationsCount == lastIter) {
                Assert.assertEquals(lastEval + 1, evaluationscount);
            } else {
                Assert.assertEquals(lastIter + 1, iterationsCount);
            }
            lastIter = iterationsCount;
            lastEval = evaluationscount;
            Assert.assertEquals(measurements.size(), evaluationsProvider.getNumber());
            try {
                evaluationsProvider.getEstimatedMeasurement(-1);
                Assert.fail("an exception should have been thrown");
            } catch (OrekitException oe) {
                Assert.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
            }
            try {
                evaluationsProvider.getEstimatedMeasurement(measurements.size());
                Assert.fail("an exception should have been thrown");
            } catch (OrekitException oe) {
                Assert.assertEquals(LocalizedCoreFormats.OUT_OF_RANGE_SIMPLE, oe.getSpecifier());
            }
            AbsoluteDate previous = AbsoluteDate.PAST_INFINITY;
            for (int i = 0; i < evaluationsProvider.getNumber(); ++i) {
                AbsoluteDate current = evaluationsProvider.getEstimatedMeasurement(i).getDate();
                Assert.assertTrue(current.compareTo(previous) >= 0);
                previous = current;
            }
        }
    });
    ParameterDriver aDriver = estimator.getOrbitalParametersDrivers(true).getDrivers().get(0);
    Assert.assertEquals("a", aDriver.getName());
    aDriver.setValue(aDriver.getValue() + 1.2);
    aDriver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
    EstimationTestUtils.checkFit(context, estimator, 2, 3, 0.0, 1.1e-6, 0.0, 2.8e-6, 0.0, 4.0e-7, 0.0, 2.2e-10);
    // got a default one
    for (final ParameterDriver driver : estimator.getOrbitalParametersDrivers(true).getDrivers()) {
        if ("a".equals(driver.getName())) {
            // user-specified reference date
            Assert.assertEquals(0, driver.getReferenceDate().durationFrom(AbsoluteDate.GALILEO_EPOCH), 1.0e-15);
        } else {
            // default reference date
            Assert.assertEquals(0, driver.getReferenceDate().durationFrom(propagatorBuilder.getInitialOrbitDate()), 1.0e-15);
        }
    }
}
Also used : Context(org.orekit.estimation.Context) Evaluation(org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem.Evaluation) Orbit(org.orekit.orbits.Orbit) CartesianOrbit(org.orekit.orbits.CartesianOrbit) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) ParameterDriver(org.orekit.utils.ParameterDriver) AbsoluteDate(org.orekit.time.AbsoluteDate) LevenbergMarquardtOptimizer(org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer) ParameterDriversList(org.orekit.utils.ParameterDriversList) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) BoundedPropagator(org.orekit.propagation.BoundedPropagator) Propagator(org.orekit.propagation.Propagator) OrekitException(org.orekit.errors.OrekitException) RangeMeasurementCreator(org.orekit.estimation.measurements.RangeMeasurementCreator) InterSatellitesRangeMeasurementCreator(org.orekit.estimation.measurements.InterSatellitesRangeMeasurementCreator) ObservedMeasurement(org.orekit.estimation.measurements.ObservedMeasurement) EstimationsProvider(org.orekit.estimation.measurements.EstimationsProvider) Test(org.junit.Test)

Aggregations

EstimationsProvider (org.orekit.estimation.measurements.EstimationsProvider)7 ObservedMeasurement (org.orekit.estimation.measurements.ObservedMeasurement)7 CartesianOrbit (org.orekit.orbits.CartesianOrbit)7 KeplerianOrbit (org.orekit.orbits.KeplerianOrbit)7 Orbit (org.orekit.orbits.Orbit)7 NumericalPropagatorBuilder (org.orekit.propagation.conversion.NumericalPropagatorBuilder)7 ParameterDriversList (org.orekit.utils.ParameterDriversList)7 Evaluation (org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem.Evaluation)5 LevenbergMarquardtOptimizer (org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer)5 Test (org.junit.Test)5 Context (org.orekit.estimation.Context)5 InterSatellitesRangeMeasurementCreator (org.orekit.estimation.measurements.InterSatellitesRangeMeasurementCreator)5 RangeMeasurementCreator (org.orekit.estimation.measurements.RangeMeasurementCreator)5 BoundedPropagator (org.orekit.propagation.BoundedPropagator)5 Propagator (org.orekit.propagation.Propagator)5 ParameterDriver (org.orekit.utils.ParameterDriver)5 OrekitException (org.orekit.errors.OrekitException)4 Vector3D (org.hipparchus.geometry.euclidean.threed.Vector3D)3 Range (org.orekit.estimation.measurements.Range)3 AbsoluteDate (org.orekit.time.AbsoluteDate)3