Examples with NumericalPropagatorBuilder org.orekit.propagation.conversion.NumericalPropagatorBuilder used on opensource projects

Example 16 with NumericalPropagatorBuilder

use of org.orekit.propagation.conversion.NumericalPropagatorBuilder in project Orekit by CS-SI.

the class ModelTest method testBackwardPropagation.

@Test
public void testBackwardPropagation() throws OrekitException {
    final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 0.001);
    final NumericalPropagatorBuilder[] builders = { propagatorBuilder };
    // create perfect PV measurements
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new PVMeasurementCreator(), 0.0, -1.0, 300.0);
    final ParameterDriversList estimatedMeasurementsParameters = new ParameterDriversList();
    for (ObservedMeasurement<?> measurement : measurements) {
        for (final ParameterDriver driver : measurement.getParametersDrivers()) {
            if (driver.isSelected()) {
                estimatedMeasurementsParameters.add(driver);
            }
        }
    }
    // create model
    final ModelObserver modelObserver = new ModelObserver() {

        /**
         * {@inheritDoc}
         */
        @Override
        public void modelCalled(final Orbit[] newOrbits, final Map<ObservedMeasurement<?>, EstimatedMeasurement<?>> newEvaluations) {
        // Do nothing here
        }
    };
    final Model model = new Model(builders, measurements, estimatedMeasurementsParameters, modelObserver);
    // Test forward propagation flag to false
    assertEquals(false, model.isForwardPropagation());
}

Example 17 with NumericalPropagatorBuilder

use of org.orekit.propagation.conversion.NumericalPropagatorBuilder 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 18 with NumericalPropagatorBuilder

use of org.orekit.propagation.conversion.NumericalPropagatorBuilder in project Orekit by CS-SI.

the class AngularRaDecTest method testStateDerivatives.

@Test
public void testStateDerivatives() throws OrekitException {
    Context context = EstimationTestUtils.geoStationnaryContext("regular-data:potential:tides");
    final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.EQUINOCTIAL, PositionAngle.TRUE, false, 1.0e-6, 60.0, 0.001);
    // create perfect azimuth-elevation measurements
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new AngularRaDecMeasurementCreator(context), 0.25, 3.0, 600.0);
    propagator.setSlaveMode();
    // Compute measurements.
    double[] RaerrorsP = new double[3 * measurements.size()];
    double[] RaerrorsV = new double[3 * measurements.size()];
    double[] DecerrorsP = new double[3 * measurements.size()];
    double[] DecerrorsV = new double[3 * measurements.size()];
    int RaindexP = 0;
    int RaindexV = 0;
    int DecindexP = 0;
    int DecindexV = 0;
    for (final ObservedMeasurement<?> measurement : measurements) {
        // parameter corresponding to station position offset
        final GroundStation stationParameter = ((AngularRaDec) measurement).getStation();
        // We intentionally propagate to a date which is close to the
        // real spacecraft state but is *not* the accurate date, by
        // compensating only part of the downlink delay. This is done
        // in order to validate the partial derivatives with respect
        // to velocity. If we had chosen the proper state date, the
        // angular would have depended only on the current position but
        // not on the current velocity.
        final AbsoluteDate datemeas = measurement.getDate();
        SpacecraftState state = propagator.propagate(datemeas);
        final Vector3D stationP = stationParameter.getOffsetToInertial(state.getFrame(), datemeas).transformPosition(Vector3D.ZERO);
        final double meanDelay = AbstractMeasurement.signalTimeOfFlight(state.getPVCoordinates(), stationP, datemeas);
        final AbsoluteDate date = measurement.getDate().shiftedBy(-0.75 * meanDelay);
        state = propagator.propagate(date);
        final EstimatedMeasurement<?> estimated = measurement.estimate(0, 0, new SpacecraftState[] { state });
        Assert.assertEquals(2, estimated.getParticipants().length);
        final double[][] jacobian = estimated.getStateDerivatives(0);
        // compute a reference value using finite differences
        final double[][] finiteDifferencesJacobian = Differentiation.differentiate(new StateFunction() {

            public double[] value(final SpacecraftState state) throws OrekitException {
                return measurement.estimate(0, 0, new SpacecraftState[] { state }).getEstimatedValue();
            }
        }, measurement.getDimension(), propagator.getAttitudeProvider(), OrbitType.CARTESIAN, PositionAngle.TRUE, 250.0, 4).value(state);
        Assert.assertEquals(finiteDifferencesJacobian.length, jacobian.length);
        Assert.assertEquals(finiteDifferencesJacobian[0].length, jacobian[0].length);
        final double smallest = FastMath.ulp((double) 1.0);
        for (int i = 0; i < jacobian.length; ++i) {
            for (int j = 0; j < jacobian[i].length; ++j) {
                double relativeError = FastMath.abs((finiteDifferencesJacobian[i][j] - jacobian[i][j]) / finiteDifferencesJacobian[i][j]);
                if ((FastMath.sqrt(finiteDifferencesJacobian[i][j]) < smallest) && (FastMath.sqrt(jacobian[i][j]) < smallest)) {
                    relativeError = 0.0;
                }
                if (j < 3) {
                    if (i == 0) {
                        RaerrorsP[RaindexP++] = relativeError;
                    } else {
                        DecerrorsP[DecindexP++] = relativeError;
                    }
                } else {
                    if (i == 0) {
                        RaerrorsV[RaindexV++] = relativeError;
                    } else {
                        DecerrorsV[DecindexV++] = relativeError;
                    }
                }
            }
        }
    }
    // median errors on Azimuth
    Assert.assertEquals(0.0, new Median().evaluate(RaerrorsP), 4.8e-11);
    Assert.assertEquals(0.0, new Median().evaluate(RaerrorsV), 2.2e-5);
    // median errors on Elevation
    Assert.assertEquals(0.0, new Median().evaluate(DecerrorsP), 1.5e-11);
    Assert.assertEquals(0.0, new Median().evaluate(DecerrorsV), 5.4e-6);
}

