use of org.orekit.orbits.KeplerianOrbit in project Orekit by CS-SI.
the class ThirdBodyAttractionTest method testGlobalStateJacobian.
@Test
public void testGlobalStateJacobian() 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);
}
use of org.orekit.orbits.KeplerianOrbit in project Orekit by CS-SI.
the class ConstantThrustManeuverInitializationTest method setUp.
@Before
public void setUp() throws OrekitException {
startDate = new AbsoluteDate();
double a = Constants.EGM96_EARTH_EQUATORIAL_RADIUS + 400e3;
double e = 0.001;
double i = (Math.PI / 4);
double pa = 0.0;
double raan = 0.0;
double anomaly = 0.0;
PositionAngle type = PositionAngle.MEAN;
Frame frame = FramesFactory.getEME2000();
double mu = Constants.EGM96_EARTH_MU;
Orbit orbit = new KeplerianOrbit(a, e, i, pa, raan, anomaly, type, frame, startDate, mu);
initialState = new SpacecraftState(orbit, mass);
// Numerical Propagator
double minStep = 0.001;
double maxStep = 1000.0;
double positionTolerance = 10.;
OrbitType propagationType = OrbitType.KEPLERIAN;
double[][] tolerances = NumericalPropagator.tolerances(positionTolerance, orbit, propagationType);
AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tolerances[0], tolerances[1]);
// Set up propagator
propagator = new NumericalPropagator(integrator);
propagator.setOrbitType(propagationType);
// Control deltaVs and mass changes
double flowRate = -thrust / (Constants.G0_STANDARD_GRAVITY * isp);
massControlFullForward = mass + (flowRate * duration);
deltaVControlFullForward = isp * Constants.G0_STANDARD_GRAVITY * FastMath.log(mass / massControlFullForward);
massControlHalfForward = mass + (flowRate * duration / 2);
massControlFullReverse = mass - (flowRate * duration);
deltaVControlFullReverse = isp * Constants.G0_STANDARD_GRAVITY * FastMath.log(massControlFullReverse / mass);
massControlHalfReverse = mass - (flowRate * duration / 2);
deltaVControlHalfReverse = isp * Constants.G0_STANDARD_GRAVITY * FastMath.log(massControlHalfReverse / mass);
}
use of org.orekit.orbits.KeplerianOrbit in project Orekit by CS-SI.
the class ConstantThrustManeuverTest 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, 5);
DerivativeStructure a_0 = factory.variable(0, 7e7);
DerivativeStructure e_0 = factory.variable(1, 0.4);
DerivativeStructure i_0 = factory.variable(2, 85 * FastMath.PI / 180);
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();
final OrbitType type = OrbitType.KEPLERIAN;
double[][] tolerance = NumericalPropagator.tolerances(10.0, FKO.toOrbit(), type);
AdaptiveStepsizeFieldIntegrator<DerivativeStructure> integrator = new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
integrator.setInitialStepSize(zero.add(60));
AdaptiveStepsizeIntegrator RIntegrator = new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
RIntegrator.setInitialStepSize(60);
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 ConstantThrustManeuver forceModel = new ConstantThrustManeuver(J2000.toAbsoluteDate().shiftedBy(100), 100.0, 400.0, 300.0, Vector3D.PLUS_K);
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.01, 0.01, 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-8);
Assert.assertEquals(y_DS, y, FastMath.abs(y - pos_DS.getY().getReal()) * 1e-8);
Assert.assertEquals(z_DS, z, FastMath.abs(z - pos_DS.getZ().getReal()) * 1e-8);
// 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-9);
Assert.assertEquals(vy_DS, vy, FastMath.abs(vy) * 1e-9);
Assert.assertEquals(vz_DS, vz, FastMath.abs(vz) * 1e-9);
// 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-8);
Assert.assertEquals(ay_DS, ay, FastMath.abs(ay) * 1e-8);
Assert.assertEquals(az_DS, az, FastMath.abs(az) * 1e-8);
}
}
use of org.orekit.orbits.KeplerianOrbit in project Orekit by CS-SI.
the class ConstantThrustManeuverTest method testForwardAndBackward.
@Test
public void testForwardAndBackward() 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);
Assert.assertEquals(f, maneuver.getThrust(), 1.0e-10);
Assert.assertEquals(isp, maneuver.getISP(), 1.0e-10);
double[][] tol = NumericalPropagator.tolerances(1.0, orbit, OrbitType.KEPLERIAN);
AdaptiveStepsizeIntegrator integrator1 = new DormandPrince853Integrator(0.001, 1000, tol[0], tol[1]);
integrator1.setInitialStepSize(60);
final NumericalPropagator propagator1 = new NumericalPropagator(integrator1);
propagator1.setInitialState(initialState);
propagator1.setAttitudeProvider(law);
propagator1.addForceModel(maneuver);
final SpacecraftState finalState = propagator1.propagate(fireDate.shiftedBy(3800));
AdaptiveStepsizeIntegrator integrator2 = new DormandPrince853Integrator(0.001, 1000, tol[0], tol[1]);
integrator2.setInitialStepSize(60);
final NumericalPropagator propagator2 = new NumericalPropagator(integrator2);
propagator2.setInitialState(finalState);
propagator2.setAttitudeProvider(law);
propagator2.addForceModel(maneuver);
final SpacecraftState recoveredState = propagator2.propagate(orbit.getDate());
final Vector3D refPosition = initialState.getPVCoordinates().getPosition();
final Vector3D recoveredPosition = recoveredState.getPVCoordinates().getPosition();
Assert.assertEquals(0.0, Vector3D.distance(refPosition, recoveredPosition), 30.0);
Assert.assertEquals(initialState.getMass(), recoveredState.getMass(), 1.5e-10);
}
use of org.orekit.orbits.KeplerianOrbit in project Orekit by CS-SI.
the class ConstantThrustManeuverTest method testJacobianVs80Implementation.
@Test
public void testJacobianVs80Implementation() 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 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 AttitudeProvider law = new LofOffset(orbit.getFrame(), LOFType.VVLH);
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);
// before maneuver (Jacobian wrt. state is zero)
checkStateJacobianVs80Implementation(initialState, maneuver, law, 1.0e-50, false);
// in maneuver
SpacecraftState startState = initialState.shiftedBy(fireDate.durationFrom(initDate));
maneuver.getEventsDetectors().findFirst().get().eventOccurred(startState, true);
SpacecraftState midState = startState.shiftedBy(duration / 2.0);
checkStateJacobianVs80Implementation(midState, maneuver, law, 1.0e-20, false);
}
Aggregations