Examples with DSSTForceModel org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel used on opensource projects

Example 11 with DSSTForceModel

use of org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel in project Orekit by CS-SI.

the class DSSTPropagator method beforeIntegration.

/**
 * Method called just before integration.
 * <p>
 * The default implementation does nothing, it may be specialized in subclasses.
 * </p>
 * @param initialState initial state
 * @param tEnd target date at which state should be propagated
 * @exception OrekitException if hook cannot be run
 */
@Override
protected void beforeIntegration(final SpacecraftState initialState, final AbsoluteDate tEnd) throws OrekitException {
    // compute common auxiliary elements
    final AuxiliaryElements aux = new AuxiliaryElements(initialState.getOrbit(), I);
    // check if only mean elements must be used
    final boolean meanOnly = isMeanOrbit();
    // initialize all perturbing forces
    final List<ShortPeriodTerms> shortPeriodTerms = new ArrayList<ShortPeriodTerms>();
    for (final DSSTForceModel force : forceModels) {
        shortPeriodTerms.addAll(force.initialize(aux, meanOnly));
    }
    mapper.setShortPeriodTerms(shortPeriodTerms);
    // if required, insert the special short periodics step handler
    if (!meanOnly) {
        final ShortPeriodicsHandler spHandler = new ShortPeriodicsHandler(forceModels);
        final Collection<ODEStepHandler> stepHandlers = new ArrayList<ODEStepHandler>();
        stepHandlers.add(spHandler);
        final ODEIntegrator integrator = getIntegrator();
        final Collection<ODEStepHandler> existing = integrator.getStepHandlers();
        stepHandlers.addAll(existing);
        integrator.clearStepHandlers();
        // add back the existing handlers after the short periodics one
        for (final ODEStepHandler sp : stepHandlers) {
            integrator.addStepHandler(sp);
        }
    }
}

Example 12 with DSSTForceModel

use of org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel in project Orekit by CS-SI.

the class DSSTPropagator method computeMeanOrbit.

/**
 * Compute mean state from osculating state.
 * <p>
 * Compute in a DSST sense the mean state corresponding to the input osculating state.
 * </p><p>
 * The computing is done through a fixed-point iteration process.
 * </p>
 * @param osculating initial osculating state
 * @param attitudeProvider attitude provider (may be null if there are no Gaussian force models
 * like atmospheric drag, radiation pressure or specific user-defined models)
 * @param forceModels force models
 * @return mean state
 * @throws OrekitException if the underlying computation of short periodic variation fails
 */
