Search in sources :

Example 36 with ParameterDriver

use of org.orekit.utils.ParameterDriver in project Orekit by CS-SI.

the class BatchLSEstimatorTest method testKeplerRange.

/**
 * Perfect range measurements with a biased start
 * @throws OrekitException
 */
@Test
public void testKeplerRange() throws OrekitException {
    Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 1.0);
    // create perfect range measurements
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new RangeMeasurementCreator(context), 1.0, 3.0, 300.0);
    // create orbit estimator
    final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(), propagatorBuilder);
    for (final ObservedMeasurement<?> range : measurements) {
        estimator.addMeasurement(range);
    }
    estimator.setParametersConvergenceThreshold(1.0e-2);
    estimator.setMaxIterations(10);
    estimator.setMaxEvaluations(20);
    estimator.setObserver(new BatchLSObserver() {

        int lastIter = 0;

        int lastEval = 0;

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

Example 37 with ParameterDriver

use of org.orekit.utils.ParameterDriver in project Orekit by CS-SI.

the class ModelTest method testPerfectValue.

@Test
public void testPerfectValue() 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) {
            Assert.assertEquals(1, newOrbits.length);
            Assert.assertEquals(0, context.initialOrbit.getDate().durationFrom(newOrbits[0].getDate()), 1.0e-15);
            Assert.assertEquals(0, Vector3D.distance(context.initialOrbit.getPVCoordinates().getPosition(), newOrbits[0].getPVCoordinates().getPosition()), 1.0e-15);
            Assert.assertEquals(measurements.size(), newEvaluations.size());
        }
    };
    final Model model = new Model(builders, measurements, estimatedMeasurementsParameters, modelObserver);
    model.setIterationsCounter(new Incrementor(100));
    model.setEvaluationsCounter(new Incrementor(100));
    // Test forward propagation flag to true
    assertEquals(true, model.isForwardPropagation());
    // evaluate model on perfect start point
    final double[] normalizedProp = propagatorBuilder.getSelectedNormalizedParameters();
    final double[] normalized = new double[normalizedProp.length + estimatedMeasurementsParameters.getNbParams()];
    System.arraycopy(normalizedProp, 0, normalized, 0, normalizedProp.length);
    int i = normalizedProp.length;
    for (final ParameterDriver driver : estimatedMeasurementsParameters.getDrivers()) {
        normalized[i++] = driver.getNormalizedValue();
    }
    Pair<RealVector, RealMatrix> value = model.value(new ArrayRealVector(normalized));
    int index = 0;
    for (ObservedMeasurement<?> measurement : measurements) {
        for (int k = 0; k < measurement.getDimension(); ++k) {
            // the value is already a weighted residual
            Assert.assertEquals(0.0, value.getFirst().getEntry(index++), 1.6e-7);
        }
    }
    Assert.assertEquals(index, value.getFirst().getDimension());
}
Also used : Context(org.orekit.estimation.Context) ArrayRealVector(org.hipparchus.linear.ArrayRealVector) Incrementor(org.hipparchus.util.Incrementor) ParameterDriver(org.orekit.utils.ParameterDriver) RealMatrix(org.hipparchus.linear.RealMatrix) ParameterDriversList(org.orekit.utils.ParameterDriversList) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) ArrayRealVector(org.hipparchus.linear.ArrayRealVector) RealVector(org.hipparchus.linear.RealVector) Propagator(org.orekit.propagation.Propagator) Map(java.util.Map) ObservedMeasurement(org.orekit.estimation.measurements.ObservedMeasurement) PVMeasurementCreator(org.orekit.estimation.measurements.PVMeasurementCreator) Test(org.junit.Test)

Example 38 with ParameterDriver

use of org.orekit.utils.ParameterDriver in project Orekit by CS-SI.

the class OrbitDeterminationTest method createPropagatorBuilder.

