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

Example 21 with HolmesFeatherstoneAttractionModel

use of org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel 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);
}

Example 22 with HolmesFeatherstoneAttractionModel

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

the class NumericalConverterTest method checkFit.

protected void checkFit(final Orbit orbit, final double duration, final double stepSize, final double threshold, final double expectedRMS, final String... freeParameters) throws OrekitException, IOException, ParseException {
    NumericalPropagatorBuilder builder = new NumericalPropagatorBuilder(OrbitType.CARTESIAN.convertType(orbit), new DormandPrince853IntegratorBuilder(minStep, maxStep, dP), PositionAngle.TRUE, dP);
    ForceModel guessedDrag = drag;
    ForceModel guessedGravity = gravity;
    for (String param : freeParameters) {
        if (DragSensitive.DRAG_COEFFICIENT.equals(param)) {
            // we want to adjust drag coefficient, we need to start from a wrong value
            ParameterDriver driver = drag.getParameterDriver(param);
            double coeff = driver.getReferenceValue() - driver.getScale();
            guessedDrag = new DragForce(atmosphere, new IsotropicDrag(crossSection, coeff));
        } else if (NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT.equals(param)) {
            // we want to adjust mu, we need to start from  a wrong value
            guessedGravity = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), GravityFieldFactory.getNormalizedProvider(2, 0));
            ParameterDriver driver = guessedGravity.getParameterDriver(param);
            driver.setValue(driver.getReferenceValue() + driver.getScale());
        }
    }
    builder.addForceModel(guessedDrag);
    builder.addForceModel(guessedGravity);
    JacobianPropagatorConverter fitter = new JacobianPropagatorConverter(builder, threshold, 5000);
    fitter.convert(propagator, duration, 1 + (int) (duration / stepSize), freeParameters);
    NumericalPropagator prop = (NumericalPropagator) fitter.getAdaptedPropagator();
    Orbit fitted = prop.getInitialState().getOrbit();
    for (String param : freeParameters) {
        for (ForceModel force : propagator.getAllForceModels()) {
            if (force.isSupported(param)) {
                for (ForceModel model : prop.getAllForceModels()) {
                    if (model.isSupported(param)) {
                        Assert.assertEquals(force.getParameterDriver(param).getValue(), model.getParameterDriver(param).getValue(), 3.0e-4 * FastMath.abs(force.getParameterDriver(param).getValue()));
                    }
                }
            }
        }
    }
    Assert.assertEquals(expectedRMS, fitter.getRMS(), 0.01 * expectedRMS);
    Assert.assertEquals(orbit.getPVCoordinates().getPosition().getX(), fitted.getPVCoordinates().getPosition().getX(), 1.1);
    Assert.assertEquals(orbit.getPVCoordinates().getPosition().getY(), fitted.getPVCoordinates().getPosition().getY(), 1.1);
    Assert.assertEquals(orbit.getPVCoordinates().getPosition().getZ(), fitted.getPVCoordinates().getPosition().getZ(), 1.1);
    Assert.assertEquals(orbit.getPVCoordinates().getVelocity().getX(), fitted.getPVCoordinates().getVelocity().getX(), 0.0005);
    Assert.assertEquals(orbit.getPVCoordinates().getVelocity().getY(), fitted.getPVCoordinates().getVelocity().getY(), 0.0005);
    Assert.assertEquals(orbit.getPVCoordinates().getVelocity().getZ(), fitted.getPVCoordinates().getVelocity().getZ(), 0.0005);
}

Example 23 with HolmesFeatherstoneAttractionModel

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

the class NumericalConverterTest method setUp.

