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

Example 6 with ThirdBodyAttraction

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

the class OrbitDeterminationTest method createPropagatorBuilder.

/**
 * Create a propagator builder from input parameters
 * @param parser input file parser
 * @param conventions IERS conventions to use
 * @param gravityField gravity field
 * @param body central body
 * @param orbit first orbit estimate
 * @return propagator builder
 * @throws NoSuchElementException if input parameters are missing
 * @throws OrekitException if body frame cannot be created
 */
private NumericalPropagatorBuilder createPropagatorBuilder(final KeyValueFileParser<ParameterKey> parser, final IERSConventions conventions, final NormalizedSphericalHarmonicsProvider gravityField, final OneAxisEllipsoid body, final Orbit orbit) throws NoSuchElementException, OrekitException {
    final double minStep;
    if (!parser.containsKey(ParameterKey.PROPAGATOR_MIN_STEP)) {
        minStep = 0.001;
    } else {
        minStep = parser.getDouble(ParameterKey.PROPAGATOR_MIN_STEP);
    }
    final double maxStep;
    if (!parser.containsKey(ParameterKey.PROPAGATOR_MAX_STEP)) {
        maxStep = 300;
    } else {
        maxStep = parser.getDouble(ParameterKey.PROPAGATOR_MAX_STEP);
    }
    final double dP;
    if (!parser.containsKey(ParameterKey.PROPAGATOR_POSITION_ERROR)) {
        dP = 10.0;
    } else {
        dP = parser.getDouble(ParameterKey.PROPAGATOR_POSITION_ERROR);
    }
    final double positionScale;
    if (!parser.containsKey(ParameterKey.ESTIMATOR_ORBITAL_PARAMETERS_POSITION_SCALE)) {
        positionScale = dP;
    } else {
        positionScale = parser.getDouble(ParameterKey.ESTIMATOR_ORBITAL_PARAMETERS_POSITION_SCALE);
    }
    final NumericalPropagatorBuilder propagatorBuilder = new NumericalPropagatorBuilder(orbit, new DormandPrince853IntegratorBuilder(minStep, maxStep, dP), PositionAngle.MEAN, positionScale);
    // initial mass
    final double mass;
    if (!parser.containsKey(ParameterKey.MASS)) {
        mass = 1000.0;
    } else {
        mass = parser.getDouble(ParameterKey.MASS);
    }
    propagatorBuilder.setMass(mass);
    // gravity field force model
    propagatorBuilder.addForceModel(new HolmesFeatherstoneAttractionModel(body.getBodyFrame(), gravityField));
    // ocean tides force model
    if (parser.containsKey(ParameterKey.OCEAN_TIDES_DEGREE) && parser.containsKey(ParameterKey.OCEAN_TIDES_ORDER)) {
        final int degree = parser.getInt(ParameterKey.OCEAN_TIDES_DEGREE);
        final int order = parser.getInt(ParameterKey.OCEAN_TIDES_ORDER);
        if (degree > 0 && order > 0) {
            propagatorBuilder.addForceModel(new OceanTides(body.getBodyFrame(), gravityField.getAe(), gravityField.getMu(), degree, order, conventions, TimeScalesFactory.getUT1(conventions, true)));
        }
    }
    // solid tides force model
    List<CelestialBody> solidTidesBodies = new ArrayList<CelestialBody>();
    if (parser.containsKey(ParameterKey.SOLID_TIDES_SUN) && parser.getBoolean(ParameterKey.SOLID_TIDES_SUN)) {
        solidTidesBodies.add(CelestialBodyFactory.getSun());
    }
    if (parser.containsKey(ParameterKey.SOLID_TIDES_MOON) && parser.getBoolean(ParameterKey.SOLID_TIDES_MOON)) {
        solidTidesBodies.add(CelestialBodyFactory.getMoon());
    }
    if (!solidTidesBodies.isEmpty()) {
        propagatorBuilder.addForceModel(new SolidTides(body.getBodyFrame(), gravityField.getAe(), gravityField.getMu(), gravityField.getTideSystem(), conventions, TimeScalesFactory.getUT1(conventions, true), solidTidesBodies.toArray(new CelestialBody[solidTidesBodies.size()])));
    }
    // third body attraction
    if (parser.containsKey(ParameterKey.THIRD_BODY_SUN) && parser.getBoolean(ParameterKey.THIRD_BODY_SUN)) {
        propagatorBuilder.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
    }
    if (parser.containsKey(ParameterKey.THIRD_BODY_MOON) && parser.getBoolean(ParameterKey.THIRD_BODY_MOON)) {
        propagatorBuilder.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
    }
    // drag
    if (parser.containsKey(ParameterKey.DRAG) && parser.getBoolean(ParameterKey.DRAG)) {
        final double cd = parser.getDouble(ParameterKey.DRAG_CD);
        final double area = parser.getDouble(ParameterKey.DRAG_AREA);
        final boolean cdEstimated = parser.getBoolean(ParameterKey.DRAG_CD_ESTIMATED);
        MarshallSolarActivityFutureEstimation msafe = new MarshallSolarActivityFutureEstimation("(?:Jan|Feb|Mar|Apr|May|Jun|Jul|Aug|Sep|Oct|Nov|Dec)\\p{Digit}\\p{Digit}\\p{Digit}\\p{Digit}F10\\.(?:txt|TXT)", MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
        DataProvidersManager manager = DataProvidersManager.getInstance();
        manager.feed(msafe.getSupportedNames(), msafe);
        Atmosphere atmosphere = new DTM2000(msafe, CelestialBodyFactory.getSun(), body);
        propagatorBuilder.addForceModel(new DragForce(atmosphere, new IsotropicDrag(area, cd)));
        if (cdEstimated) {
            for (final ParameterDriver driver : propagatorBuilder.getPropagationParametersDrivers().getDrivers()) {
                if (driver.getName().equals(DragSensitive.DRAG_COEFFICIENT)) {
                    driver.setSelected(true);
                }
            }
        }
    }
    // solar radiation pressure
    if (parser.containsKey(ParameterKey.SOLAR_RADIATION_PRESSURE) && parser.getBoolean(ParameterKey.SOLAR_RADIATION_PRESSURE)) {
        final double cr = parser.getDouble(ParameterKey.SOLAR_RADIATION_PRESSURE_CR);
        final double area = parser.getDouble(ParameterKey.SOLAR_RADIATION_PRESSURE_AREA);
        final boolean cREstimated = parser.getBoolean(ParameterKey.SOLAR_RADIATION_PRESSURE_CR_ESTIMATED);
        propagatorBuilder.addForceModel(new SolarRadiationPressure(CelestialBodyFactory.getSun(), body.getEquatorialRadius(), new IsotropicRadiationSingleCoefficient(area, cr)));
        if (cREstimated) {
            for (final ParameterDriver driver : propagatorBuilder.getPropagationParametersDrivers().getDrivers()) {
                if (driver.getName().equals(RadiationSensitive.REFLECTION_COEFFICIENT)) {
                    driver.setSelected(true);
                }
            }
        }
    }
    // post-Newtonian correction force due to general relativity
    if (parser.containsKey(ParameterKey.GENERAL_RELATIVITY) && parser.getBoolean(ParameterKey.GENERAL_RELATIVITY)) {
        propagatorBuilder.addForceModel(new Relativity(gravityField.getMu()));
    }
    // extra polynomial accelerations
    if (parser.containsKey(ParameterKey.POLYNOMIAL_ACCELERATION_NAME)) {
        final String[] names = parser.getStringArray(ParameterKey.POLYNOMIAL_ACCELERATION_NAME);
        final Vector3D[] directions = parser.getVectorArray(ParameterKey.POLYNOMIAL_ACCELERATION_DIRECTION_X, ParameterKey.POLYNOMIAL_ACCELERATION_DIRECTION_Y, ParameterKey.POLYNOMIAL_ACCELERATION_DIRECTION_Z);
        final List<String>[] coefficients = parser.getStringsListArray(ParameterKey.POLYNOMIAL_ACCELERATION_COEFFICIENTS, ',');
        final boolean[] estimated = parser.getBooleanArray(ParameterKey.POLYNOMIAL_ACCELERATION_ESTIMATED);
        for (int i = 0; i < names.length; ++i) {
            final PolynomialParametricAcceleration ppa = new PolynomialParametricAcceleration(directions[i], true, names[i], null, coefficients[i].size() - 1);
            for (int k = 0; k < coefficients[i].size(); ++k) {
                final ParameterDriver driver = ppa.getParameterDriver(names[i] + "[" + k + "]");
                driver.setValue(Double.parseDouble(coefficients[i].get(k)));
                driver.setSelected(estimated[i]);
            }
            propagatorBuilder.addForceModel(ppa);
        }
    }
    return propagatorBuilder;
}

Example 7 with ThirdBodyAttraction

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

the class NumericalPropagatorTest method createPropagator.

private static synchronized NumericalPropagator createPropagator(SpacecraftState spacecraftState, OrbitType orbitType, PositionAngle angleType) throws OrekitException {
    final double minStep = 0.001;
    final double maxStep = 120.0;
    final double positionTolerance = 0.1;
    final int degree = 20;
    final int order = 20;
    final double spacecraftArea = 1.0;
    final double spacecraftDragCoefficient = 2.0;
    final double spacecraftReflectionCoefficient = 2.0;
    // propagator main configuration
    final double[][] tol = NumericalPropagator.tolerances(positionTolerance, spacecraftState.getOrbit(), orbitType);
    final ODEIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
    final NumericalPropagator np = new NumericalPropagator(integrator);
    np.setOrbitType(orbitType);
    np.setPositionAngleType(angleType);
    np.setInitialState(spacecraftState);
    // Earth gravity field
    final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
    final NormalizedSphericalHarmonicsProvider harmonicsGravityProvider = GravityFieldFactory.getNormalizedProvider(degree, order);
    np.addForceModel(new HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), harmonicsGravityProvider));
    // Sun and Moon attraction
    np.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
    np.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
    // atmospheric drag
    MarshallSolarActivityFutureEstimation msafe = new MarshallSolarActivityFutureEstimation("Jan2000F10-edited-data\\.txt", MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
    DataProvidersManager.getInstance().feed(msafe.getSupportedNames(), msafe);
    DTM2000 atmosphere = new DTM2000(msafe, CelestialBodyFactory.getSun(), earth);
    np.addForceModel(new DragForce(atmosphere, new IsotropicDrag(spacecraftArea, spacecraftDragCoefficient)));
    // solar radiation pressure
    np.addForceModel(new SolarRadiationPressure(CelestialBodyFactory.getSun(), earth.getEquatorialRadius(), new IsotropicRadiationSingleCoefficient(spacecraftArea, spacecraftReflectionCoefficient)));
    return np;
}