/**
 * Create a propagator builder from input parameters
 * @param parser input file parser
 * @param conventions IERS conventions to use
 * @param gravityField gravity field
 * @param body central body
 * @param orbit first orbit estimate
 * @return propagator builder
 * @throws NoSuchElementException if input parameters are missing
 * @throws OrekitException if body frame cannot be created
 */
private NumericalPropagatorBuilder createPropagatorBuilder(final KeyValueFileParser<ParameterKey> parser, final IERSConventions conventions, final NormalizedSphericalHarmonicsProvider gravityField, final OneAxisEllipsoid body, final Orbit orbit) throws NoSuchElementException, OrekitException {
    final double minStep;
    if (!parser.containsKey(ParameterKey.PROPAGATOR_MIN_STEP)) {
        minStep = 0.001;
    } else {
        minStep = parser.getDouble(ParameterKey.PROPAGATOR_MIN_STEP);
    }
    final double maxStep;
    if (!parser.containsKey(ParameterKey.PROPAGATOR_MAX_STEP)) {
        maxStep = 300;
    } else {
        maxStep = parser.getDouble(ParameterKey.PROPAGATOR_MAX_STEP);
    }
    final double dP;
    if (!parser.containsKey(ParameterKey.PROPAGATOR_POSITION_ERROR)) {
        dP = 10.0;
    } else {
        dP = parser.getDouble(ParameterKey.PROPAGATOR_POSITION_ERROR);
    }
    final double positionScale;
    if (!parser.containsKey(ParameterKey.ESTIMATOR_ORBITAL_PARAMETERS_POSITION_SCALE)) {
        positionScale = dP;
    } else {
        positionScale = parser.getDouble(ParameterKey.ESTIMATOR_ORBITAL_PARAMETERS_POSITION_SCALE);
    }
    final NumericalPropagatorBuilder propagatorBuilder = new NumericalPropagatorBuilder(orbit, new DormandPrince853IntegratorBuilder(minStep, maxStep, dP), PositionAngle.MEAN, positionScale);
    // initial mass
    final double mass;
    if (!parser.containsKey(ParameterKey.MASS)) {
        mass = 1000.0;
    } else {
        mass = parser.getDouble(ParameterKey.MASS);
    }
    propagatorBuilder.setMass(mass);
    // gravity field force model
    propagatorBuilder.addForceModel(new HolmesFeatherstoneAttractionModel(body.getBodyFrame(), gravityField));
    // ocean tides force model
    if (parser.containsKey(ParameterKey.OCEAN_TIDES_DEGREE) && parser.containsKey(ParameterKey.OCEAN_TIDES_ORDER)) {
        final int degree = parser.getInt(ParameterKey.OCEAN_TIDES_DEGREE);
        final int order = parser.getInt(ParameterKey.OCEAN_TIDES_ORDER);
        if (degree > 0 && order > 0) {
            propagatorBuilder.addForceModel(new OceanTides(body.getBodyFrame(), gravityField.getAe(), gravityField.getMu(), degree, order, conventions, TimeScalesFactory.getUT1(conventions, true)));
        }
    }
    // solid tides force model
    List<CelestialBody> solidTidesBodies = new ArrayList<CelestialBody>();
    if (parser.containsKey(ParameterKey.SOLID_TIDES_SUN) && parser.getBoolean(ParameterKey.SOLID_TIDES_SUN)) {
        solidTidesBodies.add(CelestialBodyFactory.getSun());
    }
    if (parser.containsKey(ParameterKey.SOLID_TIDES_MOON) && parser.getBoolean(ParameterKey.SOLID_TIDES_MOON)) {
        solidTidesBodies.add(CelestialBodyFactory.getMoon());
    }
    if (!solidTidesBodies.isEmpty()) {
        propagatorBuilder.addForceModel(new SolidTides(body.getBodyFrame(), gravityField.getAe(), gravityField.getMu(), gravityField.getTideSystem(), conventions, TimeScalesFactory.getUT1(conventions, true), solidTidesBodies.toArray(new CelestialBody[solidTidesBodies.size()])));
    }
    // third body attraction
    if (parser.containsKey(ParameterKey.THIRD_BODY_SUN) && parser.getBoolean(ParameterKey.THIRD_BODY_SUN)) {
        propagatorBuilder.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
    }
    if (parser.containsKey(ParameterKey.THIRD_BODY_MOON) && parser.getBoolean(ParameterKey.THIRD_BODY_MOON)) {
        propagatorBuilder.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
    }
    // drag
    if (parser.containsKey(ParameterKey.DRAG) && parser.getBoolean(ParameterKey.DRAG)) {
        final double cd = parser.getDouble(ParameterKey.DRAG_CD);
        final double area = parser.getDouble(ParameterKey.DRAG_AREA);
        final boolean cdEstimated = parser.getBoolean(ParameterKey.DRAG_CD_ESTIMATED);
        MarshallSolarActivityFutureEstimation msafe = new MarshallSolarActivityFutureEstimation("(?:Jan|Feb|Mar|Apr|May|Jun|Jul|Aug|Sep|Oct|Nov|Dec)\\p{Digit}\\p{Digit}\\p{Digit}\\p{Digit}F10\\.(?:txt|TXT)", MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
        DataProvidersManager manager = DataProvidersManager.getInstance();
        manager.feed(msafe.getSupportedNames(), msafe);
        Atmosphere atmosphere = new DTM2000(msafe, CelestialBodyFactory.getSun(), body);
        propagatorBuilder.addForceModel(new DragForce(atmosphere, new IsotropicDrag(area, cd)));
        if (cdEstimated) {
            for (final ParameterDriver driver : propagatorBuilder.getPropagationParametersDrivers().getDrivers()) {
                if (driver.getName().equals(DragSensitive.DRAG_COEFFICIENT)) {
                    driver.setSelected(true);
                }
            }
        }
    }
    // solar radiation pressure
    if (parser.containsKey(ParameterKey.SOLAR_RADIATION_PRESSURE) && parser.getBoolean(ParameterKey.SOLAR_RADIATION_PRESSURE)) {
        final double cr = parser.getDouble(ParameterKey.SOLAR_RADIATION_PRESSURE_CR);
        final double area = parser.getDouble(ParameterKey.SOLAR_RADIATION_PRESSURE_AREA);
        final boolean cREstimated = parser.getBoolean(ParameterKey.SOLAR_RADIATION_PRESSURE_CR_ESTIMATED);
        propagatorBuilder.addForceModel(new SolarRadiationPressure(CelestialBodyFactory.getSun(), body.getEquatorialRadius(), new IsotropicRadiationSingleCoefficient(area, cr)));
        if (cREstimated) {
            for (final ParameterDriver driver : propagatorBuilder.getPropagationParametersDrivers().getDrivers()) {
                if (driver.getName().equals(RadiationSensitive.REFLECTION_COEFFICIENT)) {
                    driver.setSelected(true);
                }
            }
        }
    }
    // post-Newtonian correction force due to general relativity
    if (parser.containsKey(ParameterKey.GENERAL_RELATIVITY) && parser.getBoolean(ParameterKey.GENERAL_RELATIVITY)) {
        propagatorBuilder.addForceModel(new Relativity(gravityField.getMu()));
    }
    // extra polynomial accelerations
    if (parser.containsKey(ParameterKey.POLYNOMIAL_ACCELERATION_NAME)) {
        final String[] names = parser.getStringArray(ParameterKey.POLYNOMIAL_ACCELERATION_NAME);
        final Vector3D[] directions = parser.getVectorArray(ParameterKey.POLYNOMIAL_ACCELERATION_DIRECTION_X, ParameterKey.POLYNOMIAL_ACCELERATION_DIRECTION_Y, ParameterKey.POLYNOMIAL_ACCELERATION_DIRECTION_Z);
        final List<String>[] coefficients = parser.getStringsListArray(ParameterKey.POLYNOMIAL_ACCELERATION_COEFFICIENTS, ',');
        final boolean[] estimated = parser.getBooleanArray(ParameterKey.POLYNOMIAL_ACCELERATION_ESTIMATED);
        for (int i = 0; i < names.length; ++i) {
            final PolynomialParametricAcceleration ppa = new PolynomialParametricAcceleration(directions[i], true, names[i], null, coefficients[i].size() - 1);
            for (int k = 0; k < coefficients[i].size(); ++k) {
                final ParameterDriver driver = ppa.getParameterDriver(names[i] + "[" + k + "]");
                driver.setValue(Double.parseDouble(coefficients[i].get(k)));
                driver.setSelected(estimated[i]);
            }
            propagatorBuilder.addForceModel(ppa);
        }
    }
    return propagatorBuilder;
}
Also used : IsotropicDrag(org.orekit.forces.drag.IsotropicDrag) PolynomialParametricAcceleration(org.orekit.forces.PolynomialParametricAcceleration) OceanTides(org.orekit.forces.gravity.OceanTides) Relativity(org.orekit.forces.gravity.Relativity) ArrayList(java.util.ArrayList) SolarRadiationPressure(org.orekit.forces.radiation.SolarRadiationPressure) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) CelestialBody(org.orekit.bodies.CelestialBody) ArrayList(java.util.ArrayList) ParameterDriversList(org.orekit.utils.ParameterDriversList) List(java.util.List) IsotropicRadiationSingleCoefficient(org.orekit.forces.radiation.IsotropicRadiationSingleCoefficient) DTM2000(org.orekit.forces.drag.atmosphere.DTM2000) SolidTides(org.orekit.forces.gravity.SolidTides) ParameterDriver(org.orekit.utils.ParameterDriver) GeodeticPoint(org.orekit.bodies.GeodeticPoint) MarshallSolarActivityFutureEstimation(org.orekit.forces.drag.atmosphere.data.MarshallSolarActivityFutureEstimation) ThirdBodyAttraction(org.orekit.forces.gravity.ThirdBodyAttraction) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) Atmosphere(org.orekit.forces.drag.atmosphere.Atmosphere) DragForce(org.orekit.forces.drag.DragForce) DormandPrince853IntegratorBuilder(org.orekit.propagation.conversion.DormandPrince853IntegratorBuilder) DataProvidersManager(org.orekit.data.DataProvidersManager) HolmesFeatherstoneAttractionModel(org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel)