@Before
public void setUp() throws OrekitException, IOException, ParseException {
    Utils.setDataRoot("regular-data:potential/shm-format");
    gravity = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), GravityFieldFactory.getNormalizedProvider(2, 0));
    mu = gravity.getParameterDriver(NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT).getValue();
    minStep = 1.0;
    maxStep = 600.0;
    dP = 10.0;
    // use a orbit that comes close to Earth so the drag coefficient has an effect
    final Vector3D position = new Vector3D(7.0e6, 1.0e6, 4.0e6).normalize().scalarMultiply(Constants.WGS84_EARTH_EQUATORIAL_RADIUS + 300e3);
    final Vector3D velocity = new Vector3D(-500.0, 8000.0, 1000.0);
    final AbsoluteDate initDate = new AbsoluteDate(2010, 10, 10, 10, 10, 10.0, TimeScalesFactory.getUTC());
    orbit = new EquinoctialOrbit(new PVCoordinates(position, velocity), FramesFactory.getEME2000(), initDate, mu);
    final double[][] tol = NumericalPropagator.tolerances(dP, orbit, OrbitType.CARTESIAN);
    propagator = new NumericalPropagator(new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]));
    propagator.setInitialState(new SpacecraftState(orbit));
    propagator.setOrbitType(OrbitType.CARTESIAN);
    final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
    earth.setAngularThreshold(1.e-7);
    atmosphere = new SimpleExponentialAtmosphere(earth, 0.0004, 42000.0, 7500.0);
    final double dragCoef = 2.0;
    crossSection = 10.0;
    drag = new DragForce(atmosphere, new IsotropicDrag(crossSection, dragCoef));
    propagator.addForceModel(gravity);
    propagator.addForceModel(drag);
}

Example 24 with HolmesFeatherstoneAttractionModel

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

the class AttitudesSequenceTest method testResetDuringTransitionBackward.

@Test
public void testResetDuringTransitionBackward() throws OrekitException {
    // Initial state definition : date, orbit
    final AbsoluteDate initialDate = new AbsoluteDate(2004, 01, 01, 23, 30, 00.000, TimeScalesFactory.getUTC());
    final Vector3D position = new Vector3D(-6142438.668, 3492467.560, -25767.25680);
    final Vector3D velocity = new Vector3D(505.8479685, 942.7809215, 7435.922231);
    final Orbit initialOrbit = new KeplerianOrbit(new PVCoordinates(position, velocity), FramesFactory.getEME2000(), initialDate, Constants.EIGEN5C_EARTH_MU);
    final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
    final TopocentricFrame volgograd = new TopocentricFrame(earth, new GeodeticPoint(FastMath.toRadians(48.7), FastMath.toRadians(44.5), 24.0), "Волгоград");
    final AttitudesSequence attitudesSequence = new AttitudesSequence();
    final double transitionTime = 250.0;
    final AttitudeProvider nadirPointing = new NadirPointing(initialOrbit.getFrame(), earth);
    final AttitudeProvider targetPointing = new TargetPointing(initialOrbit.getFrame(), volgograd.getPoint(), earth);
    final ElevationDetector eventDetector = new ElevationDetector(volgograd).withConstantElevation(FastMath.toRadians(5.0)).withHandler(new ContinueOnEvent<>());
    final List<AbsoluteDate> nadirToTarget = new ArrayList<>();
    attitudesSequence.addSwitchingCondition(nadirPointing, targetPointing, eventDetector, true, false, transitionTime, AngularDerivativesFilter.USE_RR, (previous, next, state) -> nadirToTarget.add(state.getDate()));
    final double[][] tolerance = NumericalPropagator.tolerances(10.0, initialOrbit, initialOrbit.getType());
    final AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 300.0, tolerance[0], tolerance[1]);
    final NumericalPropagator propagator = new NumericalPropagator(integrator);
    GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("g007_eigen_05c_coef", false));
    propagator.addForceModel(new HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), GravityFieldFactory.getNormalizedProvider(8, 8)));
    propagator.setInitialState(new SpacecraftState(initialOrbit, nadirPointing.getAttitude(initialOrbit, initialOrbit.getDate(), initialOrbit.getFrame())));
    propagator.setAttitudeProvider(attitudesSequence);
    attitudesSequence.registerSwitchEvents(propagator);
    propagator.propagate(initialDate.shiftedBy(6000));
    // check that if we restart a backward propagation from an intermediate state
    // we properly get an interpolated attitude despite we missed the event trigger
    final AbsoluteDate midTransition = nadirToTarget.get(0).shiftedBy(0.5 * transitionTime);
    SpacecraftState state = propagator.propagate(midTransition.shiftedBy(+60), midTransition);
    Rotation nadirR = nadirPointing.getAttitude(state.getOrbit(), state.getDate(), state.getFrame()).getRotation();
    Rotation targetR = targetPointing.getAttitude(state.getOrbit(), state.getDate(), state.getFrame()).getRotation();
    final double reorientationAngle = Rotation.distance(nadirR, targetR);
    Assert.assertEquals(0.5 * reorientationAngle, Rotation.distance(state.getAttitude().getRotation(), targetR), 0.03 * reorientationAngle);
}

