Search in sources :

Example 16 with FieldPVCoordinates

use of org.orekit.utils.FieldPVCoordinates in project Orekit by CS-SI.

the class FieldTransformTest method doTestJacobianP.

private <T extends RealFieldElement<T>> void doTestJacobianP(Field<T> field) {
    // base directions for finite differences
    @SuppressWarnings("unchecked") FieldPVCoordinates<T>[] directions = (FieldPVCoordinates<T>[]) Array.newInstance(FieldPVCoordinates.class, 3);
    directions[0] = new FieldPVCoordinates<>(FieldVector3D.getPlusI(field), FieldVector3D.getZero(field), FieldVector3D.getZero(field));
    directions[1] = new FieldPVCoordinates<>(FieldVector3D.getPlusJ(field), FieldVector3D.getZero(field), FieldVector3D.getZero(field));
    directions[2] = new FieldPVCoordinates<>(FieldVector3D.getPlusK(field), FieldVector3D.getZero(field), FieldVector3D.getZero(field));
    double h = 0.01;
    RandomGenerator random = new Well19937a(0x47fd0d6809f4b173l);
    for (int i = 0; i < 20; ++i) {
        // generate a random transform
        FieldTransform<T> combined = randomTransform(field, random);
        // compute Jacobian
        T[][] jacobian = MathArrays.buildArray(field, 9, 9);
        for (int l = 0; l < jacobian.length; ++l) {
            for (int c = 0; c < jacobian[l].length; ++c) {
                jacobian[l][c] = field.getZero().add(l + 0.1 * c);
            }
        }
        combined.getJacobian(CartesianDerivativesFilter.USE_P, jacobian);
        for (int j = 0; j < 100; ++j) {
            FieldPVCoordinates<T> pv0 = new FieldPVCoordinates<>(randomVector(field, 1e3, random), randomVector(field, 1.0, random), randomVector(field, 1.0e-3, random));
            double epsilonP = 2.0e-12 * pv0.getPosition().getNorm().getReal();
            for (int c = 0; c < directions.length; ++c) {
                // eight points finite differences estimation of a Jacobian column
                FieldPVCoordinates<T> pvm4h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, -4 * h, directions[c]));
                FieldPVCoordinates<T> pvm3h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, -3 * h, directions[c]));
                FieldPVCoordinates<T> pvm2h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, -2 * h, directions[c]));
                FieldPVCoordinates<T> pvm1h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, -1 * h, directions[c]));
                FieldPVCoordinates<T> pvp1h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, +1 * h, directions[c]));
                FieldPVCoordinates<T> pvp2h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, +2 * h, directions[c]));
                FieldPVCoordinates<T> pvp3h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, +3 * h, directions[c]));
                FieldPVCoordinates<T> pvp4h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, +4 * h, directions[c]));
                FieldPVCoordinates<T> d4 = new FieldPVCoordinates<>(pvm4h, pvp4h);
                FieldPVCoordinates<T> d3 = new FieldPVCoordinates<>(pvm3h, pvp3h);
                FieldPVCoordinates<T> d2 = new FieldPVCoordinates<>(pvm2h, pvp2h);
                FieldPVCoordinates<T> d1 = new FieldPVCoordinates<>(pvm1h, pvp1h);
                double d = 1.0 / (840 * h);
                FieldPVCoordinates<T> estimatedColumn = new FieldPVCoordinates<>(-3 * d, d4, 32 * d, d3, -168 * d, d2, 672 * d, d1);
                // check analytical Jacobian against finite difference reference
                Assert.assertEquals(estimatedColumn.getPosition().getX().getReal(), jacobian[0][c].getReal(), epsilonP);
                Assert.assertEquals(estimatedColumn.getPosition().getY().getReal(), jacobian[1][c].getReal(), epsilonP);
                Assert.assertEquals(estimatedColumn.getPosition().getZ().getReal(), jacobian[2][c].getReal(), epsilonP);
                // check the rest of the matrix remains untouched
                for (int l = 3; l < jacobian.length; ++l) {
                    Assert.assertEquals(l + 0.1 * c, jacobian[l][c].getReal(), 1.0e-15);
                }
            }
            // check the rest of the matrix remains untouched
            for (int c = directions.length; c < jacobian[0].length; ++c) {
                for (int l = 0; l < jacobian.length; ++l) {
                    Assert.assertEquals(l + 0.1 * c, jacobian[l][c].getReal(), 1.0e-15);
                }
            }
        }
    }
}
Also used : FieldPVCoordinates(org.orekit.utils.FieldPVCoordinates) TimeStampedFieldPVCoordinates(org.orekit.utils.TimeStampedFieldPVCoordinates) Well19937a(org.hipparchus.random.Well19937a) RandomGenerator(org.hipparchus.random.RandomGenerator)

