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

Example 56 with NumericalPropagator

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

the class KalmanEstimatorTest method testKeplerianRange.

/**
 * Perfect range measurements with a biased start
 * Keplerian formalism
 * @throws OrekitException
 */
@Test
public void testKeplerianRange() throws OrekitException {
    // Create context
    Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    // Create initial orbit and propagator builder
    final OrbitType orbitType = OrbitType.KEPLERIAN;
    final PositionAngle positionAngle = PositionAngle.TRUE;
    final boolean perfectStart = true;
    final double minStep = 1.e-6;
    final double maxStep = 60.;
    final double dP = 1.;
    final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(orbitType, positionAngle, perfectStart, minStep, maxStep, dP);
    // Create perfect range measurements
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new RangeMeasurementCreator(context), 1.0, 4.0, 60.0);
    // Reference propagator for estimation performances
    final NumericalPropagator referencePropagator = propagatorBuilder.buildPropagator(propagatorBuilder.getSelectedNormalizedParameters());
    // Reference position/velocity at last measurement date
    final Orbit refOrbit = referencePropagator.propagate(measurements.get(measurements.size() - 1).getDate()).getOrbit();
    // Change semi-major axis of 1.2m as in the batch test
    ParameterDriver aDriver = propagatorBuilder.getOrbitalParametersDrivers().getDrivers().get(0);
    aDriver.setValue(aDriver.getValue() + 1.2);
    aDriver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
    // Cartesian covariance matrix initialization
    // 100m on position / 1e-2m/s on velocity
    final RealMatrix cartesianP = MatrixUtils.createRealDiagonalMatrix(new double[] { 100., 100., 100., 1e-2, 1e-2, 1e-2 });
    // Jacobian of the orbital parameters w/r to Cartesian
    final Orbit initialOrbit = orbitType.convertType(context.initialOrbit);
    final double[][] dYdC = new double[6][6];
    initialOrbit.getJacobianWrtCartesian(PositionAngle.TRUE, dYdC);
    final RealMatrix Jac = MatrixUtils.createRealMatrix(dYdC);
    // Keplerian initial covariance matrix
    final RealMatrix initialP = Jac.multiply(cartesianP.multiply(Jac.transpose()));
    // Process noise matrix is set to 0 here
    RealMatrix Q = MatrixUtils.createRealMatrix(6, 6);
    // Build the Kalman filter
    final KalmanEstimatorBuilder kalmanBuilder = new KalmanEstimatorBuilder();
    kalmanBuilder.builder(propagatorBuilder);
    kalmanBuilder.estimatedMeasurementsParameters(new ParameterDriversList());
    kalmanBuilder.initialCovarianceMatrix(initialP);
    kalmanBuilder.processNoiseMatrixProvider(new ConstantProcessNoise(Q));
    final KalmanEstimator kalman = kalmanBuilder.build();
    // Filter the measurements and check the results
    final double expectedDeltaPos = 0.;
    final double posEps = 1.77e-4;
    final double expectedDeltaVel = 0.;
    final double velEps = 7.93e-8;
    final double[] expectedSigmasPos = { 0.742488, 0.281910, 0.563217 };
    final double sigmaPosEps = 1e-6;
    final double[] expectedSigmasVel = { 2.206622e-4, 1.306669e-4, 1.293996e-4 };
    final double sigmaVelEps = 1e-10;
    EstimationTestUtils.checkKalmanFit(context, kalman, measurements, refOrbit, positionAngle, expectedDeltaPos, posEps, expectedDeltaVel, velEps, expectedSigmasPos, sigmaPosEps, expectedSigmasVel, sigmaVelEps);
}

Example 57 with NumericalPropagator

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

the class KalmanOrbitDeterminationTest method runReference.

/**
 * Use the physical models in the input file
 * Incorporate the initial reference values
 * And run the propagation until the last measurement to get the reference orbit at the same date
 * as the Kalman filter
 * @param input Input configuration file
 * @param orbitType Orbit type to use (calculation and display)
 * @param refPosition Initial reference position
 * @param refVelocity Initial reference velocity
 * @param refPropagationParameters Reference propagation parameters
 * @param kalmanFinalDate The final date of the Kalman filter
 * @return The reference orbit at the same date as the Kalman filter
 * @throws IOException Input file cannot be opened
 * @throws IllegalArgumentException Issue in key/value reading of input file
 * @throws OrekitException An Orekit exception... should be explicit
 * @throws ParseException Parsing of the input file or measurement file failed
 */
