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

Example 16 with DormandPrince853Integrator

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

the class DSSTPropagatorTest method setDSSTProp.

private void setDSSTProp(SpacecraftState initialState) throws OrekitException {
    initialState.getDate();
    final double minStep = initialState.getKeplerianPeriod();
    final double maxStep = 100. * minStep;
    final double[][] tol = DSSTPropagator.tolerances(1.0, initialState.getOrbit());
    AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
    dsstProp = new DSSTPropagator(integrator);
    dsstProp.setInitialState(initialState, false);
}

Example 17 with DormandPrince853Integrator

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

the class DSSTPropagatorTest method testEphemerisGeneration.

@Test
public void testEphemerisGeneration() throws OrekitException {
    Utils.setDataRoot("regular-data:potential/icgem-format");
    GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("^eigen-6s-truncated$", false));
    UnnormalizedSphericalHarmonicsProvider nshp = GravityFieldFactory.getUnnormalizedProvider(8, 8);
    Orbit orbit = new KeplerianOrbit(13378000, 0.05, 0, 0, FastMath.PI, 0, PositionAngle.MEAN, FramesFactory.getTOD(false), new AbsoluteDate(2003, 5, 6, TimeScalesFactory.getUTC()), nshp.getMu());
    double period = orbit.getKeplerianPeriod();
    double[][] tolerance = DSSTPropagator.tolerances(1.0, orbit);
    AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(period / 100, period * 100, tolerance[0], tolerance[1]);
    integrator.setInitialStepSize(10 * period);
    DSSTPropagator propagator = new DSSTPropagator(integrator, false);
    OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getGTOD(false));
    CelestialBody sun = CelestialBodyFactory.getSun();
    CelestialBody moon = CelestialBodyFactory.getMoon();
    propagator.addForceModel(new DSSTZonal(nshp, 8, 7, 17));
    propagator.addForceModel(new DSSTTesseral(earth.getBodyFrame(), Constants.WGS84_EARTH_ANGULAR_VELOCITY, nshp, 8, 8, 4, 12, 8, 8, 4));
    propagator.addForceModel(new DSSTThirdBody(sun));
    propagator.addForceModel(new DSSTThirdBody(moon));
    propagator.addForceModel(new DSSTAtmosphericDrag(new HarrisPriester(sun, earth), 2.1, 180));
    propagator.addForceModel(new DSSTSolarRadiationPressure(1.2, 180, sun, earth.getEquatorialRadius()));
    propagator.setInterpolationGridToMaxTimeGap(0.5 * Constants.JULIAN_DAY);
    // direct generation of states
    propagator.setInitialState(new SpacecraftState(orbit, 45.0), false);
    final List<SpacecraftState> states = new ArrayList<SpacecraftState>();
    propagator.setMasterMode(600, (currentState, isLast) -> states.add(currentState));
    propagator.propagate(orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY));
    // ephemeris generation
    propagator.setInitialState(new SpacecraftState(orbit, 45.0), false);
    propagator.setEphemerisMode();
    propagator.propagate(orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY));
    BoundedPropagator ephemeris = propagator.getGeneratedEphemeris();
    double maxError = 0;
    for (final SpacecraftState state : states) {
        final SpacecraftState fromEphemeris = ephemeris.propagate(state.getDate());
        final double error = Vector3D.distance(state.getPVCoordinates().getPosition(), fromEphemeris.getPVCoordinates().getPosition());
        maxError = FastMath.max(maxError, error);
    }
    Assert.assertEquals(0.0, maxError, 1.0e-10);
}

Example 18 with DormandPrince853Integrator

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

the class SecularAndHarmonicTest method createPropagator.

private NumericalPropagator createPropagator(CircularOrbit orbit) throws OrekitException {
    OrbitType type = OrbitType.CIRCULAR;
    double[][] tolerances = NumericalPropagator.tolerances(0.1, orbit, type);
    DormandPrince853Integrator integrator = new DormandPrince853Integrator(1.0, 600, tolerances[0], tolerances[1]);
    integrator.setInitialStepSize(60.0);
    NumericalPropagator propagator = new NumericalPropagator(integrator);
    propagator.addForceModel(new HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravityField));
    propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
    propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
    propagator.setOrbitType(type);
    propagator.resetInitialState(new SpacecraftState(orbit));
    return propagator;
}

Example 19 with DormandPrince853Integrator

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

the class TimeStampedAngularCoordinatesTest method interpolationErrors.

