Search in sources :

Example 11 with Max

use of org.hipparchus.stat.descriptive.rank.Max in project Orekit by CS-SI.

the class TurnAroundRangeAnalyticTest method genericTestStateDerivatives.

/**
 * Generic test function for derivatives with respect to state
 * @param isModifier Use of atmospheric modifiers
 * @param isFiniteDifferences Finite differences reference calculation if true, TurnAroundRange class otherwise
 * @param printResults Print the results ?
 * @throws OrekitException
 */
void genericTestStateDerivatives(final boolean isModifier, final boolean isFiniteDifferences, final boolean printResults, final double refErrorsPMedian, final double refErrorsPMean, final double refErrorsPMax, final double refErrorsVMedian, final double refErrorsVMean, final double refErrorsVMax) throws OrekitException {
    Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    // Context context = EstimationTestUtils.geoStationnaryContext();
    final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 0.001);
    // create perfect range2 measurements
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new TurnAroundRangeMeasurementCreator(context), 1.0, 3.0, 300.0);
    propagator.setSlaveMode();
    double[] errorsP = new double[3 * measurements.size()];
    double[] errorsV = new double[3 * measurements.size()];
    int indexP = 0;
    int indexV = 0;
    // Print the results ? Header
    if (printResults) {
        System.out.format(Locale.US, "%-15s  %-15s  %-23s  %-23s  " + "%10s  %10s  %10s  " + "%10s  %10s  %10s  " + "%10s  %10s  %10s  " + "%10s  %10s  %10s%n", "Master Station", "Slave Station", "Measurement Date", "State Date", "ΔdPx", "ΔdPy", "ΔdPz", "ΔdVx", "ΔdVy", "ΔdVz", "rel ΔdPx", "rel ΔdPy", "rel ΔdPz", "rel ΔdVx", "rel ΔdVy", "rel ΔdVz");
    }
    // Loop on the measurements
    for (final ObservedMeasurement<?> measurement : measurements) {
        // Add modifiers if test implies it
        final TurnAroundRangeTroposphericDelayModifier modifier = new TurnAroundRangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
        if (isModifier) {
            ((TurnAroundRange) measurement).addModifier(modifier);
        }
        // We intentionally propagate to a date which is close to the
        // real spacecraft state but is *not* the accurate date, by
        // compensating only part of the downlink delay. This is done
        // in order to validate the partial derivatives with respect
        // to velocity. If we had chosen the proper state date, the
        // range would have depended only on the current position but
        // not on the current velocity.
        final double meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
        final AbsoluteDate date = measurement.getDate().shiftedBy(-0.75 * meanDelay);
        final SpacecraftState state = propagator.propagate(date);
        final EstimatedMeasurement<TurnAroundRange> TAR = new TurnAroundRangeAnalytic((TurnAroundRange) measurement).theoreticalEvaluationAnalytic(0, 0, propagator.getInitialState(), state);
        if (isModifier) {
            modifier.modify(TAR);
        }
        final double[][] jacobian = TAR.getStateDerivatives(0);
        // Jacobian reference value
        final double[][] jacobianRef;
        if (isFiniteDifferences) {
            // Compute a reference value using finite differences
            jacobianRef = Differentiation.differentiate(new StateFunction() {

                public double[] value(final SpacecraftState state) throws OrekitException {
                    return measurement.estimate(0, 0, new SpacecraftState[] { state }).getEstimatedValue();
                }
            }, measurement.getDimension(), propagator.getAttitudeProvider(), OrbitType.CARTESIAN, PositionAngle.TRUE, 2.0, 3).value(state);
        } else {
            // Compute a reference value using TurnAroundRange class function
            jacobianRef = ((TurnAroundRange) measurement).theoreticalEvaluation(0, 0, new SpacecraftState[] { state }).getStateDerivatives(0);
        }
        // //Test: Test point by point with the debugger
        // if (!isFiniteDifferences && !isModifier) {
        // final EstimatedMeasurement<TurnAroundRange> test =
        // new TurnAroundRangeAnalytic((TurnAroundRange)measurement).theoreticalEvaluationValidation(0, 0, state);
        // }
        // //Test
        Assert.assertEquals(jacobianRef.length, jacobian.length);
        Assert.assertEquals(jacobianRef[0].length, jacobian[0].length);
        double[][] dJacobian = new double[jacobian.length][jacobian[0].length];
        double[][] dJacobianRelative = new double[jacobian.length][jacobian[0].length];
        for (int i = 0; i < jacobian.length; ++i) {
            for (int j = 0; j < jacobian[i].length; ++j) {
                dJacobian[i][j] = jacobian[i][j] - jacobianRef[i][j];
                dJacobianRelative[i][j] = FastMath.abs(dJacobian[i][j] / jacobianRef[i][j]);
                if (j < 3) {
                    errorsP[indexP++] = dJacobianRelative[i][j];
                } else {
                    errorsV[indexV++] = dJacobianRelative[i][j];
                }
            }
        }
        // Print results on the console ? Print the Jacobian
        if (printResults) {
            String masterStationName = ((TurnAroundRange) measurement).getMasterStation().getBaseFrame().getName();
            String slaveStationName = ((TurnAroundRange) measurement).getSlaveStation().getBaseFrame().getName();
            System.out.format(Locale.US, "%-15s  %-15s  %-23s  %-23s  " + "%10.3e  %10.3e  %10.3e  " + "%10.3e  %10.3e  %10.3e  " + "%10.3e  %10.3e  %10.3e  " + "%10.3e  %10.3e  %10.3e%n", masterStationName, slaveStationName, measurement.getDate(), date, dJacobian[0][0], dJacobian[0][1], dJacobian[0][2], dJacobian[0][3], dJacobian[0][4], dJacobian[0][5], dJacobianRelative[0][0], dJacobianRelative[0][1], dJacobianRelative[0][2], dJacobianRelative[0][3], dJacobianRelative[0][4], dJacobianRelative[0][5]);
        }
    }
    // End loop on the measurements
    // Compute some statistics
    final double errorsPMedian = new Median().evaluate(errorsP);
    final double errorsPMean = new Mean().evaluate(errorsP);
    final double errorsPMax = new Max().evaluate(errorsP);
    final double errorsVMedian = new Median().evaluate(errorsV);
    final double errorsVMean = new Mean().evaluate(errorsV);
    final double errorsVMax = new Max().evaluate(errorsV);
    // Print the results on console ? Final results
    if (printResults) {
        System.out.println();
        System.out.format(Locale.US, "Relative errors dR/dP -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n", errorsPMedian, errorsPMean, errorsPMax);
        System.out.format(Locale.US, "Relative errors dR/dV -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n", errorsVMedian, errorsVMean, errorsVMax);
    }
    // Assert the results / max values depend on the test
    Assert.assertEquals(0.0, errorsPMedian, refErrorsPMedian);
    Assert.assertEquals(0.0, errorsPMean, refErrorsPMean);
    Assert.assertEquals(0.0, errorsPMax, refErrorsPMax);
    Assert.assertEquals(0.0, errorsVMedian, refErrorsVMedian);
    Assert.assertEquals(0.0, errorsVMean, refErrorsVMean);
    Assert.assertEquals(0.0, errorsVMax, refErrorsVMax);
}
Also used : Context(org.orekit.estimation.Context) Mean(org.hipparchus.stat.descriptive.moment.Mean) Max(org.hipparchus.stat.descriptive.rank.Max) Median(org.hipparchus.stat.descriptive.rank.Median) AbsoluteDate(org.orekit.time.AbsoluteDate) SpacecraftState(org.orekit.propagation.SpacecraftState) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) Propagator(org.orekit.propagation.Propagator) TurnAroundRangeTroposphericDelayModifier(org.orekit.estimation.measurements.modifiers.TurnAroundRangeTroposphericDelayModifier) OrekitException(org.orekit.errors.OrekitException) StateFunction(org.orekit.utils.StateFunction)