Example 25 with HolmesFeatherstoneAttractionModel

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

the class AdapterPropagatorTest method getEphemeris.

private BoundedPropagator getEphemeris(final Orbit orbit, final double mass, final int nbOrbits, final AttitudeProvider law, final AbsoluteDate t0, final Vector3D dV, final double f, final double isp, final boolean sunAttraction, final boolean moonAttraction, final NormalizedSphericalHarmonicsProvider gravityField) throws OrekitException, ParseException, IOException {
    SpacecraftState initialState = new SpacecraftState(orbit, law.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), mass);
    // add some dummy additional states
    initialState = initialState.addAdditionalState("dummy 1", 1.25, 2.5);
    initialState = initialState.addAdditionalState("dummy 2", 5.0);
    // set up numerical propagator
    final double dP = 1.0;
    double[][] tolerances = NumericalPropagator.tolerances(dP, orbit, orbit.getType());
    AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 1000, tolerances[0], tolerances[1]);
    integrator.setInitialStepSize(orbit.getKeplerianPeriod() / 100.0);
    final NumericalPropagator propagator = new NumericalPropagator(integrator);
    propagator.addAdditionalStateProvider(new AdditionalStateProvider() {

        public String getName() {
            return "dummy 2";
        }

        public double[] getAdditionalState(SpacecraftState state) {
            return new double[] { 5.0 };
        }
    });
    propagator.setInitialState(initialState);
    propagator.setAttitudeProvider(law);
    if (dV.getNorm() > 1.0e-6) {
        // set up maneuver
        final double vExhaust = Constants.G0_STANDARD_GRAVITY * isp;
        final double dt = -(mass * vExhaust / f) * FastMath.expm1(-dV.getNorm() / vExhaust);
        final ConstantThrustManeuver maneuver = new ConstantThrustManeuver(t0.shiftedBy(-0.5 * dt), dt, f, isp, dV.normalize());
        propagator.addForceModel(maneuver);
    }
    if (sunAttraction) {
        propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
    }
    if (moonAttraction) {
        propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
    }
    if (gravityField != null) {
        propagator.addForceModel(new HolmesFeatherstoneAttractionModel(FramesFactory.getGTOD(false), gravityField));
    }
    propagator.setEphemerisMode();
    propagator.propagate(t0.shiftedBy(nbOrbits * orbit.getKeplerianPeriod()));
    final BoundedPropagator ephemeris = propagator.getGeneratedEphemeris();
    // both the initial propagator and generated ephemeris manage one of the two
    // additional states, but they also contain unmanaged copies of the other one
    Assert.assertFalse(propagator.isAdditionalStateManaged("dummy 1"));
    Assert.assertTrue(propagator.isAdditionalStateManaged("dummy 2"));
    Assert.assertFalse(ephemeris.isAdditionalStateManaged("dummy 1"));
    Assert.assertTrue(ephemeris.isAdditionalStateManaged("dummy 2"));
    Assert.assertEquals(2, ephemeris.getInitialState().getAdditionalState("dummy 1").length);
    Assert.assertEquals(1, ephemeris.getInitialState().getAdditionalState("dummy 2").length);
    return ephemeris;
}