private Orbit runReference(final File input, final OrbitType orbitType, final Vector3D refPosition, final Vector3D refVelocity, final ParameterDriversList refPropagationParameters, final AbsoluteDate kalmanFinalDate) throws IOException, IllegalArgumentException, OrekitException, ParseException {
    // Read input parameters
    KeyValueFileParser<ParameterKey> parser = new KeyValueFileParser<ParameterKey>(ParameterKey.class);
    parser.parseInput(input.getAbsolutePath(), new FileInputStream(input));
    // Gravity field
    GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-5c.gfc", true));
    final NormalizedSphericalHarmonicsProvider gravityField = createGravityField(parser);
    // Orbit initial guess
    Orbit initialRefOrbit = new CartesianOrbit(new PVCoordinates(refPosition, refVelocity), parser.getInertialFrame(ParameterKey.INERTIAL_FRAME), parser.getDate(ParameterKey.ORBIT_DATE, TimeScalesFactory.getUTC()), gravityField.getMu());
    // Convert to desired orbit type
    initialRefOrbit = orbitType.convertType(initialRefOrbit);
    // 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, initialRefOrbit);
    // Force the selected propagation parameters to their reference values
    if (refPropagationParameters != null) {
        for (DelegatingDriver refDriver : refPropagationParameters.getDrivers()) {
            for (DelegatingDriver driver : propagatorBuilder.getPropagationParametersDrivers().getDrivers()) {
                if (driver.getName().equals(refDriver.getName())) {
                    driver.setValue(refDriver.getValue());
                }
            }
        }
    }
    // Build the reference propagator
    final NumericalPropagator propagator = propagatorBuilder.buildPropagator(propagatorBuilder.getSelectedNormalizedParameters());
    // Propagate until last date and return the orbit
    return propagator.propagate(kalmanFinalDate).getOrbit();
}

Example 58 with NumericalPropagator

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

the class ConstantThrustManeuverTest method testRoughBehaviour.

@Test
public void testRoughBehaviour() throws OrekitException {
    final double isp = 318;
    final double mass = 2500;
    final double a = 24396159;
    final double e = 0.72831215;
    final double i = FastMath.toRadians(7);
    final double omega = FastMath.toRadians(180);
    final double OMEGA = FastMath.toRadians(261);
    final double lv = 0;
    final double duration = 3653.99;
    final double f = 420;
    final double delta = FastMath.toRadians(-7.4978);
    final double alpha = FastMath.toRadians(351);
    final AttitudeProvider law = new InertialProvider(new Rotation(new Vector3D(alpha, delta), Vector3D.PLUS_I));
    final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2004, 01, 01), new TimeComponents(23, 30, 00.000), TimeScalesFactory.getUTC());
    final Orbit orbit = new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngle.TRUE, FramesFactory.getEME2000(), initDate, mu);
    final SpacecraftState initialState = new SpacecraftState(orbit, law.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), mass);
    final AbsoluteDate fireDate = new AbsoluteDate(new DateComponents(2004, 01, 02), new TimeComponents(04, 15, 34.080), TimeScalesFactory.getUTC());
    final ConstantThrustManeuver maneuver = new ConstantThrustManeuver(fireDate, duration, f, isp, Vector3D.PLUS_I);
    Assert.assertEquals(f, maneuver.getThrust(), 1.0e-10);
    Assert.assertEquals(isp, maneuver.getISP(), 1.0e-10);
    double[] absTolerance = { 0.001, 1.0e-9, 1.0e-9, 1.0e-6, 1.0e-6, 1.0e-6, 0.001 };
    double[] relTolerance = { 1.0e-7, 1.0e-4, 1.0e-4, 1.0e-7, 1.0e-7, 1.0e-7, 1.0e-7 };
    AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 1000, absTolerance, relTolerance);
    integrator.setInitialStepSize(60);
    final NumericalPropagator propagator = new NumericalPropagator(integrator);
    propagator.setInitialState(initialState);
    propagator.setAttitudeProvider(law);
    propagator.addForceModel(maneuver);
    final SpacecraftState finalorb = propagator.propagate(fireDate.shiftedBy(3800));
    final double massTolerance = FastMath.abs(maneuver.getFlowRate()) * maneuver.getEventsDetectors().findFirst().get().getThreshold();
    Assert.assertEquals(2007.8824544261233, finalorb.getMass(), massTolerance);
    Assert.assertEquals(2.6872, FastMath.toDegrees(MathUtils.normalizeAngle(finalorb.getI(), FastMath.PI)), 1e-4);
    Assert.assertEquals(28970, finalorb.getA() / 1000, 1);
}

Example 59 with NumericalPropagator

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

the class ConstantThrustManeuverTest method testInertialManeuver.

