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

Example 16 with UnnormalizedSphericalHarmonicsProvider

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