Examples with FieldPVCoordinates org.orekit.utils.FieldPVCoordinates used on opensource projects

Example 31 with FieldPVCoordinates

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

the class ThirdBodyAttractionTest method RealFieldTest.

/**
 *Testing if the propagation between the FieldPropagation and the propagation
 * is equivalent.
 * Also testing if propagating X+dX with the propagation is equivalent to
 * propagation X with the FieldPropagation and then applying the taylor
 * expansion of dX to the result.
 */
@Test
public void RealFieldTest() throws OrekitException {
    DSFactory factory = new DSFactory(6, 5);
    DerivativeStructure a_0 = factory.variable(0, 7e7);
    DerivativeStructure e_0 = factory.variable(1, 0.4);
    DerivativeStructure i_0 = factory.variable(2, 85 * FastMath.PI / 180);
    DerivativeStructure R_0 = factory.variable(3, 0.7);
    DerivativeStructure O_0 = factory.variable(4, 0.5);
    DerivativeStructure n_0 = factory.variable(5, 0.1);
    Field<DerivativeStructure> field = a_0.getField();
    DerivativeStructure zero = field.getZero();
    FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
    Frame EME = FramesFactory.getEME2000();
    FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0, PositionAngle.MEAN, EME, J2000, Constants.EIGEN5C_EARTH_MU);
    FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
    SpacecraftState iSR = initialState.toSpacecraftState();
    OrbitType type = OrbitType.KEPLERIAN;
    double[][] tolerance = NumericalPropagator.tolerances(10.0, FKO.toOrbit(), type);
    AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator = new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
    integrator.setInitialStepSize(zero.add(60));
    AdaptiveStepsizeIntegrator RIntegrator = new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
    RIntegrator.setInitialStepSize(60);
    FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
    FNP.setOrbitType(type);
    FNP.setInitialState(initialState);
    NumericalPropagator NP = new NumericalPropagator(RIntegrator);
    NP.setOrbitType(type);
    NP.setInitialState(iSR);
    final ThirdBodyAttraction forceModel = new ThirdBodyAttraction(CelestialBodyFactory.getSun());
    FNP.addForceModel(forceModel);
    NP.addForceModel(forceModel);
    FieldAbsoluteDate<DerivativeStructure> target = J2000.shiftedBy(1000.);
    FieldSpacecraftState<DerivativeStructure> finalState_DS = FNP.propagate(target);
    SpacecraftState finalState_R = NP.propagate(target.toAbsoluteDate());
    FieldPVCoordinates<DerivativeStructure> finPVC_DS = finalState_DS.getPVCoordinates();
    PVCoordinates finPVC_R = finalState_R.getPVCoordinates();
    Assert.assertEquals(finPVC_DS.toPVCoordinates().getPosition().getX(), finPVC_R.getPosition().getX(), FastMath.abs(finPVC_R.getPosition().getX()) * 1e-11);
    Assert.assertEquals(finPVC_DS.toPVCoordinates().getPosition().getY(), finPVC_R.getPosition().getY(), FastMath.abs(finPVC_R.getPosition().getY()) * 1e-11);
    Assert.assertEquals(finPVC_DS.toPVCoordinates().getPosition().getZ(), finPVC_R.getPosition().getZ(), FastMath.abs(finPVC_R.getPosition().getZ()) * 1e-11);
    long number = 23091991;
    RandomGenerator RG = new Well19937a(number);
    GaussianRandomGenerator NGG = new GaussianRandomGenerator(RG);
    UncorrelatedRandomVectorGenerator URVG = new UncorrelatedRandomVectorGenerator(new double[] { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, new double[] { 1e3, 0.01, 0.01, 0.01, 0.01, 0.01 }, NGG);
    double a_R = a_0.getReal();
    double e_R = e_0.getReal();
    double i_R = i_0.getReal();
    double R_R = R_0.getReal();
    double O_R = O_0.getReal();
    double n_R = n_0.getReal();
    double maxP = 0;
    double maxV = 0;
    double maxA = 0;
    for (int ii = 0; ii < 1; ii++) {
        double[] rand_next = URVG.nextVector();
        double a_shift = a_R + rand_next[0];
        double e_shift = e_R + rand_next[1];
        double i_shift = i_R + rand_next[2];
        double R_shift = R_R + rand_next[3];
        double O_shift = O_R + rand_next[4];
        double n_shift = n_R + rand_next[5];
        KeplerianOrbit shiftedOrb = new KeplerianOrbit(a_shift, e_shift, i_shift, R_shift, O_shift, n_shift, PositionAngle.MEAN, EME, J2000.toAbsoluteDate(), Constants.EIGEN5C_EARTH_MU);
        SpacecraftState shift_iSR = new SpacecraftState(shiftedOrb);
        NumericalPropagator shift_NP = new NumericalPropagator(RIntegrator);
        shift_NP.setInitialState(shift_iSR);
        shift_NP.addForceModel(forceModel);
        SpacecraftState finalState_shift = shift_NP.propagate(target.toAbsoluteDate());
        PVCoordinates finPVC_shift = finalState_shift.getPVCoordinates();
        // position check
        FieldVector3D<DerivativeStructure> pos_DS = finPVC_DS.getPosition();
        double x_DS = pos_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double y_DS = pos_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double z_DS = pos_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double x = finPVC_shift.getPosition().getX();
        double y = finPVC_shift.getPosition().getY();
        double z = finPVC_shift.getPosition().getZ();
        maxP = FastMath.max(maxP, FastMath.abs((x_DS - x) / (x - pos_DS.getX().getReal())));
        maxP = FastMath.max(maxP, FastMath.abs((y_DS - y) / (y - pos_DS.getY().getReal())));
        maxP = FastMath.max(maxP, FastMath.abs((z_DS - z) / (z - pos_DS.getZ().getReal())));
        // velocity check
        FieldVector3D<DerivativeStructure> vel_DS = finPVC_DS.getVelocity();
        double vx_DS = vel_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double vy_DS = vel_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double vz_DS = vel_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double vx = finPVC_shift.getVelocity().getX();
        double vy = finPVC_shift.getVelocity().getY();
        double vz = finPVC_shift.getVelocity().getZ();
        maxV = FastMath.max(maxV, FastMath.abs((vx_DS - vx) / vx));
        maxV = FastMath.max(maxV, FastMath.abs((vy_DS - vy) / vy));
        maxV = FastMath.max(maxV, FastMath.abs((vz_DS - vz) / vz));
        // acceleration check
        FieldVector3D<DerivativeStructure> acc_DS = finPVC_DS.getAcceleration();
        double ax_DS = acc_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double ay_DS = acc_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double az_DS = acc_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
        double ax = finPVC_shift.getAcceleration().getX();
        double ay = finPVC_shift.getAcceleration().getY();
        double az = finPVC_shift.getAcceleration().getZ();
        maxA = FastMath.max(maxA, FastMath.abs((ax_DS - ax) / ax));
        maxA = FastMath.max(maxA, FastMath.abs((ay_DS - ay) / ay));
        maxA = FastMath.max(maxA, FastMath.abs((az_DS - az) / az));
    }
    Assert.assertEquals(0, maxP, 5.0e-9);
    Assert.assertEquals(0, maxV, 3.0e-10);
    Assert.assertEquals(0, maxA, 8.0e-8);
}

