Examples with ParameterDriver org.orekit.utils.ParameterDriver used on opensource projects

Example 51 with ParameterDriver

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

the class RangeRateTest method testParameterDerivativesWithModifier.

@Test
public void testParameterDerivativesWithModifier() 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, 0.001);
    // create perfect range rate 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 RangeRateMeasurementCreator(context, false), 1.0, 3.0, 300.0);
    propagator.setSlaveMode();
    double maxRelativeError = 0;
    for (final ObservedMeasurement<?> measurement : measurements) {
        final RangeRateTroposphericDelayModifier modifier = new RangeRateTroposphericDelayModifier(SaastamoinenModel.getStandardModel(), true);
        ((RangeRate) measurement).addModifier(modifier);
        // parameter corresponding to station position offset
        final GroundStation stationParameter = ((RangeRate) 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
        // range would have depended only on the current position but
        // not on the current velocity.
        final double meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
        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(1, measurement.getDimension());
            Assert.assertEquals(1, gradient.length);
            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()[0];
                }
            }, drivers[i], 3, 20.0);
            final double ref = dMkdP.value(drivers[i]);
            maxRelativeError = FastMath.max(maxRelativeError, FastMath.abs((ref - gradient[0]) / ref));
        }
    }
    Assert.assertEquals(0, maxRelativeError, 1.2e-6);
}

Example 52 with ParameterDriver

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

the class PartialDerivativesTest method doTestParametersDerivatives.

