Examples with ObservedMeasurement org.orekit.estimation.measurements.ObservedMeasurement used on opensource projects

Example 16 with ObservedMeasurement

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

the class KalmanEstimatorTest method testKeplerianRangeAndRangeRate.

/**
 * Perfect range and range rate measurements with a perfect start
 * @throws OrekitException
 */
@Test
public void testKeplerianRangeAndRangeRate() 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 & range rate measurements
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurementsRange = EstimationTestUtils.createMeasurements(propagator, new RangeMeasurementCreator(context), 1.0, 3.0, 300.0);
    final List<ObservedMeasurement<?>> measurementsRangeRate = EstimationTestUtils.createMeasurements(propagator, new RangeRateMeasurementCreator(context, false), 1.0, 3.0, 300.0);
    // Concatenate measurements
    final List<ObservedMeasurement<?>> measurements = new ArrayList<ObservedMeasurement<?>>();
    measurements.addAll(measurementsRange);
    measurements.addAll(measurementsRangeRate);
    // 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-2, 1e-2, 1e-2, 1e-8, 1e-8, 1e-8 });
    // 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
    final RealMatrix cartesianQ = MatrixUtils.createRealDiagonalMatrix(new double[] { 1.e-4, 1.e-4, 1.e-4, 1.e-10, 1.e-10, 1.e-10 });
    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 = 5.96e-3;
    final double expectedDeltaVel = 0.;
    final double velEps = 2.06e-6;
    final double[] expectedSigmasPos = { 0.341538, 8.175281, 4.634384 };
    final double sigmaPosEps = 1e-6;
    final double[] expectedSigmasVel = { 1.167838e-3, 1.036437e-3, 2.834385e-3 };
    final double sigmaVelEps = 1e-9;
    EstimationTestUtils.checkKalmanFit(context, kalman, measurements, refOrbit, positionAngle, expectedDeltaPos, posEps, expectedDeltaVel, velEps, expectedSigmasPos, sigmaPosEps, expectedSigmasVel, sigmaVelEps);
}

Example 17 with ObservedMeasurement

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

the class KalmanEstimatorTest method testKeplerianPV.

/**
 * Perfect PV measurements with a perfect start
 * Keplerian formalism
 * @throws OrekitException
 */
@Test
public void testKeplerianPV() 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 PV measurements
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new PVMeasurementCreator(), 0.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();
    // Covariance matrix initialization
    final RealMatrix initialP = MatrixUtils.createRealDiagonalMatrix(new double[] { 1e-2, 1e-2, 1e-2, 1e-5, 1e-5, 1e-5 });
    // Process noise matrix
    RealMatrix Q = MatrixUtils.createRealDiagonalMatrix(new double[] { 1.e-8, 1.e-8, 1.e-8, 1.e-8, 1.e-8, 1.e-8 });
    // 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 = 5.80e-8;
    final double expectedDeltaVel = 0.;
    final double velEps = 2.28e-11;
    final double[] expectedsigmasPos = { 0.998872, 0.933655, 0.997516 };
    final double sigmaPosEps = 1e-6;
    final double[] expectedSigmasVel = { 9.478853e-4, 9.910788e-4, 5.0438709e-4 };
    final double sigmaVelEps = 1e-10;
    EstimationTestUtils.checkKalmanFit(context, kalman, measurements, refOrbit, positionAngle, expectedDeltaPos, posEps, expectedDeltaVel, velEps, expectedsigmasPos, sigmaPosEps, expectedSigmasVel, sigmaVelEps);
}

Example 18 with ObservedMeasurement

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

the class KalmanOrbitDeterminationTest method readMeasurements.

/**
 * Read a measurements file.
 * @param file measurements file
 * @param stations name to stations data map
 * @param pvData PV measurements data
 * @param satRangeBias range bias due to transponder delay
 * @param weights base weights for measurements
 * @param rangeOutliersManager manager for range measurements outliers (null if none configured)
 * @param rangeRateOutliersManager manager for range-rate measurements outliers (null if none configured)
 * @param azElOutliersManager manager for azimuth-elevation measurements outliers (null if none configured)
 * @param pvOutliersManager manager for PV measurements outliers (null if none configured)
 * @return measurements list
 */
