Examples with EcksteinHechlerPropagator org.orekit.propagation.analytical.EcksteinHechlerPropagator used on opensource projects

Example 26 with EcksteinHechlerPropagator

use of org.orekit.propagation.analytical.EcksteinHechlerPropagator in project Orekit by CS-SI.

the class AttitudeTest method testInterpolation.

@Test
public void testInterpolation() throws OrekitException {
    Utils.setDataRoot("regular-data");
    final double ehMu = 3.9860047e14;
    final double ae = 6.378137e6;
    final double c20 = -1.08263e-3;
    final double c30 = 2.54e-6;
    final double c40 = 1.62e-6;
    final double c50 = 2.3e-7;
    final double c60 = -5.5e-7;
    final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
    final Vector3D position = new Vector3D(3220103., 69623., 6449822.);
    final Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
    final CircularOrbit initialOrbit = new CircularOrbit(new PVCoordinates(position, velocity), FramesFactory.getEME2000(), date, ehMu);
    EcksteinHechlerPropagator propagator = new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
    OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
    propagator.setAttitudeProvider(new BodyCenterPointing(initialOrbit.getFrame(), earth));
    final Attitude initialAttitude = propagator.propagate(initialOrbit.getDate()).getAttitude();
    // set up a 5 points sample
    List<Attitude> sample = new ArrayList<Attitude>();
    for (double dt = 0; dt < 251.0; dt += 60.0) {
        sample.add(propagator.propagate(date.shiftedBy(dt)).getAttitude());
    }
    // well inside the sample, interpolation should be better than quadratic shift
    double maxShiftAngleError = 0;
    double maxInterpolationAngleError = 0;
    double maxShiftRateError = 0;
    double maxInterpolationRateError = 0;
    for (double dt = 0; dt < 240.0; dt += 1.0) {
        AbsoluteDate t = initialOrbit.getDate().shiftedBy(dt);
        Attitude propagated = propagator.propagate(t).getAttitude();
        double shiftAngleError = Rotation.distance(propagated.getRotation(), initialAttitude.shiftedBy(dt).getRotation());
        double interpolationAngleError = Rotation.distance(propagated.getRotation(), initialAttitude.interpolate(t, sample).getRotation());
        double shiftRateError = Vector3D.distance(propagated.getSpin(), initialAttitude.shiftedBy(dt).getSpin());
        double interpolationRateError = Vector3D.distance(propagated.getSpin(), initialAttitude.interpolate(t, sample).getSpin());
        maxShiftAngleError = FastMath.max(maxShiftAngleError, shiftAngleError);
        maxInterpolationAngleError = FastMath.max(maxInterpolationAngleError, interpolationAngleError);
        maxShiftRateError = FastMath.max(maxShiftRateError, shiftRateError);
        maxInterpolationRateError = FastMath.max(maxInterpolationRateError, interpolationRateError);
    }
    Assert.assertTrue(maxShiftAngleError > 4.0e-6);
    Assert.assertTrue(maxInterpolationAngleError < 1.5e-13);
    Assert.assertTrue(maxShiftRateError > 6.0e-8);
    Assert.assertTrue(maxInterpolationRateError < 2.5e-14);
    // past sample end, interpolation error should increase, but still be far better than quadratic shift
    maxShiftAngleError = 0;
    maxInterpolationAngleError = 0;
    maxShiftRateError = 0;
    maxInterpolationRateError = 0;
    for (double dt = 250.0; dt < 300.0; dt += 1.0) {
        AbsoluteDate t = initialOrbit.getDate().shiftedBy(dt);
        Attitude propagated = propagator.propagate(t).getAttitude();
        double shiftAngleError = Rotation.distance(propagated.getRotation(), initialAttitude.shiftedBy(dt).getRotation());
        double interpolationAngleError = Rotation.distance(propagated.getRotation(), initialAttitude.interpolate(t, sample).getRotation());
        double shiftRateError = Vector3D.distance(propagated.getSpin(), initialAttitude.shiftedBy(dt).getSpin());
        double interpolationRateError = Vector3D.distance(propagated.getSpin(), initialAttitude.interpolate(t, sample).getSpin());
        maxShiftAngleError = FastMath.max(maxShiftAngleError, shiftAngleError);
        maxInterpolationAngleError = FastMath.max(maxInterpolationAngleError, interpolationAngleError);
        maxShiftRateError = FastMath.max(maxShiftRateError, shiftRateError);
        maxInterpolationRateError = FastMath.max(maxInterpolationRateError, interpolationRateError);
    }
    Assert.assertTrue(maxShiftAngleError > 9.0e-6);
    Assert.assertTrue(maxInterpolationAngleError < 6.0e-11);
    Assert.assertTrue(maxShiftRateError > 9.0e-8);
    Assert.assertTrue(maxInterpolationRateError < 4.0e-12);
}

