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.propagation.numerical.NumericalPropagatorTest.java
@Test public void testEventDetectionBug() throws OrekitException, IOException, ParseException { TimeScale utc = TimeScalesFactory.getUTC(); AbsoluteDate initialDate = new AbsoluteDate(2005, 1, 1, 0, 0, 0.0, utc); double duration = 100000.; AbsoluteDate endDate = new AbsoluteDate(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))); Orbit geo = new CartesianOrbit( new PVCoordinates(new Vector3D(rApogee, 0., 0.), new Vector3D(0., vApogee, 0.)), EME2000, initialDate, mu);// www . j a v a 2 s .c om duration = geo.getKeplerianPeriod(); endDate = new AbsoluteDate(initialDate, duration); // Numerical Integration final double minStep = 0.001; final double maxStep = 1000; final double initStep = 60; 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 }; AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, absTolerance, relTolerance); integrator.setInitialStepSize(initStep); // Numerical propagator based on the integrator propagator = new NumericalPropagator(integrator); double mass = 1000.; SpacecraftState initialState = new SpacecraftState(geo, mass); propagator.setInitialState(initialState); propagator.setOrbitType(OrbitType.CARTESIAN); // Set the events Detectors ApsideDetector event1 = new ApsideDetector(geo); propagator.addEventDetector(event1); // Set the propagation mode propagator.setSlaveMode(); // Propagate SpacecraftState finalState = propagator.propagate(endDate); // we should stop long before endDate Assert.assertTrue(endDate.durationFrom(finalState.getDate()) > 40000.0); }
From source file:org.orekit.propagation.numerical.NumericalPropagatorTest.java
@Before public void setUp() throws OrekitException { Utils.setDataRoot("regular-data:potential/shm-format"); GravityFieldFactory.addPotentialCoefficientsReader(new SHMFormatReader("^eigen_cg03c_coef$", false)); mu = GravityFieldFactory.getUnnormalizedProvider(0, 0).getMu(); final Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6); final Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0); initDate = AbsoluteDate.J2000_EPOCH; final Orbit orbit = new EquinoctialOrbit(new PVCoordinates(position, velocity), FramesFactory.getEME2000(), initDate, mu);// www . j a va2s . c o m initialState = new SpacecraftState(orbit); double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit, OrbitType.EQUINOCTIAL); AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]); integrator.setInitialStepSize(60); propagator = new NumericalPropagator(integrator); propagator.setInitialState(initialState); }
From source file:org.orekit.propagation.numerical.PartialDerivativesTest.java
@Test public void testJacobianIssue18() throws OrekitException { // Body mu/*from ww w . jav a2s .c o m*/ final double mu = 3.9860047e14; 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); 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.setAttitudeProvider(law); propagator.addForceModel(maneuver); propagator.setOrbitType(OrbitType.CARTESIAN); PartialDerivativesEquations PDE = new PartialDerivativesEquations("derivatives", propagator); PDE.selectParamAndStep("thrust", Double.NaN); Assert.assertEquals(3, PDE.getAvailableParameters().size()); Assert.assertEquals("central attraction coefficient", PDE.getAvailableParameters().get(0)); Assert.assertEquals("thrust", PDE.getAvailableParameters().get(1)); Assert.assertEquals("flow rate", PDE.getAvailableParameters().get(2)); propagator.setInitialState(PDE.setInitialJacobians(initialState, 7, 1)); final AbsoluteDate finalDate = fireDate.shiftedBy(3800); final SpacecraftState finalorb = propagator.propagate(finalDate); Assert.assertEquals(0, finalDate.durationFrom(finalorb.getDate()), 1.0e-11); }
From source file:org.orekit.propagation.numerical.PartialDerivativesTest.java
private NumericalPropagator setUpPropagator(Orbit orbit, double dP, OrbitType orbitType, PositionAngle angleType, ForceModel... models) throws OrekitException { final double minStep = 0.001; final double maxStep = 1000; double[][] tol = NumericalPropagator.tolerances(dP, orbit, orbitType); NumericalPropagator propagator = new NumericalPropagator( new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1])); propagator.setOrbitType(orbitType);/* ww w . j av a2 s. c o m*/ propagator.setPositionAngleType(angleType); for (ForceModel model : models) { propagator.addForceModel(model); } return propagator; }
From source file:org.orekit.propagation.semianalytical.dsst.DSSTPropagatorTest.java
@Test public void testEphemerisDates() throws OrekitException { //setup/*from ww w .java 2s . co m*/ TimeScale tai = TimeScalesFactory.getTAI(); AbsoluteDate initialDate = new AbsoluteDate("2015-07-01", tai); AbsoluteDate startDate = new AbsoluteDate("2015-07-03", tai).shiftedBy(-0.1); AbsoluteDate endDate = new AbsoluteDate("2015-07-04", tai); Frame eci = FramesFactory.getGCRF(); KeplerianOrbit orbit = new KeplerianOrbit(600e3 + Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 0, 0, 0, 0, 0, PositionAngle.TRUE, eci, initialDate, Constants.EIGEN5C_EARTH_MU); double[][] tol = DSSTPropagator.tolerances(1, orbit); Propagator prop = new DSSTPropagator(new DormandPrince853Integrator(0.1, 500, tol[0], tol[1])); prop.resetInitialState(new SpacecraftState(new CartesianOrbit(orbit))); //action prop.setEphemerisMode(); prop.propagate(startDate, endDate); BoundedPropagator ephemeris = prop.getGeneratedEphemeris(); //verify TimeStampedPVCoordinates actualPV = ephemeris.getPVCoordinates(startDate, eci); TimeStampedPVCoordinates expectedPV = orbit.getPVCoordinates(startDate, eci); MatcherAssert.assertThat(actualPV.getPosition(), OrekitMatchers.vectorCloseTo(expectedPV.getPosition(), 1.0)); MatcherAssert.assertThat(actualPV.getVelocity(), OrekitMatchers.vectorCloseTo(expectedPV.getVelocity(), 1.0)); MatcherAssert.assertThat(ephemeris.getMinDate().durationFrom(startDate), OrekitMatchers.closeTo(0, 0)); MatcherAssert.assertThat(ephemeris.getMaxDate().durationFrom(endDate), OrekitMatchers.closeTo(0, 0)); //test date AbsoluteDate date = endDate.shiftedBy(-0.11); Assert.assertEquals(ephemeris.propagate(date).getDate().durationFrom(date), 0, 0); }
From source file:org.orekit.propagation.semianalytical.dsst.DSSTPropagatorTest.java
@Test public void testIssue157() 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());/*w ww .j a va 2 s . c o m*/ 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, true); 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 DSSTCentralBody(earth.getBodyFrame(), Constants.WGS84_EARTH_ANGULAR_VELOCITY, nshp)); 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.resetInitialState(new SpacecraftState(orbit, 45.0)); final SpacecraftState finalState = propagator .propagate(orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY)); Assert.assertEquals(2189.4, orbit.getA() - finalState.getA(), 10.0); }
From source file:org.orekit.propagation.semianalytical.dsst.DSSTPropagatorTest.java
@Test public void testDSSTrestart() throws OrekitException { DSSTPropagator dsstProp;/*w w w. j a v a2s.com*/ // build force model geopotential 8x8 Utils.setDataRoot("regular-data:potential/icgem-format"); GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("^eigen-6s-truncated$", false)); final UnnormalizedSphericalHarmonicsProvider gravityProvider = GravityFieldFactory .getUnnormalizedProvider(8, 8); final Frame rotatingFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame(); DSSTForceModel gravityForceModel = new DSSTCentralBody(rotatingFrame, Constants.WGS84_EARTH_ANGULAR_VELOCITY, gravityProvider); // build initial state final AbsoluteDate epochDate = new AbsoluteDate(2014, 01, 01, 0, 0, 0, TimeScalesFactory.getUTC()); final KeplerianOrbit initialOrbit = new KeplerianOrbit(26562000.0, 0.72, FastMath.toRadians(63.435), FastMath.toRadians(270.0), 0.0, 0.0, PositionAngle.MEAN, FramesFactory.getEME2000(), epochDate, gravityProvider.getMu()); final SpacecraftState initialState = new SpacecraftState(new EquinoctialOrbit(initialOrbit)); // build integrator final double minStep = initialState.getKeplerianPeriod() * 0.1; final double maxStep = initialState.getKeplerianPeriod() * 10.0; final double[][] tol = DSSTPropagator.tolerances(0.1, initialState.getOrbit()); AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]); dsstProp = new DSSTPropagator(integrator); // add force model dsstProp.addForceModel(gravityForceModel); // DSST Propagation (first propagation without timing, for warm-up purposes) dsstProp.setInitialState(initialState, false); dsstProp.propagate(epochDate.shiftedBy(100.0 * 86400.0)); double refExecTime = 0; for (int i = 0; i < 5; i++) { dsstProp.setInitialState(initialState, false); long propStart = System.currentTimeMillis(); dsstProp.propagate(epochDate.shiftedBy(100.0 * 86400.0)); long propEnd = System.currentTimeMillis(); double execTime = 0.001 * (propEnd - propStart); if (refExecTime <= 0) { refExecTime = execTime; } else { Assert.assertTrue(execTime <= refExecTime * 1.5); } } }
From source file:org.orekit.propagation.semianalytical.dsst.DSSTPropagatorTest.java
@Test public void testGetInitialOsculatingState() throws IllegalArgumentException, OrekitException { final SpacecraftState initialState = getGEOrbit(); // build integrator final double minStep = initialState.getKeplerianPeriod() * 0.1; final double maxStep = initialState.getKeplerianPeriod() * 10.0; final double[][] tol = DSSTPropagator.tolerances(0.1, initialState.getOrbit()); AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]); DSSTPropagator prop = new DSSTPropagator(integrator, false); final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(4, 0); final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame(); DSSTForceModel force = new DSSTCentralBody(earthFrame, Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider); prop.addForceModel(force);// w w w. j av a2s. c o m prop.setInitialState(initialState, false); prop.getInitialState(); }
From source file:org.orekit.propagation.semianalytical.dsst.DSSTPropagatorTest.java
private void setDSSTProp(SpacecraftState initialState) throws OrekitException { initialState.getDate();/*from w w w . j a va 2s.c o m*/ 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); }
From source file:org.orekit.utils.AngularCoordinatesTest.java
@Test public void testShiftWithAcceleration() throws OrekitException { double rate = 2 * FastMath.PI / (12 * 60); double acc = 0.001; double dt = 1.0; int n = 2000; final AngularCoordinates quadratic = new AngularCoordinates(Rotation.IDENTITY, new Vector3D(rate, Vector3D.PLUS_K), new Vector3D(acc, Vector3D.PLUS_J)); final AngularCoordinates linear = new AngularCoordinates(quadratic.getRotation(), quadratic.getRotationRate(), Vector3D.ZERO); final FirstOrderDifferentialEquations ode = new FirstOrderDifferentialEquations() { public int getDimension() { return 4; }/*w w w. j a v a 2 s . c o m*/ public void computeDerivatives(final double t, final double[] q, final double[] qDot) { final double omegaX = quadratic.getRotationRate().getX() + t * quadratic.getRotationAcceleration().getX(); final double omegaY = quadratic.getRotationRate().getY() + t * quadratic.getRotationAcceleration().getY(); final double omegaZ = quadratic.getRotationRate().getZ() + t * quadratic.getRotationAcceleration().getZ(); qDot[0] = 0.5 * MathArrays.linearCombination(-q[1], omegaX, -q[2], omegaY, -q[3], omegaZ); qDot[1] = 0.5 * MathArrays.linearCombination(q[0], omegaX, -q[3], omegaY, q[2], omegaZ); qDot[2] = 0.5 * MathArrays.linearCombination(q[3], omegaX, q[0], omegaY, -q[1], omegaZ); qDot[3] = 0.5 * MathArrays.linearCombination(-q[2], omegaX, q[1], omegaY, q[0], omegaZ); } }; FirstOrderIntegrator integrator = new DormandPrince853Integrator(1.0e-6, 1.0, 1.0e-12, 1.0e-12); integrator.addStepHandler(new StepNormalizer(dt / n, new FixedStepHandler() { public void init(double t0, double[] y0, double t) { } public void handleStep(double t, double[] y, double[] yDot, boolean isLast) { Rotation reference = new Rotation(y[0], y[1], y[2], y[3], true); // the error in shiftedBy taking acceleration into account is cubic double expectedCubicError = 1.4544e-6 * t * t * t; Assert.assertEquals(expectedCubicError, Rotation.distance(reference, quadratic.shiftedBy(t).getRotation()), 0.0001 * expectedCubicError); // the error in shiftedBy not taking acceleration into account is quadratic double expectedQuadraticError = 5.0e-4 * t * t; Assert.assertEquals(expectedQuadraticError, Rotation.distance(reference, linear.shiftedBy(t).getRotation()), 0.00001 * expectedQuadraticError); } })); double[] y = new double[] { quadratic.getRotation().getQ0(), quadratic.getRotation().getQ1(), quadratic.getRotation().getQ2(), quadratic.getRotation().getQ3() }; integrator.integrate(ode, 0, y, dt, y); }