Example 8 with ThirdBodyAttraction

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

Example 9 with ThirdBodyAttraction

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

Example 10 with ThirdBodyAttraction

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

the class Phasing method run.

private void run(final File input) throws IOException, IllegalArgumentException, ParseException, OrekitException {
    // read input parameters
    KeyValueFileParser<ParameterKey> parser = new KeyValueFileParser<ParameterKey>(ParameterKey.class);
    try (final FileInputStream fis = new FileInputStream(input)) {
        parser.parseInput(input.getAbsolutePath(), fis);
    }
    TimeScale utc = TimeScalesFactory.getUTC();
    // simulation properties
    AbsoluteDate date = parser.getDate(ParameterKey.ORBIT_DATE, utc);
    int nbOrbits = parser.getInt(ParameterKey.PHASING_ORBITS_NUMBER);
    int nbDays = parser.getInt(ParameterKey.PHASING_DAYS_NUMBER);
    double latitude = parser.getAngle(ParameterKey.SUN_SYNCHRONOUS_REFERENCE_LATITUDE);
    boolean ascending = parser.getBoolean(ParameterKey.SUN_SYNCHRONOUS_REFERENCE_ASCENDING);
    double mst = parser.getTime(ParameterKey.SUN_SYNCHRONOUS_MEAN_SOLAR_TIME).getSecondsInUTCDay() / 3600;
    int degree = parser.getInt(ParameterKey.GRAVITY_FIELD_DEGREE);
    int order = parser.getInt(ParameterKey.GRAVITY_FIELD_ORDER);
    String gridOutput = parser.getString(ParameterKey.GRID_OUTPUT);
    double[] gridLatitudes = new double[] { parser.getAngle(ParameterKey.GRID_LATITUDE_1), parser.getAngle(ParameterKey.GRID_LATITUDE_2), parser.getAngle(ParameterKey.GRID_LATITUDE_3), parser.getAngle(ParameterKey.GRID_LATITUDE_4), parser.getAngle(ParameterKey.GRID_LATITUDE_5) };
    boolean[] gridAscending = new boolean[] { parser.getBoolean(ParameterKey.GRID_ASCENDING_1), parser.getBoolean(ParameterKey.GRID_ASCENDING_2), parser.getBoolean(ParameterKey.GRID_ASCENDING_3), parser.getBoolean(ParameterKey.GRID_ASCENDING_4), parser.getBoolean(ParameterKey.GRID_ASCENDING_5) };
    gravityField = GravityFieldFactory.getNormalizedProvider(degree, order);
    // initial guess for orbit
    CircularOrbit orbit = guessOrbit(date, FramesFactory.getEME2000(), nbOrbits, nbDays, latitude, ascending, mst);
    System.out.println("initial orbit: " + orbit);
    System.out.println("please wait while orbit is adjusted...");
    System.out.println();
    // numerical model for improving orbit
    double[][] tolerances = NumericalPropagator.tolerances(0.1, orbit, OrbitType.CIRCULAR);
    DormandPrince853Integrator integrator = new DormandPrince853Integrator(1.0e-4 * orbit.getKeplerianPeriod(), 1.0e-1 * orbit.getKeplerianPeriod(), tolerances[0], tolerances[1]);
    integrator.setInitialStepSize(1.0e-2 * orbit.getKeplerianPeriod());
    NumericalPropagator propagator = new NumericalPropagator(integrator);
    propagator.addForceModel(new HolmesFeatherstoneAttractionModel(FramesFactory.getGTOD(IERSConventions.IERS_2010, true), gravityField));
    propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
    propagator.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
    double deltaP = Double.POSITIVE_INFINITY;
    double deltaV = Double.POSITIVE_INFINITY;
    int counter = 0;
    DecimalFormat f = new DecimalFormat("0.000E00", new DecimalFormatSymbols(Locale.US));
    while (deltaP > 3.0e-1 || deltaV > 3.0e-4) {
        CircularOrbit previous = orbit;
        CircularOrbit tmp1 = improveEarthPhasing(previous, nbOrbits, nbDays, propagator);
        CircularOrbit tmp2 = improveSunSynchronization(tmp1, nbOrbits * tmp1.getKeplerianPeriod(), latitude, ascending, mst, propagator);
        orbit = improveFrozenEccentricity(tmp2, nbOrbits * tmp2.getKeplerianPeriod(), propagator);
        double da = orbit.getA() - previous.getA();
        double dex = orbit.getCircularEx() - previous.getCircularEx();
        double dey = orbit.getCircularEy() - previous.getCircularEy();
        double di = FastMath.toDegrees(orbit.getI() - previous.getI());
        double dr = FastMath.toDegrees(orbit.getRightAscensionOfAscendingNode() - previous.getRightAscensionOfAscendingNode());
        System.out.println(" iteration " + (++counter) + ": deltaA = " + f.format(da) + " m, deltaEx = " + f.format(dex) + ", deltaEy = " + f.format(dey) + ", deltaI = " + f.format(di) + " deg, deltaRAAN = " + f.format(dr) + " deg");
        PVCoordinates delta = new PVCoordinates(previous.getPVCoordinates(), orbit.getPVCoordinates());
        deltaP = delta.getPosition().getNorm();
        deltaV = delta.getVelocity().getNorm();
    }
    // final orbit
    System.out.println();
    System.out.println("final orbit (osculating): " + orbit);
    // generate the ground track grid file
    try (PrintStream output = new PrintStream(new File(input.getParent(), gridOutput), "UTF-8")) {
        for (int i = 0; i < gridLatitudes.length; ++i) {
            printGridPoints(output, gridLatitudes[i], gridAscending[i], orbit, propagator, nbOrbits);
        }
    }
}