Example 12 with Max

use of org.hipparchus.stat.descriptive.rank.Max in project Orekit by CS-SI.

the class TurnAroundRangeTest method genericTestParameterDerivatives.

void genericTestParameterDerivatives(final boolean isModifier, final boolean printResults, final double refErrorQMMedian, final double refErrorQMMean, final double refErrorQMMax, final double refErrorQSMedian, final double refErrorQSMean, final double refErrorQSMax) throws OrekitException {
    Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 0.001);
    // Create perfect TAR measurements
    for (Map.Entry<GroundStation, GroundStation> entry : context.TARstations.entrySet()) {
        final GroundStation masterStation = entry.getKey();
        final GroundStation slaveStation = entry.getValue();
        masterStation.getEastOffsetDriver().setSelected(true);
        masterStation.getNorthOffsetDriver().setSelected(true);
        masterStation.getZenithOffsetDriver().setSelected(true);
        slaveStation.getEastOffsetDriver().setSelected(true);
        slaveStation.getNorthOffsetDriver().setSelected(true);
        slaveStation.getZenithOffsetDriver().setSelected(true);
    }
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new TurnAroundRangeMeasurementCreator(context), 1.0, 3.0, 300.0);
    propagator.setSlaveMode();
    // Print results on console ? Header
    if (printResults) {
        System.out.format(Locale.US, "%-15s %-15s %-23s  %-23s  " + "%10s  %10s  %10s  " + "%10s  %10s  %10s  " + "%10s  %10s  %10s  " + "%10s  %10s  %10s%n", "Master Station", "Slave Station", "Measurement Date", "State Date", "ΔdQMx", "rel ΔdQMx", "ΔdQMy", "rel ΔdQMy", "ΔdQMz", "rel ΔdQMz", "ΔdQSx", "rel ΔdQSx", "ΔdQSy", "rel ΔdQSy", "ΔdQSz", "rel ΔdQSz");
    }
    // List to store the results for master and slave station
    final List<Double> relErrorQMList = new ArrayList<Double>();
    final List<Double> relErrorQSList = new ArrayList<Double>();
    // Loop on the measurements
    for (final ObservedMeasurement<?> measurement : measurements) {
        // Add modifiers if test implies it
        final TurnAroundRangeTroposphericDelayModifier modifier = new TurnAroundRangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
        if (isModifier) {
            ((TurnAroundRange) measurement).addModifier(modifier);
        }
        // parameter corresponding to station position offset
        final GroundStation masterStationParameter = ((TurnAroundRange) measurement).getMasterStation();
        final GroundStation slaveStationParameter = ((TurnAroundRange) measurement).getSlaveStation();
        // We intentionally propagate to a date which is close to the
        // real spacecraft state but is *not* the accurate date, by
        // compensating only part of the downlink delay. This is done
        // in order to validate the partial derivatives with respect
        // to velocity. If we had chosen the proper state date, the
        // range would have depended only on the current position but
        // not on the current velocity.
        final double meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
        final AbsoluteDate date = measurement.getDate().shiftedBy(-0.75 * meanDelay);
        final SpacecraftState state = propagator.propagate(date);
        final ParameterDriver[] drivers = new ParameterDriver[] { masterStationParameter.getEastOffsetDriver(), masterStationParameter.getNorthOffsetDriver(), masterStationParameter.getZenithOffsetDriver(), slaveStationParameter.getEastOffsetDriver(), slaveStationParameter.getNorthOffsetDriver(), slaveStationParameter.getZenithOffsetDriver() };
        // Print results on console ? Stations' names
        if (printResults) {
            String masterStationName = masterStationParameter.getBaseFrame().getName();
            String slaveStationName = slaveStationParameter.getBaseFrame().getName();
            System.out.format(Locale.US, "%-15s %-15s %-23s  %-23s  ", masterStationName, slaveStationName, measurement.getDate(), date);
        }
        // Loop on the parameters
        for (int i = 0; i < 6; ++i) {
            final double[] gradient = measurement.estimate(0, 0, new SpacecraftState[] { state }).getParameterDerivatives(drivers[i]);
            Assert.assertEquals(1, measurement.getDimension());
            Assert.assertEquals(1, gradient.length);
            // Compute a reference value using finite differences
            final ParameterFunction dMkdP = Differentiation.differentiate(new ParameterFunction() {

                /**
                 * {@inheritDoc}
                 */
                @Override
                public double value(final ParameterDriver parameterDriver) throws OrekitException {
                    return measurement.estimate(0, 0, new SpacecraftState[] { state }).getEstimatedValue()[0];
                }
            }, drivers[i], 3, 20.0);
            final double ref = dMkdP.value(drivers[i]);
            // Deltas
            double dGradient = gradient[0] - ref;
            double dGradientRelative = FastMath.abs(dGradient / ref);
            // Print results on console ? Gradient difference
            if (printResults) {
                System.out.format(Locale.US, "%10.3e  %10.3e  ", dGradient, dGradientRelative);
            }
            // Add relative error to the list
            if (i < 3) {
                relErrorQMList.add(dGradientRelative);
            } else {
                relErrorQSList.add(dGradientRelative);
            }
        }
        // End for loop on the parameters
        if (printResults) {
            System.out.format(Locale.US, "%n");
        }
    }
    // End for loop on the measurements
    // Convert error list to double[]
    final double[] relErrorQM = relErrorQMList.stream().mapToDouble(Double::doubleValue).toArray();
    final double[] relErrorQS = relErrorQSList.stream().mapToDouble(Double::doubleValue).toArray();
    // Compute statistics
    final double relErrorsQMMedian = new Median().evaluate(relErrorQM);
    final double relErrorsQMMean = new Mean().evaluate(relErrorQM);
    final double relErrorsQMMax = new Max().evaluate(relErrorQM);
    final double relErrorsQSMedian = new Median().evaluate(relErrorQS);
    final double relErrorsQSMean = new Mean().evaluate(relErrorQS);
    final double relErrorsQSMax = new Max().evaluate(relErrorQS);
    // Print the results on console ?
    if (printResults) {
        System.out.println();
        System.out.format(Locale.US, "Relative errors dR/dQ master station -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n", relErrorsQMMedian, relErrorsQMMean, relErrorsQMMax);
        System.out.format(Locale.US, "Relative errors dR/dQ slave station  -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n", relErrorsQSMedian, relErrorsQSMean, relErrorsQSMax);
    }
    // Check values
    Assert.assertEquals(0.0, relErrorsQMMedian, refErrorQMMedian);
    Assert.assertEquals(0.0, relErrorsQMMean, refErrorQMMean);
    Assert.assertEquals(0.0, relErrorsQMMax, refErrorQMMax);
    Assert.assertEquals(0.0, relErrorsQSMedian, refErrorQSMedian);
    Assert.assertEquals(0.0, relErrorsQSMean, refErrorQSMean);
    Assert.assertEquals(0.0, relErrorsQSMax, refErrorQSMax);
}
Also used : Mean(org.hipparchus.stat.descriptive.moment.Mean) Max(org.hipparchus.stat.descriptive.rank.Max) ArrayList(java.util.ArrayList) Median(org.hipparchus.stat.descriptive.rank.Median) AbsoluteDate(org.orekit.time.AbsoluteDate) SpacecraftState(org.orekit.propagation.SpacecraftState) Propagator(org.orekit.propagation.Propagator) TurnAroundRangeTroposphericDelayModifier(org.orekit.estimation.measurements.modifiers.TurnAroundRangeTroposphericDelayModifier) OrekitException(org.orekit.errors.OrekitException) Context(org.orekit.estimation.Context) ParameterDriver(org.orekit.utils.ParameterDriver) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) ParameterFunction(org.orekit.utils.ParameterFunction) Map(java.util.Map)