Example 27 with EcksteinHechlerPropagator

use of org.orekit.propagation.analytical.EcksteinHechlerPropagator in project Orekit by CS-SI.

the class BodyCenterPointingTest method testQDot.

@Test
public void testQDot() throws OrekitException {
    Utils.setDataRoot("regular-data");
    final double ehMu = 3.9860047e14;
    final double ae = 6.378137e6;
    final double c20 = -1.08263e-3;
    final double c30 = 2.54e-6;
    final double c40 = 1.62e-6;
    final double c50 = 2.3e-7;
    final double c60 = -5.5e-7;
    final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
    final Vector3D position = new Vector3D(3220103., 69623., 6449822.);
    final Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
    final CircularOrbit initialOrbit = new CircularOrbit(new PVCoordinates(position, velocity), FramesFactory.getEME2000(), date, ehMu);
    EcksteinHechlerPropagator propagator = new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
    propagator.setAttitudeProvider(earthCenterAttitudeLaw);
    List<WeightedObservedPoint> w0 = new ArrayList<WeightedObservedPoint>();
    List<WeightedObservedPoint> w1 = new ArrayList<WeightedObservedPoint>();
    List<WeightedObservedPoint> w2 = new ArrayList<WeightedObservedPoint>();
    List<WeightedObservedPoint> w3 = new ArrayList<WeightedObservedPoint>();
    for (double dt = -1; dt < 1; dt += 0.01) {
        Rotation rP = propagator.propagate(date.shiftedBy(dt)).getAttitude().getRotation();
        w0.add(new WeightedObservedPoint(1, dt, rP.getQ0()));
        w1.add(new WeightedObservedPoint(1, dt, rP.getQ1()));
        w2.add(new WeightedObservedPoint(1, dt, rP.getQ2()));
        w3.add(new WeightedObservedPoint(1, dt, rP.getQ3()));
    }
    double q0DotRef = PolynomialCurveFitter.create(2).fit(w0)[1];
    double q1DotRef = PolynomialCurveFitter.create(2).fit(w1)[1];
    double q2DotRef = PolynomialCurveFitter.create(2).fit(w2)[1];
    double q3DotRef = PolynomialCurveFitter.create(2).fit(w3)[1];
    Attitude a0 = propagator.propagate(date).getAttitude();
    double q0 = a0.getRotation().getQ0();
    double q1 = a0.getRotation().getQ1();
    double q2 = a0.getRotation().getQ2();
    double q3 = a0.getRotation().getQ3();
    double oX = a0.getSpin().getX();
    double oY = a0.getSpin().getY();
    double oZ = a0.getSpin().getZ();
    // first time-derivatives of the quaternion
    double q0Dot = 0.5 * MathArrays.linearCombination(-q1, oX, -q2, oY, -q3, oZ);
    double q1Dot = 0.5 * MathArrays.linearCombination(q0, oX, -q3, oY, q2, oZ);
    double q2Dot = 0.5 * MathArrays.linearCombination(q3, oX, q0, oY, -q1, oZ);
    double q3Dot = 0.5 * MathArrays.linearCombination(-q2, oX, q1, oY, q0, oZ);
    Assert.assertEquals(q0DotRef, q0Dot, 5.0e-9);
    Assert.assertEquals(q1DotRef, q1Dot, 5.0e-9);
    Assert.assertEquals(q2DotRef, q2Dot, 5.0e-9);
    Assert.assertEquals(q3DotRef, q3Dot, 5.0e-9);
}

Example 28 with EcksteinHechlerPropagator

use of org.orekit.propagation.analytical.EcksteinHechlerPropagator in project Orekit by CS-SI.

the class BodyCenterPointingTest method testSpin.

