Examples with NormalizedSphericalHarmonicsProvider org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider used on opensource projects

Example 31 with NormalizedSphericalHarmonicsProvider

use of org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider 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 32 with NormalizedSphericalHarmonicsProvider

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

the class DSSTPropagation method run.

private void run(final File input, final File output) throws IOException, IllegalArgumentException, OrekitException, ParseException {
    // read input parameters
    KeyValueFileParser<ParameterKey> parser = new KeyValueFileParser<ParameterKey>(ParameterKey.class);
    try (final FileInputStream fis = new FileInputStream(input)) {
        parser.parseInput(input.getAbsolutePath(), fis);
    }
    // check mandatory input parameters
    if (!parser.containsKey(ParameterKey.ORBIT_DATE)) {
        throw new IOException("Orbit date is not defined.");
    }
    if (!parser.containsKey(ParameterKey.DURATION) && !parser.containsKey(ParameterKey.DURATION_IN_DAYS)) {
        throw new IOException("Propagation duration is not defined.");
    }
    // All dates in UTC
    final TimeScale utc = TimeScalesFactory.getUTC();
    final double rotationRate;
    if (!parser.containsKey(ParameterKey.CENTRAL_BODY_ROTATION_RATE)) {
        rotationRate = Constants.WGS84_EARTH_ANGULAR_VELOCITY;
    } else {
        rotationRate = parser.getDouble(ParameterKey.CENTRAL_BODY_ROTATION_RATE);
    }
    final int degree = parser.getInt(ParameterKey.CENTRAL_BODY_DEGREE);
    final int order = FastMath.min(degree, parser.getInt(ParameterKey.CENTRAL_BODY_ORDER));
    // Potential coefficients providers
    final UnnormalizedSphericalHarmonicsProvider unnormalized = GravityFieldFactory.getConstantUnnormalizedProvider(degree, order);
    final NormalizedSphericalHarmonicsProvider normalized = GravityFieldFactory.getConstantNormalizedProvider(degree, order);
    // Central body attraction coefficient (m³/s²)
    final double mu = unnormalized.getMu();
    // Earth frame definition
    final Frame earthFrame;
    if (!parser.containsKey(ParameterKey.BODY_FRAME)) {
        earthFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
    } else {
        earthFrame = parser.getEarthFrame(ParameterKey.BODY_FRAME);
    }
    // Orbit definition
    final Orbit orbit = createOrbit(parser, utc, mu);
    // DSST propagator definition
    double mass = 1000.0;
    if (parser.containsKey(ParameterKey.MASS)) {
        mass = parser.getDouble(ParameterKey.MASS);
    }
    boolean initialIsOsculating = false;
    if (parser.containsKey(ParameterKey.INITIAL_ORBIT_IS_OSCULATING)) {
        initialIsOsculating = parser.getBoolean(ParameterKey.INITIAL_ORBIT_IS_OSCULATING);
    }
    boolean outputIsOsculating = initialIsOsculating;
    if (parser.containsKey(ParameterKey.OUTPUT_ORBIT_IS_OSCULATING)) {
        outputIsOsculating = parser.getBoolean(ParameterKey.OUTPUT_ORBIT_IS_OSCULATING);
    }
    List<String> shortPeriodCoefficients = null;
    if (parser.containsKey(ParameterKey.OUTPUT_SHORT_PERIOD_COEFFICIENTS)) {
        shortPeriodCoefficients = parser.getStringsList(ParameterKey.OUTPUT_SHORT_PERIOD_COEFFICIENTS, ',');
        if (shortPeriodCoefficients.size() == 1 && shortPeriodCoefficients.get(0).equalsIgnoreCase("all")) {
            // special case, we use the empty list to represent all possible (unknown) keys
            // we don't use Collections.emptyList() because we want the list to be populated later on
            shortPeriodCoefficients = new ArrayList<String>();
        } else if (shortPeriodCoefficients.size() == 1 && shortPeriodCoefficients.get(0).equalsIgnoreCase("none")) {
            // special case, we use null to select no coefficients at all
            shortPeriodCoefficients = null;
        } else {
            // general case, we have an explicit list of coefficients names
            Collections.sort(shortPeriodCoefficients);
        }
        if (shortPeriodCoefficients != null && !outputIsOsculating) {
            System.out.println("\nWARNING:");
            System.out.println("Short periodic coefficients can be output only if output orbit is osculating.");
            System.out.println("No coefficients will be computed here.\n");
        }
    }
    double fixedStepSize = -1.;
    double minStep = 6000.0;
    double maxStep = 86400.0;
    double dP = 1.0;
    if (parser.containsKey(ParameterKey.FIXED_INTEGRATION_STEP)) {
        fixedStepSize = parser.getDouble(ParameterKey.FIXED_INTEGRATION_STEP);
    } else {
        if (parser.containsKey(ParameterKey.MIN_VARIABLE_INTEGRATION_STEP)) {
            minStep = parser.getDouble(ParameterKey.MIN_VARIABLE_INTEGRATION_STEP);
        }
        if (parser.containsKey(ParameterKey.MAX_VARIABLE_INTEGRATION_STEP)) {
            maxStep = parser.getDouble(ParameterKey.MAX_VARIABLE_INTEGRATION_STEP);
        }
        if (parser.containsKey(ParameterKey.POSITION_TOLERANCE_VARIABLE_INTEGRATION_STEP)) {
            dP = parser.getDouble(ParameterKey.POSITION_TOLERANCE_VARIABLE_INTEGRATION_STEP);
        }
    }
    final DSSTPropagator dsstProp = createDSSTProp(orbit, mass, initialIsOsculating, outputIsOsculating, fixedStepSize, minStep, maxStep, dP, shortPeriodCoefficients);
    if (parser.containsKey(ParameterKey.FIXED_NUMBER_OF_INTERPOLATION_POINTS)) {
        if (parser.containsKey(ParameterKey.MAX_TIME_GAP_BETWEEN_INTERPOLATION_POINTS)) {
            throw new OrekitException(LocalizedCoreFormats.SIMPLE_MESSAGE, "cannot specify both fixed.number.of.interpolation.points" + " and max.time.gap.between.interpolation.points");
        }
        dsstProp.setInterpolationGridToFixedNumberOfPoints(parser.getInt(ParameterKey.FIXED_NUMBER_OF_INTERPOLATION_POINTS));
    } else if (parser.containsKey(ParameterKey.MAX_TIME_GAP_BETWEEN_INTERPOLATION_POINTS)) {
        dsstProp.setInterpolationGridToMaxTimeGap(parser.getDouble(ParameterKey.MAX_TIME_GAP_BETWEEN_INTERPOLATION_POINTS));
    } else {
        dsstProp.setInterpolationGridToFixedNumberOfPoints(3);
    }
    // Set Force models
    setForceModel(parser, unnormalized, earthFrame, rotationRate, dsstProp);
    // Simulation properties
    AbsoluteDate start;
    if (parser.containsKey(ParameterKey.START_DATE)) {
        start = parser.getDate(ParameterKey.START_DATE, utc);
    } else {
        start = parser.getDate(ParameterKey.ORBIT_DATE, utc);
    }
    double duration = 0.;
    if (parser.containsKey(ParameterKey.DURATION)) {
        duration = parser.getDouble(ParameterKey.DURATION);
    }
    if (parser.containsKey(ParameterKey.DURATION_IN_DAYS)) {
        duration = parser.getDouble(ParameterKey.DURATION_IN_DAYS) * Constants.JULIAN_DAY;
    }
    double outStep = parser.getDouble(ParameterKey.OUTPUT_STEP);
    boolean displayKeplerian = true;
    if (parser.containsKey(ParameterKey.OUTPUT_KEPLERIAN)) {
        displayKeplerian = parser.getBoolean(ParameterKey.OUTPUT_KEPLERIAN);
    }
    boolean displayEquinoctial = true;
    if (parser.containsKey(ParameterKey.OUTPUT_EQUINOCTIAL)) {
        displayEquinoctial = parser.getBoolean(ParameterKey.OUTPUT_EQUINOCTIAL);
    }
    boolean displayCartesian = true;
    if (parser.containsKey(ParameterKey.OUTPUT_CARTESIAN)) {
        displayCartesian = parser.getBoolean(ParameterKey.OUTPUT_CARTESIAN);
    }
    // DSST Propagation
    dsstProp.setEphemerisMode();
    final double dsstOn = System.currentTimeMillis();
    dsstProp.propagate(start, start.shiftedBy(duration));
    final double dsstOff = System.currentTimeMillis();
    System.out.println("DSST execution time (without large file write) : " + (dsstOff - dsstOn) / 1000.);
    System.out.println("writing file...");
    final BoundedPropagator dsstEphem = dsstProp.getGeneratedEphemeris();
    dsstEphem.setMasterMode(outStep, new OutputHandler(output, displayKeplerian, displayEquinoctial, displayCartesian, shortPeriodCoefficients));
    dsstEphem.propagate(start, start.shiftedBy(duration));
    System.out.println("DSST results saved as file " + output);
    // Check if we want to compare numerical to DSST propagator (default is false)
    if (parser.containsKey(ParameterKey.NUMERICAL_COMPARISON) && parser.getBoolean(ParameterKey.NUMERICAL_COMPARISON)) {
        if (!outputIsOsculating) {
            System.out.println("\nWARNING:");
            System.out.println("The DSST propagator considers a mean orbit while the numerical will consider an osculating one.");
            System.out.println("The comparison will be meaningless.\n");
        }
        // Numerical propagator definition
        final NumericalPropagator numProp = createNumProp(orbit, mass);
        // Set Force models
        setForceModel(parser, normalized, earthFrame, numProp);
        // Numerical Propagation without output
        numProp.setEphemerisMode();
        final double numOn = System.currentTimeMillis();
        numProp.propagate(start, start.shiftedBy(duration));
        final double numOff = System.currentTimeMillis();
        System.out.println("Numerical execution time (including output): " + (numOff - numOn) / 1000.);
        // Add output
        final BoundedPropagator numEphemeris = numProp.getGeneratedEphemeris();
        File numOutput = new File(input.getParentFile(), "numerical-propagation.out");
        numEphemeris.setMasterMode(outStep, new OutputHandler(numOutput, displayKeplerian, displayEquinoctial, displayCartesian, null));
        System.out.println("Writing file, this may take some time ...");
        numEphemeris.propagate(numEphemeris.getMaxDate());
        System.out.println("Numerical results saved as file " + numOutput);
    }
}

