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

Example 26 with HolmesFeatherstoneAttractionModel

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

Example 27 with HolmesFeatherstoneAttractionModel

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

the class OrbitDetermination 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 28 with HolmesFeatherstoneAttractionModel

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

the class PropagatorConversion method main.

/**
 * Program entry point.
 * @param args program arguments (unused here)
 */
public static void main(String[] args) {
    try {
        // configure Orekit
        File home = new File(System.getProperty("user.home"));
        File orekitData = new File(home, "orekit-data");
        if (!orekitData.exists()) {
            System.err.format(Locale.US, "Failed to find %s folder%n", orekitData.getAbsolutePath());
            System.err.format(Locale.US, "You need to download %s from the %s page and unzip it in %s for this tutorial to work%n", "orekit-data.zip", "https://www.orekit.org/forge/projects/orekit/files", home.getAbsolutePath());
            System.exit(1);
        }
        DataProvidersManager manager = DataProvidersManager.getInstance();
        manager.addProvider(new DirectoryCrawler(orekitData));
        // gravity field
        NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(2, 0);
        double mu = provider.getMu();
        // inertial frame
        Frame inertialFrame = FramesFactory.getEME2000();
        // Initial date
        AbsoluteDate initialDate = new AbsoluteDate(2004, 01, 01, 23, 30, 00.000, TimeScalesFactory.getUTC());
        // Initial orbit (GTO)
        // semi major axis in meters
        final double a = 24396159;
        // eccentricity
        final double e = 0.72831215;
        // inclination
        final double i = FastMath.toRadians(7);
        // perigee argument
        final double omega = FastMath.toRadians(180);
        // right ascention of ascending node
        final double raan = FastMath.toRadians(261);
        // mean anomaly
        final double lM = 0;
        Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.MEAN, inertialFrame, initialDate, mu);
        final double period = initialOrbit.getKeplerianPeriod();
        // Initial state definition
        final SpacecraftState initialState = new SpacecraftState(initialOrbit);
        // Adaptive step integrator with a minimum step of 0.001 and a maximum step of 1000
        final double minStep = 0.001;
        final double maxStep = 1000.;
        final double dP = 1.e-2;
        final OrbitType orbType = OrbitType.CARTESIAN;
        final double[][] tol = NumericalPropagator.tolerances(dP, initialOrbit, orbType);
        final AbstractIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tol[0], tol[1]);
        // Propagator
        NumericalPropagator numProp = new NumericalPropagator(integrator);
        numProp.setInitialState(initialState);
        numProp.setOrbitType(orbType);
        // Force Models:
        // 1 - Perturbing gravity field (only J2 is considered here)
        ForceModel gravity = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
        // Add force models to the propagator
        numProp.addForceModel(gravity);
        // Propagator factory
        PropagatorBuilder builder = new KeplerianPropagatorBuilder(initialOrbit, PositionAngle.TRUE, dP);
        // Propagator converter
        PropagatorConverter fitter = new FiniteDifferencePropagatorConverter(builder, 1.e-6, 5000);
        // Resulting propagator
        KeplerianPropagator kepProp = (KeplerianPropagator) fitter.convert(numProp, 2 * period, 251);
        // Step handlers
        StatesHandler numStepHandler = new StatesHandler();
        StatesHandler kepStepHandler = new StatesHandler();
        // Set up operating mode for the propagator as master mode
        // with fixed step and specialized step handler
        numProp.setMasterMode(60., numStepHandler);
        kepProp.setMasterMode(60., kepStepHandler);
        // Extrapolate from the initial to the final date
        numProp.propagate(initialDate.shiftedBy(10. * period));
        kepProp.propagate(initialDate.shiftedBy(10. * period));
        // retrieve the states
        List<SpacecraftState> numStates = numStepHandler.getStates();
        List<SpacecraftState> kepStates = kepStepHandler.getStates();
        // Print the results on the output file
        File output = new File(new File(System.getProperty("user.home")), "elements.dat");
        try (final PrintStream stream = new PrintStream(output, "UTF-8")) {
            stream.println("# date Anum Akep Enum Ekep Inum Ikep LMnum LMkep");
            for (SpacecraftState numState : numStates) {
                for (SpacecraftState kepState : kepStates) {
                    if (numState.getDate().compareTo(kepState.getDate()) == 0) {
                        stream.println(numState.getDate() + " " + numState.getA() + " " + kepState.getA() + " " + numState.getE() + " " + kepState.getE() + " " + FastMath.toDegrees(numState.getI()) + " " + FastMath.toDegrees(kepState.getI()) + " " + FastMath.toDegrees(MathUtils.normalizeAngle(numState.getLM(), FastMath.PI)) + " " + FastMath.toDegrees(MathUtils.normalizeAngle(kepState.getLM(), FastMath.PI)));
                        break;
                    }
                }
            }
        }
        System.out.println("Results saved as file " + output);
        File output1 = new File(new File(System.getProperty("user.home")), "elts_pv.dat");
        try (final PrintStream stream = new PrintStream(output1, "UTF-8")) {
            stream.println("# date pxn pyn pzn vxn vyn vzn pxk pyk pzk vxk vyk vzk");
            for (SpacecraftState numState : numStates) {
                for (SpacecraftState kepState : kepStates) {
                    if (numState.getDate().compareTo(kepState.getDate()) == 0) {
                        final double pxn = numState.getPVCoordinates().getPosition().getX();
                        final double pyn = numState.getPVCoordinates().getPosition().getY();
                        final double pzn = numState.getPVCoordinates().getPosition().getZ();
                        final double vxn = numState.getPVCoordinates().getVelocity().getX();
                        final double vyn = numState.getPVCoordinates().getVelocity().getY();
                        final double vzn = numState.getPVCoordinates().getVelocity().getZ();
                        final double pxk = kepState.getPVCoordinates().getPosition().getX();
                        final double pyk = kepState.getPVCoordinates().getPosition().getY();
                        final double pzk = kepState.getPVCoordinates().getPosition().getZ();
                        final double vxk = kepState.getPVCoordinates().getVelocity().getX();
                        final double vyk = kepState.getPVCoordinates().getVelocity().getY();
                        final double vzk = kepState.getPVCoordinates().getVelocity().getZ();
                        stream.println(numState.getDate() + " " + pxn + " " + pyn + " " + pzn + " " + vxn + " " + vyn + " " + vzn + " " + pxk + " " + pyk + " " + pzk + " " + vxk + " " + vyk + " " + vzk);
                        break;
                    }
                }
            }
        }
        System.out.println("Results saved as file " + output1);
    } catch (OrekitException oe) {
        System.err.println(oe.getLocalizedMessage());
        System.exit(1);
    } catch (IOException ioe) {
        System.err.println(ioe.getLocalizedMessage());
        System.exit(1);
    }
}

