Search in sources :

Example 31 with AdaptiveStepsizeIntegrator

use of org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator in project Orekit by CS-SI.

the class HarmonicParametricAccelerationTest method doTestEquivalentManeuver.

private void doTestEquivalentManeuver(final double mass, final AttitudeProvider maneuverLaw, final ConstantThrustManeuver maneuver, final AttitudeProvider accelerationLaw, final HarmonicParametricAcceleration parametricAcceleration, final double positionTolerance) throws OrekitException {
    SpacecraftState initialState = new SpacecraftState(initialOrbit, maneuverLaw.getAttitude(initialOrbit, initialOrbit.getDate(), initialOrbit.getFrame()), mass);
    double[][] tolerance = NumericalPropagator.tolerances(10, initialOrbit, initialOrbit.getType());
    // propagator 0 uses a maneuver that is so efficient it does not consume any fuel
    // (hence mass remains constant)
    AdaptiveStepsizeIntegrator integrator0 = new DormandPrince853Integrator(0.001, 100, tolerance[0], tolerance[1]);
    integrator0.setInitialStepSize(60);
    final NumericalPropagator propagator0 = new NumericalPropagator(integrator0);
    propagator0.setInitialState(initialState);
    propagator0.setAttitudeProvider(maneuverLaw);
    propagator0.addForceModel(maneuver);
    // propagator 1 uses a constant acceleration
    AdaptiveStepsizeIntegrator integrator1 = new DormandPrince853Integrator(0.001, 100, tolerance[0], tolerance[1]);
    integrator1.setInitialStepSize(60);
    final NumericalPropagator propagator1 = new NumericalPropagator(integrator1);
    propagator1.setInitialState(initialState);
    propagator1.setAttitudeProvider(accelerationLaw);
    propagator1.addForceModel(parametricAcceleration);
    MultiSatStepHandler handler = (interpolators, isLast) -> {
        Vector3D p0 = interpolators.get(0).getCurrentState().getPVCoordinates().getPosition();
        Vector3D p1 = interpolators.get(1).getCurrentState().getPVCoordinates().getPosition();
        Assert.assertEquals(0.0, Vector3D.distance(p0, p1), positionTolerance);
    };
    PropagatorsParallelizer parallelizer = new PropagatorsParallelizer(Arrays.asList(propagator0, propagator1), handler);
    parallelizer.propagate(initialOrbit.getDate(), initialOrbit.getDate().shiftedBy(1000.0));
}
Also used : ParameterDriver(org.orekit.utils.ParameterDriver) Arrays(java.util.Arrays) AdaptiveStepsizeFieldIntegrator(org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator) LOFType(org.orekit.frames.LOFType) Rotation(org.hipparchus.geometry.euclidean.threed.Rotation) AdaptiveStepsizeIntegrator(org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator) AttitudeProvider(org.orekit.attitudes.AttitudeProvider) PVCoordinates(org.orekit.utils.PVCoordinates) SpacecraftState(org.orekit.propagation.SpacecraftState) LevenbergMarquardtOptimizer(org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer) PositionAngle(org.orekit.orbits.PositionAngle) Utils(org.orekit.Utils) CartesianOrbit(org.orekit.orbits.CartesianOrbit) PropagatorsParallelizer(org.orekit.propagation.PropagatorsParallelizer) FramesFactory(org.orekit.frames.FramesFactory) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) ObservedMeasurement(org.orekit.estimation.measurements.ObservedMeasurement) LofOffset(org.orekit.attitudes.LofOffset) List(java.util.List) InertialProvider(org.orekit.attitudes.InertialProvider) RealFieldElement(org.hipparchus.RealFieldElement) CelestialBodyPointed(org.orekit.attitudes.CelestialBodyPointed) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) ConstantThrustManeuver(org.orekit.forces.maneuvers.ConstantThrustManeuver) PV(org.orekit.estimation.measurements.PV) DormandPrince853Integrator(org.hipparchus.ode.nonstiff.DormandPrince853Integrator) MultiSatStepHandler(org.orekit.propagation.sampling.MultiSatStepHandler) FieldNumericalPropagator(org.orekit.propagation.numerical.FieldNumericalPropagator) DormandPrince853FieldIntegrator(org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator) Orbit(org.orekit.orbits.Orbit) FieldSpacecraftState(org.orekit.propagation.FieldSpacecraftState) ArrayList(java.util.ArrayList) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) DormandPrince853IntegratorBuilder(org.orekit.propagation.conversion.DormandPrince853IntegratorBuilder) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) FastMath(org.hipparchus.util.FastMath) FieldBoundedPropagator(org.orekit.propagation.FieldBoundedPropagator) Before(org.junit.Before) CircularOrbit(org.orekit.orbits.CircularOrbit) Constants(org.orekit.utils.Constants) DateComponents(org.orekit.time.DateComponents) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) Test(org.junit.Test) Field(org.hipparchus.Field) OrekitException(org.orekit.errors.OrekitException) CelestialBodyFactory(org.orekit.bodies.CelestialBodyFactory) TimeScalesFactory(org.orekit.time.TimeScalesFactory) Decimal64Field(org.hipparchus.util.Decimal64Field) TimeComponents(org.orekit.time.TimeComponents) Assert(org.junit.Assert) BatchLSEstimator(org.orekit.estimation.leastsquares.BatchLSEstimator) AbsoluteDate(org.orekit.time.AbsoluteDate) SpacecraftState(org.orekit.propagation.SpacecraftState) FieldSpacecraftState(org.orekit.propagation.FieldSpacecraftState) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) FieldNumericalPropagator(org.orekit.propagation.numerical.FieldNumericalPropagator) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) AdaptiveStepsizeIntegrator(org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator) PropagatorsParallelizer(org.orekit.propagation.PropagatorsParallelizer) DormandPrince853Integrator(org.hipparchus.ode.nonstiff.DormandPrince853Integrator) MultiSatStepHandler(org.orekit.propagation.sampling.MultiSatStepHandler)