private void doTestParametersDerivatives(String parameterName, double tolerance, OrbitType... orbitTypes) throws OrekitException {
    OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
    ForceModel drag = new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(), earth), new IsotropicDrag(2.5, 1.2));
    NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(5, 5);
    ForceModel gravityField = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
    Orbit baseOrbit = new KeplerianOrbit(7000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.TRUE, FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, provider.getMu());
    double dt = 900;
    double dP = 1.0;
    for (OrbitType orbitType : orbitTypes) {
        final Orbit initialOrbit = orbitType.convertType(baseOrbit);
        for (PositionAngle angleType : PositionAngle.values()) {
            NumericalPropagator propagator = setUpPropagator(initialOrbit, dP, orbitType, angleType, drag, gravityField);
            propagator.setMu(provider.getMu());
            for (final ForceModel forceModel : propagator.getAllForceModels()) {
                for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
                    driver.setValue(driver.getReferenceValue());
                    driver.setSelected(driver.getName().equals(parameterName));
                }
            }
            PartialDerivativesEquations partials = new PartialDerivativesEquations("partials", propagator);
            final SpacecraftState initialState = partials.setInitialJacobians(new SpacecraftState(initialOrbit));
            propagator.setInitialState(initialState);
            final JacobiansMapper mapper = partials.getMapper();
            PickUpHandler pickUp = new PickUpHandler(mapper, null);
            propagator.setMasterMode(pickUp);
            propagator.propagate(initialState.getDate().shiftedBy(dt));
            double[][] dYdP = pickUp.getdYdP();
            // compute reference Jacobian using finite differences
            double[][] dYdPRef = new double[6][1];
            NumericalPropagator propagator2 = setUpPropagator(initialOrbit, dP, orbitType, angleType, drag, gravityField);
            propagator2.setMu(provider.getMu());
            ParameterDriversList bound = new ParameterDriversList();
            for (final ForceModel forceModel : propagator2.getAllForceModels()) {
                for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
                    if (driver.getName().equals(parameterName)) {
                        driver.setSelected(true);
                        bound.add(driver);
                    } else {
                        driver.setSelected(false);
                    }
                }
            }
            ParameterDriver selected = bound.getDrivers().get(0);
            double p0 = selected.getReferenceValue();
            double h = selected.getScale();
            selected.setValue(p0 - 4 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sM4h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            selected.setValue(p0 - 3 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sM3h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            selected.setValue(p0 - 2 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sM2h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            selected.setValue(p0 - 1 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sM1h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            selected.setValue(p0 + 1 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sP1h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            selected.setValue(p0 + 2 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sP2h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            selected.setValue(p0 + 3 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sP3h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            selected.setValue(p0 + 4 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sP4h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            fillJacobianColumn(dYdPRef, 0, orbitType, angleType, h, sM4h, sM3h, sM2h, sM1h, sP1h, sP2h, sP3h, sP4h);
            for (int i = 0; i < 6; ++i) {
                Assert.assertEquals(dYdPRef[i][0], dYdP[i][0], FastMath.abs(dYdPRef[i][0] * tolerance));
            }
        }
    }
}

Example 53 with ParameterDriver

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

the class JacobianPropagatorConverterTest method doTestDerivatives.

private void doTestDerivatives(double tolP, double tolV, String... names) throws OrekitException {
    // we use a fixed step integrator on purpose
    // as the test is based on external differentiation using finite differences,
    // an adaptive step size integrator would introduce *lots* of numerical noise
    NumericalPropagatorBuilder builder = new NumericalPropagatorBuilder(OrbitType.CARTESIAN.convertType(orbit), new LutherIntegratorBuilder(10.0), PositionAngle.TRUE, dP);
    builder.setMass(200.0);
    builder.addForceModel(drag);
    builder.addForceModel(gravity);
    // retrieve a state slightly different from the initial state,
    // using normalized values different from 0.0 for the sake of generality
    RandomGenerator random = new Well19937a(0xe67f19c1a678d037l);
    List<ParameterDriver> all = new ArrayList<ParameterDriver>();
    for (final ParameterDriver driver : builder.getOrbitalParametersDrivers().getDrivers()) {
        all.add(driver);
    }
    for (final ParameterDriver driver : builder.getPropagationParametersDrivers().getDrivers()) {
        all.add(driver);
    }
    double[] normalized = new double[names.length];
    List<ParameterDriver> selected = new ArrayList<ParameterDriver>(names.length);
    int index = 0;
    for (final ParameterDriver driver : all) {
        boolean found = false;
        for (final String name : names) {
            if (name.equals(driver.getName())) {
                found = true;
                normalized[index++] = driver.getNormalizedValue() + (2 * random.nextDouble() - 1);
                selected.add(driver);
            }
        }
        driver.setSelected(found);
    }
    // create a one hour sample that starts 10 minutes after initial state
    // the 10 minutes offset implies even the first point is influenced by model parameters
    final List<SpacecraftState> sample = new ArrayList<SpacecraftState>();
    Propagator propagator = builder.buildPropagator(normalized);
    propagator.setMasterMode(60.0, new OrekitFixedStepHandler() {

        @Override
        public void handleStep(SpacecraftState currentState, boolean isLast) {
            sample.add(currentState);
        }
    });
    propagator.propagate(orbit.getDate().shiftedBy(600.0), orbit.getDate().shiftedBy(4200.0));
    JacobianPropagatorConverter fitter = new JacobianPropagatorConverter(builder, 1.0e-3, 5000);
    try {
        Method setSample = AbstractPropagatorConverter.class.getDeclaredMethod("setSample", List.class);
        setSample.setAccessible(true);
        setSample.invoke(fitter, sample);
    } catch (NoSuchMethodException | SecurityException | IllegalAccessException | IllegalArgumentException | InvocationTargetException e) {
        Assert.fail(e.getLocalizedMessage());
    }
    MultivariateVectorFunction f = fitter.getObjectiveFunction();
    Pair<RealVector, RealMatrix> p = fitter.getModel().value(new ArrayRealVector(normalized));
    // check derivatives
    // a h offset on normalized parameter represents a physical offset of h * scale
    RealMatrix m = p.getSecond();
    double h = 10.0;
    double[] shifted = normalized.clone();
    double maxErrorP = 0;
    double maxErrorV = 0;
    for (int j = 0; j < selected.size(); ++j) {
        shifted[j] = normalized[j] + 2.0 * h;
        double[] valueP2 = f.value(shifted);
        shifted[j] = normalized[j] + 1.0 * h;
        double[] valueP1 = f.value(shifted);
        shifted[j] = normalized[j] - 1.0 * h;
        double[] valueM1 = f.value(shifted);
        shifted[j] = normalized[j] - 2.0 * h;
        double[] valueM2 = f.value(shifted);
        shifted[j] = normalized[j];
        for (int i = 0; i < valueP2.length; ++i) {
            double d = (8 * (valueP1[i] - valueM1[i]) - (valueP2[i] - valueM2[i])) / (12 * h);
            if (i % 6 < 3) {
                // position
                maxErrorP = FastMath.max(maxErrorP, FastMath.abs(m.getEntry(i, j) - d));
            } else {
                // velocity
                maxErrorV = FastMath.max(maxErrorV, FastMath.abs(m.getEntry(i, j) - d));
            }
        }
    }
    Assert.assertEquals(0.0, maxErrorP, tolP);
    Assert.assertEquals(0.0, maxErrorV, tolV);
}

Example 54 with ParameterDriver

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

the class NumericalConverterTest method checkFit.

protected void checkFit(final Orbit orbit, final double duration, final double stepSize, final double threshold, final double expectedRMS, final String... freeParameters) throws OrekitException, IOException, ParseException {
    NumericalPropagatorBuilder builder = new NumericalPropagatorBuilder(OrbitType.CARTESIAN.convertType(orbit), new DormandPrince853IntegratorBuilder(minStep, maxStep, dP), PositionAngle.TRUE, dP);
    ForceModel guessedDrag = drag;
    ForceModel guessedGravity = gravity;
    for (String param : freeParameters) {
        if (DragSensitive.DRAG_COEFFICIENT.equals(param)) {
            // we want to adjust drag coefficient, we need to start from a wrong value
            ParameterDriver driver = drag.getParameterDriver(param);
            double coeff = driver.getReferenceValue() - driver.getScale();
            guessedDrag = new DragForce(atmosphere, new IsotropicDrag(crossSection, coeff));
        } else if (NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT.equals(param)) {
            // we want to adjust mu, we need to start from  a wrong value
            guessedGravity = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), GravityFieldFactory.getNormalizedProvider(2, 0));
            ParameterDriver driver = guessedGravity.getParameterDriver(param);
            driver.setValue(driver.getReferenceValue() + driver.getScale());
        }
    }
    builder.addForceModel(guessedDrag);
    builder.addForceModel(guessedGravity);
    JacobianPropagatorConverter fitter = new JacobianPropagatorConverter(builder, threshold, 5000);
    fitter.convert(propagator, duration, 1 + (int) (duration / stepSize), freeParameters);
    NumericalPropagator prop = (NumericalPropagator) fitter.getAdaptedPropagator();
    Orbit fitted = prop.getInitialState().getOrbit();
    for (String param : freeParameters) {
        for (ForceModel force : propagator.getAllForceModels()) {
            if (force.isSupported(param)) {
                for (ForceModel model : prop.getAllForceModels()) {
                    if (model.isSupported(param)) {
                        Assert.assertEquals(force.getParameterDriver(param).getValue(), model.getParameterDriver(param).getValue(), 3.0e-4 * FastMath.abs(force.getParameterDriver(param).getValue()));
                    }
                }
            }
        }
    }
    Assert.assertEquals(expectedRMS, fitter.getRMS(), 0.01 * expectedRMS);
    Assert.assertEquals(orbit.getPVCoordinates().getPosition().getX(), fitted.getPVCoordinates().getPosition().getX(), 1.1);
    Assert.assertEquals(orbit.getPVCoordinates().getPosition().getY(), fitted.getPVCoordinates().getPosition().getY(), 1.1);
    Assert.assertEquals(orbit.getPVCoordinates().getPosition().getZ(), fitted.getPVCoordinates().getPosition().getZ(), 1.1);
    Assert.assertEquals(orbit.getPVCoordinates().getVelocity().getX(), fitted.getPVCoordinates().getVelocity().getX(), 0.0005);
    Assert.assertEquals(orbit.getPVCoordinates().getVelocity().getY(), fitted.getPVCoordinates().getVelocity().getY(), 0.0005);
    Assert.assertEquals(orbit.getPVCoordinates().getVelocity().getZ(), fitted.getPVCoordinates().getVelocity().getZ(), 0.0005);
}

Example 55 with ParameterDriver

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

the class Model method getSelectedPropagationDriversForBuilder.

/**
 * Get the selected propagation drivers for a propagatorBuilder.
 * @param iBuilder index of the builder in the builders' array
 * @return the list of selected propagation drivers for propagatorBuilder of index iBuilder
 * @exception OrekitException if orbit cannot be created with the current point
 */
public ParameterDriversList getSelectedPropagationDriversForBuilder(final int iBuilder) throws OrekitException {
    // Lazy evaluation, create the list only if it hasn't been created yet
    if (estimatedPropagationParameters[iBuilder] == null) {
        // Gather the drivers
        final ParameterDriversList selectedPropagationDrivers = new ParameterDriversList();
        for (final DelegatingDriver delegating : builders[iBuilder].getPropagationParametersDrivers().getDrivers()) {
            if (delegating.isSelected()) {
                for (final ParameterDriver driver : delegating.getRawDrivers()) {
                    selectedPropagationDrivers.add(driver);
                }
            }
        }
        // List of propagation drivers are sorted in the BatchLSEstimator class.
        // Hence we need to sort this list so the parameters' indexes match
        selectedPropagationDrivers.sort();
        // Add the list of selected propagation drivers to the array
        estimatedPropagationParameters[iBuilder] = selectedPropagationDrivers;
    }
    return estimatedPropagationParameters[iBuilder];
}