@Test
public void testSpin() throws OrekitException {
    Utils.setDataRoot("regular-data");
    final double ehMu = 3.9860047e14;
    final double ae = 6.378137e6;
    final double c20 = -1.08263e-3;
    final double c30 = 2.54e-6;
    final double c40 = 1.62e-6;
    final double c50 = 2.3e-7;
    final double c60 = -5.5e-7;
    final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
    final Vector3D position = new Vector3D(3220103., 69623., 6449822.);
    final Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
    final CircularOrbit initialOrbit = new CircularOrbit(new PVCoordinates(position, velocity), FramesFactory.getEME2000(), date, ehMu);
    EcksteinHechlerPropagator propagator = new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
    propagator.setAttitudeProvider(earthCenterAttitudeLaw);
    double h = 0.01;
    SpacecraftState s0 = propagator.propagate(date);
    SpacecraftState sMinus = propagator.propagate(date.shiftedBy(-h));
    SpacecraftState sPlus = propagator.propagate(date.shiftedBy(h));
    // check spin is consistent with attitude evolution
    double errorAngleMinus = Rotation.distance(sMinus.shiftedBy(h).getAttitude().getRotation(), s0.getAttitude().getRotation());
    double evolutionAngleMinus = Rotation.distance(sMinus.getAttitude().getRotation(), s0.getAttitude().getRotation());
    Assert.assertEquals(0.0, errorAngleMinus, 1.0e-6 * evolutionAngleMinus);
    double errorAnglePlus = Rotation.distance(s0.getAttitude().getRotation(), sPlus.shiftedBy(-h).getAttitude().getRotation());
    double evolutionAnglePlus = Rotation.distance(s0.getAttitude().getRotation(), sPlus.getAttitude().getRotation());
    Assert.assertEquals(0.0, errorAnglePlus, 1.0e-6 * evolutionAnglePlus);
    Vector3D spin0 = s0.getAttitude().getSpin();
    Vector3D reference = AngularCoordinates.estimateRate(sMinus.getAttitude().getRotation(), sPlus.getAttitude().getRotation(), 2 * h);
    Assert.assertTrue(spin0.getNorm() > 1.0e-3);
    Assert.assertEquals(0.0, spin0.subtract(reference).getNorm(), 1.0e-13);
}

Example 29 with EcksteinHechlerPropagator

use of org.orekit.propagation.analytical.EcksteinHechlerPropagator in project Orekit by CS-SI.

the class AttitudesSequenceTest method testDayNightSwitch.

