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

Example 1 with PV

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

the class HarmonicParametricAccelerationTest method testCoefficientsDetermination.

@Test
public void testCoefficientsDetermination() throws OrekitException {
    final double mass = 2500;
    final Orbit orbit = new CircularOrbit(7500000.0, 1.0e-4, 1.0e-3, 1.7, 0.3, 0.5, PositionAngle.TRUE, FramesFactory.getEME2000(), new AbsoluteDate(new DateComponents(2004, 2, 3), TimeComponents.H00, TimeScalesFactory.getUTC()), Constants.EIGEN5C_EARTH_MU);
    final double period = orbit.getKeplerianPeriod();
    AttitudeProvider maneuverLaw = new LofOffset(orbit.getFrame(), LOFType.VNC);
    SpacecraftState initialState = new SpacecraftState(orbit, maneuverLaw.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), mass);
    double dP = 10.0;
    double minStep = 0.001;
    double maxStep = 100;
    double[][] tolerance = NumericalPropagator.tolerances(dP, orbit, orbit.getType());
    // generate PV measurements corresponding to a tangential maneuver
    AdaptiveStepsizeIntegrator integrator0 = new DormandPrince853Integrator(minStep, maxStep, tolerance[0], tolerance[1]);
    integrator0.setInitialStepSize(60);
    final NumericalPropagator propagator0 = new NumericalPropagator(integrator0);
    propagator0.setInitialState(initialState);
    propagator0.setAttitudeProvider(maneuverLaw);
    ForceModel hpaRefX1 = new HarmonicParametricAcceleration(Vector3D.PLUS_I, true, "refX1", null, period, 1);
    ForceModel hpaRefY1 = new HarmonicParametricAcceleration(Vector3D.PLUS_J, true, "refY1", null, period, 1);
    ForceModel hpaRefZ2 = new HarmonicParametricAcceleration(Vector3D.PLUS_K, true, "refZ2", null, period, 2);
    hpaRefX1.getParametersDrivers()[0].setValue(2.4e-2);
    hpaRefX1.getParametersDrivers()[1].setValue(3.1);
    hpaRefY1.getParametersDrivers()[0].setValue(4.0e-2);
    hpaRefY1.getParametersDrivers()[1].setValue(0.3);
    hpaRefZ2.getParametersDrivers()[0].setValue(1.0e-2);
    hpaRefZ2.getParametersDrivers()[1].setValue(1.8);
    propagator0.addForceModel(hpaRefX1);
    propagator0.addForceModel(hpaRefY1);
    propagator0.addForceModel(hpaRefZ2);
    final List<ObservedMeasurement<?>> measurements = new ArrayList<>();
    propagator0.setMasterMode(10.0, (state, isLast) -> measurements.add(new PV(state.getDate(), state.getPVCoordinates().getPosition(), state.getPVCoordinates().getVelocity(), 1.0e-3, 1.0e-6, 1.0)));
    propagator0.propagate(orbit.getDate().shiftedBy(900));
    // set up an estimator to retrieve the maneuver as several harmonic accelerations in inertial frame
    final NumericalPropagatorBuilder propagatorBuilder = new NumericalPropagatorBuilder(orbit, new DormandPrince853IntegratorBuilder(minStep, maxStep, dP), PositionAngle.TRUE, dP);
    propagatorBuilder.addForceModel(new HarmonicParametricAcceleration(Vector3D.PLUS_I, true, "X1", null, period, 1));
    propagatorBuilder.addForceModel(new HarmonicParametricAcceleration(Vector3D.PLUS_J, true, "Y1", null, period, 1));
    propagatorBuilder.addForceModel(new HarmonicParametricAcceleration(Vector3D.PLUS_K, true, "Z2", null, period, 2));
    final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(), propagatorBuilder);
    estimator.setParametersConvergenceThreshold(1.0e-2);
    estimator.setMaxIterations(20);
    estimator.setMaxEvaluations(100);
    for (final ObservedMeasurement<?> measurement : measurements) {
        estimator.addMeasurement(measurement);
    }
    // we will estimate only the force model parameters, not the orbit
    for (final ParameterDriver d : estimator.getOrbitalParametersDrivers(false).getDrivers()) {
        d.setSelected(false);
    }
    setParameter(estimator, "X1 γ", 1.0e-2);
    setParameter(estimator, "X1 φ", 4.0);
    setParameter(estimator, "Y1 γ", 1.0e-2);
    setParameter(estimator, "Y1 φ", 0.0);
    setParameter(estimator, "Z2 γ", 1.0e-2);
    setParameter(estimator, "Z2 φ", 1.0);
    estimator.estimate();
    Assert.assertTrue(estimator.getIterationsCount() < 15);
    Assert.assertTrue(estimator.getEvaluationsCount() < 15);
    Assert.assertEquals(0.0, estimator.getOptimum().getRMS(), 1.0e-5);
    Assert.assertEquals(hpaRefX1.getParametersDrivers()[0].getValue(), getParameter(estimator, "X1 γ"), 1.e-12);
    Assert.assertEquals(hpaRefX1.getParametersDrivers()[1].getValue(), getParameter(estimator, "X1 φ"), 1.e-12);
    Assert.assertEquals(hpaRefY1.getParametersDrivers()[0].getValue(), getParameter(estimator, "Y1 γ"), 1.e-12);
    Assert.assertEquals(hpaRefY1.getParametersDrivers()[1].getValue(), getParameter(estimator, "Y1 φ"), 1.e-12);
    Assert.assertEquals(hpaRefZ2.getParametersDrivers()[0].getValue(), getParameter(estimator, "Z2 γ"), 1.e-12);
    Assert.assertEquals(hpaRefZ2.getParametersDrivers()[1].getValue(), getParameter(estimator, "Z2 φ"), 1.e-12);
}

Example 2 with PV

use of org.orekit.estimation.measurements.PV 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));
}

Example 3 with PV

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

use of org.orekit.estimation.measurements.PV 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());
}

Example 5 with PV

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

the class OrbitDetermination 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 (rangeRateOutliersManager != 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;
}