Example 13 with Max

use of org.hipparchus.stat.descriptive.rank.Max in project Orekit by CS-SI.

the class RangeAnalyticTest method genericTestValues.

/**
 * Generic test function for values of the range
 * @param printResults Print the results ?
 * @throws OrekitException
 */
void genericTestValues(final boolean printResults) throws OrekitException {
    Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 0.001);
    // Create perfect range measurements
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new RangeMeasurementCreator(context), 1.0, 3.0, 300.0);
    // Lists for results' storage - Used only for derivatives with respect to state
    // "final" value to be seen by "handleStep" function of the propagator
    final List<Double> absoluteErrors = new ArrayList<Double>();
    final List<Double> relativeErrors = new ArrayList<Double>();
    // Set master mode
    // Use a lambda function to implement "handleStep" function
    propagator.setMasterMode((OrekitStepInterpolator interpolator, boolean isLast) -> {
        for (final ObservedMeasurement<?> measurement : measurements) {
            // Play test if the measurement date is between interpolator previous and current date
            if ((measurement.getDate().durationFrom(interpolator.getPreviousState().getDate()) > 0.) && (measurement.getDate().durationFrom(interpolator.getCurrentState().getDate()) <= 0.)) {
                // We intentionally propagate to a date which is close to the
                // real spacecraft state but is *not* the accurate date, by
                // compensating only part of the downlink delay. This is done
                // in order to validate the partial derivatives with respect
                // to velocity. If we had chosen the proper state date, the
                // range would have depended only on the current position but
                // not on the current velocity.
                final double meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
                final AbsoluteDate date = measurement.getDate().shiftedBy(-0.75 * meanDelay);
                final SpacecraftState state = interpolator.getInterpolatedState(date);
                // Values of the RangeAnalytic & errors
                final double RangeObserved = measurement.getObservedValue()[0];
                final double RangeEstimated = new RangeAnalytic((Range) measurement).theoreticalEvaluationAnalytic(0, 0, state).getEstimatedValue()[0];
                final double absoluteError = RangeEstimated - RangeObserved;
                absoluteErrors.add(absoluteError);
                relativeErrors.add(FastMath.abs(absoluteError) / FastMath.abs(RangeObserved));
                // Print results on console ?
                if (printResults) {
                    final AbsoluteDate measurementDate = measurement.getDate();
                    String stationName = ((Range) measurement).getStation().getBaseFrame().getName();
                    System.out.format(Locale.US, "%-15s  %-23s  %-23s  %19.6f  %19.6f  %13.6e  %13.6e%n", stationName, measurementDate, date, RangeObserved, RangeEstimated, FastMath.abs(RangeEstimated - RangeObserved), FastMath.abs((RangeEstimated - RangeObserved) / RangeObserved));
                }
            }
        // End if measurement date between previous and current interpolator step
        }
    // End for loop on the measurements
    });
    // Print results on console ? Header
    if (printResults) {
        System.out.format(Locale.US, "%-15s  %-23s  %-23s  %19s  %19s  %13s  %13s%n", "Station", "Measurement Date", "State Date", "Range observed [m]", "Range estimated [m]", "ΔRange [m]", "rel ΔRange");
    }
    // Rewind the propagator to initial date
    propagator.propagate(context.initialOrbit.getDate());
    // Sort measurements chronologically
    measurements.sort(new ChronologicalComparator());
    // Propagate to final measurement's date
    propagator.propagate(measurements.get(measurements.size() - 1).getDate());
    // Convert lists to double array
    final double[] absErrors = absoluteErrors.stream().mapToDouble(Double::doubleValue).toArray();
    final double[] relErrors = relativeErrors.stream().mapToDouble(Double::doubleValue).toArray();
    // Statistics' assertion
    final double absErrorsMedian = new Median().evaluate(absErrors);
    final double absErrorsMin = new Min().evaluate(absErrors);
    final double absErrorsMax = new Max().evaluate(absErrors);
    final double relErrorsMedian = new Median().evaluate(relErrors);
    final double relErrorsMax = new Max().evaluate(relErrors);
    // Print the results on console ? Final results
    if (printResults) {
        System.out.println();
        System.out.println("Absolute errors median: " + absErrorsMedian);
        System.out.println("Absolute errors min   : " + absErrorsMin);
        System.out.println("Absolute errors max   : " + absErrorsMax);
        System.out.println("Relative errors median: " + relErrorsMedian);
        System.out.println("Relative errors max   : " + relErrorsMax);
    }
    Assert.assertEquals(0.0, absErrorsMedian, 3.8e-08);
    Assert.assertEquals(0.0, absErrorsMin, 2.0e-07);
    Assert.assertEquals(0.0, absErrorsMax, 2.3e-07);
    Assert.assertEquals(0.0, relErrorsMedian, 6.5e-15);
    Assert.assertEquals(0.0, relErrorsMax, 2.4e-14);
}
Also used : Context(org.orekit.estimation.Context) Max(org.hipparchus.stat.descriptive.rank.Max) ArrayList(java.util.ArrayList) Median(org.hipparchus.stat.descriptive.rank.Median) AbsoluteDate(org.orekit.time.AbsoluteDate) OrekitStepInterpolator(org.orekit.propagation.sampling.OrekitStepInterpolator) SpacecraftState(org.orekit.propagation.SpacecraftState) Min(org.hipparchus.stat.descriptive.rank.Min) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) Propagator(org.orekit.propagation.Propagator) ChronologicalComparator(org.orekit.time.ChronologicalComparator)