@Test
public void testDayNightSwitch() throws OrekitException {
    // Initial state definition : date, orbit
    final AbsoluteDate initialDate = new AbsoluteDate(2004, 01, 01, 23, 30, 00.000, TimeScalesFactory.getUTC());
    final Vector3D position = new Vector3D(-6142438.668, 3492467.560, -25767.25680);
    final Vector3D velocity = new Vector3D(505.8479685, 942.7809215, 7435.922231);
    final Orbit initialOrbit = new KeplerianOrbit(new PVCoordinates(position, velocity), FramesFactory.getEME2000(), initialDate, Constants.EIGEN5C_EARTH_MU);
    final EventsLogger // Attitudes sequence definition
    logger = new EventsLogger();
    final AttitudesSequence attitudesSequence = new AttitudesSequence();
    final AttitudeProvider dayObservationLaw = new LofOffset(initialOrbit.getFrame(), LOFType.VVLH, RotationOrder.XYZ, FastMath.toRadians(20), FastMath.toRadians(40), 0);
    final AttitudeProvider nightRestingLaw = new LofOffset(initialOrbit.getFrame(), LOFType.VVLH);
    final PVCoordinatesProvider sun = CelestialBodyFactory.getSun();
    final PVCoordinatesProvider earth = CelestialBodyFactory.getEarth();
    final EclipseDetector ed = new EclipseDetector(sun, 696000000., earth, Constants.WGS84_EARTH_EQUATORIAL_RADIUS).withHandler(new ContinueOnEvent<EclipseDetector>() {

        private static final long serialVersionUID = 1L;

        public EventHandler.Action eventOccurred(final SpacecraftState s, final EclipseDetector d, final boolean increasing) {
            setInEclipse(s.getDate(), !increasing);
            return EventHandler.Action.RESET_STATE;
        }
    });
    final EventDetector monitored = logger.monitorDetector(ed);
    final Handler dayToNightHandler = new Handler(dayObservationLaw, nightRestingLaw);
    final Handler nightToDayHandler = new Handler(nightRestingLaw, dayObservationLaw);
    attitudesSequence.addSwitchingCondition(dayObservationLaw, nightRestingLaw, monitored, false, true, 300.0, AngularDerivativesFilter.USE_RRA, dayToNightHandler);
    attitudesSequence.addSwitchingCondition(nightRestingLaw, dayObservationLaw, monitored, true, false, 300.0, AngularDerivativesFilter.USE_RRA, nightToDayHandler);
    SpacecraftState initialState = new SpacecraftState(initialOrbit);
    initialState = initialState.addAdditionalState("fortyTwo", 42.0);
    if (ed.g(initialState) >= 0) {
        // initial position is in daytime
        setInEclipse(initialDate, false);
        attitudesSequence.resetActiveProvider(dayObservationLaw);
    } else {
        // initial position is in nighttime
        setInEclipse(initialDate, true);
        attitudesSequence.resetActiveProvider(nightRestingLaw);
    }
    // Propagator : consider the analytical Eckstein-Hechler model
    final Propagator propagator = new EcksteinHechlerPropagator(initialOrbit, attitudesSequence, Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS, Constants.EIGEN5C_EARTH_MU, Constants.EIGEN5C_EARTH_C20, Constants.EIGEN5C_EARTH_C30, Constants.EIGEN5C_EARTH_C40, Constants.EIGEN5C_EARTH_C50, Constants.EIGEN5C_EARTH_C60);
    // Register the switching events to the propagator
    attitudesSequence.registerSwitchEvents(propagator);
    propagator.setMasterMode(60.0, new OrekitFixedStepHandler() {

        public void handleStep(SpacecraftState currentState, boolean isLast) throws OrekitException {
            // the Earth position in spacecraft frame should be along spacecraft Z axis
            // during night time and away from it during day time due to roll and pitch offsets
            final Vector3D earth = currentState.toTransform().transformPosition(Vector3D.ZERO);
            final double pointingOffset = Vector3D.angle(earth, Vector3D.PLUS_K);
            // the g function is the eclipse indicator, its an angle between Sun and Earth limb,
            // positive when Sun is outside of Earth limb, negative when Sun is hidden by Earth limb
            final double eclipseAngle = ed.g(currentState);
            if (currentState.getDate().durationFrom(lastChange) > 300) {
                if (inEclipse) {
                    Assert.assertTrue(eclipseAngle <= 0);
                    Assert.assertEquals(0.0, pointingOffset, 1.0e-6);
                } else {
                    Assert.assertTrue(eclipseAngle >= 0);
                    Assert.assertEquals(0.767215, pointingOffset, 1.0e-6);
                }
            } else {
                // we are in transition
                Assert.assertTrue(pointingOffset + " " + (0.767215 - pointingOffset), pointingOffset <= 0.7672155);
            }
        }
    });
    // Propagate from the initial date for the fixed duration
    propagator.propagate(initialDate.shiftedBy(12600.));
    // as we have 2 switch events (even if they share the same underlying event detector),
    // and these events are triggered at both eclipse entry and exit, we get 8
    // raw events on 2 orbits
    Assert.assertEquals(8, logger.getLoggedEvents().size());
    // we have 4 attitudes switch on 2 orbits, 2 of each type
    Assert.assertEquals(2, dayToNightHandler.dates.size());
    Assert.assertEquals(2, nightToDayHandler.dates.size());
}

Example 30 with EcksteinHechlerPropagator

use of org.orekit.propagation.analytical.EcksteinHechlerPropagator in project Orekit by CS-SI.

the class CartesianOrbitTest method doTestInterpolation.

