Search in sources :

Example 41 with Orbit

use of org.orekit.orbits.Orbit in project Orekit by CS-SI.

the class OrbitDetermination method run.

private void run(final File input, final File home) 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);
    }
    // log file
    final String baseName;
    final PrintStream logStream;
    if (parser.containsKey(ParameterKey.OUTPUT_BASE_NAME) && parser.getString(ParameterKey.OUTPUT_BASE_NAME).length() > 0) {
        baseName = parser.getString(ParameterKey.OUTPUT_BASE_NAME);
        logStream = new PrintStream(new File(home, baseName + "-log.out"), "UTF-8");
    } else {
        baseName = null;
        logStream = null;
    }
    final RangeLog rangeLog = new RangeLog(home, baseName);
    final RangeRateLog rangeRateLog = new RangeRateLog(home, baseName);
    final AzimuthLog azimuthLog = new AzimuthLog(home, baseName);
    final ElevationLog elevationLog = new ElevationLog(home, baseName);
    final PositionLog positionLog = new PositionLog(home, baseName);
    final VelocityLog velocityLog = new VelocityLog(home, baseName);
    try {
        // gravity field
        final NormalizedSphericalHarmonicsProvider gravityField = createGravityField(parser);
        // Orbit initial guess
        final Orbit initialGuess = createOrbit(parser, gravityField.getMu());
        // IERS conventions
        final IERSConventions conventions;
        if (!parser.containsKey(ParameterKey.IERS_CONVENTIONS)) {
            conventions = IERSConventions.IERS_2010;
        } else {
            conventions = IERSConventions.valueOf("IERS_" + parser.getInt(ParameterKey.IERS_CONVENTIONS));
        }
        // central body
        final OneAxisEllipsoid body = createBody(parser);
        // propagator builder
        final NumericalPropagatorBuilder propagatorBuilder = createPropagatorBuilder(parser, conventions, gravityField, body, initialGuess);
        // estimator
        final BatchLSEstimator estimator = createEstimator(parser, propagatorBuilder);
        // measurements
        final List<ObservedMeasurement<?>> measurements = new ArrayList<ObservedMeasurement<?>>();
        for (final String fileName : parser.getStringsList(ParameterKey.MEASUREMENTS_FILES, ',')) {
            measurements.addAll(readMeasurements(new File(input.getParentFile(), fileName), createStationsData(parser, body), createPVData(parser), createSatRangeBias(parser), createWeights(parser), createRangeOutliersManager(parser), createRangeRateOutliersManager(parser), createAzElOutliersManager(parser), createPVOutliersManager(parser)));
        }
        for (ObservedMeasurement<?> measurement : measurements) {
            estimator.addMeasurement(measurement);
        }
        // estimate orbit
        estimator.setObserver(new BatchLSObserver() {

            private PVCoordinates previousPV;

            {
                previousPV = initialGuess.getPVCoordinates();
                final String header = "iteration evaluations      ΔP(m)        ΔV(m/s)           RMS          nb Range    nb Range-rate  nb Angular     nb PV%n";
                System.out.format(Locale.US, header);
                if (logStream != null) {
                    logStream.format(Locale.US, header);
                }
            }

            /**
             * {@inheritDoc}
             */
            @Override
            public void evaluationPerformed(final int iterationsCount, final int evaluationsCount, final Orbit[] orbits, final ParameterDriversList estimatedOrbitalParameters, final ParameterDriversList estimatedPropagatorParameters, final ParameterDriversList estimatedMeasurementsParameters, final EstimationsProvider evaluationsProvider, final LeastSquaresProblem.Evaluation lspEvaluation) {
                PVCoordinates currentPV = orbits[0].getPVCoordinates();
                final String format0 = "    %2d         %2d                                 %16.12f     %s       %s     %s     %s%n";
                final String format = "    %2d         %2d      %13.6f %12.9f %16.12f     %s       %s     %s     %s%n";
                final EvaluationCounter<Range> rangeCounter = new EvaluationCounter<Range>();
                final EvaluationCounter<RangeRate> rangeRateCounter = new EvaluationCounter<RangeRate>();
                final EvaluationCounter<AngularAzEl> angularCounter = new EvaluationCounter<AngularAzEl>();
                final EvaluationCounter<PV> pvCounter = new EvaluationCounter<PV>();
                for (final Map.Entry<ObservedMeasurement<?>, EstimatedMeasurement<?>> entry : estimator.getLastEstimations().entrySet()) {
                    if (entry.getKey() instanceof Range) {
                        @SuppressWarnings("unchecked") EstimatedMeasurement<Range> evaluation = (EstimatedMeasurement<Range>) entry.getValue();
                        rangeCounter.add(evaluation);
                    } else if (entry.getKey() instanceof RangeRate) {
                        @SuppressWarnings("unchecked") EstimatedMeasurement<RangeRate> evaluation = (EstimatedMeasurement<RangeRate>) entry.getValue();
                        rangeRateCounter.add(evaluation);
                    } else if (entry.getKey() instanceof AngularAzEl) {
                        @SuppressWarnings("unchecked") EstimatedMeasurement<AngularAzEl> evaluation = (EstimatedMeasurement<AngularAzEl>) entry.getValue();
                        angularCounter.add(evaluation);
                    } else if (entry.getKey() instanceof PV) {
                        @SuppressWarnings("unchecked") EstimatedMeasurement<PV> evaluation = (EstimatedMeasurement<PV>) entry.getValue();
                        pvCounter.add(evaluation);
                    }
                }
                if (evaluationsCount == 1) {
                    System.out.format(Locale.US, format0, iterationsCount, evaluationsCount, lspEvaluation.getRMS(), rangeCounter.format(8), rangeRateCounter.format(8), angularCounter.format(8), pvCounter.format(8));
                    if (logStream != null) {
                        logStream.format(Locale.US, format0, iterationsCount, evaluationsCount, lspEvaluation.getRMS(), rangeCounter.format(8), rangeRateCounter.format(8), angularCounter.format(8), pvCounter.format(8));
                    }
                } else {
                    System.out.format(Locale.US, format, iterationsCount, evaluationsCount, Vector3D.distance(previousPV.getPosition(), currentPV.getPosition()), Vector3D.distance(previousPV.getVelocity(), currentPV.getVelocity()), lspEvaluation.getRMS(), rangeCounter.format(8), rangeRateCounter.format(8), angularCounter.format(8), pvCounter.format(8));
                    if (logStream != null) {
                        logStream.format(Locale.US, format, iterationsCount, evaluationsCount, Vector3D.distance(previousPV.getPosition(), currentPV.getPosition()), Vector3D.distance(previousPV.getVelocity(), currentPV.getVelocity()), lspEvaluation.getRMS(), rangeCounter.format(8), rangeRateCounter.format(8), angularCounter.format(8), pvCounter.format(8));
                    }
                }
                previousPV = currentPV;
            }
        });
        Orbit estimated = estimator.estimate()[0].getInitialState().getOrbit();
        // compute some statistics
        for (final Map.Entry<ObservedMeasurement<?>, EstimatedMeasurement<?>> entry : estimator.getLastEstimations().entrySet()) {
            if (entry.getKey() instanceof Range) {
                @SuppressWarnings("unchecked") EstimatedMeasurement<Range> evaluation = (EstimatedMeasurement<Range>) entry.getValue();
                rangeLog.add(evaluation);
            } else if (entry.getKey() instanceof RangeRate) {
                @SuppressWarnings("unchecked") EstimatedMeasurement<RangeRate> evaluation = (EstimatedMeasurement<RangeRate>) entry.getValue();
                rangeRateLog.add(evaluation);
            } else if (entry.getKey() instanceof AngularAzEl) {
                @SuppressWarnings("unchecked") EstimatedMeasurement<AngularAzEl> evaluation = (EstimatedMeasurement<AngularAzEl>) entry.getValue();
                azimuthLog.add(evaluation);
                elevationLog.add(evaluation);
            } else if (entry.getKey() instanceof PV) {
                @SuppressWarnings("unchecked") EstimatedMeasurement<PV> evaluation = (EstimatedMeasurement<PV>) entry.getValue();
                positionLog.add(evaluation);
                velocityLog.add(evaluation);
            }
        }
        System.out.println("Estimated orbit: " + estimated);
        if (logStream != null) {
            logStream.println("Estimated orbit: " + estimated);
        }
        final ParameterDriversList orbitalParameters = estimator.getOrbitalParametersDrivers(true);
        final ParameterDriversList propagatorParameters = estimator.getPropagatorParametersDrivers(true);
        final ParameterDriversList measurementsParameters = estimator.getMeasurementsParametersDrivers(true);
        int length = 0;
        for (final ParameterDriver parameterDriver : orbitalParameters.getDrivers()) {
            length = FastMath.max(length, parameterDriver.getName().length());
        }
        for (final ParameterDriver parameterDriver : propagatorParameters.getDrivers()) {
            length = FastMath.max(length, parameterDriver.getName().length());
        }
        for (final ParameterDriver parameterDriver : measurementsParameters.getDrivers()) {
            length = FastMath.max(length, parameterDriver.getName().length());
        }
        displayParametersChanges(System.out, "Estimated orbital parameters changes: ", false, length, orbitalParameters);
        if (logStream != null) {
            displayParametersChanges(logStream, "Estimated orbital parameters changes: ", false, length, orbitalParameters);
        }
        displayParametersChanges(System.out, "Estimated propagator parameters changes: ", true, length, propagatorParameters);
        if (logStream != null) {
            displayParametersChanges(logStream, "Estimated propagator parameters changes: ", true, length, propagatorParameters);
        }
        displayParametersChanges(System.out, "Estimated measurements parameters changes: ", true, length, measurementsParameters);
        if (logStream != null) {
            displayParametersChanges(logStream, "Estimated measurements parameters changes: ", true, length, measurementsParameters);
        }
        System.out.println("Number of iterations: " + estimator.getIterationsCount());
        System.out.println("Number of evaluations: " + estimator.getEvaluationsCount());
        rangeLog.displaySummary(System.out);
        rangeRateLog.displaySummary(System.out);
        azimuthLog.displaySummary(System.out);
        elevationLog.displaySummary(System.out);
        positionLog.displaySummary(System.out);
        velocityLog.displaySummary(System.out);
        if (logStream != null) {
            logStream.println("Number of iterations: " + estimator.getIterationsCount());
            logStream.println("Number of evaluations: " + estimator.getEvaluationsCount());
            rangeLog.displaySummary(logStream);
            rangeRateLog.displaySummary(logStream);
            azimuthLog.displaySummary(logStream);
            elevationLog.displaySummary(logStream);
            positionLog.displaySummary(logStream);
            velocityLog.displaySummary(logStream);
        }
        rangeLog.displayResiduals();
        rangeRateLog.displayResiduals();
        azimuthLog.displayResiduals();
        elevationLog.displayResiduals();
        positionLog.displayResiduals();
        velocityLog.displayResiduals();
    } finally {
        if (logStream != null) {
            logStream.close();
        }
        rangeLog.close();
        rangeRateLog.close();
        azimuthLog.close();
        elevationLog.close();
        positionLog.close();
        velocityLog.close();
    }
}
Also used : OneAxisEllipsoid(org.orekit.bodies.OneAxisEllipsoid) PV(org.orekit.estimation.measurements.PV) ArrayList(java.util.ArrayList) PVCoordinates(org.orekit.utils.PVCoordinates) ObservedMeasurement(org.orekit.estimation.measurements.ObservedMeasurement) EstimatedMeasurement(org.orekit.estimation.measurements.EstimatedMeasurement) BatchLSObserver(org.orekit.estimation.leastsquares.BatchLSObserver) EquinoctialOrbit(org.orekit.orbits.EquinoctialOrbit) CartesianOrbit(org.orekit.orbits.CartesianOrbit) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) Orbit(org.orekit.orbits.Orbit) CircularOrbit(org.orekit.orbits.CircularOrbit) Range(org.orekit.estimation.measurements.Range) FileInputStream(java.io.FileInputStream) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) RangeRate(org.orekit.estimation.measurements.RangeRate) File(java.io.File) AngularAzEl(org.orekit.estimation.measurements.AngularAzEl) Map(java.util.Map) HashMap(java.util.HashMap) BatchLSEstimator(org.orekit.estimation.leastsquares.BatchLSEstimator) ParameterDriversList(org.orekit.utils.ParameterDriversList) LeastSquaresProblem(org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem) NormalizedSphericalHarmonicsProvider(org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider) EstimationsProvider(org.orekit.estimation.measurements.EstimationsProvider) PrintStream(java.io.PrintStream) KeyValueFileParser(fr.cs.examples.KeyValueFileParser) IERSConventions(org.orekit.utils.IERSConventions) ParameterDriver(org.orekit.utils.ParameterDriver) GeodeticPoint(org.orekit.bodies.GeodeticPoint)

