Search in sources :

Example 26 with ObservedMeasurement

use of org.orekit.estimation.measurements.ObservedMeasurement 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 27 with ObservedMeasurement

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

the class KalmanEstimatorTest method testCartesianRangeRate.

/**
 * Perfect range rate measurements with a perfect start
 * Cartesian formalism
 * @throws OrekitException
 */
@Test
public void testCartesianRangeRate() throws OrekitException {
    // Create context
    Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    // Create initial orbit and propagator builder
    final OrbitType orbitType = OrbitType.CARTESIAN;
    final PositionAngle positionAngle = PositionAngle.TRUE;
    final boolean perfectStart = true;
    final double minStep = 1.e-6;
    final double maxStep = 60.;
    final double dP = 1.;
    final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(orbitType, positionAngle, perfectStart, minStep, maxStep, dP);
    // Create perfect range measurements
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new RangeRateMeasurementCreator(context, false), 1.0, 3.0, 300.0);
    // Reference propagator for estimation performances
    final NumericalPropagator referencePropagator = propagatorBuilder.buildPropagator(propagatorBuilder.getSelectedNormalizedParameters());
    // Reference position/velocity at last measurement date
    final Orbit refOrbit = referencePropagator.propagate(measurements.get(measurements.size() - 1).getDate()).getOrbit();
    // Cartesian covariance matrix initialization
    // 100m on position / 1e-2m/s on velocity
    final RealMatrix cartesianP = MatrixUtils.createRealDiagonalMatrix(new double[] { 1e-4, 1e-4, 1e-4, 1e-10, 1e-10, 1e-10 });
    // Jacobian of the orbital parameters w/r to Cartesian
    final Orbit initialOrbit = orbitType.convertType(context.initialOrbit);
    final double[][] dYdC = new double[6][6];
    initialOrbit.getJacobianWrtCartesian(PositionAngle.TRUE, dYdC);
    final RealMatrix Jac = MatrixUtils.createRealMatrix(dYdC);
    // Initial covariance matrix
    final RealMatrix initialP = Jac.multiply(cartesianP.multiply(Jac.transpose()));
    // Process noise matrix
    final RealMatrix cartesianQ = MatrixUtils.createRealDiagonalMatrix(new double[] { 1.e-6, 1.e-6, 1.e-6, 1.e-12, 1.e-12, 1.e-12 });
    final RealMatrix Q = Jac.multiply(cartesianQ.multiply(Jac.transpose()));
    // Build the Kalman filter
    final KalmanEstimatorBuilder kalmanBuilder = new KalmanEstimatorBuilder();
    kalmanBuilder.builder(propagatorBuilder);
    kalmanBuilder.estimatedMeasurementsParameters(new ParameterDriversList());
    kalmanBuilder.initialCovarianceMatrix(initialP);
    kalmanBuilder.processNoiseMatrixProvider(new ConstantProcessNoise(Q));
    final KalmanEstimator kalman = kalmanBuilder.build();
    // Filter the measurements and check the results
    final double expectedDeltaPos = 0.;
    final double posEps = 9.50e-4;
    final double expectedDeltaVel = 0.;
    final double velEps = 3.49e-7;
    final double[] expectedSigmasPos = { 0.324398, 1.347031, 1.743310 };
    final double sigmaPosEps = 1e-6;
    final double[] expectedSigmasVel = { 2.856883e-4, 5.765844e-4, 5.056186e-4 };
    final double sigmaVelEps = 1e-10;
    EstimationTestUtils.checkKalmanFit(context, kalman, measurements, refOrbit, positionAngle, expectedDeltaPos, posEps, expectedDeltaVel, velEps, expectedSigmasPos, sigmaPosEps, expectedSigmasVel, sigmaVelEps);
}
Also used : Context(org.orekit.estimation.Context) Orbit(org.orekit.orbits.Orbit) PositionAngle(org.orekit.orbits.PositionAngle) RealMatrix(org.hipparchus.linear.RealMatrix) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) ParameterDriversList(org.orekit.utils.ParameterDriversList) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) Propagator(org.orekit.propagation.Propagator) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) OrbitType(org.orekit.orbits.OrbitType) ObservedMeasurement(org.orekit.estimation.measurements.ObservedMeasurement) RangeRateMeasurementCreator(org.orekit.estimation.measurements.RangeRateMeasurementCreator) Test(org.junit.Test)