Example 19 with NumericalPropagatorBuilder

use of org.orekit.propagation.conversion.NumericalPropagatorBuilder in project Orekit by CS-SI.

the class AngularRaDecTest method testParameterDerivatives.

@Test
public void testParameterDerivatives() throws OrekitException {
    Context context = EstimationTestUtils.geoStationnaryContext("regular-data:potential:tides");
    final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.EQUINOCTIAL, PositionAngle.TRUE, false, 1.0e-6, 60.0, 0.001);
    // create perfect azimuth-elevation measurements
    for (final GroundStation station : context.stations) {
        station.getEastOffsetDriver().setSelected(true);
        station.getNorthOffsetDriver().setSelected(true);
        station.getZenithOffsetDriver().setSelected(true);
    }
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new AngularRaDecMeasurementCreator(context), 0.25, 3.0, 600.0);
    propagator.setSlaveMode();
    for (final ObservedMeasurement<?> measurement : measurements) {
        // parameter corresponding to station position offset
        final GroundStation stationParameter = ((AngularRaDec) measurement).getStation();
        // We intentionally propagate to a date which is close to the
        // real spacecraft state but is *not* the accurate date, by
        // compensating only part of the downlink delay. This is done
        // in order to validate the partial derivatives with respect
        // to velocity. If we had chosen the proper state date, the
        // angular would have depended only on the current position but
        // not on the current velocity.
        final AbsoluteDate datemeas = measurement.getDate();
        final SpacecraftState stateini = propagator.propagate(datemeas);
        final Vector3D stationP = stationParameter.getOffsetToInertial(stateini.getFrame(), datemeas).transformPosition(Vector3D.ZERO);
        final double meanDelay = AbstractMeasurement.signalTimeOfFlight(stateini.getPVCoordinates(), stationP, datemeas);
        final AbsoluteDate date = measurement.getDate().shiftedBy(-0.75 * meanDelay);
        final SpacecraftState state = propagator.propagate(date);
        final ParameterDriver[] drivers = new ParameterDriver[] { stationParameter.getEastOffsetDriver(), stationParameter.getNorthOffsetDriver(), stationParameter.getZenithOffsetDriver() };
        for (int i = 0; i < 3; ++i) {
            final double[] gradient = measurement.estimate(0, 0, new SpacecraftState[] { state }).getParameterDerivatives(drivers[i]);
            Assert.assertEquals(2, measurement.getDimension());
            Assert.assertEquals(2, gradient.length);
            for (final int k : new int[] { 0, 1 }) {
                final ParameterFunction dMkdP = Differentiation.differentiate(new ParameterFunction() {

                    /**
                     * {@inheritDoc}
                     */
                    @Override
                    public double value(final ParameterDriver parameterDriver) throws OrekitException {
                        return measurement.estimate(0, 0, new SpacecraftState[] { state }).getEstimatedValue()[k];
                    }
                }, drivers[i], 3, 50.0);
                final double ref = dMkdP.value(drivers[i]);
                if (ref > 1.e-12) {
                    Assert.assertEquals(ref, gradient[k], 3e-9 * FastMath.abs(ref));
                }
            }
        }
    }
}

Example 20 with NumericalPropagatorBuilder

use of org.orekit.propagation.conversion.NumericalPropagatorBuilder in project Orekit by CS-SI.

the class GroundStationTest method testEstimateStationPosition.

