List of usage examples for org.apache.commons.math3.ode.nonstiff DormandPrince853Integrator DormandPrince853Integrator
public DormandPrince853Integrator(final double minStep, final double maxStep, final double[] vecAbsoluteTolerance, final double[] vecRelativeTolerance)
From source file:org.orekit.forces.gravity.ThirdBodyAttractionTest.java
@Test public void testStateJacobian() throws OrekitException { // initialization AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01), new TimeComponents(13, 59, 27.816), TimeScalesFactory.getUTC()); double i = FastMath.toRadians(98.7); double omega = FastMath.toRadians(93.0); double OMEGA = FastMath.toRadians(15.0 * 22.5); Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i, omega, OMEGA, 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date, Constants.EIGEN5C_EARTH_MU); OrbitType integrationType = OrbitType.CARTESIAN; double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType); NumericalPropagator propagator = new NumericalPropagator( new DormandPrince853Integrator(1.0e-3, 120, tolerances[0], tolerances[1])); propagator.setOrbitType(integrationType); final CelestialBody moon = CelestialBodyFactory.getMoon(); final ThirdBodyAttraction forceModel = new ThirdBodyAttraction(moon); propagator.addForceModel(forceModel); SpacecraftState state0 = new SpacecraftState(orbit); checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0), 1e4, tolerances[0], 2.0e-9); }
From source file:org.orekit.forces.maneuvers.ConstantThrustManeuverTest.java
@Test public void testRoughBehaviour() throws OrekitException { final double isp = 318; final double mass = 2500; final double a = 24396159; final double e = 0.72831215; final double i = FastMath.toRadians(7); final double omega = FastMath.toRadians(180); final double OMEGA = FastMath.toRadians(261); final double lv = 0; final double duration = 3653.99; final double f = 420; final double delta = FastMath.toRadians(-7.4978); final double alpha = FastMath.toRadians(351); final AttitudeProvider law = new InertialProvider( new Rotation(new Vector3D(alpha, delta), Vector3D.PLUS_I)); final AbsoluteDate initDate = new AbsoluteDate(new DateComponents(2004, 01, 01), new TimeComponents(23, 30, 00.000), TimeScalesFactory.getUTC()); final Orbit orbit = new KeplerianOrbit(a, e, i, omega, OMEGA, lv, PositionAngle.TRUE, FramesFactory.getEME2000(), initDate, mu); final SpacecraftState initialState = new SpacecraftState(orbit, law.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), mass); final AbsoluteDate fireDate = new AbsoluteDate(new DateComponents(2004, 01, 02), new TimeComponents(04, 15, 34.080), TimeScalesFactory.getUTC()); final ConstantThrustManeuver maneuver = new ConstantThrustManeuver(fireDate, duration, f, isp, Vector3D.PLUS_I);/*from www.jav a 2 s.c o m*/ Assert.assertEquals(f, maneuver.getThrust(), 1.0e-10); Assert.assertEquals(isp, maneuver.getISP(), 1.0e-10); double[] absTolerance = { 0.001, 1.0e-9, 1.0e-9, 1.0e-6, 1.0e-6, 1.0e-6, 0.001 }; double[] relTolerance = { 1.0e-7, 1.0e-4, 1.0e-4, 1.0e-7, 1.0e-7, 1.0e-7, 1.0e-7 }; AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 1000, absTolerance, relTolerance); integrator.setInitialStepSize(60); final NumericalPropagator propagator = new NumericalPropagator(integrator); propagator.setInitialState(initialState); propagator.setAttitudeProvider(law); propagator.addForceModel(maneuver); final SpacecraftState finalorb = propagator.propagate(fireDate.shiftedBy(3800)); final double massTolerance = FastMath.abs(maneuver.getFlowRate()) * maneuver.getEventsDetectors()[0].getThreshold(); Assert.assertEquals(2007.8824544261233, finalorb.getMass(), massTolerance); Assert.assertEquals(2.6872, FastMath.toDegrees(MathUtils.normalizeAngle(finalorb.getI(), FastMath.PI)), 1e-4); Assert.assertEquals(28970, finalorb.getA() / 1000, 1); }
From source file:org.orekit.forces.maneuvers.SmallManeuverAnalyticalModelTest.java
private BoundedPropagator getEphemeris(final Orbit orbit, final double mass, final AbsoluteDate t0, final Vector3D dV, final double f, final double isp) throws OrekitException { AttitudeProvider law = new LofOffset(orbit.getFrame(), LOFType.LVLH); final SpacecraftState initialState = new SpacecraftState(orbit, law.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), mass); // set up numerical propagator final double dP = 1.0; double[][] tolerances = NumericalPropagator.tolerances(dP, orbit, orbit.getType()); AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 1000, tolerances[0], tolerances[1]);/*from ww w . ja v a2 s. co m*/ integrator.setInitialStepSize(orbit.getKeplerianPeriod() / 100.0); final NumericalPropagator propagator = new NumericalPropagator(integrator); propagator.setOrbitType(orbit.getType()); propagator.setInitialState(initialState); propagator.setAttitudeProvider(law); if (dV.getNorm() > 1.0e-6) { // set up maneuver final double vExhaust = Constants.G0_STANDARD_GRAVITY * isp; final double dt = -(mass * vExhaust / f) * FastMath.expm1(-dV.getNorm() / vExhaust); final ConstantThrustManeuver maneuver = new ConstantThrustManeuver(t0, dt, f, isp, dV.normalize()); propagator.addForceModel(maneuver); } propagator.setEphemerisMode(); propagator.propagate(t0.shiftedBy(5 * orbit.getKeplerianPeriod())); return propagator.getGeneratedEphemeris(); }
From source file:org.orekit.forces.radiation.SolarRadiationPressureTest.java
@Test public void testStateJacobianSphere() throws OrekitException { // initialization AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01), new TimeComponents(13, 59, 27.816), TimeScalesFactory.getUTC()); double i = FastMath.toRadians(98.7); double omega = FastMath.toRadians(93.0); double OMEGA = FastMath.toRadians(15.0 * 22.5); Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i, omega, OMEGA, 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date, Constants.EIGEN5C_EARTH_MU); OrbitType integrationType = OrbitType.CARTESIAN; double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType); NumericalPropagator propagator = new NumericalPropagator( new DormandPrince853Integrator(1.0e-3, 120, tolerances[0], tolerances[1])); propagator.setOrbitType(integrationType); SolarRadiationPressure forceModel = new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS, new SphericalSpacecraft(2.5, 1.2, 0.7, 0.2)); propagator.addForceModel(forceModel); SpacecraftState state0 = new SpacecraftState(orbit); checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0), 1e3, tolerances[0], 2.0e-6); }
From source file:org.orekit.forces.radiation.SolarRadiationPressureTest.java
@Test public void testStateJacobianBox() throws OrekitException { // initialization AbsoluteDate date = new AbsoluteDate(new DateComponents(2003, 03, 01), new TimeComponents(13, 59, 27.816), TimeScalesFactory.getUTC()); double i = FastMath.toRadians(98.7); double omega = FastMath.toRadians(93.0); double OMEGA = FastMath.toRadians(15.0 * 22.5); Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i, omega, OMEGA, 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date, Constants.EIGEN5C_EARTH_MU); OrbitType integrationType = OrbitType.CARTESIAN; double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, integrationType); NumericalPropagator propagator = new NumericalPropagator( new DormandPrince853Integrator(1.0e-3, 120, tolerances[0], tolerances[1])); propagator.setOrbitType(integrationType); SolarRadiationPressure forceModel = new SolarRadiationPressure(CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS, new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0, Vector3D.PLUS_J, 1.2, 0.7, 0.2)); propagator.addForceModel(forceModel); SpacecraftState state0 = new SpacecraftState(orbit); checkStateJacobian(propagator, state0, date.shiftedBy(3.5 * 3600.0), 1e3, tolerances[0], 2.0e-5); }
From source file:org.orekit.forces.radiation.SolarRadiationPressureTest.java
@Test public void testRoughOrbitalModifs() throws ParseException, OrekitException, FileNotFoundException { // initialization AbsoluteDate date = new AbsoluteDate(new DateComponents(1970, 7, 1), new TimeComponents(13, 59, 27.816), TimeScalesFactory.getUTC()); Orbit orbit = new EquinoctialOrbit(42164000, 10e-3, 10e-3, FastMath.tan(0.001745329) * FastMath.cos(2 * FastMath.PI / 3), FastMath.tan(0.001745329) * FastMath.sin(2 * FastMath.PI / 3), 0.1, PositionAngle.TRUE, FramesFactory.getEME2000(), date, mu); final double period = orbit.getKeplerianPeriod(); Assert.assertEquals(86164, period, 1); PVCoordinatesProvider sun = CelestialBodyFactory.getSun(); // creation of the force model OneAxisEllipsoid earth = new OneAxisEllipsoid(6378136.46, 1.0 / 298.25765, FramesFactory.getITRF(IERSConventions.IERS_2010, true)); SolarRadiationPressure SRP = new SolarRadiationPressure(sun, earth.getEquatorialRadius(), (RadiationSensitive) new SphericalSpacecraft(500.0, 0.7, 0.7, 0.7)); // creation of the propagator double[] absTolerance = { 0.1, 1.0e-9, 1.0e-9, 1.0e-5, 1.0e-5, 1.0e-5, 0.001 }; double[] relTolerance = { 1.0e-4, 1.0e-4, 1.0e-4, 1.0e-6, 1.0e-6, 1.0e-6, 1.0e-7 }; AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(900.0, 60000, absTolerance, relTolerance);/*from w w w . j av a2 s . c om*/ integrator.setInitialStepSize(3600); final NumericalPropagator calc = new NumericalPropagator(integrator); calc.addForceModel(SRP); // Step Handler calc.setMasterMode(FastMath.floor(period), new SolarStepHandler()); AbsoluteDate finalDate = date.shiftedBy(10 * period); calc.setInitialState(new SpacecraftState(orbit, 1500.0)); calc.propagate(finalDate); Assert.assertTrue(calc.getCalls() < 7100); }
From source file:org.orekit.propagation.analytical.AdapterPropagatorTest.java
private BoundedPropagator getEphemeris(final Orbit orbit, final double mass, final int nbOrbits, final AttitudeProvider law, final AbsoluteDate t0, final Vector3D dV, final double f, final double isp, final boolean sunAttraction, final boolean moonAttraction, final NormalizedSphericalHarmonicsProvider gravityField) throws OrekitException, ParseException, IOException { SpacecraftState initialState = new SpacecraftState(orbit, law.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), mass); // add some dummy additional states initialState = initialState.addAdditionalState("dummy 1", 1.25, 2.5); initialState = initialState.addAdditionalState("dummy 2", 5.0); // set up numerical propagator final double dP = 1.0; double[][] tolerances = NumericalPropagator.tolerances(dP, orbit, orbit.getType()); AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 1000, tolerances[0], tolerances[1]);// w w w . ja v a2 s.com integrator.setInitialStepSize(orbit.getKeplerianPeriod() / 100.0); final NumericalPropagator propagator = new NumericalPropagator(integrator); propagator.addAdditionalStateProvider(new AdditionalStateProvider() { public String getName() { return "dummy 2"; } public double[] getAdditionalState(SpacecraftState state) { return new double[] { 5.0 }; } }); propagator.setInitialState(initialState); propagator.setAttitudeProvider(law); if (dV.getNorm() > 1.0e-6) { // set up maneuver final double vExhaust = Constants.G0_STANDARD_GRAVITY * isp; final double dt = -(mass * vExhaust / f) * FastMath.expm1(-dV.getNorm() / vExhaust); final ConstantThrustManeuver maneuver = new ConstantThrustManeuver(t0.shiftedBy(-0.5 * dt), dt, f, isp, dV.normalize()); propagator.addForceModel(maneuver); } if (sunAttraction) { propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun())); } if (moonAttraction) { propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon())); } if (gravityField != null) { propagator.addForceModel( new HolmesFeatherstoneAttractionModel(FramesFactory.getGTOD(false), gravityField)); } propagator.setEphemerisMode(); propagator.propagate(t0.shiftedBy(nbOrbits * orbit.getKeplerianPeriod())); final BoundedPropagator ephemeris = propagator.getGeneratedEphemeris(); // both the initial propagator and generated ephemeris manage one of the two // additional states, but they also contain unmanaged copies of the other one Assert.assertFalse(propagator.isAdditionalStateManaged("dummy 1")); Assert.assertTrue(propagator.isAdditionalStateManaged("dummy 2")); Assert.assertFalse(ephemeris.isAdditionalStateManaged("dummy 1")); Assert.assertTrue(ephemeris.isAdditionalStateManaged("dummy 2")); Assert.assertEquals(2, ephemeris.getInitialState().getAdditionalState("dummy 1").length); Assert.assertEquals(1, ephemeris.getInitialState().getAdditionalState("dummy 2").length); return ephemeris; }
From source file:org.orekit.propagation.analytical.EcksteinHechlerPropagatorTest.java
@Test public void testInitializationCorrectness() throws OrekitException, IOException { // Definition of initial conditions AbsoluteDate date = AbsoluteDate.J2000_EPOCH.shiftedBy(154.); Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true); Frame eme2000 = FramesFactory.getEME2000(); Vector3D pole = itrf.getTransformTo(eme2000, date).transformVector(Vector3D.PLUS_K); Frame poleAligned = new Frame(FramesFactory.getEME2000(), new Transform(date, new Rotation(pole, Vector3D.PLUS_K)), "pole aligned", true); CircularOrbit initial = new CircularOrbit(7208669.8179538045, 1.3740461966386876E-4, -3.2364250248363356E-5, FastMath.toRadians(97.40236024565775), FastMath.toRadians(166.15873160992115), FastMath.toRadians(90.1282370098961), PositionAngle.MEAN, poleAligned, date, provider.getMu()); // find the default Eckstein-Hechler propagator initialized from the initial orbit EcksteinHechlerPropagator defaultEH = new EcksteinHechlerPropagator(initial, provider); // the osculating parameters recomputed by the default Eckstein-Hechler propagator are quite different // from initial orbit CircularOrbit defaultOrbit = (CircularOrbit) OrbitType.CIRCULAR .convertType(defaultEH.propagateOrbit(initial.getDate())); Assert.assertEquals(267.4, defaultOrbit.getA() - initial.getA(), 0.1); // the position on the other hand match perfectly Assert.assertEquals(0.0, Vector3D.distance(defaultOrbit.getPVCoordinates().getPosition(), initial.getPVCoordinates().getPosition()), 1.0e-8); // set up a reference numerical propagator starting for the specified start orbit // using the same force models (i.e. the first few zonal terms) double[][] tol = NumericalPropagator.tolerances(0.1, initial, OrbitType.CIRCULAR); AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 1000, tol[0], tol[1]); integrator.setInitialStepSize(60);//from www .jav a 2 s . c o m NumericalPropagator num = new NumericalPropagator(integrator); num.addForceModel( new HolmesFeatherstoneAttractionModel(itrf, GravityFieldFactory.getNormalizedProvider(provider))); num.setInitialState(new SpacecraftState(initial)); num.setOrbitType(OrbitType.CIRCULAR); // find the best Eckstein-Hechler propagator that match the orbit evolution PropagatorConverter converter = new FiniteDifferencePropagatorConverter( new EcksteinHechlerPropagatorBuilder(poleAligned, provider, OrbitType.CIRCULAR, PositionAngle.TRUE), 1.0e-6, 100); EcksteinHechlerPropagator fittedEH = (EcksteinHechlerPropagator) converter.convert(num, 3 * initial.getKeplerianPeriod(), 300); // the default Eckstein-Hechler propagator did however quite a good job, as it found // an orbit close to the best fitting CircularOrbit fittedOrbit = (CircularOrbit) OrbitType.CIRCULAR .convertType(fittedEH.propagateOrbit(initial.getDate())); Assert.assertEquals(0.623, defaultOrbit.getA() - fittedOrbit.getA(), 0.1); // the position on the other hand are slightly different // because the fitted orbit minimizes the residuals over a complete time span, // not on a single point Assert.assertEquals(58.0, Vector3D.distance(defaultOrbit.getPVCoordinates().getPosition(), fittedOrbit.getPVCoordinates().getPosition()), 0.1); }
From source file:org.orekit.propagation.conversion.DormandPrince853IntegratorBuilder.java
/** {@inheritDoc} */ public AbstractIntegrator buildIntegrator(final Orbit orbit) throws PropagationException { final double[][] tol = NumericalPropagator.tolerances(dP, orbit, OrbitType.CARTESIAN); return new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]); }
From source file:org.orekit.propagation.conversion.NumericalConverterTest.java
@Before public void setUp() throws OrekitException, IOException, ParseException { Utils.setDataRoot("regular-data:potential/shm-format"); gravity = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), GravityFieldFactory.getNormalizedProvider(2, 0)); mu = gravity.getParameter(NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT); minStep = 0.001;/*from w w w .j a v a 2 s. c o m*/ maxStep = 200.0; dP = 0.01; //use a orbit that comes close to Earth so the drag coefficient has an effect final Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6).normalize() .scalarMultiply(Constants.WGS84_EARTH_EQUATORIAL_RADIUS + 300e3); final Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0); final AbsoluteDate initDate = new AbsoluteDate(2010, 10, 10, 10, 10, 10.0, TimeScalesFactory.getUTC()); orbit = new EquinoctialOrbit(new PVCoordinates(position, velocity), FramesFactory.getEME2000(), initDate, mu); final double[][] tol = NumericalPropagator.tolerances(dP, orbit, OrbitType.CARTESIAN); propagator = new NumericalPropagator(new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1])); propagator.setInitialState(new SpacecraftState(orbit)); propagator.setOrbitType(OrbitType.CARTESIAN); final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true)); earth.setAngularThreshold(1.e-7); final Atmosphere atmosphere = new SimpleExponentialAtmosphere(earth, 0.0004, 42000.0, 7500.0); final double dragCoef = 2.0; final SphericalSpacecraft ss = new SphericalSpacecraft(10., dragCoef, 0., 0.); drag = new DragForce(atmosphere, ss); propagator.addForceModel(gravity); propagator.addForceModel(drag); }