Example 39 with ParameterDriver

use of org.orekit.utils.ParameterDriver in project Orekit by CS-SI.

the class AngularAzElTest 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 AngularAzElMeasurementCreator(context), 0.25, 3.0, 600.0);
    propagator.setSlaveMode();
    for (final ObservedMeasurement<?> measurement : measurements) {
        // parameter corresponding to station position offset
        final GroundStation stationParameter = ((AngularAzEl) 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-10 * FastMath.abs(ref));
                }
            }
        }
    }
}
Also used : Context(org.orekit.estimation.Context) ParameterDriver(org.orekit.utils.ParameterDriver) AbsoluteDate(org.orekit.time.AbsoluteDate) SpacecraftState(org.orekit.propagation.SpacecraftState) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) ParameterFunction(org.orekit.utils.ParameterFunction) Propagator(org.orekit.propagation.Propagator) OrekitException(org.orekit.errors.OrekitException) Test(org.junit.Test)

Example 40 with ParameterDriver

use of org.orekit.utils.ParameterDriver in project Orekit by CS-SI.

the class GroundStationTest method doTestCartesianDerivatives.

private void doTestCartesianDerivatives(double latitude, double longitude, double altitude, double stepFactor, double relativeTolerancePositionValue, double relativeTolerancePositionDerivative, double relativeToleranceVelocityValue, double relativeToleranceVelocityDerivative, String... parameterPattern) throws OrekitException {
    Utils.setDataRoot("regular-data");
    final Frame eme2000 = FramesFactory.getEME2000();
    final AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
    final AbsoluteDate date0 = date.shiftedBy(50000);
    final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
    final GroundStation station = new GroundStation(new TopocentricFrame(earth, new GeodeticPoint(latitude, longitude, altitude), "dummy"));
    final DSFactory factory = new DSFactory(parameterPattern.length, 1);
    final FieldAbsoluteDate<DerivativeStructure> dateDS = new FieldAbsoluteDate<>(factory.getDerivativeField(), date);
    ParameterDriver[] selectedDrivers = new ParameterDriver[parameterPattern.length];
    UnivariateDifferentiableVectorFunction[] dFCartesian = new UnivariateDifferentiableVectorFunction[parameterPattern.length];
    final ParameterDriver[] allDrivers = selectAllDrivers(station);
    for (ParameterDriver driver : allDrivers) {
        driver.setReferenceDate(date0);
    }
    Map<String, Integer> indices = new HashMap<>();
    for (int k = 0; k < dFCartesian.length; ++k) {
        for (int i = 0; i < allDrivers.length; ++i) {
            if (allDrivers[i].getName().matches(parameterPattern[k])) {
                selectedDrivers[k] = allDrivers[i];
                dFCartesian[k] = differentiatedStationPV(station, eme2000, date, selectedDrivers[k], stepFactor);
                indices.put(selectedDrivers[k].getName(), k);
            }
        }
    }
    ;
    DSFactory factory11 = new DSFactory(1, 1);
    RandomGenerator generator = new Well19937a(0x084d58a19c498a54l);
    double maxPositionValueRelativeError = 0;
    double maxPositionDerivativeRelativeError = 0;
    double maxVelocityValueRelativeError = 0;
    double maxVelocityDerivativeRelativeError = 0;
    for (int i = 0; i < 1000; ++i) {
        // randomly change one parameter
        ParameterDriver changed = allDrivers[generator.nextInt(allDrivers.length)];
        changed.setNormalizedValue(2 * generator.nextDouble() - 1);
        // transform to check
        FieldTransform<DerivativeStructure> t = station.getOffsetToInertial(eme2000, dateDS, factory, indices);
        FieldPVCoordinates<DerivativeStructure> pv = t.transformPVCoordinates(FieldPVCoordinates.getZero(factory.getDerivativeField()));
        for (int k = 0; k < dFCartesian.length; ++k) {
            // reference values and derivatives computed using finite differences
            DerivativeStructure[] refCartesian = dFCartesian[k].value(factory11.variable(0, selectedDrivers[k].getValue()));
            // position
            final Vector3D refP = new Vector3D(refCartesian[0].getValue(), refCartesian[1].getValue(), refCartesian[2].getValue());
            final Vector3D resP = new Vector3D(pv.getPosition().getX().getValue(), pv.getPosition().getY().getValue(), pv.getPosition().getZ().getValue());
            maxPositionValueRelativeError = FastMath.max(maxPositionValueRelativeError, Vector3D.distance(refP, resP) / refP.getNorm());
            final Vector3D refPD = new Vector3D(refCartesian[0].getPartialDerivative(1), refCartesian[1].getPartialDerivative(1), refCartesian[2].getPartialDerivative(1));
            final Vector3D resPD = new Vector3D(pv.getPosition().getX().getAllDerivatives()[k + 1], pv.getPosition().getY().getAllDerivatives()[k + 1], pv.getPosition().getZ().getAllDerivatives()[k + 1]);
            maxPositionDerivativeRelativeError = FastMath.max(maxPositionDerivativeRelativeError, Vector3D.distance(refPD, resPD) / refPD.getNorm());
            // velocity
            final Vector3D refV = new Vector3D(refCartesian[3].getValue(), refCartesian[4].getValue(), refCartesian[5].getValue());
            final Vector3D resV = new Vector3D(pv.getVelocity().getX().getValue(), pv.getVelocity().getY().getValue(), pv.getVelocity().getZ().getValue());
            maxVelocityValueRelativeError = FastMath.max(maxVelocityValueRelativeError, Vector3D.distance(refV, resV) / refV.getNorm());
            final Vector3D refVD = new Vector3D(refCartesian[3].getPartialDerivative(1), refCartesian[4].getPartialDerivative(1), refCartesian[5].getPartialDerivative(1));
            final Vector3D resVD = new Vector3D(pv.getVelocity().getX().getAllDerivatives()[k + 1], pv.getVelocity().getY().getAllDerivatives()[k + 1], pv.getVelocity().getZ().getAllDerivatives()[k + 1]);
            maxVelocityDerivativeRelativeError = FastMath.max(maxVelocityDerivativeRelativeError, Vector3D.distance(refVD, resVD) / refVD.getNorm());
        }
    }
    Assert.assertEquals(0.0, maxPositionValueRelativeError, relativeTolerancePositionValue);
    Assert.assertEquals(0.0, maxPositionDerivativeRelativeError, relativeTolerancePositionDerivative);
    Assert.assertEquals(0.0, maxVelocityValueRelativeError, relativeToleranceVelocityValue);
    Assert.assertEquals(0.0, maxVelocityDerivativeRelativeError, relativeToleranceVelocityDerivative);
}
Also used : Frame(org.orekit.frames.Frame) TopocentricFrame(org.orekit.frames.TopocentricFrame) OneAxisEllipsoid(org.orekit.bodies.OneAxisEllipsoid) HashMap(java.util.HashMap) TopocentricFrame(org.orekit.frames.TopocentricFrame) Well19937a(org.hipparchus.random.Well19937a) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) AbsoluteDate(org.orekit.time.AbsoluteDate) RandomGenerator(org.hipparchus.random.RandomGenerator) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) GeodeticPoint(org.orekit.bodies.GeodeticPoint) DerivativeStructure(org.hipparchus.analysis.differentiation.DerivativeStructure) DSFactory(org.hipparchus.analysis.differentiation.DSFactory) ParameterDriver(org.orekit.utils.ParameterDriver) GeodeticPoint(org.orekit.bodies.GeodeticPoint) UnivariateDifferentiableVectorFunction(org.hipparchus.analysis.differentiation.UnivariateDifferentiableVectorFunction) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate)