private void doTestInterpolation(boolean useDerivatives, double shiftPositionErrorWithin, double interpolationPositionErrorWithin, double shiftVelocityErrorWithin, double interpolationVelocityErrorWithin, double shiftPositionErrorFarPast, double interpolationPositionErrorFarPast, double shiftVelocityErrorFarPast, double interpolationVelocityErrorFarPast) throws OrekitException {
    final double ehMu = 3.9860047e14;
    final double ae = 6.378137e6;
    final double c20 = -1.08263e-3;
    final double c30 = 2.54e-6;
    final double c40 = 1.62e-6;
    final double c50 = 2.3e-7;
    final double c60 = -5.5e-7;
    final AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
    final Vector3D position = new Vector3D(3220103., 69623., 6449822.);
    final Vector3D velocity = new Vector3D(6414.7, -2006., -3180.);
    final CartesianOrbit initialOrbit = new CartesianOrbit(new PVCoordinates(position, velocity), FramesFactory.getEME2000(), date, ehMu);
    EcksteinHechlerPropagator propagator = new EcksteinHechlerPropagator(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
    // set up a 5 points sample
    List<Orbit> sample = new ArrayList<Orbit>();
    for (double dt = 0; dt < 251.0; dt += 60.0) {
        Orbit orbit = propagator.propagate(date.shiftedBy(dt)).getOrbit();
        if (!useDerivatives) {
            // remove derivatives
            double[] stateVector = new double[6];
            orbit.getType().mapOrbitToArray(orbit, PositionAngle.TRUE, stateVector, null);
            orbit = orbit.getType().mapArrayToOrbit(stateVector, null, PositionAngle.TRUE, orbit.getDate(), orbit.getMu(), orbit.getFrame());
        }
        sample.add(orbit);
    }
    // well inside the sample, interpolation should be much better than Keplerian shift
    // this is because we take the full non-Keplerian acceleration into account in
    // the Cartesian parameters, which in this case is preserved by the
    // Eckstein-Hechler propagator
    double maxShiftPError = 0;
    double maxInterpolationPError = 0;
    double maxShiftVError = 0;
    double maxInterpolationVError = 0;
    for (double dt = 0; dt < 240.0; dt += 1.0) {
        AbsoluteDate t = initialOrbit.getDate().shiftedBy(dt);
        PVCoordinates propagated = propagator.propagate(t).getPVCoordinates();
        PVCoordinates shiftError = new PVCoordinates(propagated, initialOrbit.shiftedBy(dt).getPVCoordinates());
        PVCoordinates interpolationError = new PVCoordinates(propagated, initialOrbit.interpolate(t, sample).getPVCoordinates());
        maxShiftPError = FastMath.max(maxShiftPError, shiftError.getPosition().getNorm());
        maxInterpolationPError = FastMath.max(maxInterpolationPError, interpolationError.getPosition().getNorm());
        maxShiftVError = FastMath.max(maxShiftVError, shiftError.getVelocity().getNorm());
        maxInterpolationVError = FastMath.max(maxInterpolationVError, interpolationError.getVelocity().getNorm());
    }
    Assert.assertEquals(shiftPositionErrorWithin, maxShiftPError, 0.01 * shiftPositionErrorWithin);
    Assert.assertEquals(interpolationPositionErrorWithin, maxInterpolationPError, 0.01 * interpolationPositionErrorWithin);
    Assert.assertEquals(shiftVelocityErrorWithin, maxShiftVError, 0.01 * shiftVelocityErrorWithin);
    Assert.assertEquals(interpolationVelocityErrorWithin, maxInterpolationVError, 0.01 * interpolationVelocityErrorWithin);
    // if we go far past sample end, interpolation becomes worse than Keplerian shift
    maxShiftPError = 0;
    maxInterpolationPError = 0;
    maxShiftVError = 0;
    maxInterpolationVError = 0;
    for (double dt = 500.0; dt < 650.0; dt += 1.0) {
        AbsoluteDate t = initialOrbit.getDate().shiftedBy(dt);
        PVCoordinates propagated = propagator.propagate(t).getPVCoordinates();
        PVCoordinates shiftError = new PVCoordinates(propagated, initialOrbit.shiftedBy(dt).getPVCoordinates());
        PVCoordinates interpolationError = new PVCoordinates(propagated, initialOrbit.interpolate(t, sample).getPVCoordinates());
        maxShiftPError = FastMath.max(maxShiftPError, shiftError.getPosition().getNorm());
        maxInterpolationPError = FastMath.max(maxInterpolationPError, interpolationError.getPosition().getNorm());
        maxShiftVError = FastMath.max(maxShiftVError, shiftError.getVelocity().getNorm());
        maxInterpolationVError = FastMath.max(maxInterpolationVError, interpolationError.getVelocity().getNorm());
    }
    Assert.assertEquals(shiftPositionErrorFarPast, maxShiftPError, 0.01 * shiftPositionErrorFarPast);
    Assert.assertEquals(interpolationPositionErrorFarPast, maxInterpolationPError, 0.01 * interpolationPositionErrorFarPast);
    Assert.assertEquals(shiftVelocityErrorFarPast, maxShiftVError, 0.01 * shiftVelocityErrorFarPast);
    Assert.assertEquals(interpolationVelocityErrorFarPast, maxInterpolationVError, 0.01 * interpolationVelocityErrorFarPast);
}