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

Example 1 with FieldEcksteinHechlerPropagator

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

the class FieldSpacecraftStateTest method doTestShiftVsEcksteinHechlerError.

private <T extends RealFieldElement<T>> void doTestShiftVsEcksteinHechlerError(final Field<T> field) throws OrekitException {
    T zero = field.getZero();
    T mass = zero.add(2500.);
    T a = zero.add(rOrbit.getA());
    T e = zero.add(rOrbit.getE());
    T i = zero.add(rOrbit.getI());
    T pa = zero.add(1.9674147913622104);
    T raan = zero.add(FastMath.toRadians(261));
    T lv = zero.add(0);
    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;
    // polynomial models for interpolation error in position, velocity, acceleration and attitude
    // these models grow as follows
    // interpolation time (s)    position error (m)   velocity error (m/s)   acceleration error (m/s²)  attitude error (°)
    // 60                        2                    0.07                  0.002               0.00002
    // 120                       12                    0.3                   0.005               0.00009
    // 300                      170                    1.6                   0.012               0.0009
    // 600                     1200                    5.7                   0.024               0.006
    // 900                     3600                   10.6                   0.034               0.02
    // the expected maximum residuals with respect to these models are about 0.4m, 0.5mm/s, 8μm/s² and 3e-6°
    PolynomialFunction pModel = new PolynomialFunction(new double[] { 1.5664070631933846e-01, 7.5504722733047560e-03, -8.2460562451009510e-05, 6.9546332080305580e-06, -1.7045365367533077e-09, -4.2187860791066264e-13 });
    PolynomialFunction vModel = new PolynomialFunction(new double[] { -3.5472364019908720e-04, 1.6568103861124980e-05, 1.9637913327830596e-05, -3.4248792843039766e-09, -5.6565135131014254e-12, 1.4730170946808630e-15 });
    PolynomialFunction aModel = new PolynomialFunction(new double[] { 3.0731707577766896e-06, 3.9770746399850350e-05, 1.9779039254538660e-09, 8.0263328220724900e-12, -1.5600835252366078e-14, 1.1785257001549687e-18 });
    PolynomialFunction rModel = new PolynomialFunction(new double[] { -2.7689062063188115e-06, 1.7406542538258334e-07, 2.5109795349592287e-09, 2.0399322661074575e-11, 9.9126348912426750e-15, -3.5015638905729510e-18 });
    FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, new DateComponents(2004, 01, 01), TimeComponents.H00, TimeScalesFactory.getUTC());
    FieldKeplerianOrbit<T> orbit = new FieldKeplerianOrbit<>(a, e, i, pa, raan, lv, PositionAngle.TRUE, FramesFactory.getEME2000(), date, mu);
    BodyCenterPointing attitudeLaw = new BodyCenterPointing(orbit.getFrame(), earth);
    FieldPropagator<T> propagator = new FieldEcksteinHechlerPropagator<>(orbit, attitudeLaw, mass, ae, mu, c20, c30, c40, c50, c60);
    FieldAbsoluteDate<T> centerDate = orbit.getDate().shiftedBy(100.0);
    FieldSpacecraftState<T> centerState = propagator.propagate(centerDate);
    double maxResidualP = 0;
    double maxResidualV = 0;
    double maxResidualA = 0;
    double maxResidualR = 0;
    for (T dt = field.getZero(); dt.getReal() < 900.0; dt = dt.add(5)) {
        FieldSpacecraftState<T> shifted = centerState.shiftedBy(dt);
        FieldSpacecraftState<T> propagated = propagator.propagate(centerDate.shiftedBy(dt));
        FieldPVCoordinates<T> dpv = new FieldPVCoordinates<>(propagated.getPVCoordinates(), shifted.getPVCoordinates());
        double residualP = pModel.value(dt.getReal()) - dpv.getPosition().getNorm().getReal();
        double residualV = vModel.value(dt.getReal()) - dpv.getVelocity().getNorm().getReal();
        double residualA = aModel.value(dt.getReal()) - dpv.getAcceleration().getNorm().getReal();
        double residualR = rModel.value(dt.getReal()) - FastMath.toDegrees(FieldRotation.distance(shifted.getAttitude().getRotation(), propagated.getAttitude().getRotation()).getReal());
        maxResidualP = FastMath.max(maxResidualP, FastMath.abs(residualP));
        maxResidualV = FastMath.max(maxResidualV, FastMath.abs(residualV));
        maxResidualA = FastMath.max(maxResidualA, FastMath.abs(residualA));
        maxResidualR = FastMath.max(maxResidualR, FastMath.abs(residualR));
    }
    Assert.assertEquals(0.40, maxResidualP, 0.01);
    Assert.assertEquals(4.9e-4, maxResidualV, 1.0e-5);
    Assert.assertEquals(2.8e-6, maxResidualR, 1.0e-1);
}