Example 14 with Max

use of org.hipparchus.stat.descriptive.rank.Max in project Orekit by CS-SI.

the class RangeTest method genericTestStateDerivatives.

void genericTestStateDerivatives(final boolean isModifier, final boolean printResults, final double refErrorsPMedian, final double refErrorsPMean, final double refErrorsPMax, final double refErrorsVMedian, final double refErrorsVMean, final double refErrorsVMax) throws OrekitException {
    Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 0.001);
    // Create perfect range measurements
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new RangeMeasurementCreator(context), 1.0, 3.0, 300.0);
    // Lists for results' storage - Used only for derivatives with respect to state
    // "final" value to be seen by "handleStep" function of the propagator
    final List<Double> errorsP = new ArrayList<Double>();
    final List<Double> errorsV = new ArrayList<Double>();
    // Set master mode
    // Use a lambda function to implement "handleStep" function
    propagator.setMasterMode((OrekitStepInterpolator interpolator, boolean isLast) -> {
        for (final ObservedMeasurement<?> measurement : measurements) {
            // Play test if the measurement date is between interpolator previous and current date
            if ((measurement.getDate().durationFrom(interpolator.getPreviousState().getDate()) > 0.) && (measurement.getDate().durationFrom(interpolator.getCurrentState().getDate()) <= 0.)) {
                // Add modifiers if test implies it
                final RangeTroposphericDelayModifier modifier = new RangeTroposphericDelayModifier(SaastamoinenModel.getStandardModel());
                if (isModifier) {
                    ((Range) measurement).addModifier(modifier);
                }
                // We intentionally propagate to a date which is close to the
                // real spacecraft state but is *not* the accurate date, by
                // compensating only part of the downlink delay. This is done
                // in order to validate the partial derivatives with respect
                // to velocity. If we had chosen the proper state date, the
                // range would have depended only on the current position but
                // not on the current velocity.
                final double meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
                final AbsoluteDate date = measurement.getDate().shiftedBy(-0.75 * meanDelay);
                final SpacecraftState state = interpolator.getInterpolatedState(date);
                final double[][] jacobian = measurement.estimate(0, 0, new SpacecraftState[] { state }).getStateDerivatives(0);
                // Jacobian reference value
                final double[][] jacobianRef;
                // Compute a reference value using finite differences
                jacobianRef = Differentiation.differentiate(new StateFunction() {

                    public double[] value(final SpacecraftState state) throws OrekitException {
                        return measurement.estimate(0, 0, new SpacecraftState[] { state }).getEstimatedValue();
                    }
                }, measurement.getDimension(), propagator.getAttitudeProvider(), OrbitType.CARTESIAN, PositionAngle.TRUE, 2.0, 3).value(state);
                Assert.assertEquals(jacobianRef.length, jacobian.length);
                Assert.assertEquals(jacobianRef[0].length, jacobian[0].length);
                // Errors & relative errors on the Jacobian
                double[][] dJacobian = new double[jacobian.length][jacobian[0].length];
                double[][] dJacobianRelative = new double[jacobian.length][jacobian[0].length];
                for (int i = 0; i < jacobian.length; ++i) {
                    for (int j = 0; j < jacobian[i].length; ++j) {
                        dJacobian[i][j] = jacobian[i][j] - jacobianRef[i][j];
                        dJacobianRelative[i][j] = FastMath.abs(dJacobian[i][j] / jacobianRef[i][j]);
                        if (j < 3) {
                            errorsP.add(dJacobianRelative[i][j]);
                        } else {
                            errorsV.add(dJacobianRelative[i][j]);
                        }
                    }
                }
                // Print values in console ?
                if (printResults) {
                    String stationName = ((Range) measurement).getStation().getBaseFrame().getName();
                    System.out.format(Locale.US, "%-15s  %-23s  %-23s  " + "%10.3e  %10.3e  %10.3e  " + "%10.3e  %10.3e  %10.3e  " + "%10.3e  %10.3e  %10.3e  " + "%10.3e  %10.3e  %10.3e%n", stationName, measurement.getDate(), date, dJacobian[0][0], dJacobian[0][1], dJacobian[0][2], dJacobian[0][3], dJacobian[0][4], dJacobian[0][5], dJacobianRelative[0][0], dJacobianRelative[0][1], dJacobianRelative[0][2], dJacobianRelative[0][3], dJacobianRelative[0][4], dJacobianRelative[0][5]);
                }
            }
        // End if measurement date between previous and current interpolator step
        }
    // End for loop on the measurements
    });
    // Print results on console ?
    if (printResults) {
        System.out.format(Locale.US, "%-15s  %-23s  %-23s  " + "%10s  %10s  %10s  " + "%10s  %10s  %10s  " + "%10s  %10s  %10s  " + "%10s  %10s  %10s%n", "Station", "Measurement Date", "State Date", "ΔdPx", "ΔdPy", "ΔdPz", "ΔdVx", "ΔdVy", "ΔdVz", "rel ΔdPx", "rel ΔdPy", "rel ΔdPz", "rel ΔdVx", "rel ΔdVy", "rel ΔdVz");
    }
    // Rewind the propagator to initial date
    propagator.propagate(context.initialOrbit.getDate());
    // Sort measurements chronologically
    measurements.sort(new ChronologicalComparator());
    // Propagate to final measurement's date
    propagator.propagate(measurements.get(measurements.size() - 1).getDate());
    // Convert lists to double[] and evaluate some statistics
    final double[] relErrorsP = errorsP.stream().mapToDouble(Double::doubleValue).toArray();
    final double[] relErrorsV = errorsV.stream().mapToDouble(Double::doubleValue).toArray();
    final double errorsPMedian = new Median().evaluate(relErrorsP);
    final double errorsPMean = new Mean().evaluate(relErrorsP);
    final double errorsPMax = new Max().evaluate(relErrorsP);
    final double errorsVMedian = new Median().evaluate(relErrorsV);
    final double errorsVMean = new Mean().evaluate(relErrorsV);
    final double errorsVMax = new Max().evaluate(relErrorsV);
    // Print the results on console ?
    if (printResults) {
        System.out.println();
        System.out.format(Locale.US, "Relative errors dR/dP -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n", errorsPMedian, errorsPMean, errorsPMax);
        System.out.format(Locale.US, "Relative errors dR/dV -> Median: %6.3e / Mean: %6.3e / Max: %6.3e%n", errorsVMedian, errorsVMean, errorsVMax);
    }
    Assert.assertEquals(0.0, errorsPMedian, refErrorsPMedian);
    Assert.assertEquals(0.0, errorsPMean, refErrorsPMean);
    Assert.assertEquals(0.0, errorsPMax, refErrorsPMax);
    Assert.assertEquals(0.0, errorsVMedian, refErrorsVMedian);
    Assert.assertEquals(0.0, errorsVMean, refErrorsVMean);
    Assert.assertEquals(0.0, errorsVMax, refErrorsVMax);
}
Also used : Mean(org.hipparchus.stat.descriptive.moment.Mean) Max(org.hipparchus.stat.descriptive.rank.Max) ArrayList(java.util.ArrayList) Median(org.hipparchus.stat.descriptive.rank.Median) AbsoluteDate(org.orekit.time.AbsoluteDate) SpacecraftState(org.orekit.propagation.SpacecraftState) Propagator(org.orekit.propagation.Propagator) OrekitException(org.orekit.errors.OrekitException) Context(org.orekit.estimation.Context) RangeTroposphericDelayModifier(org.orekit.estimation.measurements.modifiers.RangeTroposphericDelayModifier) OrekitStepInterpolator(org.orekit.propagation.sampling.OrekitStepInterpolator) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) StateFunction(org.orekit.utils.StateFunction) ChronologicalComparator(org.orekit.time.ChronologicalComparator)

