use of org.hipparchus.ode.nonstiff.ClassicalRungeKuttaIntegrator in project Orekit by CS-SI.
the class DragForceTest method RealFieldTest.
/**
*Testing if the propagation between the FieldPropagation and the propagation
* is equivalent.
* Also testing if propagating X+dX with the propagation is equivalent to
* propagation X with the FieldPropagation and then applying the taylor
* expansion of dX to the result.
*/
@Test
public void RealFieldTest() throws OrekitException {
DSFactory factory = new DSFactory(6, 4);
DerivativeStructure a_0 = factory.variable(0, 7e6);
DerivativeStructure e_0 = factory.variable(1, 0.01);
DerivativeStructure i_0 = factory.variable(2, 1.2);
DerivativeStructure R_0 = factory.variable(3, 0.7);
DerivativeStructure O_0 = factory.variable(4, 0.5);
DerivativeStructure n_0 = factory.variable(5, 0.1);
Field<DerivativeStructure> field = a_0.getField();
DerivativeStructure zero = field.getZero();
FieldAbsoluteDate<DerivativeStructure> J2000 = new FieldAbsoluteDate<>(field);
Frame EME = FramesFactory.getEME2000();
FieldKeplerianOrbit<DerivativeStructure> FKO = new FieldKeplerianOrbit<>(a_0, e_0, i_0, R_0, O_0, n_0, PositionAngle.MEAN, EME, J2000, Constants.EIGEN5C_EARTH_MU);
FieldSpacecraftState<DerivativeStructure> initialState = new FieldSpacecraftState<>(FKO);
SpacecraftState iSR = initialState.toSpacecraftState();
ClassicalRungeKuttaFieldIntegrator<DerivativeStructure> integrator = new ClassicalRungeKuttaFieldIntegrator<>(field, zero.add(6));
ClassicalRungeKuttaIntegrator RIntegrator = new ClassicalRungeKuttaIntegrator(6);
OrbitType type = OrbitType.EQUINOCTIAL;
FieldNumericalPropagator<DerivativeStructure> FNP = new FieldNumericalPropagator<>(field, integrator);
FNP.setOrbitType(type);
FNP.setInitialState(initialState);
NumericalPropagator NP = new NumericalPropagator(RIntegrator);
NP.setOrbitType(type);
NP.setInitialState(iSR);
final DragForce forceModel = new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(), new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true))), new BoxAndSolarArraySpacecraft(1.5, 2.0, 1.8, CelestialBodyFactory.getSun(), 20.0, Vector3D.PLUS_J, 1.2, 0.7, 0.2));
FNP.addForceModel(forceModel);
NP.addForceModel(forceModel);
FieldAbsoluteDate<DerivativeStructure> target = J2000.shiftedBy(1000.);
FieldSpacecraftState<DerivativeStructure> finalState_DS = FNP.propagate(target);
SpacecraftState finalState_R = NP.propagate(target.toAbsoluteDate());
FieldPVCoordinates<DerivativeStructure> finPVC_DS = finalState_DS.getPVCoordinates();
PVCoordinates finPVC_R = finalState_R.getPVCoordinates();
Assert.assertEquals(finPVC_DS.toPVCoordinates().getPosition().getX(), finPVC_R.getPosition().getX(), FastMath.abs(finPVC_R.getPosition().getX()) * 1e-11);
Assert.assertEquals(finPVC_DS.toPVCoordinates().getPosition().getY(), finPVC_R.getPosition().getY(), FastMath.abs(finPVC_R.getPosition().getY()) * 1e-11);
Assert.assertEquals(finPVC_DS.toPVCoordinates().getPosition().getZ(), finPVC_R.getPosition().getZ(), FastMath.abs(finPVC_R.getPosition().getZ()) * 1e-11);
long number = 23091991;
RandomGenerator RG = new Well19937a(number);
GaussianRandomGenerator NGG = new GaussianRandomGenerator(RG);
UncorrelatedRandomVectorGenerator URVG = new UncorrelatedRandomVectorGenerator(new double[] { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }, new double[] { 1e3, 0.005, 0.005, 0.01, 0.01, 0.01 }, NGG);
double a_R = a_0.getReal();
double e_R = e_0.getReal();
double i_R = i_0.getReal();
double R_R = R_0.getReal();
double O_R = O_0.getReal();
double n_R = n_0.getReal();
for (int ii = 0; ii < 1; ii++) {
double[] rand_next = URVG.nextVector();
double a_shift = a_R + rand_next[0];
double e_shift = e_R + rand_next[1];
double i_shift = i_R + rand_next[2];
double R_shift = R_R + rand_next[3];
double O_shift = O_R + rand_next[4];
double n_shift = n_R + rand_next[5];
KeplerianOrbit shiftedOrb = new KeplerianOrbit(a_shift, e_shift, i_shift, R_shift, O_shift, n_shift, PositionAngle.MEAN, EME, J2000.toAbsoluteDate(), Constants.EIGEN5C_EARTH_MU);
SpacecraftState shift_iSR = new SpacecraftState(shiftedOrb);
NumericalPropagator shift_NP = new NumericalPropagator(RIntegrator);
shift_NP.setInitialState(shift_iSR);
shift_NP.addForceModel(forceModel);
SpacecraftState finalState_shift = shift_NP.propagate(target.toAbsoluteDate());
PVCoordinates finPVC_shift = finalState_shift.getPVCoordinates();
// position check
FieldVector3D<DerivativeStructure> pos_DS = finPVC_DS.getPosition();
double x_DS = pos_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
double y_DS = pos_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
double z_DS = pos_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
// System.out.println(pos_DS.getX().getPartialDerivative(1));
double x = finPVC_shift.getPosition().getX();
double y = finPVC_shift.getPosition().getY();
double z = finPVC_shift.getPosition().getZ();
Assert.assertEquals(x_DS, x, FastMath.abs(x - pos_DS.getX().getReal()) * 1e-5);
Assert.assertEquals(y_DS, y, FastMath.abs(y - pos_DS.getY().getReal()) * 1e-5);
Assert.assertEquals(z_DS, z, FastMath.abs(z - pos_DS.getZ().getReal()) * 1e-5);
// velocity check
FieldVector3D<DerivativeStructure> vel_DS = finPVC_DS.getVelocity();
double vx_DS = vel_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
double vy_DS = vel_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
double vz_DS = vel_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
double vx = finPVC_shift.getVelocity().getX();
double vy = finPVC_shift.getVelocity().getY();
double vz = finPVC_shift.getVelocity().getZ();
Assert.assertEquals(vx_DS, vx, FastMath.abs(vx) * 1e-7);
Assert.assertEquals(vy_DS, vy, FastMath.abs(vy) * 1e-7);
Assert.assertEquals(vz_DS, vz, FastMath.abs(vz) * 1e-7);
// acceleration check
FieldVector3D<DerivativeStructure> acc_DS = finPVC_DS.getAcceleration();
double ax_DS = acc_DS.getX().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
double ay_DS = acc_DS.getY().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
double az_DS = acc_DS.getZ().taylor(rand_next[0], rand_next[1], rand_next[2], rand_next[3], rand_next[4], rand_next[5]);
double ax = finPVC_shift.getAcceleration().getX();
double ay = finPVC_shift.getAcceleration().getY();
double az = finPVC_shift.getAcceleration().getZ();
Assert.assertEquals(ax_DS, ax, FastMath.abs(ax) * 1e-5);
Assert.assertEquals(ay_DS, ay, FastMath.abs(ay) * 1e-5);
Assert.assertEquals(az_DS, az, FastMath.abs(az) * 1e-5);
}
}
use of org.hipparchus.ode.nonstiff.ClassicalRungeKuttaIntegrator in project Orekit by CS-SI.
the class NumericalPropagatorTest method testNotInitialised2.
@Test(expected = OrekitException.class)
public void testNotInitialised2() throws OrekitException {
final AbstractIntegratedPropagator notInitialised = new NumericalPropagator(new ClassicalRungeKuttaIntegrator(10.0));
notInitialised.propagate(AbsoluteDate.J2000_EPOCH, AbsoluteDate.J2000_EPOCH.shiftedBy(3600));
}
use of org.hipparchus.ode.nonstiff.ClassicalRungeKuttaIntegrator in project Orekit by CS-SI.
the class NumericalPropagatorTest method testNotInitialised1.
@Test(expected = OrekitException.class)
public void testNotInitialised1() throws OrekitException {
final AbstractIntegratedPropagator notInitialised = new NumericalPropagator(new ClassicalRungeKuttaIntegrator(10.0));
notInitialised.propagate(AbsoluteDate.J2000_EPOCH);
}
use of org.hipparchus.ode.nonstiff.ClassicalRungeKuttaIntegrator in project Orekit by CS-SI.
the class DSSTPropagatorTest method testPropagationWithSolarRadiationPressure.
@Test
public void testPropagationWithSolarRadiationPressure() throws OrekitException {
// Central Body geopotential 2x0
final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(2, 0);
DSSTForceModel zonal = new DSSTZonal(provider, 2, 1, 5);
DSSTForceModel tesseral = new DSSTTesseral(CelestialBodyFactory.getEarth().getBodyOrientedFrame(), Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider, 2, 0, 0, 2, 2, 0, 0);
// SRP Force Model
DSSTForceModel srp = new DSSTSolarRadiationPressure(1.2, 100., CelestialBodyFactory.getSun(), Constants.WGS84_EARTH_EQUATORIAL_RADIUS);
// GEO Orbit
final AbsoluteDate initDate = new AbsoluteDate(2003, 9, 16, 0, 0, 00.000, TimeScalesFactory.getUTC());
final Orbit orbit = new KeplerianOrbit(42166258., 0.0001, FastMath.toRadians(0.001), FastMath.toRadians(315.4985), FastMath.toRadians(130.7562), FastMath.toRadians(44.2377), PositionAngle.MEAN, FramesFactory.getGCRF(), initDate, provider.getMu());
// Set propagator with state and force model
dsstProp = new DSSTPropagator(new ClassicalRungeKuttaIntegrator(86400.));
dsstProp.setInitialState(new SpacecraftState(orbit), false);
dsstProp.addForceModel(zonal);
dsstProp.addForceModel(tesseral);
dsstProp.addForceModel(srp);
// 10 days propagation
final SpacecraftState state = dsstProp.propagate(initDate.shiftedBy(10. * 86400.));
// Ref Standalone_DSST:
// a = 42166257.99807995 m
// h/ey = -0.1191876027555493D-03
// k/ex = -0.1781865038201885D-05
// p/hy = 0.6618387121369373D-05
// q/hx = -0.5624363171289686D-05
// lM = 140°3496229467104
Assert.assertEquals(42166257.99807995, state.getA(), 0.8);
Assert.assertEquals(-0.1781865038201885e-05, state.getEquinoctialEx(), 3.e-7);
Assert.assertEquals(-0.1191876027555493e-03, state.getEquinoctialEy(), 4.e-6);
Assert.assertEquals(-0.5624363171289686e-05, state.getHx(), 4.e-9);
Assert.assertEquals(0.6618387121369373e-05, state.getHy(), 3.e-10);
Assert.assertEquals(140.3496229467104, FastMath.toDegrees(MathUtils.normalizeAngle(state.getLM(), FastMath.PI)), 2.e-4);
}
use of org.hipparchus.ode.nonstiff.ClassicalRungeKuttaIntegrator in project Orekit by CS-SI.
the class OrekitStepHandlerTest method testIsInterpolated.
/**
* Check {@link OrekitStepInterpolator#isPreviousStateInterpolated()} and {@link
* OrekitStepInterpolator#isCurrentStateInterpolated()}.
*
* @throws OrekitException on error.
*/
@Test
public void testIsInterpolated() throws OrekitException {
// setup
NumericalPropagator propagator = new NumericalPropagator(new ClassicalRungeKuttaIntegrator(60));
AbsoluteDate date = AbsoluteDate.J2000_EPOCH;
Frame eci = FramesFactory.getGCRF();
SpacecraftState ic = new SpacecraftState(new KeplerianOrbit(6378137 + 500e3, 1e-3, 0, 0, 0, 0, PositionAngle.TRUE, eci, date, Constants.EIGEN5C_EARTH_MU));
propagator.setInitialState(ic);
propagator.setOrbitType(OrbitType.CARTESIAN);
// detector triggers half way through second step
DateDetector detector = new DateDetector(date.shiftedBy(90)).withHandler(new ContinueOnEvent<>());
propagator.addEventDetector(detector);
// action and verify
Queue<Boolean> expected = new ArrayDeque<>(Arrays.asList(false, false, false, true, true, false));
propagator.setMasterMode(new OrekitStepHandler() {
@Override
public void handleStep(OrekitStepInterpolator interpolator, boolean isLast) {
assertEquals(expected.poll(), interpolator.isPreviousStateInterpolated());
assertEquals(expected.poll(), interpolator.isCurrentStateInterpolated());
}
});
final AbsoluteDate end = date.shiftedBy(120);
assertEquals(end, propagator.propagate(end).getDate());
}
Aggregations