use of org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator in project Orekit by CS-SI.
the class NumericalPropagatorTest method propagateInType.
private PVCoordinates propagateInType(SpacecraftState state, double dP, OrbitType type, PositionAngle angle) throws OrekitException {
final double dt = 3200;
final double minStep = 0.001;
final double maxStep = 1000;
double[][] tol = NumericalPropagator.tolerances(dP, state.getOrbit(), type);
AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
NumericalPropagator newPropagator = new NumericalPropagator(integrator);
newPropagator.setOrbitType(type);
newPropagator.setPositionAngleType(angle);
newPropagator.setInitialState(state);
for (ForceModel force : propagator.getAllForceModels()) {
newPropagator.addForceModel(force);
}
return newPropagator.propagate(state.getDate().shiftedBy(dt)).getPVCoordinates();
}
use of org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator in project Orekit by CS-SI.
the class NumericalPropagatorTest method testPropagationTypesElliptical.
@Test
public void testPropagationTypesElliptical() throws OrekitException, ParseException, IOException {
// setup
AbsoluteDate initDate = new AbsoluteDate();
SpacecraftState initialState;
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);
initialState = new SpacecraftState(orbit);
OrbitType type = OrbitType.EQUINOCTIAL;
double[][] tolerance = NumericalPropagator.tolerances(0.001, orbit, type);
AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 200, tolerance[0], tolerance[1]);
integrator.setInitialStepSize(60);
propagator = new NumericalPropagator(integrator);
propagator.setOrbitType(type);
propagator.setInitialState(initialState);
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 PVCoordinates pv = initialState.getPVCoordinates();
final double dP = 0.001;
final double dV = initialState.getMu() * dP / (pv.getPosition().getNormSq() * pv.getVelocity().getNorm());
final PVCoordinates pvcM = propagateInType(initialState, dP, OrbitType.CARTESIAN, PositionAngle.MEAN);
final PVCoordinates pviM = propagateInType(initialState, dP, OrbitType.CIRCULAR, PositionAngle.MEAN);
final PVCoordinates pveM = propagateInType(initialState, dP, OrbitType.EQUINOCTIAL, PositionAngle.MEAN);
final PVCoordinates pvkM = propagateInType(initialState, dP, OrbitType.KEPLERIAN, PositionAngle.MEAN);
final PVCoordinates pvcE = propagateInType(initialState, dP, OrbitType.CARTESIAN, PositionAngle.ECCENTRIC);
final PVCoordinates pviE = propagateInType(initialState, dP, OrbitType.CIRCULAR, PositionAngle.ECCENTRIC);
final PVCoordinates pveE = propagateInType(initialState, dP, OrbitType.EQUINOCTIAL, PositionAngle.ECCENTRIC);
final PVCoordinates pvkE = propagateInType(initialState, dP, OrbitType.KEPLERIAN, PositionAngle.ECCENTRIC);
final PVCoordinates pvcT = propagateInType(initialState, dP, OrbitType.CARTESIAN, PositionAngle.TRUE);
final PVCoordinates pviT = propagateInType(initialState, dP, OrbitType.CIRCULAR, PositionAngle.TRUE);
final PVCoordinates pveT = propagateInType(initialState, dP, OrbitType.EQUINOCTIAL, PositionAngle.TRUE);
final PVCoordinates pvkT = propagateInType(initialState, dP, OrbitType.KEPLERIAN, PositionAngle.TRUE);
Assert.assertEquals(0, pvcM.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 3.0);
Assert.assertEquals(0, pvcM.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 2.0);
Assert.assertEquals(0, pviM.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 0.6);
Assert.assertEquals(0, pviM.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 0.4);
Assert.assertEquals(0, pvkM.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 0.5);
Assert.assertEquals(0, pvkM.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 0.3);
Assert.assertEquals(0, pveM.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 0.2);
Assert.assertEquals(0, pveM.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 0.2);
Assert.assertEquals(0, pvcE.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 3.0);
Assert.assertEquals(0, pvcE.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 2.0);
Assert.assertEquals(0, pviE.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 0.03);
Assert.assertEquals(0, pviE.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 0.04);
Assert.assertEquals(0, pvkE.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 0.4);
Assert.assertEquals(0, pvkE.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 0.3);
Assert.assertEquals(0, pveE.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 0.2);
Assert.assertEquals(0, pveE.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 0.07);
Assert.assertEquals(0, pvcT.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 3.0);
Assert.assertEquals(0, pvcT.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 2.0);
Assert.assertEquals(0, pviT.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 0.3);
Assert.assertEquals(0, pviT.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 0.2);
Assert.assertEquals(0, pvkT.getPosition().subtract(pveT.getPosition()).getNorm() / dP, 0.4);
Assert.assertEquals(0, pvkT.getVelocity().subtract(pveT.getVelocity()).getNorm() / dV, 0.2);
}
use of org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator in project Orekit by CS-SI.
the class DSSTPropagatorTest method setDSSTProp.
private void setDSSTProp(SpacecraftState initialState) throws OrekitException {
initialState.getDate();
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);
}
use of org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator in project Orekit by CS-SI.
the class DSSTPropagatorTest method testEphemerisGeneration.
@Test
public void testEphemerisGeneration() 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());
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, false);
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 DSSTZonal(nshp, 8, 7, 17));
propagator.addForceModel(new DSSTTesseral(earth.getBodyFrame(), Constants.WGS84_EARTH_ANGULAR_VELOCITY, nshp, 8, 8, 4, 12, 8, 8, 4));
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.setInterpolationGridToMaxTimeGap(0.5 * Constants.JULIAN_DAY);
// direct generation of states
propagator.setInitialState(new SpacecraftState(orbit, 45.0), false);
final List<SpacecraftState> states = new ArrayList<SpacecraftState>();
propagator.setMasterMode(600, (currentState, isLast) -> states.add(currentState));
propagator.propagate(orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY));
// ephemeris generation
propagator.setInitialState(new SpacecraftState(orbit, 45.0), false);
propagator.setEphemerisMode();
propagator.propagate(orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY));
BoundedPropagator ephemeris = propagator.getGeneratedEphemeris();
double maxError = 0;
for (final SpacecraftState state : states) {
final SpacecraftState fromEphemeris = ephemeris.propagate(state.getDate());
final double error = Vector3D.distance(state.getPVCoordinates().getPosition(), fromEphemeris.getPVCoordinates().getPosition());
maxError = FastMath.max(maxError, error);
}
Assert.assertEquals(0.0, maxError, 1.0e-10);
}
use of org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator in project Orekit by CS-SI.
the class RelativityTest 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, 0);
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 Relativity forceModel = new Relativity(Constants.EIGEN5C_EARTH_MU);
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.assertEquals(0, Vector3D.distance(finPVC_DS.toPVCoordinates().getPosition(), finPVC_R.getPosition()), 8.0e-13 * finPVC_R.getPosition().getNorm());
}
Aggregations