Example 42 with Orbit

use of org.orekit.orbits.Orbit in project Orekit by CS-SI.

the class TrackCorridor method createPropagator.

/**
 * Create an orbit propagator for a circular orbit
 * @param a  semi-major axis (m)
 * @param ex e cos(ω), first component of circular eccentricity vector
 * @param ey e sin(ω), second component of circular eccentricity vector
 * @param i inclination (rad)
 * @param raan right ascension of ascending node (Ω, rad)
 * @param alpha  an + ω, mean latitude argument (rad)
 * @param date date of the orbital parameters
 * @return an orbit propagator
 * @exception OrekitException if propagator cannot be built
 */
private Propagator createPropagator(final AbsoluteDate date, final double a, final double ex, final double ey, final double i, final double raan, final double alpha) throws OrekitException {
    // create orbit
    Orbit initialOrbit = new CircularOrbit(a, ex, ey, i, raan, alpha, PositionAngle.MEAN, FramesFactory.getEME2000(), date, Constants.EIGEN5C_EARTH_MU);
    // create propagator
    Propagator propagator = new EcksteinHechlerPropagator(initialOrbit, new LofOffset(initialOrbit.getFrame(), LOFType.TNW), Constants.EIGEN5C_EARTH_EQUATORIAL_RADIUS, Constants.EIGEN5C_EARTH_MU, Constants.EIGEN5C_EARTH_C20, Constants.EIGEN5C_EARTH_C30, Constants.EIGEN5C_EARTH_C40, Constants.EIGEN5C_EARTH_C50, Constants.EIGEN5C_EARTH_C60);
    return propagator;
}
Also used : EcksteinHechlerPropagator(org.orekit.propagation.analytical.EcksteinHechlerPropagator) Orbit(org.orekit.orbits.Orbit) CircularOrbit(org.orekit.orbits.CircularOrbit) CircularOrbit(org.orekit.orbits.CircularOrbit) TLEPropagator(org.orekit.propagation.analytical.tle.TLEPropagator) EcksteinHechlerPropagator(org.orekit.propagation.analytical.EcksteinHechlerPropagator) Propagator(org.orekit.propagation.Propagator) LofOffset(org.orekit.attitudes.LofOffset)