private List<ObservedMeasurement<?>> readMeasurements(final File file, final Map<String, StationData> stations, final PVData pvData, final Bias<Range> satRangeBias, final Weights weights, final OutlierFilter<Range> rangeOutliersManager, final OutlierFilter<RangeRate> rangeRateOutliersManager, final OutlierFilter<AngularAzEl> azElOutliersManager, final OutlierFilter<PV> pvOutliersManager) throws UnsupportedEncodingException, IOException, OrekitException {
    final List<ObservedMeasurement<?>> measurements = new ArrayList<ObservedMeasurement<?>>();
    BufferedReader br = null;
    try {
        br = new BufferedReader(new InputStreamReader(new FileInputStream(file), "UTF-8"));
        int lineNumber = 0;
        for (String line = br.readLine(); line != null; line = br.readLine()) {
            ++lineNumber;
            line = line.trim();
            if (line.length() > 0 && !line.startsWith("#")) {
                String[] fields = line.split("\\s+");
                if (fields.length < 2) {
                    throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE, lineNumber, file.getName(), line);
                }
                switch(fields[1]) {
                    case "RANGE":
                        final Range range = new RangeParser().parseFields(fields, stations, pvData, satRangeBias, weights, line, lineNumber, file.getName());
                        if (rangeOutliersManager != null) {
                            range.addModifier(rangeOutliersManager);
                        }
                        addIfNonZeroWeight(range, measurements);
                        break;
                    case "RANGE_RATE":
                        final RangeRate rangeRate = new RangeRateParser().parseFields(fields, stations, pvData, satRangeBias, weights, line, lineNumber, file.getName());
                        if (rangeOutliersManager != null) {
                            rangeRate.addModifier(rangeRateOutliersManager);
                        }
                        addIfNonZeroWeight(rangeRate, measurements);
                        break;
                    case "AZ_EL":
                        final AngularAzEl angular = new AzElParser().parseFields(fields, stations, pvData, satRangeBias, weights, line, lineNumber, file.getName());
                        if (azElOutliersManager != null) {
                            angular.addModifier(azElOutliersManager);
                        }
                        addIfNonZeroWeight(angular, measurements);
                        break;
                    case "PV":
                        final PV pv = new PVParser().parseFields(fields, stations, pvData, satRangeBias, weights, line, lineNumber, file.getName());
                        if (pvOutliersManager != null) {
                            pv.addModifier(pvOutliersManager);
                        }
                        addIfNonZeroWeight(pv, measurements);
                        break;
                    default:
                        throw new OrekitException(LocalizedCoreFormats.SIMPLE_MESSAGE, "unknown measurement type " + fields[1] + " at line " + lineNumber + " in file " + file.getName());
                }
            }
        }
    } finally {
        if (br != null) {
            br.close();
        }
    }
    if (measurements.isEmpty()) {
        throw new OrekitException(LocalizedCoreFormats.SIMPLE_MESSAGE, "not measurements read from file " + file.getAbsolutePath());
    }
    return measurements;
}

Example 19 with ObservedMeasurement

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

the class EstimationTestUtils method checkFit.

/**
 * 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
 * @throws OrekitException
 */
public static void checkFit(final Context 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) throws OrekitException {
    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));
        }
    }
    Assert.assertEquals(expectedRMS, FastMath.sqrt(sum / k), rmsEps);
    Assert.assertEquals(expectedMax, max, maxEps);
    Assert.assertEquals(expectedDeltaPos, Vector3D.distance(context.initialOrbit.getPVCoordinates().getPosition(), estimatedPosition), posEps);
    Assert.assertEquals(expectedDeltaVel, Vector3D.distance(context.initialOrbit.getPVCoordinates().getVelocity(), estimatedVelocity), velEps);
}

Example 20 with ObservedMeasurement

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

the class IodGibbsTest method testGibbs1.

@Test
public void testGibbs1() throws OrekitException {
    final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    final double mu = context.initialOrbit.getMu();
    final Frame frame = context.initialOrbit.getFrame();
    final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 0.001);
    // create perfect range measurements
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new PVMeasurementCreator(), 0.0, 1.0, 60.0);
    final Vector3D position1 = new Vector3D(measurements.get(0).getObservedValue()[0], measurements.get(0).getObservedValue()[1], measurements.get(0).getObservedValue()[2]);
    final PV pv1 = new PV(measurements.get(0).getDate(), position1, Vector3D.ZERO, 0., 0., 1.);
    final Vector3D position2 = new Vector3D(measurements.get(1).getObservedValue()[0], measurements.get(1).getObservedValue()[1], measurements.get(1).getObservedValue()[2]);
    final PV pv2 = new PV(measurements.get(1).getDate(), position2, Vector3D.ZERO, 0., 0., 1.);
    final Vector3D position3 = new Vector3D(measurements.get(2).getObservedValue()[0], measurements.get(2).getObservedValue()[1], measurements.get(2).getObservedValue()[2]);
    final PV pv3 = new PV(measurements.get(2).getDate(), position3, Vector3D.ZERO, 0., 0., 1.);
    // instantiate the IOD method
    final IodGibbs gibbs = new IodGibbs(mu);
    final KeplerianOrbit orbit = gibbs.estimate(frame, pv1, pv2, pv3);
    Assert.assertEquals(context.initialOrbit.getA(), orbit.getA(), 1.0e-9 * context.initialOrbit.getA());
    Assert.assertEquals(context.initialOrbit.getE(), orbit.getE(), 1.0e-9 * context.initialOrbit.getE());
    Assert.assertEquals(context.initialOrbit.getI(), orbit.getI(), 1.0e-9 * context.initialOrbit.getI());
}