Examples with PartialDerivativesEquations org.orekit.propagation.numerical.PartialDerivativesEquations used on opensource projects

Example 1 with PartialDerivativesEquations

use of org.orekit.propagation.numerical.PartialDerivativesEquations in project Orekit by CS-SI.

the class JacobianPropagatorConverter method getModel.

/**
 * {@inheritDoc}
 */
protected MultivariateJacobianFunction getModel() {
    return new MultivariateJacobianFunction() {

        /**
         * {@inheritDoc}
         */
        public Pair<RealVector, RealMatrix> value(final RealVector point) throws IllegalArgumentException, OrekitExceptionWrapper {
            try {
                final RealVector value = new ArrayRealVector(getTargetSize());
                final RealMatrix jacobian = MatrixUtils.createRealMatrix(getTargetSize(), point.getDimension());
                final NumericalPropagator prop = builder.buildPropagator(point.toArray());
                final int stateSize = isOnlyPosition() ? 3 : 6;
                final ParameterDriversList orbitalParameters = builder.getOrbitalParametersDrivers();
                final PartialDerivativesEquations pde = new PartialDerivativesEquations("pde", prop);
                final ParameterDriversList propagationParameters = pde.getSelectedParameters();
                prop.setInitialState(pde.setInitialJacobians(prop.getInitialState()));
                final JacobiansMapper mapper = pde.getMapper();
                final List<SpacecraftState> sample = getSample();
                for (int i = 0; i < sample.size(); ++i) {
                    final int row = i * stateSize;
                    if (prop.getInitialState().getDate().equals(sample.get(i).getDate())) {
                        // use initial state and Jacobians
                        fillRows(value, jacobian, row, prop.getInitialState(), stateSize, orbitalParameters, propagationParameters, mapper);
                    } else {
                        // use a date detector to pick up state and Jacobians
                        prop.addEventDetector(new DateDetector(sample.get(i).getDate()).withHandler(new EventHandler<DateDetector>() {

                            /**
                             * {@inheritDoc}
                             */
                            @Override
                            public Action eventOccurred(final SpacecraftState state, final DateDetector detector, final boolean increasing) throws OrekitException {
                                fillRows(value, jacobian, row, state, stateSize, orbitalParameters, propagationParameters, mapper);
                                return row + stateSize >= getTargetSize() ? Action.STOP : Action.CONTINUE;
                            }
                        }));
                    }
                }
                prop.propagate(sample.get(sample.size() - 1).getDate().shiftedBy(10.0));
                return new Pair<RealVector, RealMatrix>(value, jacobian);
            } catch (OrekitException ex) {
                throw new OrekitExceptionWrapper(ex);
            }
        }
    };
}

Example 2 with PartialDerivativesEquations

use of org.orekit.propagation.numerical.PartialDerivativesEquations in project Orekit by CS-SI.

the class Model method configureDerivatives.

/**
 * Configure the propagator to compute derivatives.
 * @param propagator {@link Propagator} to configure
 * @return mapper for this propagator
 * @exception OrekitException if orbit cannot be created with the current point
 */
private JacobiansMapper configureDerivatives(final NumericalPropagator propagator) throws OrekitException {
    final String equationName = Model.class.getName() + "-derivatives";
    final PartialDerivativesEquations partials = new PartialDerivativesEquations(equationName, propagator);
    // add the derivatives to the initial state
    final SpacecraftState rawState = propagator.getInitialState();
    final SpacecraftState stateWithDerivatives = partials.setInitialJacobians(rawState);
    propagator.resetInitialState(stateWithDerivatives);
    return partials.getMapper();
}

Example 3 with PartialDerivativesEquations

use of org.orekit.propagation.numerical.PartialDerivativesEquations in project Orekit by CS-SI.

the class Model method updateReferenceTrajectory.

/**
 * Update the reference trajectory using the propagator as input.
 * @param propagator The new propagator to use
 * @throws OrekitException if setting up the partial derivatives failed
 */
