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

Example 6 with NormalizedSphericalHarmonicsProvider

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

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

the class FieldNumericalPropagatorTest method createPropagator.

private static <T extends RealFieldElement<T>> FieldNumericalPropagator<T> createPropagator(FieldSpacecraftState<T> spacecraftState, OrbitType orbitType, PositionAngle angleType) throws OrekitException {
    final Field<T> field = spacecraftState.getDate().getField();
    final T zero = field.getZero();
    final double minStep = 0.001;
    final double maxStep = 120.0;
    final T positionTolerance = zero.add(0.1);
    final int degree = 20;
    final int order = 20;
    final double spacecraftArea = 1.0;
    final double spacecraftDragCoefficient = 2.0;
    final double spacecraftReflectionCoefficient = 2.0;
    // propagator main configuration
    final double[][] tol = FieldNumericalPropagator.tolerances(positionTolerance, spacecraftState.getOrbit(), orbitType);
    final FieldODEIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, minStep, maxStep, tol[0], tol[1]);
    final FieldNumericalPropagator<T> np = new FieldNumericalPropagator<>(field, integrator);
    np.setOrbitType(orbitType);
    np.setPositionAngleType(angleType);
    np.setInitialState(spacecraftState);
    // Earth gravity field
    final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getITRF(IERSConventions.IERS_2010, true));
    final NormalizedSphericalHarmonicsProvider harmonicsGravityProvider = GravityFieldFactory.getNormalizedProvider(degree, order);
    np.addForceModel(new HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), harmonicsGravityProvider));
    // Sun and Moon attraction
    np.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getSun()));
    np.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory.getMoon()));
    // atmospheric drag
    MarshallSolarActivityFutureEstimation msafe = new MarshallSolarActivityFutureEstimation("Jan2000F10-edited-data\\.txt", MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE);
    DataProvidersManager.getInstance().feed(msafe.getSupportedNames(), msafe);
    DTM2000 atmosphere = new DTM2000(msafe, CelestialBodyFactory.getSun(), earth);
    np.addForceModel(new DragForce(atmosphere, new IsotropicDrag(spacecraftArea, spacecraftDragCoefficient)));
    // solar radiation pressure
    np.addForceModel(new SolarRadiationPressure(CelestialBodyFactory.getSun(), earth.getEquatorialRadius(), new IsotropicRadiationSingleCoefficient(spacecraftArea, spacecraftReflectionCoefficient)));
    return np;
}

Example 8 with NormalizedSphericalHarmonicsProvider

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

Example 9 with NormalizedSphericalHarmonicsProvider

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

the class SolidTidesFieldTest method testInterpolationAccuracy.

@Test
public void testInterpolationAccuracy() throws OrekitException {
    // The shortest periods are slightly below one half day for the tidal waves
    // considered here. This implies the sampling rate should be fast enough.
    // The tuning parameters we have finally settled correspond to a two hours
    // sample containing 12 points (i.e. one new point is computed every 10 minutes).
    // The observed relative interpolation error with these settings are essentially
    // due to Runge phenomenon at points sampling rate. Plotting the errors shows
    // singular peaks pointing out of merely numerical noise.
    final IERSConventions conventions = IERSConventions.IERS_2010;
    Frame itrf = FramesFactory.getITRF(conventions, true);
    TimeScale utc = TimeScalesFactory.getUTC();
    UT1Scale ut1 = TimeScalesFactory.getUT1(conventions, true);
    NormalizedSphericalHarmonicsProvider gravityField = GravityFieldFactory.getConstantNormalizedProvider(5, 5);
    SolidTidesField raw = new SolidTidesField(conventions.getLoveNumbers(), conventions.getTideFrequencyDependenceFunction(ut1), conventions.getPermanentTide(), conventions.getSolidPoleTide(ut1.getEOPHistory()), itrf, gravityField.getAe(), gravityField.getMu(), gravityField.getTideSystem(), CelestialBodyFactory.getSun(), CelestialBodyFactory.getMoon());
    int step = 600;
    int nbPoints = 12;
    CachedNormalizedSphericalHarmonicsProvider interpolated = new CachedNormalizedSphericalHarmonicsProvider(raw, step, nbPoints, OrekitConfiguration.getCacheSlotsNumber(), 7 * Constants.JULIAN_DAY, 0.5 * Constants.JULIAN_DAY);
    // the following time range is located around the maximal observed error
    AbsoluteDate start = new AbsoluteDate(2003, 6, 12, utc);
    AbsoluteDate end = start.shiftedBy(3 * Constants.JULIAN_DAY);
    StreamingStatistics stat = new StreamingStatistics();
    for (AbsoluteDate date = start; date.compareTo(end) < 0; date = date.shiftedBy(60)) {
        NormalizedSphericalHarmonics rawHarmonics = raw.onDate(date);
        NormalizedSphericalHarmonics interpolatedHarmonics = interpolated.onDate(date);
        for (int n = 2; n < 5; ++n) {
            for (int m = 0; m <= n; ++m) {
                if (n < 4 || m < 3) {
                    double cnmRaw = rawHarmonics.getNormalizedCnm(n, m);
                    double cnmInterp = interpolatedHarmonics.getNormalizedCnm(n, m);
                    double errorC = (cnmInterp - cnmRaw) / FastMath.abs(cnmRaw);
                    stat.addValue(errorC);
                    if (m > 0) {
                        double snmRaw = rawHarmonics.getNormalizedSnm(n, m);
                        double snmInterp = interpolatedHarmonics.getNormalizedSnm(n, m);
                        double errorS = (snmInterp - snmRaw) / FastMath.abs(snmRaw);
                        stat.addValue(errorS);
                    }
                }
            }
        }
    }
    Assert.assertEquals(0.0, stat.getMean(), 2.0e-12);
    Assert.assertTrue(stat.getStandardDeviation() < 2.0e-9);
    Assert.assertTrue(stat.getMin() > -9.0e-8);
    Assert.assertTrue(stat.getMax() < 2.2e-7);
}