Example 28 with ObservedMeasurement

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

the class KalmanEstimatorTest method testCircularAzimuthElevation.

/**
 * Perfect azimuth/elevation measurements with a perfect start
 * Circular formalism
 * @throws OrekitException
 */
@Test
public void testCircularAzimuthElevation() throws OrekitException {
    // Create context
    Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    // Create initial orbit and propagator builder
    final OrbitType orbitType = OrbitType.CIRCULAR;
    final PositionAngle positionAngle = PositionAngle.TRUE;
    final boolean perfectStart = true;
    final double minStep = 1.e-6;
    final double maxStep = 60.;
    final double dP = 1.;
    final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(orbitType, positionAngle, perfectStart, minStep, maxStep, dP);
    // Create perfect range measurements
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new AngularAzElMeasurementCreator(context), 1.0, 4.0, 60.0);
    // Reference propagator for estimation performances
    final NumericalPropagator referencePropagator = propagatorBuilder.buildPropagator(propagatorBuilder.getSelectedNormalizedParameters());
    // Reference position/velocity at last measurement date
    final Orbit refOrbit = referencePropagator.propagate(measurements.get(measurements.size() - 1).getDate()).getOrbit();
    // Cartesian covariance matrix initialization
    final RealMatrix cartesianP = MatrixUtils.createRealDiagonalMatrix(new double[] { 1e-4, 1e-4, 1e-4, 1e-10, 1e-10, 1e-10 });
    // Jacobian of the orbital parameters w/r to Cartesian
    final Orbit initialOrbit = orbitType.convertType(context.initialOrbit);
    final double[][] dYdC = new double[6][6];
    initialOrbit.getJacobianWrtCartesian(PositionAngle.TRUE, dYdC);
    final RealMatrix Jac = MatrixUtils.createRealMatrix(dYdC);
    // Initial covariance matrix
    final RealMatrix initialP = Jac.multiply(cartesianP.multiply(Jac.transpose()));
    // Process noise matrix
    final RealMatrix cartesianQ = MatrixUtils.createRealDiagonalMatrix(new double[] { 1.e-6, 1.e-6, 1.e-6, 1.e-12, 1.e-12, 1.e-12 });
    final RealMatrix Q = Jac.multiply(cartesianQ.multiply(Jac.transpose()));
    // Build the Kalman filter
    final KalmanEstimatorBuilder kalmanBuilder = new KalmanEstimatorBuilder();
    kalmanBuilder.builder(propagatorBuilder);
    kalmanBuilder.estimatedMeasurementsParameters(new ParameterDriversList());
    kalmanBuilder.initialCovarianceMatrix(initialP);
    kalmanBuilder.processNoiseMatrixProvider(new ConstantProcessNoise(Q));
    final KalmanEstimator kalman = kalmanBuilder.build();
    // Filter the measurements and check the results
    final double expectedDeltaPos = 0.;
    final double posEps = 4.78e-7;
    final double expectedDeltaVel = 0.;
    final double velEps = 1.54e-10;
    final double[] expectedSigmasPos = { 0.356902, 1.297507, 1.798551 };
    final double sigmaPosEps = 1e-6;
    final double[] expectedSigmasVel = { 2.468745e-4, 5.810027e-4, 3.887394e-4 };
    final double sigmaVelEps = 1e-10;
    EstimationTestUtils.checkKalmanFit(context, kalman, measurements, refOrbit, positionAngle, expectedDeltaPos, posEps, expectedDeltaVel, velEps, expectedSigmasPos, sigmaPosEps, expectedSigmasVel, sigmaVelEps);
}
Also used : Context(org.orekit.estimation.Context) Orbit(org.orekit.orbits.Orbit) PositionAngle(org.orekit.orbits.PositionAngle) RealMatrix(org.hipparchus.linear.RealMatrix) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) ParameterDriversList(org.orekit.utils.ParameterDriversList) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) Propagator(org.orekit.propagation.Propagator) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) AngularAzElMeasurementCreator(org.orekit.estimation.measurements.AngularAzElMeasurementCreator) OrbitType(org.orekit.orbits.OrbitType) ObservedMeasurement(org.orekit.estimation.measurements.ObservedMeasurement) Test(org.junit.Test)