Example 2 with FieldEcksteinHechlerPropagator

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

the class BodyCenterPointingTest method doTestSpin.

private <T extends RealFieldElement<T>> void doTestSpin(final Field<T> field) 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;
    // Satellite position as circular parameters
    final T zero = field.getZero();
    final T a = zero.add(7178000.0);
    final T e = zero.add(7e-5);
    final T i = zero.add(FastMath.toRadians(50.));
    final T pa = zero.add(FastMath.toRadians(45.));
    final T raan = zero.add(FastMath.toRadians(270.));
    final T m = zero.add(FastMath.toRadians(5.3 - 270));
    // Computation date
    FieldAbsoluteDate<T> date_R = new FieldAbsoluteDate<>(field, new DateComponents(2008, 04, 07), TimeComponents.H00, TimeScalesFactory.getUTC());
    // Orbit
    FieldKeplerianOrbit<T> circ = new FieldKeplerianOrbit<>(a, e, i, pa, raan, m, PositionAngle.MEAN, FramesFactory.getEME2000(), date_R, ehMu);
    // WGS84 Earth model
    OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
    // Earth center pointing attitude provider
    BodyCenterPointing earthCenterAttitudeLaw = new BodyCenterPointing(circ.getFrame(), earth);
    final FieldAbsoluteDate<T> date = FieldAbsoluteDate.getJ2000Epoch(field).shiftedBy(584.);
    final FieldVector3D<T> position = new FieldVector3D<>(zero.add(3220103.), zero.add(69623.), zero.add(6449822.));
    final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(6414.7), zero.add(-2006.), zero.add(-3180.));
    final FieldCircularOrbit<T> initialOrbit = new FieldCircularOrbit<>(new FieldPVCoordinates<>(position, velocity), FramesFactory.getEME2000(), date, ehMu);
    FieldEcksteinHechlerPropagator<T> propagator = new FieldEcksteinHechlerPropagator<>(initialOrbit, ae, ehMu, c20, c30, c40, c50, c60);
    propagator.setAttitudeProvider(earthCenterAttitudeLaw);
    double h = 0.01;
    FieldSpacecraftState<T> s0 = propagator.propagate(date);
    FieldSpacecraftState<T> sMinus = propagator.propagate(date.shiftedBy(-h));
    FieldSpacecraftState<T> sPlus = propagator.propagate(date.shiftedBy(h));
    // check spin is consistent with attitude evolution
    T errorAngleMinus = FieldRotation.distance(sMinus.shiftedBy(zero.add(h)).getAttitude().getRotation(), s0.getAttitude().getRotation());
    T evolutionAngleMinus = FieldRotation.distance(sMinus.getAttitude().getRotation(), s0.getAttitude().getRotation());
    Assert.assertEquals(0.0, errorAngleMinus.getReal(), 1.0e-6 * evolutionAngleMinus.getReal());
    T errorAnglePlus = FieldRotation.distance(s0.getAttitude().getRotation(), sPlus.shiftedBy(zero.add(-h)).getAttitude().getRotation());
    T evolutionAnglePlus = FieldRotation.distance(s0.getAttitude().getRotation(), sPlus.getAttitude().getRotation());
    Assert.assertEquals(0.0, errorAnglePlus.getReal(), 1.0e-6 * evolutionAnglePlus.getReal());
    FieldVector3D<T> spin0 = s0.getAttitude().getSpin();
    FieldVector3D<T> reference = FieldAngularCoordinates.estimateRate(sMinus.getAttitude().getRotation(), sPlus.getAttitude().getRotation(), 2 * h);
    Assert.assertTrue(spin0.getNorm().getReal() > 1.0e-3);
    Assert.assertEquals(0.0, spin0.subtract(reference).getNorm().getReal(), 1.0e-13);
}

Example 3 with FieldEcksteinHechlerPropagator

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

the class FieldApsideDetectorTest method doTestSimple.

