Examples with PolynomialFunction org.hipparchus.analysis.polynomials.PolynomialFunction used on opensource projects

Example 1 with PolynomialFunction

use of org.hipparchus.analysis.polynomials.PolynomialFunction in project Orekit by CS-SI.

the class TimeStampedPVCoordinatesTest method testInterpolatePolynomialPV.

@Test
public void testInterpolatePolynomialPV() {
    Random random = new Random(0xae7771c9933407bdl);
    AbsoluteDate t0 = AbsoluteDate.J2000_EPOCH;
    for (int i = 0; i < 20; ++i) {
        PolynomialFunction px = randomPolynomial(5, random);
        PolynomialFunction py = randomPolynomial(5, random);
        PolynomialFunction pz = randomPolynomial(5, random);
        PolynomialFunction pxDot = px.polynomialDerivative();
        PolynomialFunction pyDot = py.polynomialDerivative();
        PolynomialFunction pzDot = pz.polynomialDerivative();
        PolynomialFunction pxDotDot = pxDot.polynomialDerivative();
        PolynomialFunction pyDotDot = pyDot.polynomialDerivative();
        PolynomialFunction pzDotDot = pzDot.polynomialDerivative();
        List<TimeStampedPVCoordinates> sample = new ArrayList<TimeStampedPVCoordinates>();
        for (double dt : new double[] { 0.0, 0.5, 1.0 }) {
            Vector3D position = new Vector3D(px.value(dt), py.value(dt), pz.value(dt));
            Vector3D velocity = new Vector3D(pxDot.value(dt), pyDot.value(dt), pzDot.value(dt));
            sample.add(new TimeStampedPVCoordinates(t0.shiftedBy(dt), position, velocity, Vector3D.ZERO));
        }
        for (double dt = 0; dt < 1.0; dt += 0.01) {
            TimeStampedPVCoordinates interpolated = TimeStampedPVCoordinates.interpolate(t0.shiftedBy(dt), CartesianDerivativesFilter.USE_PV, sample);
            Vector3D p = interpolated.getPosition();
            Vector3D v = interpolated.getVelocity();
            Vector3D a = interpolated.getAcceleration();
            Assert.assertEquals(px.value(dt), p.getX(), 4.0e-16 * p.getNorm());
            Assert.assertEquals(py.value(dt), p.getY(), 4.0e-16 * p.getNorm());
            Assert.assertEquals(pz.value(dt), p.getZ(), 4.0e-16 * p.getNorm());
            Assert.assertEquals(pxDot.value(dt), v.getX(), 9.0e-16 * v.getNorm());
            Assert.assertEquals(pyDot.value(dt), v.getY(), 9.0e-16 * v.getNorm());
            Assert.assertEquals(pzDot.value(dt), v.getZ(), 9.0e-16 * v.getNorm());
            Assert.assertEquals(pxDotDot.value(dt), a.getX(), 1.0e-14 * a.getNorm());
            Assert.assertEquals(pyDotDot.value(dt), a.getY(), 1.0e-14 * a.getNorm());
            Assert.assertEquals(pzDotDot.value(dt), a.getZ(), 1.0e-14 * a.getNorm());
        }
    }
}

Example 2 with PolynomialFunction

use of org.hipparchus.analysis.polynomials.PolynomialFunction in project Orekit by CS-SI.

the class TimeStampedPVCoordinatesTest method testInterpolatePolynomialPVA.