Example 33 with NormalizedSphericalHarmonicsProvider

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

Example 34 with NormalizedSphericalHarmonicsProvider

use of org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider in project SpriteOrbits by ProjectPersephone.

the class SpritePropOrig method createPropagator.

/**
 * Create a numerical propagator for a state.
 * @param state state to propagate
 * @param attitudeProvider provider for the attitude
 * @param crossSection cross section of the object
 * @param dragCoeff drag coefficient
 */
private Propagator createPropagator(final SpacecraftState state, final AttitudeProvider attitudeProvider, final double crossSection, final double dragCoeff) throws OrekitException {
    // see https://www.orekit.org/static/architecture/propagation.html
    // steps limits
    final double minStep = 0.001;
    final double maxStep = 1000;
    final double initStep = 60;
    // error control parameters (absolute and relative)
    final double positionError = 10.0;
    // we will propagate in Cartesian coordinates
    final OrbitType orbitType = OrbitType.CARTESIAN;
    final double[][] tolerances = NumericalPropagator.tolerances(positionError, state.getOrbit(), orbitType);
    // set up mathematical integrator
    AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(minStep, maxStep, tolerances[0], tolerances[1]);
    integrator.setInitialStepSize(initStep);
    // set up space dynamics propagator
    NumericalPropagator propagator = new NumericalPropagator(integrator);
    propagator.setOrbitType(orbitType);
    // add gravity field force model
    final NormalizedSphericalHarmonicsProvider gravityProvider = GravityFieldFactory.getNormalizedProvider(8, 8);
    propagator.addForceModel(new HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravityProvider));
    // add atmospheric drag force model
    propagator.addForceModel(new DragForce(new HarrisPriester(sun, earth), new SphericalSpacecraft(crossSection, dragCoeff, 0.0, 0.0)));
    // set attitude mode
    propagator.setAttitudeProvider(attitudeProvider);
    propagator.setInitialState(state);
    return propagator;
}