use of org.hipparchus.ode.ODEIntegrator in project Orekit by CS-SI.
the class NumericalPropagatorTest method createPropagator.
private static synchronized NumericalPropagator createPropagator(SpacecraftState spacecraftState, OrbitType orbitType, PositionAngle angleType) throws OrekitException {
final double minStep = 0.001;
final double maxStep = 120.0;
final double positionTolerance = 0.1;
final int degree = 20;
final int order = 20;
final double spacecraftArea = 1.0;
final double spacecraftDragCoefficient = 2.0;
final double spacecraftReflectionCoefficient = 2.0;
// propagator main configuration
final double[][] tol = NumericalPropagator.tolerances(positionTolerance, spacecraftState.getOrbit(), orbitType);
final ODEIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
final NumericalPropagator np = new NumericalPropagator(integrator);
np.setOrbitType(orbitType);
np.setPositionAngleType(angleType);
np.setInitialState(spacecraftState);
// Earth gravity field
final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
final NormalizedSphericalHarmonicsProvider harmonicsGravityProvider = GravityFieldFactory.getNormalizedProvider(degree, order);
np.addForceModel(new HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), harmonicsGravityProvider));
// Sun and Moon attraction
np.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
np.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
// atmospheric drag
MarshallSolarActivityFutureEstimation msafe = new MarshallSolarActivityFutureEstimation("Jan2000F10-edited-data\\.txt", MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
DataProvidersManager.getInstance().feed(msafe.getSupportedNames(), msafe);
DTM2000 atmosphere = new DTM2000(msafe, CelestialBodyFactory.getSun(), earth);
np.addForceModel(new DragForce(atmosphere, new IsotropicDrag(spacecraftArea, spacecraftDragCoefficient)));
// solar radiation pressure
np.addForceModel(new SolarRadiationPressure(CelestialBodyFactory.getSun(), earth.getEquatorialRadius(), new IsotropicRadiationSingleCoefficient(spacecraftArea, spacecraftReflectionCoefficient)));
return np;
}
use of org.hipparchus.ode.ODEIntegrator in project Orekit by CS-SI.
the class AngularCoordinatesTest method testShiftWithAcceleration.
@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 OrdinaryDifferentialEquation ode = new OrdinaryDifferentialEquation() {
public int getDimension() {
return 4;
}
public double[] computeDerivatives(final double t, final double[] q) {
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();
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) };
}
};
ODEIntegrator integrator = new DormandPrince853Integrator(1.0e-6, 1.0, 1.0e-12, 1.0e-12);
integrator.addStepHandler(new StepNormalizer(dt / n, new ODEFixedStepHandler() {
public void handleStep(ODEStateAndDerivative s, boolean isLast) {
final double t = s.getTime();
final double[] y = s.getPrimaryState();
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, new ODEState(0, y), dt);
}
Aggregations