private double[] interpolationErrors(final TimeStampedAngularCoordinates reference, double dt) throws OrekitException {
    final OrdinaryDifferentialEquation ode = new OrdinaryDifferentialEquation() {

        public int getDimension() {
            return 4;
        }

        public double[] computeDerivatives(final double t, final double[] q) {
            final double omegaX = reference.getRotationRate().getX() + t * reference.getRotationAcceleration().getX();
            final double omegaY = reference.getRotationRate().getY() + t * reference.getRotationAcceleration().getY();
            final double omegaZ = reference.getRotationRate().getZ() + t * reference.getRotationAcceleration().getZ();
            return new double[] { 0.5 * MathArrays.linearCombination(-q[1], omegaX, -q[2], omegaY, -q[3], omegaZ), 0.5 * MathArrays.linearCombination(q[0], omegaX, -q[3], omegaY, q[2], omegaZ), 0.5 * MathArrays.linearCombination(q[3], omegaX, q[0], omegaY, -q[1], omegaZ), 0.5 * MathArrays.linearCombination(-q[2], omegaX, q[1], omegaY, q[0], omegaZ) };
        }
    };
    final List<TimeStampedAngularCoordinates> complete = new ArrayList<TimeStampedAngularCoordinates>();
    ODEIntegrator integrator = new DormandPrince853Integrator(1.0e-6, 1.0, 1.0e-12, 1.0e-12);
    integrator.addStepHandler(new StepNormalizer(dt / 2000, new ODEFixedStepHandler() {

        public void handleStep(ODEStateAndDerivative state, boolean isLast) {
            final double t = state.getTime();
            final double[] q = state.getPrimaryState();
            complete.add(new TimeStampedAngularCoordinates(reference.getDate().shiftedBy(t), new Rotation(q[0], q[1], q[2], q[3], true), new Vector3D(1, reference.getRotationRate(), t, reference.getRotationAcceleration()), reference.getRotationAcceleration()));
        }
    }));
    double[] y = new double[] { reference.getRotation().getQ0(), reference.getRotation().getQ1(), reference.getRotation().getQ2(), reference.getRotation().getQ3() };
    integrator.integrate(ode, new ODEState(0, y), dt);
    List<TimeStampedAngularCoordinates> sample = new ArrayList<TimeStampedAngularCoordinates>();
    sample.add(complete.get(0));
    sample.add(complete.get(complete.size() / 2));
    sample.add(complete.get(complete.size() - 1));
    double maxRotationError = 0;
    double maxRateError = 0;
    double maxAccelerationError = 0;
    for (TimeStampedAngularCoordinates acRef : complete) {
        TimeStampedAngularCoordinates interpolated = TimeStampedAngularCoordinates.interpolate(acRef.getDate(), AngularDerivativesFilter.USE_RRA, sample);
        double rotationError = Rotation.distance(acRef.getRotation(), interpolated.getRotation());
        double rateError = Vector3D.distance(acRef.getRotationRate(), interpolated.getRotationRate());
        double accelerationError = Vector3D.distance(acRef.getRotationAcceleration(), interpolated.getRotationAcceleration());
        maxRotationError = FastMath.max(maxRotationError, rotationError);
        maxRateError = FastMath.max(maxRateError, rateError);
        maxAccelerationError = FastMath.max(maxAccelerationError, accelerationError);
    }
    return new double[] { maxRotationError, maxRateError, maxAccelerationError };
}

Example 20 with DormandPrince853Integrator

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

the class RelativityTest method testSmallEffectOnOrbit.

/**
 * check against example in Tapley, Schutz, and Born, p 65-66. They predict a
 * progression of perigee of 11 arcsec/year. To get the same results we must set the
 * propagation tolerances to 1e-5.
 *
 * @throws OrekitException on error
 */
@Test
public void testSmallEffectOnOrbit() throws OrekitException {
    // setup
    final double gm = Constants.EIGEN5C_EARTH_MU;
    Orbit orbit = new KeplerianOrbit(7500e3, 0.025, FastMath.toRadians(41.2), 0, 0, 0, PositionAngle.TRUE, frame, date, gm);
    double[][] tol = NumericalPropagator.tolerances(0.00001, orbit, OrbitType.CARTESIAN);
    AbstractIntegrator integrator = new DormandPrince853Integrator(1, 3600, tol[0], tol[1]);
    NumericalPropagator propagator = new NumericalPropagator(integrator);
    propagator.setOrbitType(OrbitType.CARTESIAN);
    propagator.addForceModel(new Relativity(gm));
    propagator.setInitialState(new SpacecraftState(orbit));
    // action: propagate a period
    AbsoluteDate end = orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY);
    PVCoordinates actual = propagator.getPVCoordinates(end, frame);
    // verify
    KeplerianOrbit endOrbit = new KeplerianOrbit(actual, frame, end, gm);
    KeplerianOrbit startOrbit = new KeplerianOrbit(orbit);
    double dp = endOrbit.getPerigeeArgument() - startOrbit.getPerigeeArgument();
    double dtYears = end.durationFrom(orbit.getDate()) / Constants.JULIAN_YEAR;
    double dpDeg = FastMath.toDegrees(dp);
    // change in argument of perigee in arcseconds per year
    double arcsecPerYear = dpDeg * 3600 / dtYears;
    Assert.assertEquals(11, arcsecPerYear, 0.5);
}