Example 17 with FieldPVCoordinates

use of org.orekit.utils.FieldPVCoordinates in project Orekit by CS-SI.

the class FieldEquinoctialOrbitTest method doTestNonInertialFrame.

private <T extends RealFieldElement<T>> void doTestNonInertialFrame(Field<T> field) throws IllegalArgumentException {
    T zero = field.getZero();
    FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
    FieldVector3D<T> position = new FieldVector3D<>(zero.add(4512.9), zero.add(18260.), zero.add(-5127.));
    FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(134664.6), zero.add(90066.8), zero.add(72047.6));
    FieldPVCoordinates<T> FieldPVCoordinates = new FieldPVCoordinates<>(position, velocity);
    new FieldEquinoctialOrbit<>(FieldPVCoordinates, new Frame(FramesFactory.getEME2000(), Transform.IDENTITY, "non-inertial", false), date, mu);
}
Also used : Frame(org.orekit.frames.Frame) FieldPVCoordinates(org.orekit.utils.FieldPVCoordinates) TimeStampedFieldPVCoordinates(org.orekit.utils.TimeStampedFieldPVCoordinates) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D)

Example 18 with FieldPVCoordinates

use of org.orekit.utils.FieldPVCoordinates in project Orekit by CS-SI.

the class FieldKeplerianOrbitTest method doTestVeryLargeEccentricity.

private <T extends RealFieldElement<T>> void doTestVeryLargeEccentricity(final Field<T> field) {
    FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
    final Frame eme2000 = FramesFactory.getEME2000();
    final double meanAnomaly = 1.;
    final FieldKeplerianOrbit<T> orb0 = new FieldKeplerianOrbit<>(field.getZero().add(42600e3), field.getZero().add(0.9), field.getZero().add(0.00001), field.getZero().add(0), field.getZero().add(0), field.getZero().add(FastMath.toRadians(meanAnomaly)), PositionAngle.MEAN, eme2000, date, mu);
    // big dV along Y
    final FieldVector3D<T> deltaV = new FieldVector3D<>(field.getZero().add(0.0), field.getZero().add(110000.0), field.getZero());
    final FieldPVCoordinates<T> pv1 = new FieldPVCoordinates<>(orb0.getPVCoordinates().getPosition(), orb0.getPVCoordinates().getVelocity().add(deltaV));
    final FieldKeplerianOrbit<T> orb1 = new FieldKeplerianOrbit<>(pv1, eme2000, date, mu);
    // Despite large eccentricity, the conversion of mean anomaly to hyperbolic eccentric anomaly
    // converges in less than 50 iterations (issue #114)
    final FieldPVCoordinates<T> pvTested = orb1.shiftedBy(field.getZero()).getPVCoordinates();
    final FieldVector3D<T> pTested = pvTested.getPosition();
    final FieldVector3D<T> vTested = pvTested.getVelocity();
    final FieldPVCoordinates<T> pvReference = orb1.getPVCoordinates();
    final FieldVector3D<T> pReference = pvReference.getPosition();
    final FieldVector3D<T> vReference = pvReference.getVelocity();
    final double threshold = 1.e-15;
    Assert.assertEquals(0, pTested.subtract(pReference).getNorm().getReal(), pReference.getNorm().multiply(threshold).getReal());
    Assert.assertEquals(0, vTested.subtract(vReference).getNorm().getReal(), vReference.getNorm().multiply(threshold).getReal());
}
Also used : Frame(org.orekit.frames.Frame) FieldPVCoordinates(org.orekit.utils.FieldPVCoordinates) TimeStampedFieldPVCoordinates(org.orekit.utils.TimeStampedFieldPVCoordinates) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D)

Example 19 with FieldPVCoordinates

use of org.orekit.utils.FieldPVCoordinates in project Orekit by CS-SI.

