Examples with PositionAngle org.orekit.orbits.PositionAngle used on opensource projects

Example 11 with PositionAngle

use of org.orekit.orbits.PositionAngle 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 12 with PositionAngle

use of org.orekit.orbits.PositionAngle in project Orekit by CS-SI.

the class KalmanEstimatorTest method testWrappedException.

/**
 * Test of a wrapped exception in a Kalman observer
 * @throws OrekitException
 */
@Test
public void testWrappedException() 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, 3.0, 300.0);
    // Build the Kalman filter
    final KalmanEstimatorBuilder kalmanBuilder = new KalmanEstimatorBuilder();
    kalmanBuilder.builder(propagatorBuilder);
    kalmanBuilder.estimatedMeasurementsParameters(new ParameterDriversList());
    kalmanBuilder.initialCovarianceMatrix(MatrixUtils.createRealMatrix(6, 6));
    kalmanBuilder.processNoiseMatrixProvider(new ConstantProcessNoise(MatrixUtils.createRealMatrix(6, 6)));
    final KalmanEstimator kalman = kalmanBuilder.build();
    kalman.setObserver(estimation -> {
        throw new DummyException();
    });
    try {
        // Filter the measurements and expect an exception to occur
        EstimationTestUtils.checkKalmanFit(context, kalman, measurements, context.initialOrbit, positionAngle, 0., 0., 0., 0., new double[3], 0., new double[3], 0.);
    } catch (DummyException de) {
    // expected
    }
}

Example 13 with PositionAngle

use of org.orekit.orbits.PositionAngle in project Orekit by CS-SI.

the class SmallManeuverAnalyticalModelTest method testJacobian.

@Test
public void testJacobian() throws OrekitException {
    Frame eme2000 = FramesFactory.getEME2000();
    Orbit leo = new CircularOrbit(7200000.0, -1.0e-2, 2.0e-3, FastMath.toRadians(98.0), FastMath.toRadians(123.456), 0.3, PositionAngle.MEAN, eme2000, new AbsoluteDate(new DateComponents(2004, 01, 01), new TimeComponents(23, 30, 00.000), TimeScalesFactory.getUTC()), Constants.EIGEN5C_EARTH_MU);
    double mass = 5600.0;
    AbsoluteDate t0 = leo.getDate().shiftedBy(1000.0);
    Vector3D dV0 = new Vector3D(-0.1, 0.2, 0.3);
    double f = 400.0;
    double isp = 315.0;
    for (OrbitType orbitType : OrbitType.values()) {
        for (PositionAngle positionAngle : PositionAngle.values()) {
            BoundedPropagator withoutManeuver = getEphemeris(orbitType.convertType(leo), mass, t0, Vector3D.ZERO, f, isp);
            SpacecraftState state0 = withoutManeuver.propagate(t0);
            SmallManeuverAnalyticalModel model = new SmallManeuverAnalyticalModel(state0, eme2000, dV0, isp);
            Assert.assertEquals(t0, model.getDate());
            Vector3D[] velDirs = new Vector3D[] { Vector3D.PLUS_I, Vector3D.PLUS_J, Vector3D.PLUS_K, Vector3D.ZERO };
            double[] timeDirs = new double[] { 0, 0, 0, 1 };
            double h = 1.0;
            AbsoluteDate t1 = t0.shiftedBy(20.0);
            for (int i = 0; i < 4; ++i) {
                SmallManeuverAnalyticalModel[] models = new SmallManeuverAnalyticalModel[] { new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(-4 * h * timeDirs[i])), eme2000, new Vector3D(1, dV0, -4 * h, velDirs[i]), isp), new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(-3 * h * timeDirs[i])), eme2000, new Vector3D(1, dV0, -3 * h, velDirs[i]), isp), new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(-2 * h * timeDirs[i])), eme2000, new Vector3D(1, dV0, -2 * h, velDirs[i]), isp), new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(-1 * h * timeDirs[i])), eme2000, new Vector3D(1, dV0, -1 * h, velDirs[i]), isp), new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(+1 * h * timeDirs[i])), eme2000, new Vector3D(1, dV0, +1 * h, velDirs[i]), isp), new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(+2 * h * timeDirs[i])), eme2000, new Vector3D(1, dV0, +2 * h, velDirs[i]), isp), new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(+3 * h * timeDirs[i])), eme2000, new Vector3D(1, dV0, +3 * h, velDirs[i]), isp), new SmallManeuverAnalyticalModel(withoutManeuver.propagate(t0.shiftedBy(+4 * h * timeDirs[i])), eme2000, new Vector3D(1, dV0, +4 * h, velDirs[i]), isp) };
                double[][] array = new double[models.length][6];
                Orbit orbitWithout = withoutManeuver.propagate(t1).getOrbit();
                // compute reference orbit gradient by finite differences
                double c = 1.0 / (840 * h);
                for (int j = 0; j < models.length; ++j) {
                    orbitType.mapOrbitToArray(models[j].apply(orbitWithout), positionAngle, array[j], null);
                }
                double[] orbitGradient = new double[6];
                for (int k = 0; k < orbitGradient.length; ++k) {
                    double d4 = array[7][k] - array[0][k];
                    double d3 = array[6][k] - array[1][k];
                    double d2 = array[5][k] - array[2][k];
                    double d1 = array[4][k] - array[3][k];
                    orbitGradient[k] = (-3 * d4 + 32 * d3 - 168 * d2 + 672 * d1) * c;
                }
                // analytical Jacobian to check
                double[][] jacobian = new double[6][4];
                model.getJacobian(orbitWithout, positionAngle, jacobian);
                for (int j = 0; j < orbitGradient.length; ++j) {
                    Assert.assertEquals(orbitGradient[j], jacobian[j][i], 1.6e-4 * FastMath.abs(orbitGradient[j]));
                }
            }
        }
    }
}