Example 10 with NormalizedSphericalHarmonicsProvider

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

the class SolidTidesFieldTest method testDeltaCnmSnm.

@Test
public void testDeltaCnmSnm() throws OrekitException {
    NormalizedSphericalHarmonicsProvider gravityField = GravityFieldFactory.getConstantNormalizedProvider(8, 8);
    UT1Scale ut1 = TimeScalesFactory.getUT1(IERSConventions.IERS_2010, true);
    TimeScale utc = TimeScalesFactory.getUTC();
    AbsoluteDate date = new AbsoluteDate(2003, 5, 6, 13, 43, 32.125, utc);
    SolidTidesField tidesField = new SolidTidesField(IERSConventions.IERS_2010.getLoveNumbers(), IERSConventions.IERS_2010.getTideFrequencyDependenceFunction(ut1), IERSConventions.IERS_2010.getPermanentTide(), null, FramesFactory.getITRF(IERSConventions.IERS_2010, true), gravityField.getAe(), gravityField.getMu(), TideSystem.TIDE_FREE, CelestialBodyFactory.getSun(), CelestialBodyFactory.getMoon());
    NormalizedSphericalHarmonics harmonics = tidesField.onDate(date);
    double[][] refDeltaCnm = new double[][] { { 0.0, 0.0, 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0, 0.0, 0.0 }, { -2.6732289327355114E-9, 4.9078992447259636E-9, 3.5894110538262888E-9, 0.0, 0.0 }, // { -2.6598001259383122E-9,   4.907899244804072E-9,   3.5894110542365972E-9,         0.0 ,                 0.0  },
    { -1.290639603871307E-11, -9.287425756410472E-14, 8.356574033404024E-12, -2.2644465207860626E-12, 0.0 }, { 7.888138856951149E-12, -1.4422209452877158E-11, -6.815519349970944E-12, 0.0, 0.0 } };
    double[][] refDeltaSnm = new double[][] { { 0.0, 0.0, 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0, 0.0, 0.0 }, { 0.0, 1.599927449004677E-9, 2.1815888169727694E-9, 0.0, 0.0 }, { 0.0, -4.6129961143785774E-14, 1.8097527720906976E-11, 1.633889224766215E-11, 0.0 }, { 0.0, -4.897228975221076E-12, -4.1034042689652575E-12, 0.0, 0.0 } };
    for (int n = 0; n < refDeltaCnm.length; ++n) {
        double threshold = (n == 2) ? 1.3e-17 : 1.0e-24;
        for (int m = 0; m <= n; ++m) {
            Assert.assertEquals(refDeltaCnm[n][m], harmonics.getNormalizedCnm(n, m), threshold);
            Assert.assertEquals(refDeltaSnm[n][m], harmonics.getNormalizedSnm(n, m), threshold);
        }
    }
}