@Test
public void testInertialManeuver() throws OrekitException {
    final double isp = 318;
    final double mass = 2500;
    final double a = 24396159;
    final double e = 0.72831215;
    final double i = FastMath.toRadians(7);
    final double omega = FastMath.toRadians(180);
    final double OMEGA = FastMath.toRadians(261);
    final double lv = 0;
    final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2004, 01, 01), new TimeComponents(23, 30, 00.000), TimeScalesFactory.getUTC());
    final Orbit orbit = new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngle.TRUE, FramesFactory.getEME2000(), initDate, mu);
    final double duration = 3653.99;
    final double f = 420;
    final double delta = FastMath.toRadians(-7.4978);
    final double alpha = FastMath.toRadians(351);
    final AttitudeProvider inertialLaw = new InertialProvider(new Rotation(new Vector3D(alpha, delta), Vector3D.PLUS_I));
    final AttitudeProvider lofLaw = new LofOffset(orbit.getFrame(), LOFType.VNC);
    final SpacecraftState initialState = new SpacecraftState(orbit, inertialLaw.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), mass);
    final AbsoluteDate fireDate = new AbsoluteDate(new DateComponents(2004, 01, 02), new TimeComponents(04, 15, 34.080), TimeScalesFactory.getUTC());
    final ConstantThrustManeuver maneuverWithoutOverride = new ConstantThrustManeuver(fireDate, duration, f, isp, Vector3D.PLUS_I);
    Assert.assertEquals(f, maneuverWithoutOverride.getThrust(), 1.0e-10);
    Assert.assertEquals(isp, maneuverWithoutOverride.getISP(), 1.0e-10);
    // reference propagation:
    // propagator already uses inertial law
    // maneuver does not need to override it to get an inertial maneuver
    double[][] tol = NumericalPropagator.tolerances(1.0, orbit, OrbitType.KEPLERIAN);
    AdaptiveStepsizeIntegrator integrator1 = new DormandPrince853Integrator(0.001, 1000, tol[0], tol[1]);
    integrator1.setInitialStepSize(60);
    final NumericalPropagator propagator1 = new NumericalPropagator(integrator1);
    propagator1.setInitialState(initialState);
    propagator1.setAttitudeProvider(inertialLaw);
    propagator1.addForceModel(maneuverWithoutOverride);
    final SpacecraftState finalState1 = propagator1.propagate(fireDate.shiftedBy(3800));
    // test propagation:
    // propagator uses a LOF-aligned law
    // maneuver needs to override it to get an inertial maneuver
    final ConstantThrustManeuver maneuverWithOverride = new ConstantThrustManeuver(fireDate, duration, f, isp, inertialLaw, Vector3D.PLUS_I);
    Assert.assertEquals(f, maneuverWithoutOverride.getThrust(), 1.0e-10);
    Assert.assertEquals(isp, maneuverWithoutOverride.getISP(), 1.0e-10);
    AdaptiveStepsizeIntegrator integrator2 = new DormandPrince853Integrator(0.001, 1000, tol[0], tol[1]);
    integrator2.setInitialStepSize(60);
    final NumericalPropagator propagator2 = new NumericalPropagator(integrator2);
    propagator2.setInitialState(initialState);
    propagator2.setAttitudeProvider(lofLaw);
    propagator2.addForceModel(maneuverWithOverride);
    final SpacecraftState finalState2 = propagator2.propagate(finalState1.getDate());
    Assert.assertThat(finalState2.getPVCoordinates(), OrekitMatchers.pvCloseTo(finalState1.getPVCoordinates(), 1.0e-10));
    // intentionally wrong propagation, that will produce a very different state
    // propagator uses LOF attitude,
    // maneuver forget to override it, so maneuver will be LOF-aligned in this case
    AdaptiveStepsizeIntegrator integrator3 = new DormandPrince853Integrator(0.001, 1000, tol[0], tol[1]);
    integrator3.setInitialStepSize(60);
    final NumericalPropagator propagator3 = new NumericalPropagator(integrator3);
    propagator3.setInitialState(initialState);
    propagator3.setAttitudeProvider(lofLaw);
    propagator3.addForceModel(maneuverWithoutOverride);
    final SpacecraftState finalState3 = propagator3.propagate(finalState1.getDate());
    Assert.assertEquals(345859.0, Vector3D.distance(finalState1.getPVCoordinates().getPosition(), finalState3.getPVCoordinates().getPosition()), 1.0);
}

Example 60 with NumericalPropagator

use of org.orekit.propagation.numerical.NumericalPropagator 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);
            }
        }
    }
}