Example 14 with PositionAngle

use of org.orekit.orbits.PositionAngle in project Orekit by CS-SI.

the class KalmanEstimatorTest method testEquinoctialRightAscensionDeclination.

/**
 * Perfect right-ascension/declination measurements with a perfect start
 * Equinoctial formalism
 * @throws OrekitException
 */
@Test
public void testEquinoctialRightAscensionDeclination() throws OrekitException {
    // Create context
    Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    // Create initial orbit and propagator builder
    final OrbitType orbitType = OrbitType.EQUINOCTIAL;
    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 AngularRaDecMeasurementCreator(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();
    // Cartesian covariance matrix initialization
    final RealMatrix cartesianP = MatrixUtils.createRealDiagonalMatrix(new double[] { 1e-4, 1e-4, 1e-4, 1e-10, 1e-10, 1e-10 });
    // 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, dYdC);
    final RealMatrix Jac = MatrixUtils.createRealMatrix(dYdC);
    // Keplerian initial covariance matrix
    final RealMatrix initialP = Jac.multiply(cartesianP.multiply(Jac.transpose()));
    // Process noise matrix
    final RealMatrix cartesianQ = MatrixUtils.createRealDiagonalMatrix(new double[] { 1.e-6, 1.e-6, 1.e-6, 1.e-12, 1.e-12, 1.e-12 });
    final RealMatrix Q = Jac.multiply(cartesianQ.multiply(Jac.transpose()));
    // 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.53e-5;
    final double expectedDeltaVel = 0.;
    final double velEps = 5.04e-9;
    final double[] expectedSigmasPos = { 0.356902, 1.297507, 1.798551 };
    final double sigmaPosEps = 1e-6;
    final double[] expectedSigmasVel = { 2.468745e-4, 5.810027e-4, 3.887394e-4 };
    final double sigmaVelEps = 1e-10;
    EstimationTestUtils.checkKalmanFit(context, kalman, measurements, refOrbit, positionAngle, expectedDeltaPos, posEps, expectedDeltaVel, velEps, expectedSigmasPos, sigmaPosEps, expectedSigmasVel, sigmaVelEps);
}

Example 15 with PositionAngle

use of org.orekit.orbits.PositionAngle 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));
            }
        }
    }
}