Example 32 with AdaptiveStepsizeIntegrator

use of org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator in project Orekit by CS-SI.

the class PolynomialParametricAccelerationTest method doTestEquivalentManeuver.

private void doTestEquivalentManeuver(final double mass, final AttitudeProvider maneuverLaw, final ConstantThrustManeuver maneuver, final AttitudeProvider accelerationLaw, final PolynomialParametricAcceleration parametricAcceleration, final double positionTolerance) throws OrekitException {
    SpacecraftState initialState = new SpacecraftState(initialOrbit, maneuverLaw.getAttitude(initialOrbit, initialOrbit.getDate(), initialOrbit.getFrame()), mass);
    double[][] tolerance = NumericalPropagator.tolerances(10, initialOrbit, initialOrbit.getType());
    // propagator 0 uses a maneuver that is so efficient it does not consume any fuel
    // (hence mass remains constant)
    AdaptiveStepsizeIntegrator integrator0 = new DormandPrince853Integrator(0.001, 100, tolerance[0], tolerance[1]);
    integrator0.setInitialStepSize(60);
    final NumericalPropagator propagator0 = new NumericalPropagator(integrator0);
    propagator0.setInitialState(initialState);
    propagator0.setAttitudeProvider(maneuverLaw);
    propagator0.addForceModel(maneuver);
    // propagator 1 uses a constant acceleration
    AdaptiveStepsizeIntegrator integrator1 = new DormandPrince853Integrator(0.001, 100, tolerance[0], tolerance[1]);
    integrator1.setInitialStepSize(60);
    final NumericalPropagator propagator1 = new NumericalPropagator(integrator1);
    propagator1.setInitialState(initialState);
    propagator1.setAttitudeProvider(accelerationLaw);
    propagator1.addForceModel(parametricAcceleration);
    MultiSatStepHandler handler = (interpolators, isLast) -> {
        Vector3D p0 = interpolators.get(0).getCurrentState().getPVCoordinates().getPosition();
        Vector3D p1 = interpolators.get(1).getCurrentState().getPVCoordinates().getPosition();
        Assert.assertEquals(0.0, Vector3D.distance(p0, p1), positionTolerance);
    };
    PropagatorsParallelizer parallelizer = new PropagatorsParallelizer(Arrays.asList(propagator0, propagator1), handler);
    parallelizer.propagate(initialOrbit.getDate(), initialOrbit.getDate().shiftedBy(1000.0));
}
Also used : DormandPrince853Integrator(org.hipparchus.ode.nonstiff.DormandPrince853Integrator) Arrays(java.util.Arrays) AdaptiveStepsizeFieldIntegrator(org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator) LOFType(org.orekit.frames.LOFType) MultiSatStepHandler(org.orekit.propagation.sampling.MultiSatStepHandler) Rotation(org.hipparchus.geometry.euclidean.threed.Rotation) AdaptiveStepsizeIntegrator(org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator) FieldNumericalPropagator(org.orekit.propagation.numerical.FieldNumericalPropagator) DormandPrince853FieldIntegrator(org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator) AttitudeProvider(org.orekit.attitudes.AttitudeProvider) Orbit(org.orekit.orbits.Orbit) FieldSpacecraftState(org.orekit.propagation.FieldSpacecraftState) PVCoordinates(org.orekit.utils.PVCoordinates) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) SpacecraftState(org.orekit.propagation.SpacecraftState) PositionAngle(org.orekit.orbits.PositionAngle) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) FastMath(org.hipparchus.util.FastMath) FieldBoundedPropagator(org.orekit.propagation.FieldBoundedPropagator) Utils(org.orekit.Utils) Before(org.junit.Before) CartesianOrbit(org.orekit.orbits.CartesianOrbit) Constants(org.orekit.utils.Constants) DateComponents(org.orekit.time.DateComponents) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) PropagatorsParallelizer(org.orekit.propagation.PropagatorsParallelizer) FramesFactory(org.orekit.frames.FramesFactory) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) Test(org.junit.Test) LofOffset(org.orekit.attitudes.LofOffset) Field(org.hipparchus.Field) InertialProvider(org.orekit.attitudes.InertialProvider) OrekitException(org.orekit.errors.OrekitException) RealFieldElement(org.hipparchus.RealFieldElement) CelestialBodyPointed(org.orekit.attitudes.CelestialBodyPointed) CelestialBodyFactory(org.orekit.bodies.CelestialBodyFactory) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) TimeScalesFactory(org.orekit.time.TimeScalesFactory) Decimal64Field(org.hipparchus.util.Decimal64Field) ConstantThrustManeuver(org.orekit.forces.maneuvers.ConstantThrustManeuver) TimeComponents(org.orekit.time.TimeComponents) Assert(org.junit.Assert) AbsoluteDate(org.orekit.time.AbsoluteDate) FieldSpacecraftState(org.orekit.propagation.FieldSpacecraftState) SpacecraftState(org.orekit.propagation.SpacecraftState) FieldNumericalPropagator(org.orekit.propagation.numerical.FieldNumericalPropagator) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) AdaptiveStepsizeIntegrator(org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator) PropagatorsParallelizer(org.orekit.propagation.PropagatorsParallelizer) DormandPrince853Integrator(org.hipparchus.ode.nonstiff.DormandPrince853Integrator) MultiSatStepHandler(org.orekit.propagation.sampling.MultiSatStepHandler)