Example 43 with Orbit

use of org.orekit.orbits.Orbit in project Orekit by CS-SI.

the class SpacecraftState method interpolate.

/**
 * {@inheritDoc}
 * <p>
 * The additional states that are interpolated are the ones already present
 * in the instance. The sample instances must therefore have at least the same
 * additional states has the instance. They may have more additional states,
 * but the extra ones will be ignored.
 * </p>
 * <p>
 * As this implementation of interpolation is polynomial, it should be used only
 * with small samples (about 10-20 points) in order to avoid <a
 * href="http://en.wikipedia.org/wiki/Runge%27s_phenomenon">Runge's phenomenon</a>
 * and numerical problems (including NaN appearing).
 * </p>
 */
public SpacecraftState interpolate(final AbsoluteDate date, final Stream<SpacecraftState> sample) throws OrekitException {
    // prepare interpolators
    final List<Orbit> orbits = new ArrayList<>();
    final List<Attitude> attitudes = new ArrayList<>();
    final HermiteInterpolator massInterpolator = new HermiteInterpolator();
    final Map<String, HermiteInterpolator> additionalInterpolators = new HashMap<String, HermiteInterpolator>(additional.size());
    for (final String name : additional.keySet()) {
        additionalInterpolators.put(name, new HermiteInterpolator());
    }
    // extract sample data
    try {
        sample.forEach(state -> {
            try {
                final double deltaT = state.getDate().durationFrom(date);
                orbits.add(state.getOrbit());
                attitudes.add(state.getAttitude());
                massInterpolator.addSamplePoint(deltaT, new double[] { state.getMass() });
                for (final Map.Entry<String, HermiteInterpolator> entry : additionalInterpolators.entrySet()) {
                    entry.getValue().addSamplePoint(deltaT, state.getAdditionalState(entry.getKey()));
                }
            } catch (OrekitException oe) {
                throw new OrekitExceptionWrapper(oe);
            }
        });
    } catch (OrekitExceptionWrapper oew) {
        throw oew.getException();
    }
    // perform interpolations
    final Orbit interpolatedOrbit = orbit.interpolate(date, orbits);
    final Attitude interpolatedAttitude = attitude.interpolate(date, attitudes);
    final double interpolatedMass = massInterpolator.value(0)[0];
    final Map<String, double[]> interpolatedAdditional;
    if (additional.isEmpty()) {
        interpolatedAdditional = null;
    } else {
        interpolatedAdditional = new HashMap<String, double[]>(additional.size());
        for (final Map.Entry<String, HermiteInterpolator> entry : additionalInterpolators.entrySet()) {
            interpolatedAdditional.put(entry.getKey(), entry.getValue().value(0));
        }
    }
    // create the complete interpolated state
    return new SpacecraftState(interpolatedOrbit, interpolatedAttitude, interpolatedMass, interpolatedAdditional);
}
Also used : OrekitExceptionWrapper(org.orekit.errors.OrekitExceptionWrapper) Orbit(org.orekit.orbits.Orbit) Attitude(org.orekit.attitudes.Attitude) HashMap(java.util.HashMap) ArrayList(java.util.ArrayList) HermiteInterpolator(org.hipparchus.analysis.interpolation.HermiteInterpolator) OrekitException(org.orekit.errors.OrekitException) HashMap(java.util.HashMap) Map(java.util.Map)