Example 32 with FieldPVCoordinates

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

the class ThirdBodyAttractionTest method RealFieldExpectErrorTest.

/**
 *Same test as the previous one but not adding the ForceModel to the NumericalPropagator
 *    it is a test to validate the previous test.
 *    (to test if the ForceModel it's actually
 *    doing something in the Propagator and the FieldPropagator)
 */
@Test
public void RealFieldExpectErrorTest() throws OrekitException {
    DSFactory factory = new DSFactory(6, 5);
    DerivativeStructure a_0 = factory.variable(0, 7e7);
    DerivativeStructure e_0 = factory.variable(1, 0.4);
    DerivativeStructure i_0 = factory.variable(2, 85 * FastMath.PI / 180);
    DerivativeStructure R_0 = factory.variable(3, 0.7);
    DerivativeStructure O_0 = factory.variable(4, 0.5);
    DerivativeStructure n_0 = factory.variable(5, 0.1);
    Field<DerivativeStructure> field = a_0.getField();
    DerivativeStructure zero = field.getZero();
    FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
    Frame EME = FramesFactory.getEME2000();
    FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0, PositionAngle.MEAN, EME, J2000, Constants.EIGEN5C_EARTH_MU);
    FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
    SpacecraftState iSR = initialState.toSpacecraftState();
    OrbitType type = OrbitType.KEPLERIAN;
    double[][] tolerance = NumericalPropagator.tolerances(0.001, FKO.toOrbit(), type);
    AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator = new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
    integrator.setInitialStepSize(zero.add(60));
    AdaptiveStepsizeIntegrator RIntegrator = new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
    RIntegrator.setInitialStepSize(60);
    FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
    FNP.setOrbitType(type);
    FNP.setInitialState(initialState);
    NumericalPropagator NP = new NumericalPropagator(RIntegrator);
    NP.setOrbitType(type);
    NP.setInitialState(iSR);
    final ThirdBodyAttraction forceModel = new ThirdBodyAttraction(CelestialBodyFactory.getSun());
    FNP.addForceModel(forceModel);
    // NOT ADDING THE FORCE MODEL TO THE NUMERICAL PROPAGATOR   NP.addForceModel(forceModel);
    FieldAbsoluteDate<DerivativeStructure> target = J2000.shiftedBy(1000.);
    FieldSpacecraftState<DerivativeStructure> finalState_DS = FNP.propagate(target);
    SpacecraftState finalState_R = NP.propagate(target.toAbsoluteDate());
    FieldPVCoordinates<DerivativeStructure> finPVC_DS = finalState_DS.getPVCoordinates();
    PVCoordinates finPVC_R = finalState_R.getPVCoordinates();
    Assert.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getX() - finPVC_R.getPosition().getX()) < FastMath.abs(finPVC_R.getPosition().getX()) * 1e-11);
    Assert.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getY() - finPVC_R.getPosition().getY()) < FastMath.abs(finPVC_R.getPosition().getY()) * 1e-11);
    Assert.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getZ() - finPVC_R.getPosition().getZ()) < FastMath.abs(finPVC_R.getPosition().getZ()) * 1e-11);
}