@Test
public void testEstimateStationPosition() throws OrekitException, IOException, ClassNotFoundException {
    Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    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 RangeMeasurementCreator(context), 1.0, 3.0, 300.0);
    // move one station
    final RandomGenerator random = new Well19937a(0x4adbecfc743bda60l);
    final TopocentricFrame base = context.stations.get(0).getBaseFrame();
    final BodyShape parent = base.getParentShape();
    final Vector3D baseOrigin = parent.transform(base.getPoint());
    final Vector3D deltaTopo = new Vector3D(2 * random.nextDouble() - 1, 2 * random.nextDouble() - 1, 2 * random.nextDouble() - 1);
    final Transform topoToParent = base.getTransformTo(parent.getBodyFrame(), (AbsoluteDate) null);
    final Vector3D deltaParent = topoToParent.transformVector(deltaTopo);
    final String movedSuffix = "-moved";
    final GroundStation moved = new GroundStation(new TopocentricFrame(parent, parent.transform(baseOrigin.subtract(deltaParent), parent.getBodyFrame(), null), base.getName() + movedSuffix), context.ut1.getEOPHistory(), context.stations.get(0).getDisplacements());
    // create orbit estimator
    final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(), propagatorBuilder);
    for (final ObservedMeasurement<?> measurement : measurements) {
        final Range range = (Range) measurement;
        final String name = range.getStation().getBaseFrame().getName() + movedSuffix;
        if (moved.getBaseFrame().getName().equals(name)) {
            estimator.addMeasurement(new Range(moved, range.getDate(), range.getObservedValue()[0], range.getTheoreticalStandardDeviation()[0], range.getBaseWeight()[0]));
        } else {
            estimator.addMeasurement(range);
        }
    }
    estimator.setParametersConvergenceThreshold(1.0e-3);
    estimator.setMaxIterations(100);
    estimator.setMaxEvaluations(200);
    // we want to estimate station offsets
    moved.getEastOffsetDriver().setSelected(true);
    moved.getNorthOffsetDriver().setSelected(true);
    moved.getZenithOffsetDriver().setSelected(true);
    EstimationTestUtils.checkFit(context, estimator, 2, 3, 0.0, 5.6e-7, 0.0, 1.4e-6, 0.0, 4.8e-7, 0.0, 2.6e-10);
    Assert.assertEquals(deltaTopo.getX(), moved.getEastOffsetDriver().getValue(), 4.5e-7);
    Assert.assertEquals(deltaTopo.getY(), moved.getNorthOffsetDriver().getValue(), 6.2e-7);
    Assert.assertEquals(deltaTopo.getZ(), moved.getZenithOffsetDriver().getValue(), 2.6e-7);
    GeodeticPoint result = moved.getOffsetGeodeticPoint(null);
    GeodeticPoint reference = context.stations.get(0).getBaseFrame().getPoint();
    Assert.assertEquals(reference.getLatitude(), result.getLatitude(), 1.4e-14);
    Assert.assertEquals(reference.getLongitude(), result.getLongitude(), 2.9e-14);
    Assert.assertEquals(reference.getAltitude(), result.getAltitude(), 2.6e-7);
    RealMatrix normalizedCovariances = estimator.getOptimum().getCovariances(1.0e-10);
    RealMatrix physicalCovariances = estimator.getPhysicalCovariances(1.0e-10);
    Assert.assertEquals(9, normalizedCovariances.getRowDimension());
    Assert.assertEquals(9, normalizedCovariances.getColumnDimension());
    Assert.assertEquals(9, physicalCovariances.getRowDimension());
    Assert.assertEquals(9, physicalCovariances.getColumnDimension());
    Assert.assertEquals(0.55431, physicalCovariances.getEntry(6, 6), 1.0e-5);
    Assert.assertEquals(0.22694, physicalCovariances.getEntry(7, 7), 1.0e-5);
    Assert.assertEquals(0.13106, physicalCovariances.getEntry(8, 8), 1.0e-5);
    ByteArrayOutputStream bos = new ByteArrayOutputStream();
    ObjectOutputStream oos = new ObjectOutputStream(bos);
    oos.writeObject(moved.getEstimatedEarthFrame().getTransformProvider());
    Assert.assertTrue(bos.size() > 155000);
    Assert.assertTrue(bos.size() < 160000);
    ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
    ObjectInputStream ois = new ObjectInputStream(bis);
    EstimatedEarthFrameProvider deserialized = (EstimatedEarthFrameProvider) ois.readObject();
    Assert.assertEquals(moved.getPrimeMeridianOffsetDriver().getValue(), deserialized.getPrimeMeridianOffsetDriver().getValue(), 1.0e-15);
    Assert.assertEquals(moved.getPrimeMeridianDriftDriver().getValue(), deserialized.getPrimeMeridianDriftDriver().getValue(), 1.0e-15);
    Assert.assertEquals(moved.getPolarOffsetXDriver().getValue(), deserialized.getPolarOffsetXDriver().getValue(), 1.0e-15);
    Assert.assertEquals(moved.getPolarDriftXDriver().getValue(), deserialized.getPolarDriftXDriver().getValue(), 1.0e-15);
    Assert.assertEquals(moved.getPolarOffsetYDriver().getValue(), deserialized.getPolarOffsetYDriver().getValue(), 1.0e-15);
    Assert.assertEquals(moved.getPolarDriftYDriver().getValue(), deserialized.getPolarDriftYDriver().getValue(), 1.0e-15);
}