Examples with DormandPrince853FieldIntegrator org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator used on opensource projects

Example 26 with DormandPrince853FieldIntegrator

use of org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator in project Orekit by CS-SI.

the class FieldNumericalPropagatorTest method doTestContinueEvent.

private <T extends RealFieldElement<T>> void doTestContinueEvent(Field<T> field) throws OrekitException {
    T zero = field.getZero();
    // setup
    final FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field);
    FieldSpacecraftState<T> initialState;
    FieldNumericalPropagator<T> propagator;
    final FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
    final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
    final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity), FramesFactory.getEME2000(), initDate, mu);
    initialState = new FieldSpacecraftState<>(orbit);
    OrbitType type = OrbitType.EQUINOCTIAL;
    double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit.toOrbit(), type);
    AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
    integrator.setInitialStepSize(zero.add(60));
    propagator = new FieldNumericalPropagator<>(field, integrator);
    propagator.setOrbitType(type);
    propagator.setInitialState(initialState);
    final FieldAbsoluteDate<T> resetDate = initDate.shiftedBy(1000);
    CheckingHandler<FieldDateDetector<T>, T> checking = new CheckingHandler<FieldDateDetector<T>, T>(Action.CONTINUE);
    propagator.addEventDetector(new FieldDateDetector<>(resetDate).withHandler(checking));
    final double dt = 3200;
    checking.assertEvent(false);
    Assert.assertEquals(0.0, propagator.getInitialState().getDate().durationFrom(initDate).getReal(), 1.0e-10);
    propagator.setResetAtEnd(false);
    final FieldSpacecraftState<T> finalState = propagator.propagate(initDate.shiftedBy(dt));
    Assert.assertEquals(0.0, propagator.getInitialState().getDate().durationFrom(initDate).getReal(), 1.0e-10);
    checking.assertEvent(true);
    final double n = FastMath.sqrt(initialState.getMu() / initialState.getA().getReal()) / initialState.getA().getReal();
    Assert.assertEquals(initialState.getA().getReal(), finalState.getA().getReal(), 1.0e-10);
    Assert.assertEquals(initialState.getEquinoctialEx().getReal(), finalState.getEquinoctialEx().getReal(), 1.0e-10);
    Assert.assertEquals(initialState.getEquinoctialEy().getReal(), finalState.getEquinoctialEy().getReal(), 1.0e-10);
    Assert.assertEquals(initialState.getHx().getReal(), finalState.getHx().getReal(), 1.0e-10);
    Assert.assertEquals(initialState.getHy().getReal(), finalState.getHy().getReal(), 1.0e-10);
    Assert.assertEquals(initialState.getLM().getReal() + n * dt, finalState.getLM().getReal(), 6.0e-10);
}

Example 27 with DormandPrince853FieldIntegrator

use of org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator in project Orekit by CS-SI.

the class FieldNumericalPropagatorTest method doTestCartesian.

private <T extends RealFieldElement<T>> void doTestCartesian(Field<T> field) throws OrekitException {
    T zero = field.getZero();
    // setup
    final FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field);
    FieldSpacecraftState<T> initialState;
    FieldNumericalPropagator<T> propagator;
    final FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
    final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
    final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity), FramesFactory.getEME2000(), initDate, mu);
    initialState = new FieldSpacecraftState<>(orbit);
    OrbitType type = OrbitType.EQUINOCTIAL;
    double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit.toOrbit(), type);
    AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
    integrator.setInitialStepSize(zero.add(60));
    propagator = new FieldNumericalPropagator<>(field, integrator);
    propagator.setOrbitType(type);
    propagator.setInitialState(initialState);
    // Propagation of the initial at t + dt
    final T dt = zero.add(3200);
    propagator.setOrbitType(OrbitType.CARTESIAN);
    final FieldPVCoordinates<T> finalState = propagator.propagate(initDate.shiftedBy(dt)).getPVCoordinates();
    final FieldVector3D<T> pFin = finalState.getPosition();
    final FieldVector3D<T> vFin = finalState.getVelocity();
    // Check results
    final FieldPVCoordinates<T> reference = initialState.shiftedBy(dt).getPVCoordinates();
    final FieldVector3D<T> pRef = reference.getPosition();
    final FieldVector3D<T> vRef = reference.getVelocity();
    Assert.assertEquals(0, pRef.subtract(pFin).getNorm().getReal(), 2e-4);
    Assert.assertEquals(0, vRef.subtract(vFin).getNorm().getReal(), 7e-8);
    try {
        propagator.getGeneratedEphemeris();
        Assert.fail("an exception should have been thrown");
    } catch (IllegalStateException ise) {
    // expected
    }
}

Example 28 with DormandPrince853FieldIntegrator

use of org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator in project Orekit by CS-SI.

the class FieldNumericalPropagatorTest method doTestAdditionalStateEvent.