private static Orbit computeMeanOrbit(final SpacecraftState osculating, final AttitudeProvider attitudeProvider, final Collection<DSSTForceModel> forceModels) throws OrekitException {
    // rough initialization of the mean parameters
    EquinoctialOrbit meanOrbit = (EquinoctialOrbit) OrbitType.EQUINOCTIAL.convertType(osculating.getOrbit());
    // threshold for each parameter
    final double epsilon = 1.0e-13;
    final double thresholdA = epsilon * (1 + FastMath.abs(meanOrbit.getA()));
    final double thresholdE = epsilon * (1 + meanOrbit.getE());
    final double thresholdI = epsilon * (1 + meanOrbit.getI());
    final double thresholdL = epsilon * FastMath.PI;
    // ensure all Gaussian force models can rely on attitude
    for (final DSSTForceModel force : forceModels) {
        force.registerAttitudeProvider(attitudeProvider);
    }
    int i = 0;
    while (i++ < 200) {
        final SpacecraftState meanState = new SpacecraftState(meanOrbit, osculating.getAttitude(), osculating.getMass());
        // Create the auxiliary object
        final AuxiliaryElements aux = new AuxiliaryElements(meanOrbit, I);
        // Set the force models
        final List<ShortPeriodTerms> shortPeriodTerms = new ArrayList<ShortPeriodTerms>();
        for (final DSSTForceModel force : forceModels) {
            shortPeriodTerms.addAll(force.initialize(aux, false));
            force.updateShortPeriodTerms(meanState);
        }
        // recompute the osculating parameters from the current mean parameters
        final EquinoctialOrbit rebuilt = computeOsculatingOrbit(meanState, shortPeriodTerms);
        // adapted parameters residuals
        final double deltaA = osculating.getA() - rebuilt.getA();
        final double deltaEx = osculating.getEquinoctialEx() - rebuilt.getEquinoctialEx();
        final double deltaEy = osculating.getEquinoctialEy() - rebuilt.getEquinoctialEy();
        final double deltaHx = osculating.getHx() - rebuilt.getHx();
        final double deltaHy = osculating.getHy() - rebuilt.getHy();
        final double deltaLv = MathUtils.normalizeAngle(osculating.getLv() - rebuilt.getLv(), 0.0);
        // check convergence
        if (FastMath.abs(deltaA) < thresholdA && FastMath.abs(deltaEx) < thresholdE && FastMath.abs(deltaEy) < thresholdE && FastMath.abs(deltaHx) < thresholdI && FastMath.abs(deltaHy) < thresholdI && FastMath.abs(deltaLv) < thresholdL) {
            return meanOrbit;
        }
        // update mean parameters
        meanOrbit = new EquinoctialOrbit(meanOrbit.getA() + deltaA, meanOrbit.getEquinoctialEx() + deltaEx, meanOrbit.getEquinoctialEy() + deltaEy, meanOrbit.getHx() + deltaHx, meanOrbit.getHy() + deltaHy, meanOrbit.getLv() + deltaLv, PositionAngle.TRUE, meanOrbit.getFrame(), meanOrbit.getDate(), meanOrbit.getMu());
    }
    throw new OrekitException(OrekitMessages.UNABLE_TO_COMPUTE_DSST_MEAN_PARAMETERS, i);
}

Example 13 with DSSTForceModel

use of org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel in project Orekit by CS-SI.

the class DSSTPropagatorTest method testPropagationWithDrag.

@Test
public void testPropagationWithDrag() throws OrekitException {
    // Central Body geopotential 2x0
    final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(2, 0);
    final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
    DSSTForceModel zonal = new DSSTZonal(provider, 2, 0, 5);
    DSSTForceModel tesseral = new DSSTTesseral(earthFrame, Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider, 2, 0, 0, 2, 2, 0, 0);
    // Drag Force Model
    final OneAxisEllipsoid earth = new OneAxisEllipsoid(provider.getAe(), Constants.WGS84_EARTH_FLATTENING, earthFrame);
    final Atmosphere atm = new HarrisPriester(CelestialBodyFactory.getSun(), earth, 6);
    final double cd = 2.0;
    final double area = 25.0;
    DSSTForceModel drag = new DSSTAtmosphericDrag(atm, cd, area);
    // LEO Orbit
    final AbsoluteDate initDate = new AbsoluteDate(2003, 7, 1, 0, 0, 00.000, TimeScalesFactory.getUTC());
    final Orbit orbit = new KeplerianOrbit(7204535.848109440, 0.0012402238462686, FastMath.toRadians(98.74341600466740), FastMath.toRadians(111.1990175076630), FastMath.toRadians(43.32990110790340), FastMath.toRadians(68.66852509725620), PositionAngle.MEAN, FramesFactory.getEME2000(), initDate, provider.getMu());
    // Set propagator with state and force model
    setDSSTProp(new SpacecraftState(orbit));
    dsstProp.addForceModel(zonal);
    dsstProp.addForceModel(tesseral);
    dsstProp.addForceModel(drag);
    // 5 days propagation
    final SpacecraftState state = dsstProp.propagate(initDate.shiftedBy(5. * 86400.));
    // Ref Standalone_DSST:
    // a    = 7204521.657141485 m
    // h/ey =  0.0007093755541595772
    // k/ex = -0.001016800430994036
    // p/hy =  0.8698955648709271
    // q/hx =  0.7757573478894775
    // lM   = 193°0939742953394
    Assert.assertEquals(7204521.657141485, state.getA(), 6.e-1);
    Assert.assertEquals(-0.001016800430994036, state.getEquinoctialEx(), 5.e-8);
    Assert.assertEquals(0.0007093755541595772, state.getEquinoctialEy(), 2.e-8);
    Assert.assertEquals(0.7757573478894775, state.getHx(), 5.e-8);
    Assert.assertEquals(0.8698955648709271, state.getHy(), 5.e-8);
    Assert.assertEquals(193.0939742953394, FastMath.toDegrees(MathUtils.normalizeAngle(state.getLM(), FastMath.PI)), 2.e-3);
    // Assert.assertEquals(((DSSTAtmosphericDrag)drag).getCd(), cd, 1e-9);
    // Assert.assertEquals(((DSSTAtmosphericDrag)drag).getArea(), area, 1e-9);
    Assert.assertEquals(((DSSTAtmosphericDrag) drag).getAtmosphere(), atm);
    // DSSTAtmosphericDrag.ATMOSPHERE_ALTITUDE_MAX
    final double atmosphericMaxConstant = 1000000.0;
    Assert.assertEquals(((DSSTAtmosphericDrag) drag).getRbar(), atmosphericMaxConstant + Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 1e-9);
}

