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

Example 1 with ICGEMFormatReader

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

the class HolmesFeatherstoneAttractionModelTest method testTimeDependentField.

@Test
public void testTimeDependentField() throws OrekitException {
    Utils.setDataRoot("regular-data:potential/icgem-format");
    GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-6s-truncated", true));
    final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
    final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
    final SpacecraftState spacecraftState = new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel), FramesFactory.getGCRF(), new AbsoluteDate(2005, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()), GravityFieldFactory.getUnnormalizedProvider(1, 1).getMu()));
    double dP = 0.1;
    double duration = 3 * Constants.JULIAN_DAY;
    BoundedPropagator fixedFieldEphemeris = createEphemeris(dP, spacecraftState, duration, GravityFieldFactory.getConstantNormalizedProvider(8, 8));
    BoundedPropagator varyingFieldEphemeris = createEphemeris(dP, spacecraftState, duration, GravityFieldFactory.getNormalizedProvider(8, 8));
    double step = 60.0;
    double maxDeltaT = 0;
    double maxDeltaN = 0;
    double maxDeltaW = 0;
    for (AbsoluteDate date = fixedFieldEphemeris.getMinDate(); date.compareTo(fixedFieldEphemeris.getMaxDate()) < 0; date = date.shiftedBy(step)) {
        PVCoordinates pvFixedField = fixedFieldEphemeris.getPVCoordinates(date, FramesFactory.getGCRF());
        PVCoordinates pvVaryingField = varyingFieldEphemeris.getPVCoordinates(date, FramesFactory.getGCRF());
        Vector3D t = pvFixedField.getVelocity().normalize();
        Vector3D w = pvFixedField.getMomentum().normalize();
        Vector3D n = Vector3D.crossProduct(w, t);
        Vector3D delta = pvVaryingField.getPosition().subtract(pvFixedField.getPosition());
        maxDeltaT = FastMath.max(maxDeltaT, FastMath.abs(Vector3D.dotProduct(delta, t)));
        maxDeltaN = FastMath.max(maxDeltaN, FastMath.abs(Vector3D.dotProduct(delta, n)));
        maxDeltaW = FastMath.max(maxDeltaW, FastMath.abs(Vector3D.dotProduct(delta, w)));
    }
    Assert.assertTrue(maxDeltaT > 0.15);
    Assert.assertTrue(maxDeltaT < 0.25);
    Assert.assertTrue(maxDeltaN > 0.01);
    Assert.assertTrue(maxDeltaN < 0.02);
    Assert.assertTrue(maxDeltaW > 0.05);
    Assert.assertTrue(maxDeltaW < 0.10);
}

Example 2 with ICGEMFormatReader

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

the class OrbitDeterminationTest method run.