Example 15 with Max

use of org.hipparchus.stat.descriptive.rank.Max in project Orekit by CS-SI.

the class FieldKeplerianPropagatorTest method testTuple.

@Test
public void testTuple() throws OrekitException {
    AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(584.);
    KeplerianOrbit k0 = new KeplerianOrbit(7209668.0, 0.5e-4, 1.7, 2.1, 2.9, 6.2, PositionAngle.TRUE, FramesFactory.getEME2000(), initDate, mu);
    TimeStampedPVCoordinates pv0 = k0.getPVCoordinates();
    TimeStampedPVCoordinates pv1 = new TimeStampedPVCoordinates(pv0.getDate(), pv0.getPosition(), pv0.getVelocity().add(new Vector3D(2.0, pv0.getVelocity().normalize())));
    KeplerianOrbit k1 = new KeplerianOrbit(pv1, k0.getFrame(), k0.getMu());
    FieldOrbit<Tuple> twoOrbits = new FieldKeplerianOrbit<>(new Tuple(k0.getA(), k1.getA()), new Tuple(k0.getE(), k1.getE()), new Tuple(k0.getI(), k1.getI()), new Tuple(k0.getPerigeeArgument(), k1.getPerigeeArgument()), new Tuple(k0.getRightAscensionOfAscendingNode(), k1.getRightAscensionOfAscendingNode()), new Tuple(k0.getMeanAnomaly(), k1.getMeanAnomaly()), PositionAngle.MEAN, FramesFactory.getEME2000(), new FieldAbsoluteDate<>(initDate, new Tuple(0.0, 0.0)), mu);
    Field<Tuple> field = twoOrbits.getDate().getField();
    FieldPropagator<Tuple> propagator = new FieldKeplerianPropagator<>(twoOrbits);
    Min minTangential = new Min();
    Max maxTangential = new Max();
    Min minRadial = new Min();
    Max maxRadial = new Max();
    propagator.setMasterMode(field.getZero().add(10.0), (s, isLast) -> {
        FieldVector3D<Tuple> p = s.getPVCoordinates().getPosition();
        FieldVector3D<Tuple> v = s.getPVCoordinates().getVelocity();
        Vector3D p0 = new Vector3D(p.getX().getComponent(0), p.getY().getComponent(0), p.getZ().getComponent(0));
        Vector3D v0 = new Vector3D(v.getX().getComponent(0), v.getY().getComponent(0), v.getZ().getComponent(0));
        Vector3D t = v0.normalize();
        Vector3D r = Vector3D.crossProduct(v0, Vector3D.crossProduct(p0, v0)).normalize();
        Vector3D p1 = new Vector3D(p.getX().getComponent(1), p.getY().getComponent(1), p.getZ().getComponent(1));
        double dT = Vector3D.dotProduct(p1.subtract(p0), t);
        double dR = Vector3D.dotProduct(p1.subtract(p0), r);
        minTangential.increment(dT);
        maxTangential.increment(dT);
        minRadial.increment(dR);
        maxRadial.increment(dR);
    });
    propagator.propagate(twoOrbits.getDate().shiftedBy(Constants.JULIAN_DAY / 8));
    Assert.assertEquals(-72525.685, minTangential.getResult(), 1.0e-3);
    Assert.assertEquals(926.046, maxTangential.getResult(), 1.0e-3);
    Assert.assertEquals(-92.800, minRadial.getResult(), 1.0e-3);
    Assert.assertEquals(7739.532, maxRadial.getResult(), 1.0e-3);
}
Also used : Max(org.hipparchus.stat.descriptive.rank.Max) TimeStampedPVCoordinates(org.orekit.utils.TimeStampedPVCoordinates) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) AbsoluteDate(org.orekit.time.AbsoluteDate) FieldKeplerianOrbit(org.orekit.orbits.FieldKeplerianOrbit) Min(org.hipparchus.stat.descriptive.rank.Min) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) FieldKeplerianOrbit(org.orekit.orbits.FieldKeplerianOrbit) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) Tuple(org.hipparchus.util.Tuple) Test(org.junit.Test)