Example 29 with HolmesFeatherstoneAttractionModel

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

the class DSSTPropagation method setForceModel.

/**
 * Set numerical propagator force models
 *
 *  @param parser  input file parser
 *  @param normalized spherical harmonics provider
 *  @param earthFrame Earth rotating frame
 *  @param numProp numerical propagator
 *  @throws IOException
 *  @throws OrekitException
 */
private void setForceModel(final KeyValueFileParser<ParameterKey> parser, final NormalizedSphericalHarmonicsProvider normalized, final Frame earthFrame, final NumericalPropagator numProp) throws IOException, OrekitException {
    final double ae = normalized.getAe();
    final int degree = parser.getInt(ParameterKey.CENTRAL_BODY_DEGREE);
    final int order = parser.getInt(ParameterKey.CENTRAL_BODY_ORDER);
    if (order > degree) {
        throw new IOException("Potential order cannot be higher than potential degree");
    }
    // Central Body (normalized coefficients)
    numProp.addForceModel(new HolmesFeatherstoneAttractionModel(earthFrame, normalized));
    // 3rd body (SUN)
    if (parser.containsKey(ParameterKey.THIRD_BODY_SUN) && parser.getBoolean(ParameterKey.THIRD_BODY_SUN)) {
        numProp.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
    }
    // 3rd body (MOON)
    if (parser.containsKey(ParameterKey.THIRD_BODY_MOON) && parser.getBoolean(ParameterKey.THIRD_BODY_MOON)) {
        numProp.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
    }
    // Drag
    if (parser.containsKey(ParameterKey.DRAG) && parser.getBoolean(ParameterKey.DRAG)) {
        final OneAxisEllipsoid earth = new OneAxisEllipsoid(ae, Constants.WGS84_EARTH_FLATTENING, earthFrame);
        final Atmosphere atm = new HarrisPriester(CelestialBodyFactory.getSun(), earth, 6);
        final DragSensitive ssc = new IsotropicDrag(parser.getDouble(ParameterKey.DRAG_SF), parser.getDouble(ParameterKey.DRAG_CD));
        numProp.addForceModel(new DragForce(atm, ssc));
    }
    // 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 RadiationSensitive ssc = new IsotropicRadiationSingleCoefficient(parser.getDouble(ParameterKey.SOLAR_RADIATION_PRESSURE_SF), cR);
        numProp.addForceModel(new SolarRadiationPressure(CelestialBodyFactory.getSun(), ae, ssc));
    }
}