private ResultOD run(final File input, final boolean print) throws IOException, IllegalArgumentException, OrekitException, ParseException {
    // read input parameters
    KeyValueFileParser<ParameterKey> parser = new KeyValueFileParser<ParameterKey>(ParameterKey.class);
    parser.parseInput(input.getAbsolutePath(), new FileInputStream(input));
    // log file
    final RangeLog rangeLog = new RangeLog();
    final RangeRateLog rangeRateLog = new RangeRateLog();
    final AzimuthLog azimuthLog = new AzimuthLog();
    final ElevationLog elevationLog = new ElevationLog();
    final PositionLog positionLog = new PositionLog();
    final VelocityLog velocityLog = new VelocityLog();
    // gravity field
    GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-5c.gfc", true));
    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);
    }
    if (print) {
        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);
            }

            /**
             * {@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));
                } 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));
                }
                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);
        }
    }
    // parmaters and measurements.
    final ParameterDriversList propagatorParameters = estimator.getPropagatorParametersDrivers(true);
    final ParameterDriversList measurementsParameters = estimator.getMeasurementsParametersDrivers(true);
    // instation of results
    return new ResultOD(propagatorParameters, measurementsParameters, estimator.getIterationsCount(), estimator.getEvaluationsCount(), estimated.getPVCoordinates(), rangeLog.createStatisticsSummary(), rangeRateLog.createStatisticsSummary(), azimuthLog.createStatisticsSummary(), elevationLog.createStatisticsSummary(), positionLog.createStatisticsSummary(), velocityLog.createStatisticsSummary(), estimator.getPhysicalCovariances(1.0e-10));
}

Example 3 with ICGEMFormatReader

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

the class KalmanOrbitDeterminationTest method testLageos2.

@Test
public // Orbit determination for Lageos2 based on SLR (range) measurements
void testLageos2() throws URISyntaxException, IllegalArgumentException, IOException, OrekitException, ParseException {
    // Print results on console
    final boolean print = false;
    // input in tutorial resources directory/output
    final String inputPath = KalmanOrbitDeterminationTest.class.getClassLoader().getResource("orbit-determination/Lageos2/od_test_Lageos2.in").toURI().getPath();
    final File input = new File(inputPath);
    // configure Orekit data acces
    Utils.setDataRoot("orbit-determination/Lageos2:potential/icgem-format");
    GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-6s-truncated", true));
    // Choice of an orbit type to use
    // Default for test is Cartesian
    final OrbitType orbitType = OrbitType.CARTESIAN;
    // Initial orbital Cartesian covariance matrix
    // These covariances are derived from the deltas between initial and reference orbits
    // So in a way they are "perfect"...
    // Cartesian covariance matrix initialization
    final RealMatrix cartesianOrbitalP = MatrixUtils.createRealDiagonalMatrix(new double[] { 1e4, 4e3, 1, 5e-3, 6e-5, 1e-4 });
    // Orbital Cartesian process noise matrix (Q)
    final RealMatrix cartesianOrbitalQ = MatrixUtils.createRealDiagonalMatrix(new double[] { 1.e-4, 1.e-4, 1.e-4, 1.e-10, 1.e-10, 1.e-10 });
    // Initial measurement covariance matrix and process noise matrix
    final RealMatrix measurementP = MatrixUtils.createRealDiagonalMatrix(new double[] { 1., 1., 1., 1. });
    final RealMatrix measurementQ = MatrixUtils.createRealDiagonalMatrix(new double[] { 1e-6, 1e-6, 1e-6, 1e-6 });
    // Kalman orbit determination run.
    ResultKalman kalmanLageos2 = run(input, orbitType, print, cartesianOrbitalP, cartesianOrbitalQ, null, null, measurementP, measurementQ);
    // Definition of the accuracy for the test
    final double distanceAccuracy = 0.86;
    final double velocityAccuracy = 4.12e-3;
    // Tests
    // Note: The reference initial orbit is the same as in the batch LS tests
    // -----
    // Number of measurements processed
    final int numberOfMeas = 258;
    Assert.assertEquals(numberOfMeas, kalmanLageos2.getNumberOfMeasurements());
    // Estimated position and velocity
    final Vector3D estimatedPos = kalmanLageos2.getEstimatedPV().getPosition();
    final Vector3D estimatedVel = kalmanLageos2.getEstimatedPV().getVelocity();
    // Reference position and velocity at initial date (same as in batch LS test)
    final Vector3D refPos0 = new Vector3D(-5532131.956902, 10025696.592156, -3578940.040009);
    final Vector3D refVel0 = new Vector3D(-3871.275109, -607.880985, 4280.972530);
    // Run the reference until Kalman last date
    final Orbit refOrbit = runReference(input, orbitType, refPos0, refVel0, null, kalmanLageos2.getEstimatedPV().getDate());
    final Vector3D refPos = refOrbit.getPVCoordinates().getPosition();
    final Vector3D refVel = refOrbit.getPVCoordinates().getVelocity();
    // Check distances
    final double dP = Vector3D.distance(refPos, estimatedPos);
    final double dV = Vector3D.distance(refVel, estimatedVel);
    Assert.assertEquals(0.0, dP, distanceAccuracy);
    Assert.assertEquals(0.0, dV, velocityAccuracy);
    // Print orbit deltas
    if (print) {
        System.out.println("Test performances:");
        System.out.format("\t%-30s\n", "ΔEstimated / Reference");
        System.out.format(Locale.US, "\t%-10s %20.6f\n", "ΔP [m]", dP);
        System.out.format(Locale.US, "\t%-10s %20.6f\n", "ΔV [m/s]", dV);
    }
    // Test on measurements parameters
    final List<DelegatingDriver> list = new ArrayList<DelegatingDriver>();
    list.addAll(kalmanLageos2.measurementsParameters.getDrivers());
    sortParametersChanges(list);
    // Batch LS values
    // final double[] stationOffSet = { 1.659203,  0.861250,  -0.885352 };
    // final double rangeBias = -0.286275;
    final double[] stationOffSet = { 0.298867, -0.137456, 0.013315 };
    final double rangeBias = 0.002390;
    Assert.assertEquals(stationOffSet[0], list.get(0).getValue(), distanceAccuracy);
    Assert.assertEquals(stationOffSet[1], list.get(1).getValue(), distanceAccuracy);
    Assert.assertEquals(stationOffSet[2], list.get(2).getValue(), distanceAccuracy);
    Assert.assertEquals(rangeBias, list.get(3).getValue(), distanceAccuracy);
    // test on statistic for the range residuals
    final long nbRange = 258;
    // Batch LS values
    // final double[] RefStatRange = { -2.431135, 2.218644, 0.038483, 0.982017 };
    final double[] RefStatRange = { -23.561314, 20.436464, 0.964164, 5.687187 };
    Assert.assertEquals(nbRange, kalmanLageos2.getRangeStat().getN());
    Assert.assertEquals(RefStatRange[0], kalmanLageos2.getRangeStat().getMin(), distanceAccuracy);
    Assert.assertEquals(RefStatRange[1], kalmanLageos2.getRangeStat().getMax(), distanceAccuracy);
    Assert.assertEquals(RefStatRange[2], kalmanLageos2.getRangeStat().getMean(), distanceAccuracy);
    Assert.assertEquals(RefStatRange[3], kalmanLageos2.getRangeStat().getStandardDeviation(), distanceAccuracy);
}

Example 4 with ICGEMFormatReader

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

the class DSSTPropagatorTest method testEphemerisGeneration.

@Test
public void testEphemerisGeneration() throws OrekitException {
    Utils.setDataRoot("regular-data:potential/icgem-format");
    GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("^eigen-6s-truncated$", false));
    UnnormalizedSphericalHarmonicsProvider nshp = GravityFieldFactory.getUnnormalizedProvider(8, 8);
    Orbit orbit = new KeplerianOrbit(13378000, 0.05, 0, 0, FastMath.PI, 0, PositionAngle.MEAN, FramesFactory.getTOD(false), new AbsoluteDate(2003, 5, 6, TimeScalesFactory.getUTC()), nshp.getMu());
    double period = orbit.getKeplerianPeriod();
    double[][] tolerance = DSSTPropagator.tolerances(1.0, orbit);
    AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(period / 100, period * 100, tolerance[0], tolerance[1]);
    integrator.setInitialStepSize(10 * period);
    DSSTPropagator propagator = new DSSTPropagator(integrator, false);
    OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getGTOD(false));
    CelestialBody sun = CelestialBodyFactory.getSun();
    CelestialBody moon = CelestialBodyFactory.getMoon();
    propagator.addForceModel(new DSSTZonal(nshp, 8, 7, 17));
    propagator.addForceModel(new DSSTTesseral(earth.getBodyFrame(), Constants.WGS84_EARTH_ANGULAR_VELOCITY, nshp, 8, 8, 4, 12, 8, 8, 4));
    propagator.addForceModel(new DSSTThirdBody(sun));
    propagator.addForceModel(new DSSTThirdBody(moon));
    propagator.addForceModel(new DSSTAtmosphericDrag(new HarrisPriester(sun, earth), 2.1, 180));
    propagator.addForceModel(new DSSTSolarRadiationPressure(1.2, 180, sun, earth.getEquatorialRadius()));
    propagator.setInterpolationGridToMaxTimeGap(0.5 * Constants.JULIAN_DAY);
    // direct generation of states
    propagator.setInitialState(new SpacecraftState(orbit, 45.0), false);
    final List<SpacecraftState> states = new ArrayList<SpacecraftState>();
    propagator.setMasterMode(600, (currentState, isLast) -> states.add(currentState));
    propagator.propagate(orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY));
    // ephemeris generation
    propagator.setInitialState(new SpacecraftState(orbit, 45.0), false);
    propagator.setEphemerisMode();
    propagator.propagate(orbit.getDate().shiftedBy(30 * Constants.JULIAN_DAY));
    BoundedPropagator ephemeris = propagator.getGeneratedEphemeris();
    double maxError = 0;
    for (final SpacecraftState state : states) {
        final SpacecraftState fromEphemeris = ephemeris.propagate(state.getDate());
        final double error = Vector3D.distance(state.getPVCoordinates().getPosition(), fromEphemeris.getPVCoordinates().getPosition());
        maxError = FastMath.max(maxError, error);
    }
    Assert.assertEquals(0.0, maxError, 1.0e-10);
}

Example 5 with ICGEMFormatReader

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

the class FieldIntegratedEphemerisTest method setUp.

@Before
public void setUp() {
    mu = 3.9860047e14;
    Utils.setDataRoot("regular-data:potential/icgem-format");
    GravityFieldFactory.addPotentialCoefficientsReader(new ICGEMFormatReader("eigen-6s-truncated", true));
}