Example 33 with AdaptiveStepsizeIntegrator

use of org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator in project Orekit by CS-SI.

the class HolmesFeatherstoneAttractionModelTest 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, 7201009.7124401);
    DerivativeStructure e_0 = factory.variable(1, 1e-3);
    DerivativeStructure i_0 = factory.variable(2, 98.7 * FastMath.PI / 180);
    DerivativeStructure R_0 = factory.variable(3, 15.0 * 22.5 * FastMath.PI / 180);
    DerivativeStructure O_0 = factory.variable(4, 93.0 * FastMath.PI / 180);
    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.EQUINOCTIAL;
    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);
    double[][] c = new double[3][1];
    c[0][0] = 0.0;
    c[2][0] = normalizedC20;
    double[][] s = new double[3][1];
    NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(6378136.460, mu, TideSystem.UNKNOWN, c, s);
    HolmesFeatherstoneAttractionModel forceModel = new HolmesFeatherstoneAttractionModel(itrf, provider);
    // FNP.addForceModel(forceModel);
    NP.addForceModel(forceModel);
    FieldAbsoluteDate<DerivativeStructure> target = J2000.shiftedBy(100.);
    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);
}
Also used : DormandPrince853FieldIntegrator(org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator) Frame(org.orekit.frames.Frame) FieldSpacecraftState(org.orekit.propagation.FieldSpacecraftState) AdaptiveStepsizeIntegrator(org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator) DerivativeStructure(org.hipparchus.analysis.differentiation.DerivativeStructure) DSFactory(org.hipparchus.analysis.differentiation.DSFactory) PVCoordinates(org.orekit.utils.PVCoordinates) FieldPVCoordinates(org.orekit.utils.FieldPVCoordinates) FieldKeplerianOrbit(org.orekit.orbits.FieldKeplerianOrbit) SpacecraftState(org.orekit.propagation.SpacecraftState) FieldSpacecraftState(org.orekit.propagation.FieldSpacecraftState) FieldNumericalPropagator(org.orekit.propagation.numerical.FieldNumericalPropagator) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) FieldNumericalPropagator(org.orekit.propagation.numerical.FieldNumericalPropagator) OrbitType(org.orekit.orbits.OrbitType) DormandPrince853Integrator(org.hipparchus.ode.nonstiff.DormandPrince853Integrator) NormalizedSphericalHarmonicsProvider(org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) AbstractLegacyForceModelTest(org.orekit.forces.AbstractLegacyForceModelTest) Test(org.junit.Test)

Example 34 with AdaptiveStepsizeIntegrator

use of org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator in project Orekit by CS-SI.