private <T extends RealFieldElement<T>> void doTestAdditionalStateEvent(Field<T> field) throws OrekitException {
    T zero = field.getZero();
    // setup
    final FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field);
    FieldSpacecraftState<T> initialState;
    FieldNumericalPropagator<T> propagator;
    final FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
    final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
    final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity), FramesFactory.getEME2000(), initDate, mu);
    initialState = new FieldSpacecraftState<>(orbit);
    OrbitType type = OrbitType.EQUINOCTIAL;
    double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit.toOrbit(), type);
    AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
    integrator.setInitialStepSize(zero.add(60));
    propagator = new FieldNumericalPropagator<>(field, integrator);
    propagator.setOrbitType(type);
    propagator.setInitialState(initialState);
    propagator.addAdditionalEquations(new FieldAdditionalEquations<T>() {

        public String getName() {
            return "linear";
        }

        public T[] computeDerivatives(FieldSpacecraftState<T> s, T[] pDot) {
            pDot[0] = zero.add(1.0);
            return MathArrays.buildArray(field, 7);
        }
    });
    try {
        propagator.addAdditionalEquations(new FieldAdditionalEquations<T>() {

            public String getName() {
                return "linear";
            }

            public T[] computeDerivatives(FieldSpacecraftState<T> s, T[] pDot) {
                pDot[0] = zero.add(1.0);
                return MathArrays.buildArray(field, 7);
            }
        });
        Assert.fail("an exception should have been thrown");
    } catch (OrekitException oe) {
        Assert.assertEquals(oe.getSpecifier(), OrekitMessages.ADDITIONAL_STATE_NAME_ALREADY_IN_USE);
    }
    try {
        propagator.addAdditionalStateProvider(new FieldAdditionalStateProvider<T>() {

            public String getName() {
                return "linear";
            }

            public T[] getAdditionalState(FieldSpacecraftState<T> state) {
                return null;
            }
        });
        Assert.fail("an exception should have been thrown");
    } catch (OrekitException oe) {
        Assert.assertEquals(oe.getSpecifier(), OrekitMessages.ADDITIONAL_STATE_NAME_ALREADY_IN_USE);
    }
    propagator.addAdditionalStateProvider(new FieldAdditionalStateProvider<T>() {

        public String getName() {
            return "constant";
        }

        public T[] getAdditionalState(FieldSpacecraftState<T> state) {
            T[] ret = MathArrays.buildArray(field, 1);
            ret[0] = zero.add(1.0);
            return ret;
        }
    });
    Assert.assertTrue(propagator.isAdditionalStateManaged("linear"));
    Assert.assertTrue(propagator.isAdditionalStateManaged("constant"));
    Assert.assertFalse(propagator.isAdditionalStateManaged("non-managed"));
    Assert.assertEquals(2, propagator.getManagedAdditionalStates().length);
    propagator.setInitialState(propagator.getInitialState().addAdditionalState("linear", zero.add(1.5)));
    CheckingHandler<AdditionalStateLinearDetector<T>, T> checking = new CheckingHandler<AdditionalStateLinearDetector<T>, T>(Action.STOP);
    propagator.addEventDetector(new AdditionalStateLinearDetector<T>(zero.add(10.0), zero.add(1.0e-8)).withHandler(checking));
    final double dt = 3200;
    checking.assertEvent(false);
    final FieldSpacecraftState<T> finalState = propagator.propagate(initDate.shiftedBy(dt));
    checking.assertEvent(true);
    Assert.assertEquals(3.0, finalState.getAdditionalState("linear")[0].getReal(), 1.0e-8);
    Assert.assertEquals(1.5, finalState.getDate().durationFrom(initDate).getReal(), 1.0e-8);
}

Example 29 with DormandPrince853FieldIntegrator

use of org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator in project Orekit by CS-SI.

the class FieldNumericalPropagatorTest method doTestEventDetectionBug.

private <T extends RealFieldElement<T>> void doTestEventDetectionBug(final Field<T> field) throws OrekitException {
    T zero = field.getZero();
    TimeScale utc = TimeScalesFactory.getUTC();
    FieldAbsoluteDate<T> initialDate = new FieldAbsoluteDate<>(field, 2005, 1, 1, 0, 0, 0.0, utc);
    T duration = zero.add(100000.0);
    FieldAbsoluteDate<T> endDate = new FieldAbsoluteDate<>(initialDate, duration);
    // Initialization of the frame EME2000
    Frame EME2000 = FramesFactory.getEME2000();
    // Initial orbit
    double a = 35786000. + 6378137.0;
    double e = 0.70;
    double rApogee = a * (1 + e);
    double vApogee = FastMath.sqrt(mu * (1 - e) / (a * (1 + e)));
    FieldOrbit<T> geo = new FieldCartesianOrbit<>(new FieldPVCoordinates<>(new FieldVector3D<>(zero.add(rApogee), zero, zero), new FieldVector3D<>(zero, zero.add(vApogee), zero)), EME2000, initialDate, mu);
    duration = geo.getKeplerianPeriod();
    endDate = new FieldAbsoluteDate<>(initialDate, duration);
    // Numerical Integration
    final double minStep = 0.001;
    final double maxStep = 1000;
    final double initStep = 60;
    final OrbitType type = OrbitType.EQUINOCTIAL;
    final double[] absTolerance = { 0.001, 1.0e-9, 1.0e-9, 1.0e-6, 1.0e-6, 1.0e-6, 0.001 };
    final double[] relTolerance = { 1.0e-7, 1.0e-4, 1.0e-4, 1.0e-7, 1.0e-7, 1.0e-7, 1.0e-7 };
    AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, minStep, maxStep, absTolerance, relTolerance);
    integrator.setInitialStepSize(zero.add(initStep));
    // Numerical propagator based on the integrator
    FieldNumericalPropagator<T> propagator = new FieldNumericalPropagator<>(field, integrator);
    propagator.setOrbitType(type);
    T mass = field.getZero().add(1000.0);
    FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(geo, mass);
    propagator.setInitialState(initialState);
    propagator.setOrbitType(OrbitType.CARTESIAN);
    // Set the events Detectors
    FieldApsideDetector<T> event1 = new FieldApsideDetector<>(geo);
    propagator.addEventDetector(event1);
    // Set the propagation mode
    propagator.setSlaveMode();
    // Propagate
    FieldSpacecraftState<T> finalState = propagator.propagate(endDate);
    // we should stop long before endDate
    Assert.assertTrue(endDate.durationFrom(finalState.getDate()).getReal() > 40000.0);
}