the class FieldKeplerianOrbitTest method doTestCopyNonKeplerianAcceleration.

private <T extends RealFieldElement<T>> void doTestCopyNonKeplerianAcceleration(Field<T> field) throws OrekitException {
    final Frame eme2000 = FramesFactory.getEME2000();
    // Define GEO satellite position
    final FieldVector3D<T> position = new FieldVector3D<>(field.getZero().add(42164140), field.getZero(), field.getZero());
    // Build PVCoodrinates starting from its position and computing the corresponding circular velocity
    final FieldPVCoordinates<T> pv = new FieldPVCoordinates<>(position, new FieldVector3D<>(field.getZero(), position.getNorm().reciprocal().multiply(mu).sqrt(), field.getZero()));
    // Build a KeplerianOrbit in eme2000
    final FieldOrbit<T> orbit = new FieldKeplerianOrbit<>(pv, eme2000, FieldAbsoluteDate.getJ2000Epoch(field), mu);
    // Build another KeplerianOrbit as a copy of the first one
    final FieldOrbit<T> orbitCopy = new FieldKeplerianOrbit<>(orbit);
    // Shift the orbit of a time-interval
    // This works good
    final FieldOrbit<T> shiftedOrbit = orbit.shiftedBy(10);
    // This does not work
    final FieldOrbit<T> shiftedOrbitCopy = orbitCopy.shiftedBy(10);
    Assert.assertEquals(0.0, FieldVector3D.distance(shiftedOrbit.getPVCoordinates().getPosition(), shiftedOrbitCopy.getPVCoordinates().getPosition()).getReal(), 1.0e-10);
    Assert.assertEquals(0.0, FieldVector3D.distance(shiftedOrbit.getPVCoordinates().getVelocity(), shiftedOrbitCopy.getPVCoordinates().getVelocity()).getReal(), 1.0e-10);
}
Also used : Frame(org.orekit.frames.Frame) FieldPVCoordinates(org.orekit.utils.FieldPVCoordinates) TimeStampedFieldPVCoordinates(org.orekit.utils.TimeStampedFieldPVCoordinates) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D)

Example 20 with FieldPVCoordinates

use of org.orekit.utils.FieldPVCoordinates 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);
}
Also used : PolynomialFunction(org.hipparchus.analysis.polynomials.PolynomialFunction) DateComponents(org.orekit.time.DateComponents) FieldKeplerianOrbit(org.orekit.orbits.FieldKeplerianOrbit) BodyCenterPointing(org.orekit.attitudes.BodyCenterPointing) FieldEcksteinHechlerPropagator(org.orekit.propagation.analytical.FieldEcksteinHechlerPropagator) FieldPVCoordinates(org.orekit.utils.FieldPVCoordinates) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate)

Aggregations

FieldPVCoordinates (org.orekit.utils.FieldPVCoordinates)41 Frame (org.orekit.frames.Frame)27 FieldAbsoluteDate (org.orekit.time.FieldAbsoluteDate)22 TimeStampedFieldPVCoordinates (org.orekit.utils.TimeStampedFieldPVCoordinates)19 PVCoordinates (org.orekit.utils.PVCoordinates)17 DerivativeStructure (org.hipparchus.analysis.differentiation.DerivativeStructure)16 RandomGenerator (org.hipparchus.random.RandomGenerator)16 Well19937a (org.hipparchus.random.Well19937a)16 FieldKeplerianOrbit (org.orekit.orbits.FieldKeplerianOrbit)16 FieldVector3D (org.hipparchus.geometry.euclidean.threed.FieldVector3D)15 DSFactory (org.hipparchus.analysis.differentiation.DSFactory)14 Test (org.junit.Test)14 AbstractLegacyForceModelTest (org.orekit.forces.AbstractLegacyForceModelTest)13 OrbitType (org.orekit.orbits.OrbitType)13 FieldSpacecraftState (org.orekit.propagation.FieldSpacecraftState)13 SpacecraftState (org.orekit.propagation.SpacecraftState)13 FieldNumericalPropagator (org.orekit.propagation.numerical.FieldNumericalPropagator)13 NumericalPropagator (org.orekit.propagation.numerical.NumericalPropagator)13 AdaptiveStepsizeIntegrator (org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator)12 DormandPrince853FieldIntegrator (org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator)12