the class RelativityTest 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(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);
    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);
    }
}
Also used : Frame(org.orekit.frames.Frame) GaussianRandomGenerator(org.hipparchus.random.GaussianRandomGenerator) AdaptiveStepsizeIntegrator(org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator) FieldPVCoordinates(org.orekit.utils.FieldPVCoordinates) PVCoordinates(org.orekit.utils.PVCoordinates) Well19937a(org.hipparchus.random.Well19937a) RandomGenerator(org.hipparchus.random.RandomGenerator) GaussianRandomGenerator(org.hipparchus.random.GaussianRandomGenerator) FieldKeplerianOrbit(org.orekit.orbits.FieldKeplerianOrbit) FieldSpacecraftState(org.orekit.propagation.FieldSpacecraftState) SpacecraftState(org.orekit.propagation.SpacecraftState) FieldNumericalPropagator(org.orekit.propagation.numerical.FieldNumericalPropagator) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) FieldKeplerianOrbit(org.orekit.orbits.FieldKeplerianOrbit) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) DormandPrince853Integrator(org.hipparchus.ode.nonstiff.DormandPrince853Integrator) DormandPrince853FieldIntegrator(org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator) FieldSpacecraftState(org.orekit.propagation.FieldSpacecraftState) DerivativeStructure(org.hipparchus.analysis.differentiation.DerivativeStructure) DSFactory(org.hipparchus.analysis.differentiation.DSFactory) FieldNumericalPropagator(org.orekit.propagation.numerical.FieldNumericalPropagator) OrbitType(org.orekit.orbits.OrbitType) UncorrelatedRandomVectorGenerator(org.hipparchus.random.UncorrelatedRandomVectorGenerator) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) AbstractLegacyForceModelTest(org.orekit.forces.AbstractLegacyForceModelTest) Test(org.junit.Test)

Example 35 with AdaptiveStepsizeIntegrator

use of org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator 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);
}
Also used : Frame(org.orekit.frames.Frame) GaussianRandomGenerator(org.hipparchus.random.GaussianRandomGenerator) AdaptiveStepsizeIntegrator(org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator) PVCoordinates(org.orekit.utils.PVCoordinates) FieldPVCoordinates(org.orekit.utils.FieldPVCoordinates) Well19937a(org.hipparchus.random.Well19937a) RandomGenerator(org.hipparchus.random.RandomGenerator) GaussianRandomGenerator(org.hipparchus.random.GaussianRandomGenerator) FieldKeplerianOrbit(org.orekit.orbits.FieldKeplerianOrbit) SpacecraftState(org.orekit.propagation.SpacecraftState) FieldSpacecraftState(org.orekit.propagation.FieldSpacecraftState) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) FieldNumericalPropagator(org.orekit.propagation.numerical.FieldNumericalPropagator) FieldKeplerianOrbit(org.orekit.orbits.FieldKeplerianOrbit) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) DormandPrince853Integrator(org.hipparchus.ode.nonstiff.DormandPrince853Integrator) DormandPrince853FieldIntegrator(org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator) FieldSpacecraftState(org.orekit.propagation.FieldSpacecraftState) DerivativeStructure(org.hipparchus.analysis.differentiation.DerivativeStructure) DSFactory(org.hipparchus.analysis.differentiation.DSFactory) FieldNumericalPropagator(org.orekit.propagation.numerical.FieldNumericalPropagator) OrbitType(org.orekit.orbits.OrbitType) UncorrelatedRandomVectorGenerator(org.hipparchus.random.UncorrelatedRandomVectorGenerator) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) AbstractLegacyForceModelTest(org.orekit.forces.AbstractLegacyForceModelTest) Test(org.junit.Test)

Aggregations

AdaptiveStepsizeIntegrator (org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator)51 DormandPrince853Integrator (org.hipparchus.ode.nonstiff.DormandPrince853Integrator)51 SpacecraftState (org.orekit.propagation.SpacecraftState)45 NumericalPropagator (org.orekit.propagation.numerical.NumericalPropagator)38 Test (org.junit.Test)34 AbsoluteDate (org.orekit.time.AbsoluteDate)28 KeplerianOrbit (org.orekit.orbits.KeplerianOrbit)26 PVCoordinates (org.orekit.utils.PVCoordinates)26 FieldSpacecraftState (org.orekit.propagation.FieldSpacecraftState)25 Orbit (org.orekit.orbits.Orbit)23 FieldAbsoluteDate (org.orekit.time.FieldAbsoluteDate)21 FieldNumericalPropagator (org.orekit.propagation.numerical.FieldNumericalPropagator)20 Vector3D (org.hipparchus.geometry.euclidean.threed.Vector3D)19 Frame (org.orekit.frames.Frame)19 FieldKeplerianOrbit (org.orekit.orbits.FieldKeplerianOrbit)19 AbstractLegacyForceModelTest (org.orekit.forces.AbstractLegacyForceModelTest)16 OrbitType (org.orekit.orbits.OrbitType)16 CartesianOrbit (org.orekit.orbits.CartesianOrbit)15 FieldPVCoordinates (org.orekit.utils.FieldPVCoordinates)15 DormandPrince853FieldIntegrator (org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator)14