Example 44 with Orbit

use of org.orekit.orbits.Orbit in project Orekit by CS-SI.

the class KeplerianPropagator method fixState.

/**
 * Fix state to use a specified mu and remove derivatives.
 * <p>
 * This ensures the propagation model (which is based on calling
 * {@link Orbit#shiftedBy(double)}) is Keplerian only and uses a specified mu.
 * </p>
 * @param orbit orbit to fix
 * @param attitude current attitude
 * @param mass current mass
 * @param mu gravity coefficient to use
 * @param additionalStates additional states
 * @return fixed orbit
 */
private SpacecraftState fixState(final Orbit orbit, final Attitude attitude, final double mass, final double mu, final Map<String, double[]> additionalStates) {
    final OrbitType type = orbit.getType();
    final double[] stateVector = new double[6];
    type.mapOrbitToArray(orbit, PositionAngle.TRUE, stateVector, null);
    final Orbit fixedOrbit = type.mapArrayToOrbit(stateVector, null, PositionAngle.TRUE, orbit.getDate(), mu, orbit.getFrame());
    SpacecraftState fixedState = new SpacecraftState(fixedOrbit, attitude, mass);
    for (final Map.Entry<String, double[]> entry : additionalStates.entrySet()) {
        fixedState = fixedState.addAdditionalState(entry.getKey(), entry.getValue());
    }
    return fixedState;
}
Also used : SpacecraftState(org.orekit.propagation.SpacecraftState) Orbit(org.orekit.orbits.Orbit) OrbitType(org.orekit.orbits.OrbitType) TimeSpanMap(org.orekit.utils.TimeSpanMap) Map(java.util.Map)