@Test
public void testInterpolatePolynomialPVA() {
    Random random = new Random(0xfe3945fcb8bf47cel);
    AbsoluteDate t0 = AbsoluteDate.J2000_EPOCH;
    for (int i = 0; i < 20; ++i) {
        PolynomialFunction px = randomPolynomial(5, random);
        PolynomialFunction py = randomPolynomial(5, random);
        PolynomialFunction pz = randomPolynomial(5, random);
        PolynomialFunction pxDot = px.polynomialDerivative();
        PolynomialFunction pyDot = py.polynomialDerivative();
        PolynomialFunction pzDot = pz.polynomialDerivative();
        PolynomialFunction pxDotDot = pxDot.polynomialDerivative();
        PolynomialFunction pyDotDot = pyDot.polynomialDerivative();
        PolynomialFunction pzDotDot = pzDot.polynomialDerivative();
        List<TimeStampedPVCoordinates> sample = new ArrayList<TimeStampedPVCoordinates>();
        for (double dt : new double[] { 0.0, 0.5, 1.0 }) {
            Vector3D position = new Vector3D(px.value(dt), py.value(dt), pz.value(dt));
            Vector3D velocity = new Vector3D(pxDot.value(dt), pyDot.value(dt), pzDot.value(dt));
            Vector3D acceleration = new Vector3D(pxDotDot.value(dt), pyDotDot.value(dt), pzDotDot.value(dt));
            sample.add(new TimeStampedPVCoordinates(t0.shiftedBy(dt), position, velocity, acceleration));
        }
        for (double dt = 0; dt < 1.0; dt += 0.01) {
            TimeStampedPVCoordinates interpolated = TimeStampedPVCoordinates.interpolate(t0.shiftedBy(dt), CartesianDerivativesFilter.USE_PVA, sample);
            Vector3D p = interpolated.getPosition();
            Vector3D v = interpolated.getVelocity();
            Vector3D a = interpolated.getAcceleration();
            Assert.assertEquals(px.value(dt), p.getX(), 4.0e-16 * p.getNorm());
            Assert.assertEquals(py.value(dt), p.getY(), 4.0e-16 * p.getNorm());
            Assert.assertEquals(pz.value(dt), p.getZ(), 4.0e-16 * p.getNorm());
            Assert.assertEquals(pxDot.value(dt), v.getX(), 9.0e-16 * v.getNorm());
            Assert.assertEquals(pyDot.value(dt), v.getY(), 9.0e-16 * v.getNorm());
            Assert.assertEquals(pzDot.value(dt), v.getZ(), 9.0e-16 * v.getNorm());
            Assert.assertEquals(pxDotDot.value(dt), a.getX(), 9.0e-15 * a.getNorm());
            Assert.assertEquals(pyDotDot.value(dt), a.getY(), 9.0e-15 * a.getNorm());
            Assert.assertEquals(pzDotDot.value(dt), a.getZ(), 9.0e-15 * a.getNorm());
        }
    }
}

Example 3 with PolynomialFunction

use of org.hipparchus.analysis.polynomials.PolynomialFunction in project Orekit by CS-SI.

the class CoefficientFactoryTest method getQnsPolynomialValue.

/**
 * Get the Q<sub>ns</sub> value from 2.8.1-(4) evaluated in γ This method is using the
 * Legendre polynomial to compute the Q<sub>ns</sub>'s one. This direct computation method
 * allows to store the polynomials value in a static map. If the Q<sub>ns</sub> had been
 * computed already, they just will be evaluated at γ
 *
 * @param gamma γ angle for which Q<sub>ns</sub> is evaluated
 * @param n n value
 * @param s s value
 * @return the polynomial value evaluated at γ
 */
private static double getQnsPolynomialValue(final double gamma, final int n, final int s) {
    PolynomialFunction derivative;
    if (QNS_MAP.containsKey(new NSKey(n, s))) {
        derivative = QNS_MAP.get(new NSKey(n, s));
    } else {
        final PolynomialFunction legendre = PolynomialsUtils.createLegendrePolynomial(n);
        derivative = legendre;
        for (int i = 0; i < s; i++) {
            derivative = (PolynomialFunction) derivative.polynomialDerivative();
        }
        QNS_MAP.put(new NSKey(n, s), derivative);
    }
    return derivative.value(gamma);
}

Example 4 with PolynomialFunction

use of org.hipparchus.analysis.polynomials.PolynomialFunction 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 5 with PolynomialFunction

use of org.hipparchus.analysis.polynomials.PolynomialFunction in project Orekit by CS-SI.

the class SpacecraftStateTest method testShiftVsEcksteinHechlerError.

@Test
public void testShiftVsEcksteinHechlerError() throws ParseException, OrekitException {
    // 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 });
    AbsoluteDate centerDate = orbit.getDate().shiftedBy(100.0);
    SpacecraftState centerState = propagator.propagate(centerDate);
    double maxResidualP = 0;
    double maxResidualV = 0;
    double maxResidualA = 0;
    double maxResidualR = 0;
    for (double dt = 0; dt < 900.0; dt += 5) {
        SpacecraftState shifted = centerState.shiftedBy(dt);
        SpacecraftState propagated = propagator.propagate(centerDate.shiftedBy(dt));
        PVCoordinates dpv = new PVCoordinates(propagated.getPVCoordinates(), shifted.getPVCoordinates());
        double residualP = pModel.value(dt) - dpv.getPosition().getNorm();
        double residualV = vModel.value(dt) - dpv.getVelocity().getNorm();
        double residualA = aModel.value(dt) - dpv.getAcceleration().getNorm();
        double residualR = rModel.value(dt) - FastMath.toDegrees(Rotation.distance(shifted.getAttitude().getRotation(), propagated.getAttitude().getRotation()));
        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(7.7e-6, maxResidualA, 1.0e-7);
    Assert.assertEquals(2.8e-6, maxResidualR, 1.0e-1);
}