Aggregations

ParameterDriver (org.orekit.utils.ParameterDriver)80 AbsoluteDate (org.orekit.time.AbsoluteDate)33 SpacecraftState (org.orekit.propagation.SpacecraftState)32 NumericalPropagatorBuilder (org.orekit.propagation.conversion.NumericalPropagatorBuilder)27 Test (org.junit.Test)23 Propagator (org.orekit.propagation.Propagator)23 Context (org.orekit.estimation.Context)21 ParameterDriversList (org.orekit.utils.ParameterDriversList)20 OrekitException (org.orekit.errors.OrekitException)19 Vector3D (org.hipparchus.geometry.euclidean.threed.Vector3D)16 FieldAbsoluteDate (org.orekit.time.FieldAbsoluteDate)16 ObservedMeasurement (org.orekit.estimation.measurements.ObservedMeasurement)15 Orbit (org.orekit.orbits.Orbit)15 ArrayList (java.util.ArrayList)14 DerivativeStructure (org.hipparchus.analysis.differentiation.DerivativeStructure)14 ParameterFunction (org.orekit.utils.ParameterFunction)14 DSFactory (org.hipparchus.analysis.differentiation.DSFactory)13 HashMap (java.util.HashMap)11 KeplerianOrbit (org.orekit.orbits.KeplerianOrbit)11 RealMatrix (org.hipparchus.linear.RealMatrix)10