private void updateReferenceTrajectory(final NumericalPropagator propagator) throws OrekitException {
    // Update the reference trajectory propagator
    referenceTrajectory = propagator;
    // Link the partial derivatives to this new propagator
    final String equationName = KalmanEstimator.class.getName() + "-derivatives";
    final PartialDerivativesEquations pde = new PartialDerivativesEquations(equationName, referenceTrajectory);
    // Reset the Jacobians
    final SpacecraftState rawState = referenceTrajectory.getInitialState();
    final SpacecraftState stateWithDerivatives = pde.setInitialJacobians(rawState);
    referenceTrajectory.resetInitialState(stateWithDerivatives);
    jacobiansMapper = pde.getMapper();
}

Example 4 with PartialDerivativesEquations

use of org.orekit.propagation.numerical.PartialDerivativesEquations in project Orekit by CS-SI.

the class ImpulseManeuverTest method testAdditionalStateNumerical.

@Test
public void testAdditionalStateNumerical() throws OrekitException {
    final double mu = CelestialBodyFactory.getEarth().getGM();
    final double initialX = 7100e3;
    final double initialY = 0.0;
    final double initialZ = 1300e3;
    final double initialVx = 0;
    final double initialVy = 8000;
    final double initialVz = 1000;
    final Vector3D position = new Vector3D(initialX, initialY, initialZ);
    final Vector3D velocity = new Vector3D(initialVx, initialVy, initialVz);
    final AbsoluteDate epoch = new AbsoluteDate(2010, 1, 1, 0, 0, 0, TimeScalesFactory.getUTC());
    final TimeStampedPVCoordinates pv = new TimeStampedPVCoordinates(epoch, position, velocity, Vector3D.ZERO);
    final Orbit initialOrbit = new CartesianOrbit(pv, FramesFactory.getEME2000(), mu);
    final double totalPropagationTime = 10.0;
    final double deltaX = 0.01;
    final double deltaY = 0.02;
    final double deltaZ = 0.03;
    final double isp = 300;
    final Vector3D deltaV = new Vector3D(deltaX, deltaY, deltaZ);
    final AttitudeProvider attitudeProvider = new LofOffset(initialOrbit.getFrame(), LOFType.VNC);
    final Attitude initialAttitude = attitudeProvider.getAttitude(initialOrbit, initialOrbit.getDate(), initialOrbit.getFrame());
    double[][] tolerances = NumericalPropagator.tolerances(10.0, initialOrbit, initialOrbit.getType());
    DormandPrince853Integrator integrator = new DormandPrince853Integrator(1.0e-3, 60, tolerances[0], tolerances[1]);
    NumericalPropagator propagator = new NumericalPropagator(integrator);
    propagator.setOrbitType(initialOrbit.getType());
    PartialDerivativesEquations pde = new PartialDerivativesEquations("derivatives", propagator);
    final SpacecraftState initialState = pde.setInitialJacobians(new SpacecraftState(initialOrbit, initialAttitude));
    propagator.resetInitialState(initialState);
    DateDetector dateDetector = new DateDetector(epoch.shiftedBy(0.5 * totalPropagationTime));
    InertialProvider attitudeOverride = new InertialProvider(new Rotation(RotationOrder.XYX, RotationConvention.VECTOR_OPERATOR, 0, 0, 0));
    ImpulseManeuver<DateDetector> burnAtEpoch = new ImpulseManeuver<DateDetector>(dateDetector, attitudeOverride, deltaV, isp).withThreshold(1.0e-3);
    propagator.addEventDetector(burnAtEpoch);
    SpacecraftState finalState = propagator.propagate(epoch.shiftedBy(totalPropagationTime));
    Assert.assertEquals(1, finalState.getAdditionalStates().size());
    Assert.assertEquals(36, finalState.getAdditionalState("derivatives").length);
    double[][] stateTransitionMatrix = new double[6][6];
    pde.getMapper().getStateJacobian(finalState, stateTransitionMatrix);
    for (int i = 0; i < 6; ++i) {
        for (int j = 0; j < 6; ++j) {
            double sIJ = stateTransitionMatrix[i][j];
            if (j == i) {
                // dPi/dPj and dVi/dVj are roughly 1 for small propagation times
                Assert.assertEquals(1.0, sIJ, 2.0e-4);
            } else if (j == i + 3) {
                // dVi/dPi is roughly the propagation time for small propagation times
                Assert.assertEquals(totalPropagationTime, sIJ, 4.0e-5 * totalPropagationTime);
            } else {
                // other derivatives are almost zero for small propagation times
                Assert.assertEquals(0, sIJ, 1.0e-4);
            }
        }
    }
}