Example 30 with DormandPrince853FieldIntegrator

use of org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator in project Orekit by CS-SI.

the class FieldNumericalPropagatorTest method doTestPropagationTypesHyperbolic.

private <T extends RealFieldElement<T>> void doTestPropagationTypesHyperbolic(Field<T> field) throws OrekitException {
    T zero = field.getZero();
    // setup
    final FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field);
    FieldSpacecraftState<T> initialState;
    FieldNumericalPropagator<T> propagator;
    final FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
    final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
    final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity), FramesFactory.getEME2000(), initDate, mu);
    initialState = new FieldSpacecraftState<>(orbit);
    OrbitType type = OrbitType.EQUINOCTIAL;
    double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit.toOrbit(), type);
    AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
    integrator.setInitialStepSize(zero.add(60));
    propagator = new FieldNumericalPropagator<>(field, integrator);
    propagator.setOrbitType(type);
    propagator.setInitialState(initialState);
    FieldSpacecraftState<T> state = new FieldSpacecraftState<>(new FieldKeplerianOrbit<>(zero.add(-10000000.0), zero.add(2.5), zero.add(0.3), zero, zero, zero, PositionAngle.TRUE, FramesFactory.getEME2000(), initDate, mu));
    ForceModel gravityField = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), GravityFieldFactory.getNormalizedProvider(5, 5));
    propagator.addForceModel(gravityField);
    // Propagation of the initial at t + dt
    final FieldPVCoordinates<T> pv = state.getPVCoordinates();
    final T dP = zero.add(0.001);
    final T dV = dP.multiply(state.getMu()).divide(pv.getPosition().getNormSq().multiply(pv.getVelocity().getNorm()));
    final FieldPVCoordinates<T> pvcM = propagateInType(state, dP, OrbitType.CARTESIAN, PositionAngle.MEAN, propagator);
    final FieldPVCoordinates<T> pvkM = propagateInType(state, dP, OrbitType.KEPLERIAN, PositionAngle.MEAN, propagator);
    final FieldPVCoordinates<T> pvcE = propagateInType(state, dP, OrbitType.CARTESIAN, PositionAngle.ECCENTRIC, propagator);
    final FieldPVCoordinates<T> pvkE = propagateInType(state, dP, OrbitType.KEPLERIAN, PositionAngle.ECCENTRIC, propagator);
    final FieldPVCoordinates<T> pvcT = propagateInType(state, dP, OrbitType.CARTESIAN, PositionAngle.TRUE, propagator);
    final FieldPVCoordinates<T> pvkT = propagateInType(state, dP, OrbitType.KEPLERIAN, PositionAngle.TRUE, propagator);
    Assert.assertEquals(0, pvcM.getPosition().subtract(pvkT.getPosition()).getNorm().getReal() / dP.getReal(), 0.3);
    Assert.assertEquals(0, pvcM.getVelocity().subtract(pvkT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.4);
    Assert.assertEquals(0, pvkM.getPosition().subtract(pvkT.getPosition()).getNorm().getReal() / dP.getReal(), 0.2);
    Assert.assertEquals(0, pvkM.getVelocity().subtract(pvkT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.3);
    Assert.assertEquals(0, pvcE.getPosition().subtract(pvkT.getPosition()).getNorm().getReal() / dP.getReal(), 0.3);
    Assert.assertEquals(0, pvcE.getVelocity().subtract(pvkT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.4);
    Assert.assertEquals(0, pvkE.getPosition().subtract(pvkT.getPosition()).getNorm().getReal() / dP.getReal(), 0.009);
    Assert.assertEquals(0, pvkE.getVelocity().subtract(pvkT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.006);
    Assert.assertEquals(0, pvcT.getPosition().subtract(pvkT.getPosition()).getNorm().getReal() / dP.getReal(), 0.3);
    Assert.assertEquals(0, pvcT.getVelocity().subtract(pvkT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.4);
}