Example 29 with ObservedMeasurement

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

the class KalmanEstimatorTest method testKeplerianRange.

/**
 * Perfect range measurements with a biased start
 * Keplerian formalism
 * @throws OrekitException
 */
@Test
public void testKeplerianRange() throws OrekitException {
    // Create context
    Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    // Create initial orbit and propagator builder
    final OrbitType orbitType = OrbitType.KEPLERIAN;
    final PositionAngle positionAngle = PositionAngle.TRUE;
    final boolean perfectStart = true;
    final double minStep = 1.e-6;
    final double maxStep = 60.;
    final double dP = 1.;
    final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(orbitType, positionAngle, perfectStart, minStep, maxStep, dP);
    // Create perfect range measurements
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new RangeMeasurementCreator(context), 1.0, 4.0, 60.0);
    // Reference propagator for estimation performances
    final NumericalPropagator referencePropagator = propagatorBuilder.buildPropagator(propagatorBuilder.getSelectedNormalizedParameters());
    // Reference position/velocity at last measurement date
    final Orbit refOrbit = referencePropagator.propagate(measurements.get(measurements.size() - 1).getDate()).getOrbit();
    // Change semi-major axis of 1.2m as in the batch test
    ParameterDriver aDriver = propagatorBuilder.getOrbitalParametersDrivers().getDrivers().get(0);
    aDriver.setValue(aDriver.getValue() + 1.2);
    aDriver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
    // Cartesian covariance matrix initialization
    // 100m on position / 1e-2m/s on velocity
    final RealMatrix cartesianP = MatrixUtils.createRealDiagonalMatrix(new double[] { 100., 100., 100., 1e-2, 1e-2, 1e-2 });
    // Jacobian of the orbital parameters w/r to Cartesian
    final Orbit initialOrbit = orbitType.convertType(context.initialOrbit);
    final double[][] dYdC = new double[6][6];
    initialOrbit.getJacobianWrtCartesian(PositionAngle.TRUE, dYdC);
    final RealMatrix Jac = MatrixUtils.createRealMatrix(dYdC);
    // Keplerian initial covariance matrix
    final RealMatrix initialP = Jac.multiply(cartesianP.multiply(Jac.transpose()));
    // Process noise matrix is set to 0 here
    RealMatrix Q = MatrixUtils.createRealMatrix(6, 6);
    // Build the Kalman filter
    final KalmanEstimatorBuilder kalmanBuilder = new KalmanEstimatorBuilder();
    kalmanBuilder.builder(propagatorBuilder);
    kalmanBuilder.estimatedMeasurementsParameters(new ParameterDriversList());
    kalmanBuilder.initialCovarianceMatrix(initialP);
    kalmanBuilder.processNoiseMatrixProvider(new ConstantProcessNoise(Q));
    final KalmanEstimator kalman = kalmanBuilder.build();
    // Filter the measurements and check the results
    final double expectedDeltaPos = 0.;
    final double posEps = 1.77e-4;
    final double expectedDeltaVel = 0.;
    final double velEps = 7.93e-8;
    final double[] expectedSigmasPos = { 0.742488, 0.281910, 0.563217 };
    final double sigmaPosEps = 1e-6;
    final double[] expectedSigmasVel = { 2.206622e-4, 1.306669e-4, 1.293996e-4 };
    final double sigmaVelEps = 1e-10;
    EstimationTestUtils.checkKalmanFit(context, kalman, measurements, refOrbit, positionAngle, expectedDeltaPos, posEps, expectedDeltaVel, velEps, expectedSigmasPos, sigmaPosEps, expectedSigmasVel, sigmaVelEps);
}
Also used : Context(org.orekit.estimation.Context) Orbit(org.orekit.orbits.Orbit) PositionAngle(org.orekit.orbits.PositionAngle) ParameterDriver(org.orekit.utils.ParameterDriver) RealMatrix(org.hipparchus.linear.RealMatrix) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) ParameterDriversList(org.orekit.utils.ParameterDriversList) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) Propagator(org.orekit.propagation.Propagator) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) OrbitType(org.orekit.orbits.OrbitType) RangeMeasurementCreator(org.orekit.estimation.measurements.RangeMeasurementCreator) ObservedMeasurement(org.orekit.estimation.measurements.ObservedMeasurement) Test(org.junit.Test)

