use of org.orekit.propagation.FieldSpacecraftState in project Orekit by CS-SI.
the class ThirdBodyAttractionTest 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();
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 ThirdBodyAttraction forceModel = new ThirdBodyAttraction(CelestialBodyFactory.getSun());
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();
double maxP = 0;
double maxV = 0;
double maxA = 0;
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]);
double x = finPVC_shift.getPosition().getX();
double y = finPVC_shift.getPosition().getY();
double z = finPVC_shift.getPosition().getZ();
maxP = FastMath.max(maxP, FastMath.abs((x_DS - x) / (x - pos_DS.getX().getReal())));
maxP = FastMath.max(maxP, FastMath.abs((y_DS - y) / (y - pos_DS.getY().getReal())));
maxP = FastMath.max(maxP, FastMath.abs((z_DS - z) / (z - pos_DS.getZ().getReal())));
// 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();
maxV = FastMath.max(maxV, FastMath.abs((vx_DS - vx) / vx));
maxV = FastMath.max(maxV, FastMath.abs((vy_DS - vy) / vy));
maxV = FastMath.max(maxV, FastMath.abs((vz_DS - vz) / vz));
// 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();
maxA = FastMath.max(maxA, FastMath.abs((ax_DS - ax) / ax));
maxA = FastMath.max(maxA, FastMath.abs((ay_DS - ay) / ay));
maxA = FastMath.max(maxA, FastMath.abs((az_DS - az) / az));
}
Assert.assertEquals(0, maxP, 5.0e-9);
Assert.assertEquals(0, maxV, 3.0e-10);
Assert.assertEquals(0, maxA, 8.0e-8);
}
use of org.orekit.propagation.FieldSpacecraftState in project Orekit by CS-SI.
the class ThirdBodyAttractionTest method RealFieldExpectErrorTest.
/**
*Same test as the previous one but not adding the ForceModel to the NumericalPropagator
* it is a test to validate the previous test.
* (to test if the ForceModel it's actually
* doing something in the Propagator and the FieldPropagator)
*/
@Test
public void RealFieldExpectErrorTest() 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();
OrbitType type = OrbitType.KEPLERIAN;
double[][] tolerance = NumericalPropagator.tolerances(0.001, 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 ThirdBodyAttraction forceModel = new ThirdBodyAttraction(CelestialBodyFactory.getSun());
FNP.addForceModel(forceModel);
// NOT ADDING THE FORCE MODEL TO THE NUMERICAL PROPAGATOR 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.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getX() - finPVC_R.getPosition().getX()) < FastMath.abs(finPVC_R.getPosition().getX()) * 1e-11);
Assert.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getY() - finPVC_R.getPosition().getY()) < FastMath.abs(finPVC_R.getPosition().getY()) * 1e-11);
Assert.assertFalse(FastMath.abs(finPVC_DS.toPVCoordinates().getPosition().getZ() - finPVC_R.getPosition().getZ()) < FastMath.abs(finPVC_R.getPosition().getZ()) * 1e-11);
}
use of org.orekit.propagation.FieldSpacecraftState in project Orekit by CS-SI.
the class FieldNumericalPropagatorTest method doTestResetStateEvent.
private <T extends RealFieldElement<T>> void doTestResetStateEvent(Field<T> field) throws OrekitException {
T zero = field.getZero();
// setup
final FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field);
FieldSpacecraftState<T> initialState;
FieldNumericalPropagator<T> propagator;
final FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity), FramesFactory.getEME2000(), initDate, mu);
initialState = new FieldSpacecraftState<>(orbit);
OrbitType type = OrbitType.EQUINOCTIAL;
double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit.toOrbit(), type);
AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
integrator.setInitialStepSize(zero.add(60));
propagator = new FieldNumericalPropagator<>(field, integrator);
propagator.setOrbitType(type);
propagator.setInitialState(initialState);
final FieldAbsoluteDate<T> resetDate = initDate.shiftedBy(1000);
CheckingHandler<FieldDateDetector<T>, T> checking = new CheckingHandler<FieldDateDetector<T>, T>(Action.RESET_STATE) {
public FieldSpacecraftState<T> resetState(FieldDateDetector<T> detector, FieldSpacecraftState<T> oldState) {
return new FieldSpacecraftState<>(oldState.getOrbit(), oldState.getAttitude(), oldState.getMass().subtract(200.0));
}
};
propagator.addEventDetector(new FieldDateDetector<>(resetDate).withHandler(checking));
checking.assertEvent(false);
final FieldSpacecraftState<T> finalState = propagator.propagate(initDate.shiftedBy(3200));
checking.assertEvent(true);
Assert.assertEquals(initialState.getMass().getReal() - 200, finalState.getMass().getReal(), 1.0e-10);
}
use of org.orekit.propagation.FieldSpacecraftState in project Orekit by CS-SI.
the class FieldNumericalPropagatorTest method doTestEventDetectionBug.
private <T extends RealFieldElement<T>> void doTestEventDetectionBug(final Field<T> field) throws OrekitException {
T zero = field.getZero();
TimeScale utc = TimeScalesFactory.getUTC();
FieldAbsoluteDate<T> initialDate = new FieldAbsoluteDate<>(field, 2005, 1, 1, 0, 0, 0.0, utc);
T duration = zero.add(100000.0);
FieldAbsoluteDate<T> endDate = new FieldAbsoluteDate<>(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)));
FieldOrbit<T> geo = new FieldCartesianOrbit<>(new FieldPVCoordinates<>(new FieldVector3D<>(zero.add(rApogee), zero, zero), new FieldVector3D<>(zero, zero.add(vApogee), zero)), EME2000, initialDate, mu);
duration = geo.getKeplerianPeriod();
endDate = new FieldAbsoluteDate<>(initialDate, duration);
// Numerical Integration
final double minStep = 0.001;
final double maxStep = 1000;
final double initStep = 60;
final OrbitType type = OrbitType.EQUINOCTIAL;
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 };
AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, minStep, maxStep, absTolerance, relTolerance);
integrator.setInitialStepSize(zero.add(initStep));
// Numerical propagator based on the integrator
FieldNumericalPropagator<T> propagator = new FieldNumericalPropagator<>(field, integrator);
propagator.setOrbitType(type);
T mass = field.getZero().add(1000.0);
FieldSpacecraftState<T> initialState = new FieldSpacecraftState<>(geo, mass);
propagator.setInitialState(initialState);
propagator.setOrbitType(OrbitType.CARTESIAN);
// Set the events Detectors
FieldApsideDetector<T> event1 = new FieldApsideDetector<>(geo);
propagator.addEventDetector(event1);
// Set the propagation mode
propagator.setSlaveMode();
// Propagate
FieldSpacecraftState<T> finalState = propagator.propagate(endDate);
// we should stop long before endDate
Assert.assertTrue(endDate.durationFrom(finalState.getDate()).getReal() > 40000.0);
}
use of org.orekit.propagation.FieldSpacecraftState in project Orekit by CS-SI.
the class FieldNumericalPropagatorTest method doTestPropagationTypesHyperbolic.
private <T extends RealFieldElement<T>> void doTestPropagationTypesHyperbolic(Field<T> field) throws OrekitException {
T zero = field.getZero();
// setup
final FieldAbsoluteDate<T> initDate = FieldAbsoluteDate.getJ2000Epoch(field);
FieldSpacecraftState<T> initialState;
FieldNumericalPropagator<T> propagator;
final FieldVector3D<T> position = new FieldVector3D<>(zero.add(7.0e6), zero.add(1.0e6), zero.add(4.0e6));
final FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(-500.0), zero.add(8000.0), zero.add(1000.0));
final FieldOrbit<T> orbit = new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity), FramesFactory.getEME2000(), initDate, mu);
initialState = new FieldSpacecraftState<>(orbit);
OrbitType type = OrbitType.EQUINOCTIAL;
double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit.toOrbit(), type);
AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, 0.001, 200, tolerance[0], tolerance[1]);
integrator.setInitialStepSize(zero.add(60));
propagator = new FieldNumericalPropagator<>(field, integrator);
propagator.setOrbitType(type);
propagator.setInitialState(initialState);
FieldSpacecraftState<T> state = new FieldSpacecraftState<>(new FieldKeplerianOrbit<>(zero.add(-10000000.0), zero.add(2.5), zero.add(0.3), zero, zero, zero, PositionAngle.TRUE, FramesFactory.getEME2000(), initDate, mu));
ForceModel gravityField = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), GravityFieldFactory.getNormalizedProvider(5, 5));
propagator.addForceModel(gravityField);
// Propagation of the initial at t + dt
final FieldPVCoordinates<T> pv = state.getPVCoordinates();
final T dP = zero.add(0.001);
final T dV = dP.multiply(state.getMu()).divide(pv.getPosition().getNormSq().multiply(pv.getVelocity().getNorm()));
final FieldPVCoordinates<T> pvcM = propagateInType(state, dP, OrbitType.CARTESIAN, PositionAngle.MEAN, propagator);
final FieldPVCoordinates<T> pvkM = propagateInType(state, dP, OrbitType.KEPLERIAN, PositionAngle.MEAN, propagator);
final FieldPVCoordinates<T> pvcE = propagateInType(state, dP, OrbitType.CARTESIAN, PositionAngle.ECCENTRIC, propagator);
final FieldPVCoordinates<T> pvkE = propagateInType(state, dP, OrbitType.KEPLERIAN, PositionAngle.ECCENTRIC, propagator);
final FieldPVCoordinates<T> pvcT = propagateInType(state, dP, OrbitType.CARTESIAN, PositionAngle.TRUE, propagator);
final FieldPVCoordinates<T> pvkT = propagateInType(state, dP, OrbitType.KEPLERIAN, PositionAngle.TRUE, propagator);
Assert.assertEquals(0, pvcM.getPosition().subtract(pvkT.getPosition()).getNorm().getReal() / dP.getReal(), 0.3);
Assert.assertEquals(0, pvcM.getVelocity().subtract(pvkT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.4);
Assert.assertEquals(0, pvkM.getPosition().subtract(pvkT.getPosition()).getNorm().getReal() / dP.getReal(), 0.2);
Assert.assertEquals(0, pvkM.getVelocity().subtract(pvkT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.3);
Assert.assertEquals(0, pvcE.getPosition().subtract(pvkT.getPosition()).getNorm().getReal() / dP.getReal(), 0.3);
Assert.assertEquals(0, pvcE.getVelocity().subtract(pvkT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.4);
Assert.assertEquals(0, pvkE.getPosition().subtract(pvkT.getPosition()).getNorm().getReal() / dP.getReal(), 0.009);
Assert.assertEquals(0, pvkE.getVelocity().subtract(pvkT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.006);
Assert.assertEquals(0, pvcT.getPosition().subtract(pvkT.getPosition()).getNorm().getReal() / dP.getReal(), 0.3);
Assert.assertEquals(0, pvcT.getVelocity().subtract(pvkT.getVelocity()).getNorm().getReal() / dV.getReal(), 0.4);
}
Aggregations