Aggregations

Max (org.hipparchus.stat.descriptive.rank.Max)15 AbsoluteDate (org.orekit.time.AbsoluteDate)15 Median (org.hipparchus.stat.descriptive.rank.Median)14 Context (org.orekit.estimation.Context)14 Propagator (org.orekit.propagation.Propagator)14 SpacecraftState (org.orekit.propagation.SpacecraftState)14 NumericalPropagatorBuilder (org.orekit.propagation.conversion.NumericalPropagatorBuilder)14 ArrayList (java.util.ArrayList)10 Mean (org.hipparchus.stat.descriptive.moment.Mean)9 OrekitException (org.orekit.errors.OrekitException)8 OrekitStepInterpolator (org.orekit.propagation.sampling.OrekitStepInterpolator)8 ChronologicalComparator (org.orekit.time.ChronologicalComparator)8 Min (org.hipparchus.stat.descriptive.rank.Min)6 StateFunction (org.orekit.utils.StateFunction)5 TimeStampedPVCoordinates (org.orekit.utils.TimeStampedPVCoordinates)5 RangeTroposphericDelayModifier (org.orekit.estimation.measurements.modifiers.RangeTroposphericDelayModifier)4 TurnAroundRangeTroposphericDelayModifier (org.orekit.estimation.measurements.modifiers.TurnAroundRangeTroposphericDelayModifier)4 ParameterDriver (org.orekit.utils.ParameterDriver)4 Vector3D (org.hipparchus.geometry.euclidean.threed.Vector3D)3 ParameterFunction (org.orekit.utils.ParameterFunction)3