Example 30 with ObservedMeasurement

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

the class KalmanEstimatorTest method testWrappedException.

/**
 * Test of a wrapped exception in a Kalman observer
 * @throws OrekitException
 */
@Test
public void testWrappedException() throws OrekitException {
    // Create context
    Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    // Create initial orbit and propagator builder
    final OrbitType orbitType = OrbitType.KEPLERIAN;
    final PositionAngle positionAngle = PositionAngle.TRUE;
    final boolean perfectStart = true;
    final double minStep = 1.e-6;
    final double maxStep = 60.;
    final double dP = 1.;
    final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(orbitType, positionAngle, perfectStart, minStep, maxStep, dP);
    // 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);
    // Build the Kalman filter
    final KalmanEstimatorBuilder kalmanBuilder = new KalmanEstimatorBuilder();
    kalmanBuilder.builder(propagatorBuilder);
    kalmanBuilder.estimatedMeasurementsParameters(new ParameterDriversList());
    kalmanBuilder.initialCovarianceMatrix(MatrixUtils.createRealMatrix(6, 6));
    kalmanBuilder.processNoiseMatrixProvider(new ConstantProcessNoise(MatrixUtils.createRealMatrix(6, 6)));
    final KalmanEstimator kalman = kalmanBuilder.build();
    kalman.setObserver(estimation -> {
        throw new DummyException();
    });
    try {
        // Filter the measurements and expect an exception to occur
        EstimationTestUtils.checkKalmanFit(context, kalman, measurements, context.initialOrbit, positionAngle, 0., 0., 0., 0., new double[3], 0., new double[3], 0.);
    } catch (DummyException de) {
    // expected
    }
}
Also used : Context(org.orekit.estimation.Context) PositionAngle(org.orekit.orbits.PositionAngle) ParameterDriversList(org.orekit.utils.ParameterDriversList) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) Propagator(org.orekit.propagation.Propagator) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) OrbitType(org.orekit.orbits.OrbitType) RangeMeasurementCreator(org.orekit.estimation.measurements.RangeMeasurementCreator) ObservedMeasurement(org.orekit.estimation.measurements.ObservedMeasurement) Test(org.junit.Test)

Aggregations

ObservedMeasurement (org.orekit.estimation.measurements.ObservedMeasurement)49 NumericalPropagatorBuilder (org.orekit.propagation.conversion.NumericalPropagatorBuilder)44 Test (org.junit.Test)41 Context (org.orekit.estimation.Context)40 Propagator (org.orekit.propagation.Propagator)40 Orbit (org.orekit.orbits.Orbit)20 AbsoluteDate (org.orekit.time.AbsoluteDate)19 ParameterDriversList (org.orekit.utils.ParameterDriversList)19 RangeMeasurementCreator (org.orekit.estimation.measurements.RangeMeasurementCreator)17 Vector3D (org.hipparchus.geometry.euclidean.threed.Vector3D)16 ParameterDriver (org.orekit.utils.ParameterDriver)15 Range (org.orekit.estimation.measurements.Range)14 KeplerianOrbit (org.orekit.orbits.KeplerianOrbit)13 RealMatrix (org.hipparchus.linear.RealMatrix)12 ArrayList (java.util.ArrayList)11 LevenbergMarquardtOptimizer (org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer)11 CartesianOrbit (org.orekit.orbits.CartesianOrbit)11 BoundedPropagator (org.orekit.propagation.BoundedPropagator)11 SpacecraftState (org.orekit.propagation.SpacecraftState)11 GroundStation (org.orekit.estimation.measurements.GroundStation)10