Example 5 with PartialDerivativesEquations

use of org.orekit.propagation.numerical.PartialDerivativesEquations in project Orekit by CS-SI.

the class DragForceTest method testIssue229.

@Test
public void testIssue229() throws OrekitException {
    AbsoluteDate initialDate = new AbsoluteDate(2004, 1, 1, 0, 0, 0., TimeScalesFactory.getUTC());
    Frame frame = FramesFactory.getEME2000();
    double rpe = 160.e3 + Constants.WGS84_EARTH_EQUATORIAL_RADIUS;
    double rap = 2000.e3 + Constants.WGS84_EARTH_EQUATORIAL_RADIUS;
    double inc = FastMath.toRadians(0.);
    double aop = FastMath.toRadians(0.);
    double raan = FastMath.toRadians(0.);
    double mean = FastMath.toRadians(180.);
    double mass = 100.;
    KeplerianOrbit orbit = new KeplerianOrbit(0.5 * (rpe + rap), (rap - rpe) / (rpe + rap), inc, aop, raan, mean, PositionAngle.MEAN, frame, initialDate, Constants.EIGEN5C_EARTH_MU);
    IsotropicDrag shape = new IsotropicDrag(10., 2.2);
    Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
    BodyShape earthShape = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, itrf);
    Atmosphere atmosphere = new SimpleExponentialAtmosphere(earthShape, 2.6e-10, 200000, 26000);
    double[][] tolerance = NumericalPropagator.tolerances(0.1, orbit, OrbitType.CARTESIAN);
    AbstractIntegrator integrator = new DormandPrince853Integrator(1.0e-3, 300, tolerance[0], tolerance[1]);
    NumericalPropagator propagator = new NumericalPropagator(integrator);
    propagator.setOrbitType(OrbitType.CARTESIAN);
    propagator.setMu(orbit.getMu());
    propagator.addForceModel(new DragForce(atmosphere, shape));
    PartialDerivativesEquations partials = new PartialDerivativesEquations("partials", propagator);
    propagator.setInitialState(partials.setInitialJacobians(new SpacecraftState(orbit, mass)));
    SpacecraftState state = propagator.propagate(new AbsoluteDate(2004, 1, 1, 1, 30, 0., TimeScalesFactory.getUTC()));
    double delta = 0.1;
    Orbit shifted = new CartesianOrbit(new TimeStampedPVCoordinates(orbit.getDate(), orbit.getPVCoordinates().getPosition().add(new Vector3D(delta, 0, 0)), orbit.getPVCoordinates().getVelocity()), orbit.getFrame(), orbit.getMu());
    propagator.setInitialState(partials.setInitialJacobians(new SpacecraftState(shifted, mass)));
    SpacecraftState newState = propagator.propagate(new AbsoluteDate(2004, 1, 1, 1, 30, 0., TimeScalesFactory.getUTC()));
    double[] dPVdX = new double[] { (newState.getPVCoordinates().getPosition().getX() - state.getPVCoordinates().getPosition().getX()) / delta, (newState.getPVCoordinates().getPosition().getY() - state.getPVCoordinates().getPosition().getY()) / delta, (newState.getPVCoordinates().getPosition().getZ() - state.getPVCoordinates().getPosition().getZ()) / delta, (newState.getPVCoordinates().getVelocity().getX() - state.getPVCoordinates().getVelocity().getX()) / delta, (newState.getPVCoordinates().getVelocity().getY() - state.getPVCoordinates().getVelocity().getY()) / delta, (newState.getPVCoordinates().getVelocity().getZ() - state.getPVCoordinates().getVelocity().getZ()) / delta };
    double[][] dYdY0 = new double[6][6];
    partials.getMapper().getStateJacobian(state, dYdY0);
    for (int i = 0; i < 6; ++i) {
        Assert.assertEquals(dPVdX[i], dYdY0[i][0], 6.2e-6 * FastMath.abs(dPVdX[i]));
    }
}