Example 14 with DSSTForceModel

use of org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel in project Orekit by CS-SI.

the class DSSTPropagatorTest method testMeanToOsculatingState.

@Test
public void testMeanToOsculatingState() throws IllegalArgumentException, OrekitException {
    final SpacecraftState meanState = getGEOState();
    final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(2, 0);
    final Frame earthFrame = CelestialBodyFactory.getEarth().getBodyOrientedFrame();
    final DSSTForceModel zonal = new DSSTZonal(provider, 2, 1, 5);
    final DSSTForceModel tesseral = new DSSTTesseral(earthFrame, Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider, 2, 0, 0, 2, 2, 0, 0);
    final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
    forces.add(zonal);
    forces.add(tesseral);
    final SpacecraftState osculatingState = DSSTPropagator.computeOsculatingState(meanState, null, forces);
    Assert.assertEquals(1559.1, Vector3D.distance(meanState.getPVCoordinates().getPosition(), osculatingState.getPVCoordinates().getPosition()), 1.0);
}

Example 15 with DSSTForceModel

use of org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel in project Orekit by CS-SI.

the class DSSTPropagatorTest method testHighDegreesSetting.

@Test
public void testHighDegreesSetting() throws OrekitException {
    Utils.setDataRoot("regular-data:potential/grgs-format");
    GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
    int earthDegree = 36;
    int earthOrder = 36;
    int eccPower = 4;
    final UnnormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getUnnormalizedProvider(earthDegree, earthOrder);
    final org.orekit.frames.Frame earthFrame = // terrestrial frame
    FramesFactory.getITRF(IERSConventions.IERS_2010, true);
    final DSSTForceModel force = new DSSTTesseral(earthFrame, Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider, earthDegree, earthOrder, eccPower, earthDegree + eccPower, earthDegree, earthOrder, eccPower);
    final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
    forces.add(force);
    TimeScale tai = TimeScalesFactory.getTAI();
    AbsoluteDate initialDate = new AbsoluteDate("2015-07-01", tai);
    Frame eci = FramesFactory.getGCRF();
    KeplerianOrbit orbit = new KeplerianOrbit(7120000.0, 1.0e-3, FastMath.toRadians(60.0), FastMath.toRadians(120.0), FastMath.toRadians(47.0), FastMath.toRadians(12.0), PositionAngle.TRUE, eci, initialDate, Constants.EIGEN5C_EARTH_MU);
    SpacecraftState oscuState = DSSTPropagator.computeOsculatingState(new SpacecraftState(orbit), null, forces);
    Assert.assertEquals(7119927.097122, oscuState.getA(), 0.001);
}