private <T extends RealFieldElement<T>> void doTestSimple(Field<T> field) throws OrekitException {
    final T zero = field.getZero();
    final TimeScale utc = TimeScalesFactory.getUTC();
    final FieldVector3D<T> position = new FieldVector3D<>(zero.add(-6142438.668), zero.add(3492467.56), zero.add(-25767.257));
    final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(506.0), zero.add(943.0), zero.add(7450));
    final FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, 2003, 9, 16, utc);
    final FieldOrbit<T> orbit = new FieldCartesianOrbit<>(new FieldPVCoordinates<>(position, velocity), FramesFactory.getEME2000(), date, Constants.EIGEN5C_EARTH_MU);
    FieldEcksteinHechlerPropagator<T> propagator = new FieldEcksteinHechlerPropagator<>(orbit, 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);
    FieldEventDetector<T> detector = new FieldApsideDetector<>(propagator.getInitialState().getOrbit()).withMaxCheck(zero.add(600.0)).withThreshold(zero.add(1.0e-12)).withHandler(new FieldContinueOnEvent<FieldApsideDetector<T>, T>());
    Assert.assertEquals(600.0, detector.getMaxCheckInterval().getReal(), 1.0e-15);
    Assert.assertEquals(1.0e-12, detector.getThreshold().getReal(), 1.0e-15);
    Assert.assertEquals(AbstractDetector.DEFAULT_MAX_ITER, detector.getMaxIterationCount());
    FieldEventsLogger<T> logger = new FieldEventsLogger<>();
    propagator.addEventDetector(logger.monitorDetector(detector));
    propagator.propagate(propagator.getInitialState().getOrbit().getDate().shiftedBy(Constants.JULIAN_DAY));
    Assert.assertEquals(30, logger.getLoggedEvents().size());
    for (FieldLoggedEvent<T> e : logger.getLoggedEvents()) {
        FieldKeplerianOrbit<T> o = (FieldKeplerianOrbit<T>) OrbitType.KEPLERIAN.convertType(e.getState().getOrbit());
        double expected = e.isIncreasing() ? 0.0 : FastMath.PI;
        Assert.assertEquals(expected, MathUtils.normalizeAngle(o.getMeanAnomaly().getReal(), expected), 4.0e-14);
    }
}

Example 4 with FieldEcksteinHechlerPropagator

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

the class BodyCenterPointingTest method doTestQDot.

private <T extends RealFieldElement<T>> void doTestQDot(final Field<T> field) 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;
    // Satellite position as circular parameters
    T zero = field.getZero();
    final T a = zero.add(7178000.0);
    final T e = zero.add(7e-5);
    final T i = zero.add(FastMath.toRadians(50.));
    final T pa = zero.add(FastMath.toRadians(45.));
    final T raan = zero.add(FastMath.toRadians(270.));
    final T m = zero.add(FastMath.toRadians(5.3 - 270));
    // Computation date
    FieldAbsoluteDate<T> date_comp = new FieldAbsoluteDate<>(field, new DateComponents(2008, 04, 07), TimeComponents.H00, TimeScalesFactory.getUTC());
    // Orbit
    FieldKeplerianOrbit<T> circ = new FieldKeplerianOrbit<>(a, e, i, pa, raan, m, PositionAngle.MEAN, FramesFactory.getEME2000(), date_comp, ehMu);
    // WGS84 Earth model
    OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
    // Earth center pointing attitude provider
    BodyCenterPointing earthCenterAttitudeLaw = new BodyCenterPointing(circ.getFrame(), earth);
    final FieldAbsoluteDate<T> date = FieldAbsoluteDate.getJ2000Epoch(field).shiftedBy(584.);
    final FieldVector3D<T> position = new FieldVector3D<>(zero.add(3220103.), zero.add(69623.), zero.add(6449822.));
    final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(6414.7), zero.add(-2006.), zero.add(-3180.));
    final FieldCircularOrbit<T> initialOrbit = new FieldCircularOrbit<>(new FieldPVCoordinates<>(position, velocity), FramesFactory.getEME2000(), date, ehMu);
    FieldEcksteinHechlerPropagator<T> propagator = new FieldEcksteinHechlerPropagator<>(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) {
        FieldRotation<T> rP = propagator.propagate(date.shiftedBy(dt)).getAttitude().getRotation();
        w0.add(new WeightedObservedPoint(1, dt, rP.getQ0().getReal()));
        w1.add(new WeightedObservedPoint(1, dt, rP.getQ1().getReal()));
        w2.add(new WeightedObservedPoint(1, dt, rP.getQ2().getReal()));
        w3.add(new WeightedObservedPoint(1, dt, rP.getQ3().getReal()));
    }
    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];
    FieldAttitude<T> a0 = propagator.propagate(date).getAttitude();
    T q0 = a0.getRotation().getQ0();
    T q1 = a0.getRotation().getQ1();
    T q2 = a0.getRotation().getQ2();
    T q3 = a0.getRotation().getQ3();
    T oX = a0.getSpin().getX();
    T oY = a0.getSpin().getY();
    T oZ = a0.getSpin().getZ();
    // first time-derivatives of the quaternion
    double q0Dot = 0.5 * MathArrays.linearCombination(-q1.getReal(), oX.getReal(), -q2.getReal(), oY.getReal(), -q3.getReal(), oZ.getReal());
    double q1Dot = 0.5 * MathArrays.linearCombination(q0.getReal(), oX.getReal(), -q3.getReal(), oY.getReal(), q2.getReal(), oZ.getReal());
    double q2Dot = 0.5 * MathArrays.linearCombination(q3.getReal(), oX.getReal(), q0.getReal(), oY.getReal(), -q1.getReal(), oZ.getReal());
    double q3Dot = 0.5 * MathArrays.linearCombination(-q2.getReal(), oX.getReal(), q1.getReal(), oY.getReal(), q0.getReal(), oZ.getReal());
    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 5 with FieldEcksteinHechlerPropagator

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