Example 30 with HolmesFeatherstoneAttractionModel

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

the class MasterMode method main.

/**
 * Program entry point.
 * @param args program arguments (unused here)
 */
public static void main(String[] args) {
    try {
        // configure Orekit
        File home = new File(System.getProperty("user.home"));
        File orekitData = new File(home, "orekit-data");
        if (!orekitData.exists()) {
            System.err.format(Locale.US, "Failed to find %s folder%n", orekitData.getAbsolutePath());
            System.err.format(Locale.US, "You need to download %s from the %s page and unzip it in %s for this tutorial to work%n", "orekit-data.zip", "https://www.orekit.org/forge/projects/orekit/files", home.getAbsolutePath());
            System.exit(1);
        }
        DataProvidersManager manager = DataProvidersManager.getInstance();
        manager.addProvider(new DirectoryCrawler(orekitData));
        // gravitation coefficient
        double mu = 3.986004415e+14;
        // inertial frame
        Frame inertialFrame = FramesFactory.getEME2000();
        // Initial date
        AbsoluteDate initialDate = new AbsoluteDate(2004, 01, 01, 23, 30, 00.000, TimeScalesFactory.getUTC());
        // Initial orbit
        // semi major axis in meters
        double a = 24396159;
        // eccentricity
        double e = 0.72831215;
        // inclination
        double i = FastMath.toRadians(7);
        // perigee argument
        double omega = FastMath.toRadians(180);
        // right ascention of ascending node
        double raan = FastMath.toRadians(261);
        // mean anomaly
        double lM = 0;
        Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.MEAN, inertialFrame, initialDate, mu);
        // Initial state definition
        SpacecraftState initialState = new SpacecraftState(initialOrbit);
        // Adaptive step integrator with a minimum step of 0.001 and a maximum step of 1000
        final double minStep = 0.001;
        final double maxstep = 1000.0;
        final double positionTolerance = 10.0;
        final OrbitType propagationType = OrbitType.KEPLERIAN;
        final double[][] tolerances = NumericalPropagator.tolerances(positionTolerance, initialOrbit, propagationType);
        AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]);
        // Propagator
        NumericalPropagator propagator = new NumericalPropagator(integrator);
        propagator.setOrbitType(propagationType);
        // Force Model (reduced to perturbing gravity field)
        final NormalizedSphericalHarmonicsProvider provider = GravityFieldFactory.getNormalizedProvider(10, 10);
        ForceModel holmesFeatherstone = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider);
        // Add force model to the propagator
        propagator.addForceModel(holmesFeatherstone);
        // Set up initial state in the propagator
        propagator.setInitialState(initialState);
        // Set up operating mode for the propagator as master mode
        // with fixed step and specialized step handler
        propagator.setMasterMode(60., new TutorialStepHandler());
        // Extrapolate from the initial to the final date
        SpacecraftState finalState = propagator.propagate(initialDate.shiftedBy(630.));
        KeplerianOrbit o = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(finalState.getOrbit());
        System.out.format(Locale.US, "Final state:%n%s %12.3f %10.8f %10.6f %10.6f %10.6f %10.6f%n", finalState.getDate(), o.getA(), o.getE(), FastMath.toDegrees(o.getI()), FastMath.toDegrees(o.getPerigeeArgument()), FastMath.toDegrees(o.getRightAscensionOfAscendingNode()), FastMath.toDegrees(o.getTrueAnomaly()));
    } catch (OrekitException oe) {
        System.err.println(oe.getMessage());
    }
}