Example 33 with FieldPVCoordinates

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

the class FieldEventDetectorTest method doTestWrappedException.

private <T extends RealFieldElement<T>> void doTestWrappedException(Field<T> field) throws OrekitException {
    final T zero = field.getZero();
    final Throwable dummyCause = new RuntimeException();
    try {
        // initial conditions
        Frame eme2000 = FramesFactory.getEME2000();
        TimeScale utc = TimeScalesFactory.getUTC();
        final FieldAbsoluteDate<T> initialDate = new FieldAbsoluteDate<>(field, 2011, 5, 11, utc);
        final FieldAbsoluteDate<T> exceptionDate = initialDate.shiftedBy(3600.0);
        FieldKeplerianPropagator<T> k = new FieldKeplerianPropagator<>(new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(new FieldVector3D<>(zero.add(4008462.4706055815), zero.add(-3155502.5373837613), zero.add(-5044275.9880020910)), new FieldVector3D<>(zero.add(-5012.9298276860990), zero.add(1920.3567095973078), zero.add(-5172.7403501801580))), eme2000, initialDate, Constants.WGS84_EARTH_MU));
        k.addEventDetector(new FieldDateDetector<T>(initialDate.shiftedBy(Constants.JULIAN_DAY)) {

            @Override
            public T g(final FieldSpacecraftState<T> s) throws OrekitException {
                final T dt = s.getDate().durationFrom(exceptionDate);
                if (dt.abs().getReal() < 1.0) {
                    throw new OrekitException(dummyCause, LocalizedCoreFormats.SIMPLE_MESSAGE, "dummy");
                }
                return dt;
            }
        });
        k.propagate(initialDate.shiftedBy(Constants.JULIAN_YEAR));
        Assert.fail("an exception should have been thrown");
    } catch (OrekitException oe) {
        Assert.assertSame(dummyCause, oe.getCause());
    }
}

Example 34 with FieldPVCoordinates

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

the class FieldTransformTest method doTestJacobianPV.

private <T extends RealFieldElement<T>> void doTestJacobianPV(Field<T> field) {
    // base directions for finite differences
    @SuppressWarnings("unchecked") FieldPVCoordinates<T>[] directions = (FieldPVCoordinates<T>[]) Array.newInstance(FieldPVCoordinates.class, 6);
    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));
    directions[3] = new FieldPVCoordinates<>(FieldVector3D.getZero(field), FieldVector3D.getPlusI(field), FieldVector3D.getZero(field));
    directions[4] = new FieldPVCoordinates<>(FieldVector3D.getZero(field), FieldVector3D.getPlusJ(field), FieldVector3D.getZero(field));
    directions[5] = new FieldPVCoordinates<>(FieldVector3D.getZero(field), FieldVector3D.getPlusK(field), FieldVector3D.getZero(field));
    double h = 0.01;
    RandomGenerator random = new Well19937a(0xce2bfddfbb9796bel);
    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_PV, 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();
            double epsilonV = 6.0e-11 * pv0.getVelocity().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);
                Assert.assertEquals(estimatedColumn.getVelocity().getX().getReal(), jacobian[3][c].getReal(), epsilonV);
                Assert.assertEquals(estimatedColumn.getVelocity().getY().getReal(), jacobian[4][c].getReal(), epsilonV);
                Assert.assertEquals(estimatedColumn.getVelocity().getZ().getReal(), jacobian[5][c].getReal(), epsilonV);
                // check the rest of the matrix remains untouched
                for (int l = 6; 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);
                }
            }
        }
    }
}

Example 35 with FieldPVCoordinates

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

the class FieldCartesianOrbitTest 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 FieldCartesianOrbit<>(pv, eme2000, FieldAbsoluteDate.getJ2000Epoch(field), mu);
    // Build another KeplerianOrbit as a copy of the first one
    final FieldOrbit<T> orbitCopy = new FieldCartesianOrbit<>(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);
}