the class FieldAttitudeTest method doTestInterpolation.

private <T extends RealFieldElement<T>> void doTestInterpolation(final Field<T> field) throws OrekitException {
    T zero = field.getZero();
    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 FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field).shiftedBy(584.);
    final FieldVector3D<T> position = new FieldVector3D<>(zero.add(3220103.), zero.add(69623.), zero.add(6449822.));
    final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(6414.7), zero.add(-2006.), zero.add(-3180.));
    final FieldCircularOrbit<T> initialOrbit = new FieldCircularOrbit<>(new FieldPVCoordinates<>(position, velocity), FramesFactory.getEME2000(), date, ehMu);
    FieldEcksteinHechlerPropagator<T> propagator = new FieldEcksteinHechlerPropagator<>(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));
    FieldAttitude<T> initialAttitude = propagator.propagate(initialOrbit.getDate()).getAttitude();
    // set up a 5 points sample
    List<FieldAttitude<T>> sample = new ArrayList<FieldAttitude<T>>();
    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_R = 0; dt_R < 240.0; dt_R += 1.0) {
        T dt = zero.add(dt_R);
        FieldAbsoluteDate<T> t = initialOrbit.getDate().shiftedBy(dt);
        FieldAttitude<T> propagated = propagator.propagate(t).getAttitude();
        T shiftAngleError = FieldRotation.distance(propagated.getRotation(), initialAttitude.shiftedBy(dt).getRotation());
        T interpolationAngleError = FieldRotation.distance(propagated.getRotation(), initialAttitude.interpolate(t, sample).getRotation());
        T shiftRateError = FieldVector3D.distance(propagated.getSpin(), initialAttitude.shiftedBy(dt).getSpin());
        T interpolationRateError = FieldVector3D.distance(propagated.getSpin(), initialAttitude.interpolate(t, sample).getSpin());
        maxShiftAngleError = FastMath.max(maxShiftAngleError, shiftAngleError.getReal());
        maxInterpolationAngleError = FastMath.max(maxInterpolationAngleError, interpolationAngleError.getReal());
        maxShiftRateError = FastMath.max(maxShiftRateError, shiftRateError.getReal());
        maxInterpolationRateError = FastMath.max(maxInterpolationRateError, interpolationRateError.getReal());
    }
    Assert.assertTrue(maxShiftAngleError > 6.0e-6);
    Assert.assertTrue(maxInterpolationAngleError < 6.0e-15);
    Assert.assertTrue(maxShiftRateError > 7.0e-8);
    Assert.assertTrue(maxInterpolationRateError < 2.0e-16);
    // 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_R = 250.0; dt_R < 300.0; dt_R += 1.0) {
        T dt = zero.add(dt_R);
        FieldAbsoluteDate<T> t = initialOrbit.getDate().shiftedBy(dt);
        FieldAttitude<T> propagated = propagator.propagate(t).getAttitude();
        T shiftAngleError = FieldRotation.distance(propagated.getRotation(), initialAttitude.shiftedBy(dt).getRotation());
        T interpolationAngleError = FieldRotation.distance(propagated.getRotation(), initialAttitude.interpolate(t, sample).getRotation());
        T shiftRateError = FieldVector3D.distance(propagated.getSpin(), initialAttitude.shiftedBy(dt).getSpin());
        T interpolationRateError = FieldVector3D.distance(propagated.getSpin(), initialAttitude.interpolate(t, sample).getSpin());
        maxShiftAngleError = FastMath.max(maxShiftAngleError, shiftAngleError.getReal());
        maxInterpolationAngleError = FastMath.max(maxInterpolationAngleError, interpolationAngleError.getReal());
        maxShiftRateError = FastMath.max(maxShiftRateError, shiftRateError.getReal());
        maxInterpolationRateError = FastMath.max(maxInterpolationRateError, interpolationRateError.getReal());
    }
    Assert.assertTrue(maxShiftAngleError > 1.0e-5);
    Assert.assertTrue(maxInterpolationAngleError < 8.0e-13);
    Assert.assertTrue(maxShiftRateError > 1.0e-7);
    Assert.assertTrue(maxInterpolationRateError < 6.0e-14);
}