Examples with HolmesFeatherstoneAttractionModel org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel used on opensource projects

Example 16 with HolmesFeatherstoneAttractionModel

use of org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel in project Orekit by CS-SI.

the class PartialDerivativesTest method doTestParametersDerivatives.

private void doTestParametersDerivatives(String parameterName, double tolerance, OrbitType... orbitTypes) throws OrekitException {
    OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
    ForceModel drag = new DragForce(new HarrisPriester(CelestialBodyFactory.getSun(), earth), new IsotropicDrag(2.5, 1.2));
    NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(5, 5);
    ForceModel gravityField = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
    Orbit baseOrbit = new KeplerianOrbit(7000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.TRUE, FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, provider.getMu());
    double dt = 900;
    double dP = 1.0;
    for (OrbitType orbitType : orbitTypes) {
        final Orbit initialOrbit = orbitType.convertType(baseOrbit);
        for (PositionAngle angleType : PositionAngle.values()) {
            NumericalPropagator propagator = setUpPropagator(initialOrbit, dP, orbitType, angleType, drag, gravityField);
            propagator.setMu(provider.getMu());
            for (final ForceModel forceModel : propagator.getAllForceModels()) {
                for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
                    driver.setValue(driver.getReferenceValue());
                    driver.setSelected(driver.getName().equals(parameterName));
                }
            }
            PartialDerivativesEquations partials = new PartialDerivativesEquations("partials", propagator);
            final SpacecraftState initialState = partials.setInitialJacobians(new SpacecraftState(initialOrbit));
            propagator.setInitialState(initialState);
            final JacobiansMapper mapper = partials.getMapper();
            PickUpHandler pickUp = new PickUpHandler(mapper, null);
            propagator.setMasterMode(pickUp);
            propagator.propagate(initialState.getDate().shiftedBy(dt));
            double[][] dYdP = pickUp.getdYdP();
            // compute reference Jacobian using finite differences
            double[][] dYdPRef = new double[6][1];
            NumericalPropagator propagator2 = setUpPropagator(initialOrbit, dP, orbitType, angleType, drag, gravityField);
            propagator2.setMu(provider.getMu());
            ParameterDriversList bound = new ParameterDriversList();
            for (final ForceModel forceModel : propagator2.getAllForceModels()) {
                for (final ParameterDriver driver : forceModel.getParametersDrivers()) {
                    if (driver.getName().equals(parameterName)) {
                        driver.setSelected(true);
                        bound.add(driver);
                    } else {
                        driver.setSelected(false);
                    }
                }
            }
            ParameterDriver selected = bound.getDrivers().get(0);
            double p0 = selected.getReferenceValue();
            double h = selected.getScale();
            selected.setValue(p0 - 4 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sM4h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            selected.setValue(p0 - 3 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sM3h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            selected.setValue(p0 - 2 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sM2h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            selected.setValue(p0 - 1 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sM1h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            selected.setValue(p0 + 1 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sP1h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            selected.setValue(p0 + 2 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sP2h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            selected.setValue(p0 + 3 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sP3h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            selected.setValue(p0 + 4 * h);
            propagator2.resetInitialState(arrayToState(stateToArray(initialState, orbitType, angleType, true), orbitType, angleType, initialState.getFrame(), initialState.getDate(), // the mu may have been reset above
            propagator2.getMu(), initialState.getAttitude()));
            SpacecraftState sP4h = propagator2.propagate(initialOrbit.getDate().shiftedBy(dt));
            fillJacobianColumn(dYdPRef, 0, orbitType, angleType, h, sM4h, sM3h, sM2h, sM1h, sP1h, sP2h, sP3h, sP4h);
            for (int i = 0; i < 6; ++i) {
                Assert.assertEquals(dYdPRef[i][0], dYdP[i][0], FastMath.abs(dYdPRef[i][0] * tolerance));
            }
        }
    }
}

Example 17 with HolmesFeatherstoneAttractionModel

use of org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel in project Orekit by CS-SI.

the class PartialDerivativesTest method testPropagationTypesElliptical.

@Test
public void testPropagationTypesElliptical() throws OrekitException {
    NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(5, 5);
    ForceModel gravityField = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
    Orbit initialOrbit = new KeplerianOrbit(8000000.0, 0.01, 0.1, 0.7, 0, 1.2, PositionAngle.TRUE, FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, provider.getMu());
    double dt = 900;
    double dP = 0.001;
    for (OrbitType orbitType : OrbitType.values()) {
        for (PositionAngle angleType : PositionAngle.values()) {
            // compute state Jacobian using PartialDerivatives
            NumericalPropagator propagator = setUpPropagator(initialOrbit, dP, orbitType, angleType, gravityField);
            PartialDerivativesEquations partials = new PartialDerivativesEquations("partials", propagator);
            final SpacecraftState initialState = partials.setInitialJacobians(new SpacecraftState(initialOrbit));
            propagator.setInitialState(initialState);
            final JacobiansMapper mapper = partials.getMapper();
            PickUpHandler pickUp = new PickUpHandler(mapper, null);
            propagator.setMasterMode(pickUp);
            propagator.propagate(initialState.getDate().shiftedBy(dt));
            double[][] dYdY0 = pickUp.getdYdY0();
            // compute reference state Jacobian using finite differences
            double[][] dYdY0Ref = new double[6][6];
            AbstractIntegratedPropagator propagator2 = setUpPropagator(initialOrbit, dP, orbitType, angleType, gravityField);
            double[] steps = NumericalPropagator.tolerances(1000000 * dP, initialOrbit, orbitType)[0];
            for (int i = 0; i < 6; ++i) {
                propagator2.resetInitialState(shiftState(initialState, orbitType, angleType, -4 * steps[i], i));
                SpacecraftState sM4h = propagator2.propagate(initialState.getDate().shiftedBy(dt));
                propagator2.resetInitialState(shiftState(initialState, orbitType, angleType, -3 * steps[i], i));
                SpacecraftState sM3h = propagator2.propagate(initialState.getDate().shiftedBy(dt));
                propagator2.resetInitialState(shiftState(initialState, orbitType, angleType, -2 * steps[i], i));
                SpacecraftState sM2h = propagator2.propagate(initialState.getDate().shiftedBy(dt));
                propagator2.resetInitialState(shiftState(initialState, orbitType, angleType, -1 * steps[i], i));
                SpacecraftState sM1h = propagator2.propagate(initialState.getDate().shiftedBy(dt));
                propagator2.resetInitialState(shiftState(initialState, orbitType, angleType, 1 * steps[i], i));
                SpacecraftState sP1h = propagator2.propagate(initialState.getDate().shiftedBy(dt));
                propagator2.resetInitialState(shiftState(initialState, orbitType, angleType, 2 * steps[i], i));
                SpacecraftState sP2h = propagator2.propagate(initialState.getDate().shiftedBy(dt));
                propagator2.resetInitialState(shiftState(initialState, orbitType, angleType, 3 * steps[i], i));
                SpacecraftState sP3h = propagator2.propagate(initialState.getDate().shiftedBy(dt));
                propagator2.resetInitialState(shiftState(initialState, orbitType, angleType, 4 * steps[i], i));
                SpacecraftState sP4h = propagator2.propagate(initialState.getDate().shiftedBy(dt));
                fillJacobianColumn(dYdY0Ref, i, orbitType, angleType, steps[i], sM4h, sM3h, sM2h, sM1h, sP1h, sP2h, sP3h, sP4h);
            }
            for (int i = 0; i < 6; ++i) {
                for (int j = 0; j < 6; ++j) {
                    double error = FastMath.abs((dYdY0[i][j] - dYdY0Ref[i][j]) / dYdY0Ref[i][j]);
                    Assert.assertEquals(0, error, 6.0e-2);
                }
            }
        }
    }
}

Example 18 with HolmesFeatherstoneAttractionModel

use of org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel in project Orekit by CS-SI.

the class PartialDerivativesTest method testPropagationTypesHyperbolic.

@Test
public void testPropagationTypesHyperbolic() throws OrekitException {
    NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(5, 5);
    ForceModel gravityField = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
    Orbit initialOrbit = new KeplerianOrbit(new PVCoordinates(new Vector3D(-1551946.0, 708899.0, 6788204.0), new Vector3D(-9875.0, -3941.0, -1845.0)), FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, provider.getMu());
    double dt = 900;
    double dP = 0.001;
    for (OrbitType orbitType : new OrbitType[] { OrbitType.KEPLERIAN, OrbitType.CARTESIAN }) {
        for (PositionAngle angleType : PositionAngle.values()) {
            // compute state Jacobian using PartialDerivatives
            NumericalPropagator propagator = setUpPropagator(initialOrbit, dP, orbitType, angleType, gravityField);
            PartialDerivativesEquations partials = new PartialDerivativesEquations("partials", propagator);
            final SpacecraftState initialState = partials.setInitialJacobians(new SpacecraftState(initialOrbit));
            propagator.setInitialState(initialState);
            final JacobiansMapper mapper = partials.getMapper();
            PickUpHandler pickUp = new PickUpHandler(mapper, null);
            propagator.setMasterMode(pickUp);
            propagator.propagate(initialState.getDate().shiftedBy(dt));
            double[][] dYdY0 = pickUp.getdYdY0();
            // compute reference state Jacobian using finite differences
            double[][] dYdY0Ref = new double[6][6];
            AbstractIntegratedPropagator propagator2 = setUpPropagator(initialOrbit, dP, orbitType, angleType, gravityField);
            double[] steps = NumericalPropagator.tolerances(1000000 * dP, initialOrbit, orbitType)[0];
            for (int i = 0; i < 6; ++i) {
                propagator2.resetInitialState(shiftState(initialState, orbitType, angleType, -4 * steps[i], i));
                SpacecraftState sM4h = propagator2.propagate(initialState.getDate().shiftedBy(dt));
                propagator2.resetInitialState(shiftState(initialState, orbitType, angleType, -3 * steps[i], i));
                SpacecraftState sM3h = propagator2.propagate(initialState.getDate().shiftedBy(dt));
                propagator2.resetInitialState(shiftState(initialState, orbitType, angleType, -2 * steps[i], i));
                SpacecraftState sM2h = propagator2.propagate(initialState.getDate().shiftedBy(dt));
                propagator2.resetInitialState(shiftState(initialState, orbitType, angleType, -1 * steps[i], i));
                SpacecraftState sM1h = propagator2.propagate(initialState.getDate().shiftedBy(dt));
                propagator2.resetInitialState(shiftState(initialState, orbitType, angleType, 1 * steps[i], i));
                SpacecraftState sP1h = propagator2.propagate(initialState.getDate().shiftedBy(dt));
                propagator2.resetInitialState(shiftState(initialState, orbitType, angleType, 2 * steps[i], i));
                SpacecraftState sP2h = propagator2.propagate(initialState.getDate().shiftedBy(dt));
                propagator2.resetInitialState(shiftState(initialState, orbitType, angleType, 3 * steps[i], i));
                SpacecraftState sP3h = propagator2.propagate(initialState.getDate().shiftedBy(dt));
                propagator2.resetInitialState(shiftState(initialState, orbitType, angleType, 4 * steps[i], i));
                SpacecraftState sP4h = propagator2.propagate(initialState.getDate().shiftedBy(dt));
                fillJacobianColumn(dYdY0Ref, i, orbitType, angleType, steps[i], sM4h, sM3h, sM2h, sM1h, sP1h, sP2h, sP3h, sP4h);
            }
            for (int i = 0; i < 6; ++i) {
                for (int j = 0; j < 6; ++j) {
                    double error = FastMath.abs((dYdY0[i][j] - dYdY0Ref[i][j]) / dYdY0Ref[i][j]);
                    Assert.assertEquals(0, error, 1.0e-3);
                }
            }
        }
    }
}

Example 19 with HolmesFeatherstoneAttractionModel

use of org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel in project Orekit by CS-SI.

the class DSSTPropagatorTest method getLEOStatePropagatedBy30Minutes.

private SpacecraftState getLEOStatePropagatedBy30Minutes() throws IllegalArgumentException, OrekitException {
    final Vector3D position = new Vector3D(-6142438.668, 3492467.560, -25767.25680);
    final Vector3D velocity = new Vector3D(505.8479685, 942.7809215, 7435.922231);
    // Spring equinoxe 21st mars 2003 1h00m
    final AbsoluteDate initialDate = new AbsoluteDate(new DateComponents(2003, 03, 21), new TimeComponents(1, 0, 0.), TimeScalesFactory.getUTC());
    final CartesianOrbit osculatingOrbit = new CartesianOrbit(new PVCoordinates(position, velocity), FramesFactory.getTOD(IERSConventions.IERS_1996, false), initialDate, Constants.WGS84_EARTH_MU);
    // Adaptive step integrator
    // with a minimum step of 0.001 and a maximum step of 1000
    double minStep = 0.001;
    double maxstep = 1000.0;
    double positionTolerance = 10.0;
    OrbitType propagationType = OrbitType.EQUINOCTIAL;
    double[][] tolerances = NumericalPropagator.tolerances(positionTolerance, osculatingOrbit, propagationType);
    AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
    NumericalPropagator propagator = new NumericalPropagator(integrator);
    propagator.setOrbitType(propagationType);
    NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(5, 5);
    ForceModel holmesFeatherstone = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
    propagator.addForceModel(holmesFeatherstone);
    propagator.setInitialState(new SpacecraftState(osculatingOrbit));
    return propagator.propagate(new AbsoluteDate(initialDate, 1800.));
}

Example 20 with HolmesFeatherstoneAttractionModel

use of org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel in project Orekit by CS-SI.

the class IntegratedEphemerisTest method testSerializationNumerical.

@Test
public void testSerializationNumerical() throws OrekitException, IOException, ClassNotFoundException {
    AbsoluteDate finalDate = initialOrbit.getDate().shiftedBy(Constants.JULIAN_DAY);
    numericalPropagator.setEphemerisMode();
    numericalPropagator.setInitialState(new SpacecraftState(initialOrbit));
    final Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
    final NormalizedSphericalHarmonicsProvider gravity = GravityFieldFactory.getNormalizedProvider(8, 8);
    final CelestialBody sun = CelestialBodyFactory.getSun();
    final CelestialBody moon = CelestialBodyFactory.getMoon();
    final RadiationSensitive spacecraft = new IsotropicRadiationSingleCoefficient(20.0, 2.0);
    numericalPropagator.addForceModel(new HolmesFeatherstoneAttractionModel(itrf, gravity));
    numericalPropagator.addForceModel(new ThirdBodyAttraction(sun));
    numericalPropagator.addForceModel(new ThirdBodyAttraction(moon));
    numericalPropagator.addForceModel(new SolarRadiationPressure(sun, Constants.WGS84_EARTH_EQUATORIAL_RADIUS, spacecraft));
    numericalPropagator.propagate(finalDate);
    IntegratedEphemeris ephemeris = (IntegratedEphemeris) numericalPropagator.getGeneratedEphemeris();
    ByteArrayOutputStream bos = new ByteArrayOutputStream();
    ObjectOutputStream oos = new ObjectOutputStream(bos);
    oos.writeObject(ephemeris);
    int expectedSize = 258223;
    Assert.assertTrue("size = " + bos.size(), bos.size() > 9 * expectedSize / 10);
    Assert.assertTrue("size = " + bos.size(), bos.size() < 11 * expectedSize / 10);
    Assert.assertNotNull(ephemeris.getFrame());
    Assert.assertSame(ephemeris.getFrame(), numericalPropagator.getFrame());
    ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
    ObjectInputStream ois = new ObjectInputStream(bis);
    IntegratedEphemeris deserialized = (IntegratedEphemeris) ois.readObject();
    Assert.assertEquals(deserialized.getMinDate(), deserialized.getMinDate());
    Assert.assertEquals(deserialized.getMaxDate(), deserialized.getMaxDate());
}