Example 45 with Orbit

use of org.orekit.orbits.Orbit in project Orekit by CS-SI.

the class TLEPropagatorBuilder method buildPropagator.

/**
 * {@inheritDoc}
 */
public Propagator buildPropagator(final double[] normalizedParameters) throws OrekitException {
    // create the orbit
    setParameters(normalizedParameters);
    final Orbit orbit = createInitialOrbit();
    // we really need a Keplerian orbit type
    final KeplerianOrbit kep = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(orbit);
    final TLE tle = new TLE(satelliteNumber, classification, launchYear, launchNumber, launchPiece, TLE.DEFAULT, elementNumber, orbit.getDate(), kep.getKeplerianMeanMotion(), 0.0, 0.0, kep.getE(), MathUtils.normalizeAngle(orbit.getI(), FastMath.PI), MathUtils.normalizeAngle(kep.getPerigeeArgument(), FastMath.PI), MathUtils.normalizeAngle(kep.getRightAscensionOfAscendingNode(), FastMath.PI), MathUtils.normalizeAngle(kep.getMeanAnomaly(), FastMath.PI), revolutionNumberAtEpoch, bStar);
    return TLEPropagator.selectExtrapolator(tle);
}
Also used : KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) Orbit(org.orekit.orbits.Orbit) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) TLE(org.orekit.propagation.analytical.tle.TLE)

Aggregations

Orbit (org.orekit.orbits.Orbit)211 KeplerianOrbit (org.orekit.orbits.KeplerianOrbit)161 Test (org.junit.Test)153 AbsoluteDate (org.orekit.time.AbsoluteDate)153 SpacecraftState (org.orekit.propagation.SpacecraftState)129 Vector3D (org.hipparchus.geometry.euclidean.threed.Vector3D)99 EquinoctialOrbit (org.orekit.orbits.EquinoctialOrbit)94 CartesianOrbit (org.orekit.orbits.CartesianOrbit)88 FieldAbsoluteDate (org.orekit.time.FieldAbsoluteDate)74 CircularOrbit (org.orekit.orbits.CircularOrbit)68 PVCoordinates (org.orekit.utils.PVCoordinates)66 Frame (org.orekit.frames.Frame)51 NumericalPropagator (org.orekit.propagation.numerical.NumericalPropagator)51 DateComponents (org.orekit.time.DateComponents)48 FieldSpacecraftState (org.orekit.propagation.FieldSpacecraftState)46 Propagator (org.orekit.propagation.Propagator)46 TimeComponents (org.orekit.time.TimeComponents)44 OneAxisEllipsoid (org.orekit.bodies.OneAxisEllipsoid)43 AbstractLegacyForceModelTest (org.orekit.forces.AbstractLegacyForceModelTest)41 FieldKeplerianOrbit (org.orekit.orbits.FieldKeplerianOrbit)39