Search in sources :

Example 36 with SpacecraftState

use of org.orekit.propagation.SpacecraftState in project Orekit by CS-SI.

the class OrbitDeterminationTest method createOrbit.

/**
 * Create an orbit from input parameters
 * @param parser input file parser
 * @param mu     central attraction coefficient
 * @throws NoSuchElementException if input parameters are missing
 * @throws OrekitException if inertial frame cannot be created
 */
private Orbit createOrbit(final KeyValueFileParser<ParameterKey> parser, final double mu) throws NoSuchElementException, OrekitException {
    final Frame frame;
    if (!parser.containsKey(ParameterKey.INERTIAL_FRAME)) {
        frame = FramesFactory.getEME2000();
    } else {
        frame = parser.getInertialFrame(ParameterKey.INERTIAL_FRAME);
    }
    // Orbit definition
    PositionAngle angleType = PositionAngle.MEAN;
    if (parser.containsKey(ParameterKey.ORBIT_ANGLE_TYPE)) {
        angleType = PositionAngle.valueOf(parser.getString(ParameterKey.ORBIT_ANGLE_TYPE).toUpperCase());
    }
    if (parser.containsKey(ParameterKey.ORBIT_KEPLERIAN_A)) {
        return new KeplerianOrbit(parser.getDouble(ParameterKey.ORBIT_KEPLERIAN_A), parser.getDouble(ParameterKey.ORBIT_KEPLERIAN_E), parser.getAngle(ParameterKey.ORBIT_KEPLERIAN_I), parser.getAngle(ParameterKey.ORBIT_KEPLERIAN_PA), parser.getAngle(ParameterKey.ORBIT_KEPLERIAN_RAAN), parser.getAngle(ParameterKey.ORBIT_KEPLERIAN_ANOMALY), angleType, frame, parser.getDate(ParameterKey.ORBIT_DATE, TimeScalesFactory.getUTC()), mu);
    } else if (parser.containsKey(ParameterKey.ORBIT_EQUINOCTIAL_A)) {
        return new EquinoctialOrbit(parser.getDouble(ParameterKey.ORBIT_EQUINOCTIAL_A), parser.getDouble(ParameterKey.ORBIT_EQUINOCTIAL_EX), parser.getDouble(ParameterKey.ORBIT_EQUINOCTIAL_EY), parser.getDouble(ParameterKey.ORBIT_EQUINOCTIAL_HX), parser.getDouble(ParameterKey.ORBIT_EQUINOCTIAL_HY), parser.getAngle(ParameterKey.ORBIT_EQUINOCTIAL_LAMBDA), angleType, frame, parser.getDate(ParameterKey.ORBIT_DATE, TimeScalesFactory.getUTC()), mu);
    } else if (parser.containsKey(ParameterKey.ORBIT_CIRCULAR_A)) {
        return new CircularOrbit(parser.getDouble(ParameterKey.ORBIT_CIRCULAR_A), parser.getDouble(ParameterKey.ORBIT_CIRCULAR_EX), parser.getDouble(ParameterKey.ORBIT_CIRCULAR_EY), parser.getAngle(ParameterKey.ORBIT_CIRCULAR_I), parser.getAngle(ParameterKey.ORBIT_CIRCULAR_RAAN), parser.getAngle(ParameterKey.ORBIT_CIRCULAR_ALPHA), angleType, frame, parser.getDate(ParameterKey.ORBIT_DATE, TimeScalesFactory.getUTC()), mu);
    } else if (parser.containsKey(ParameterKey.ORBIT_TLE_LINE_1)) {
        final String line1 = parser.getString(ParameterKey.ORBIT_TLE_LINE_1);
        final String line2 = parser.getString(ParameterKey.ORBIT_TLE_LINE_2);
        final TLE tle = new TLE(line1, line2);
        TLEPropagator propagator = TLEPropagator.selectExtrapolator(tle);
        // propagator.setEphemerisMode();
        AbsoluteDate initDate = tle.getDate();
        SpacecraftState initialState = propagator.getInitialState();
        // Transformation from TEME to frame.
        Transform t = FramesFactory.getTEME().getTransformTo(FramesFactory.getEME2000(), initDate.getDate());
        return new CartesianOrbit(t.transformPVCoordinates(initialState.getPVCoordinates()), frame, initDate, mu);
    } else {
        final double[] pos = { parser.getDouble(ParameterKey.ORBIT_CARTESIAN_PX), parser.getDouble(ParameterKey.ORBIT_CARTESIAN_PY), parser.getDouble(ParameterKey.ORBIT_CARTESIAN_PZ) };
        final double[] vel = { parser.getDouble(ParameterKey.ORBIT_CARTESIAN_VX), parser.getDouble(ParameterKey.ORBIT_CARTESIAN_VY), parser.getDouble(ParameterKey.ORBIT_CARTESIAN_VZ) };
        return new CartesianOrbit(new PVCoordinates(new Vector3D(pos), new Vector3D(vel)), frame, parser.getDate(ParameterKey.ORBIT_DATE, TimeScalesFactory.getUTC()), mu);
    }
}
Also used : Frame(org.orekit.frames.Frame) TopocentricFrame(org.orekit.frames.TopocentricFrame) CartesianOrbit(org.orekit.orbits.CartesianOrbit) PositionAngle(org.orekit.orbits.PositionAngle) PVCoordinates(org.orekit.utils.PVCoordinates) TimeStampedPVCoordinates(org.orekit.utils.TimeStampedPVCoordinates) TLEPropagator(org.orekit.propagation.analytical.tle.TLEPropagator) TLE(org.orekit.propagation.analytical.tle.TLE) AbsoluteDate(org.orekit.time.AbsoluteDate) SpacecraftState(org.orekit.propagation.SpacecraftState) CircularOrbit(org.orekit.orbits.CircularOrbit) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) EquinoctialOrbit(org.orekit.orbits.EquinoctialOrbit) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) Transform(org.orekit.frames.Transform)

Example 37 with SpacecraftState

use of org.orekit.propagation.SpacecraftState in project Orekit by CS-SI.

the class AngularRaDecTest method testStateDerivatives.

@Test
public void testStateDerivatives() throws OrekitException {
    Context context = EstimationTestUtils.geoStationnaryContext("regular-data:potential:tides");
    final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.EQUINOCTIAL, PositionAngle.TRUE, false, 1.0e-6, 60.0, 0.001);
    // create perfect azimuth-elevation measurements
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new AngularRaDecMeasurementCreator(context), 0.25, 3.0, 600.0);
    propagator.setSlaveMode();
    // Compute measurements.
    double[] RaerrorsP = new double[3 * measurements.size()];
    double[] RaerrorsV = new double[3 * measurements.size()];
    double[] DecerrorsP = new double[3 * measurements.size()];
    double[] DecerrorsV = new double[3 * measurements.size()];
    int RaindexP = 0;
    int RaindexV = 0;
    int DecindexP = 0;
    int DecindexV = 0;
    for (final ObservedMeasurement<?> measurement : measurements) {
        // parameter corresponding to station position offset
        final GroundStation stationParameter = ((AngularRaDec) measurement).getStation();
        // 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
        // angular would have depended only on the current position but
        // not on the current velocity.
        final AbsoluteDate datemeas = measurement.getDate();
        SpacecraftState state = propagator.propagate(datemeas);
        final Vector3D stationP = stationParameter.getOffsetToInertial(state.getFrame(), datemeas).transformPosition(Vector3D.ZERO);
        final double meanDelay = AbstractMeasurement.signalTimeOfFlight(state.getPVCoordinates(), stationP, datemeas);
        final AbsoluteDate date = measurement.getDate().shiftedBy(-0.75 * meanDelay);
        state = propagator.propagate(date);
        final EstimatedMeasurement<?> estimated = measurement.estimate(0, 0, new SpacecraftState[] { state });
        Assert.assertEquals(2, estimated.getParticipants().length);
        final double[][] jacobian = estimated.getStateDerivatives(0);
        // compute a reference value using finite differences
        final double[][] finiteDifferencesJacobian = 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, 250.0, 4).value(state);
        Assert.assertEquals(finiteDifferencesJacobian.length, jacobian.length);
        Assert.assertEquals(finiteDifferencesJacobian[0].length, jacobian[0].length);
        final double smallest = FastMath.ulp((double) 1.0);
        for (int i = 0; i < jacobian.length; ++i) {
            for (int j = 0; j < jacobian[i].length; ++j) {
                double relativeError = FastMath.abs((finiteDifferencesJacobian[i][j] - jacobian[i][j]) / finiteDifferencesJacobian[i][j]);
                if ((FastMath.sqrt(finiteDifferencesJacobian[i][j]) < smallest) && (FastMath.sqrt(jacobian[i][j]) < smallest)) {
                    relativeError = 0.0;
                }
                if (j < 3) {
                    if (i == 0) {
                        RaerrorsP[RaindexP++] = relativeError;
                    } else {
                        DecerrorsP[DecindexP++] = relativeError;
                    }
                } else {
                    if (i == 0) {
                        RaerrorsV[RaindexV++] = relativeError;
                    } else {
                        DecerrorsV[DecindexV++] = relativeError;
                    }
                }
            }
        }
    }
    // median errors on Azimuth
    Assert.assertEquals(0.0, new Median().evaluate(RaerrorsP), 4.8e-11);
    Assert.assertEquals(0.0, new Median().evaluate(RaerrorsV), 2.2e-5);
    // median errors on Elevation
    Assert.assertEquals(0.0, new Median().evaluate(DecerrorsP), 1.5e-11);
    Assert.assertEquals(0.0, new Median().evaluate(DecerrorsV), 5.4e-6);
}
Also used : Context(org.orekit.estimation.Context) Median(org.hipparchus.stat.descriptive.rank.Median) AbsoluteDate(org.orekit.time.AbsoluteDate) SpacecraftState(org.orekit.propagation.SpacecraftState) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) Propagator(org.orekit.propagation.Propagator) StateFunction(org.orekit.utils.StateFunction) Test(org.junit.Test)

Example 38 with SpacecraftState

use of org.orekit.propagation.SpacecraftState in project Orekit by CS-SI.

the class AngularRaDecTest method testParameterDerivatives.

@Test
public void testParameterDerivatives() throws OrekitException {
    Context context = EstimationTestUtils.geoStationnaryContext("regular-data:potential:tides");
    final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(OrbitType.EQUINOCTIAL, PositionAngle.TRUE, false, 1.0e-6, 60.0, 0.001);
    // create perfect azimuth-elevation measurements
    for (final GroundStation station : context.stations) {
        station.getEastOffsetDriver().setSelected(true);
        station.getNorthOffsetDriver().setSelected(true);
        station.getZenithOffsetDriver().setSelected(true);
    }
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new AngularRaDecMeasurementCreator(context), 0.25, 3.0, 600.0);
    propagator.setSlaveMode();
    for (final ObservedMeasurement<?> measurement : measurements) {
        // parameter corresponding to station position offset
        final GroundStation stationParameter = ((AngularRaDec) measurement).getStation();
        // 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
        // angular would have depended only on the current position but
        // not on the current velocity.
        final AbsoluteDate datemeas = measurement.getDate();
        final SpacecraftState stateini = propagator.propagate(datemeas);
        final Vector3D stationP = stationParameter.getOffsetToInertial(stateini.getFrame(), datemeas).transformPosition(Vector3D.ZERO);
        final double meanDelay = AbstractMeasurement.signalTimeOfFlight(stateini.getPVCoordinates(), stationP, datemeas);
        final AbsoluteDate date = measurement.getDate().shiftedBy(-0.75 * meanDelay);
        final SpacecraftState state = propagator.propagate(date);
        final ParameterDriver[] drivers = new ParameterDriver[] { stationParameter.getEastOffsetDriver(), stationParameter.getNorthOffsetDriver(), stationParameter.getZenithOffsetDriver() };
        for (int i = 0; i < 3; ++i) {
            final double[] gradient = measurement.estimate(0, 0, new SpacecraftState[] { state }).getParameterDerivatives(drivers[i]);
            Assert.assertEquals(2, measurement.getDimension());
            Assert.assertEquals(2, gradient.length);
            for (final int k : new int[] { 0, 1 }) {
                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()[k];
                    }
                }, drivers[i], 3, 50.0);
                final double ref = dMkdP.value(drivers[i]);
                if (ref > 1.e-12) {
                    Assert.assertEquals(ref, gradient[k], 3e-9 * FastMath.abs(ref));
                }
            }
        }
    }
}
Also used : Context(org.orekit.estimation.Context) ParameterDriver(org.orekit.utils.ParameterDriver) AbsoluteDate(org.orekit.time.AbsoluteDate) SpacecraftState(org.orekit.propagation.SpacecraftState) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) ParameterFunction(org.orekit.utils.ParameterFunction) Propagator(org.orekit.propagation.Propagator) OrekitException(org.orekit.errors.OrekitException) Test(org.junit.Test)

Example 39 with SpacecraftState

use of org.orekit.propagation.SpacecraftState in project Orekit by CS-SI.

the class InterSatellitesRangeTest method genericTestStateDerivatives.

void genericTestStateDerivatives(final boolean printResults, final int index, 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 inter-satellites range measurements
    final TimeStampedPVCoordinates original = context.initialOrbit.getPVCoordinates();
    final Orbit closeOrbit = new CartesianOrbit(new TimeStampedPVCoordinates(context.initialOrbit.getDate(), original.getPosition().add(new Vector3D(1000, 2000, 3000)), original.getVelocity().add(new Vector3D(-0.03, 0.01, 0.02))), context.initialOrbit.getFrame(), context.initialOrbit.getMu());
    final Propagator closePropagator = EstimationTestUtils.createPropagator(closeOrbit, propagatorBuilder);
    closePropagator.setEphemerisMode();
    closePropagator.propagate(context.initialOrbit.getDate().shiftedBy(3.5 * closeOrbit.getKeplerianPeriod()));
    final BoundedPropagator ephemeris = closePropagator.getGeneratedEphemeris();
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new InterSatellitesRangeMeasurementCreator(ephemeris), 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.)) {
                // 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.
                final double meanDelay = measurement.getObservedValue()[0] / Constants.SPEED_OF_LIGHT;
                final AbsoluteDate date = measurement.getDate().shiftedBy(-0.75 * meanDelay);
                final SpacecraftState[] states = { interpolator.getInterpolatedState(date), ephemeris.propagate(date) };
                final double[][] jacobian = measurement.estimate(0, 0, states).getStateDerivatives(index);
                // 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 {
                        final SpacecraftState[] s = states.clone();
                        s[index] = state;
                        return measurement.estimate(0, 0, s).getEstimatedValue();
                    }
                }, measurement.getDimension(), propagator.getAttitudeProvider(), OrbitType.CARTESIAN, PositionAngle.TRUE, 2.0, 3).value(states[index]);
                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) {
                    System.out.format(Locale.US, "%-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", 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, "%-23s  %-23s  " + "%10s  %10s  %10s  " + "%10s  %10s  %10s  " + "%10s  %10s  %10s  " + "%10s  %10s  %10s%n", "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) CartesianOrbit(org.orekit.orbits.CartesianOrbit) Max(org.hipparchus.stat.descriptive.rank.Max) ArrayList(java.util.ArrayList) Median(org.hipparchus.stat.descriptive.rank.Median) TimeStampedPVCoordinates(org.orekit.utils.TimeStampedPVCoordinates) AbsoluteDate(org.orekit.time.AbsoluteDate) SpacecraftState(org.orekit.propagation.SpacecraftState) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) BoundedPropagator(org.orekit.propagation.BoundedPropagator) Propagator(org.orekit.propagation.Propagator) OrekitException(org.orekit.errors.OrekitException) BoundedPropagator(org.orekit.propagation.BoundedPropagator) Context(org.orekit.estimation.Context) Orbit(org.orekit.orbits.Orbit) CartesianOrbit(org.orekit.orbits.CartesianOrbit) OrekitStepInterpolator(org.orekit.propagation.sampling.OrekitStepInterpolator) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) StateFunction(org.orekit.utils.StateFunction) ChronologicalComparator(org.orekit.time.ChronologicalComparator)

Example 40 with SpacecraftState

use of org.orekit.propagation.SpacecraftState in project Orekit by CS-SI.

the class InterSatellitesRangeTest method genericTestValues.

/**
 * Generic test function for values of the inter-satellites 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 inter-satellites range measurements
    final TimeStampedPVCoordinates original = context.initialOrbit.getPVCoordinates();
    final Orbit closeOrbit = new CartesianOrbit(new TimeStampedPVCoordinates(context.initialOrbit.getDate(), original.getPosition().add(new Vector3D(1000, 2000, 3000)), original.getVelocity().add(new Vector3D(-0.03, 0.01, 0.02))), context.initialOrbit.getFrame(), context.initialOrbit.getMu());
    final Propagator closePropagator = EstimationTestUtils.createPropagator(closeOrbit, propagatorBuilder);
    closePropagator.setEphemerisMode();
    closePropagator.propagate(context.initialOrbit.getDate().shiftedBy(3.5 * closeOrbit.getKeplerianPeriod()));
    final BoundedPropagator ephemeris = closePropagator.getGeneratedEphemeris();
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new InterSatellitesRangeMeasurementCreator(ephemeris), 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.
                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 Range & errors
                final double RangeObserved = measurement.getObservedValue()[0];
                final EstimatedMeasurement<?> estimated = measurement.estimate(0, 0, new SpacecraftState[] { state, ephemeris.propagate(state.getDate()) });
                // the real state used for estimation is adjusted according to downlink delay
                double adjustment = state.getDate().durationFrom(estimated.getStates()[0].getDate());
                Assert.assertTrue(adjustment > 0.000006);
                Assert.assertTrue(adjustment < 0.0003);
                final double RangeEstimated = estimated.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();
                    System.out.format(Locale.US, "%-23s  %-23s  %19.6f  %19.6f  %13.6e  %13.6e%n", 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, "%-23s  %-23s  %19s  %19s  %13s  %13s%n", "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, 1.3e-7);
    Assert.assertEquals(0.0, absErrorsMin, 7.3e-7);
    Assert.assertEquals(0.0, absErrorsMax, 1.8e-7);
    Assert.assertEquals(0.0, relErrorsMedian, 1.0e-12);
    Assert.assertEquals(0.0, relErrorsMax, 3.2e-12);
}
Also used : CartesianOrbit(org.orekit.orbits.CartesianOrbit) Max(org.hipparchus.stat.descriptive.rank.Max) ArrayList(java.util.ArrayList) Median(org.hipparchus.stat.descriptive.rank.Median) TimeStampedPVCoordinates(org.orekit.utils.TimeStampedPVCoordinates) AbsoluteDate(org.orekit.time.AbsoluteDate) SpacecraftState(org.orekit.propagation.SpacecraftState) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) BoundedPropagator(org.orekit.propagation.BoundedPropagator) Propagator(org.orekit.propagation.Propagator) BoundedPropagator(org.orekit.propagation.BoundedPropagator) Context(org.orekit.estimation.Context) Orbit(org.orekit.orbits.Orbit) CartesianOrbit(org.orekit.orbits.CartesianOrbit) OrekitStepInterpolator(org.orekit.propagation.sampling.OrekitStepInterpolator) Min(org.hipparchus.stat.descriptive.rank.Min) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) ChronologicalComparator(org.orekit.time.ChronologicalComparator)

Aggregations

SpacecraftState (org.orekit.propagation.SpacecraftState)470 Test (org.junit.Test)324 AbsoluteDate (org.orekit.time.AbsoluteDate)280 KeplerianOrbit (org.orekit.orbits.KeplerianOrbit)178 Vector3D (org.hipparchus.geometry.euclidean.threed.Vector3D)153 FieldAbsoluteDate (org.orekit.time.FieldAbsoluteDate)138 Orbit (org.orekit.orbits.Orbit)131 FieldSpacecraftState (org.orekit.propagation.FieldSpacecraftState)127 PVCoordinates (org.orekit.utils.PVCoordinates)98 CartesianOrbit (org.orekit.orbits.CartesianOrbit)95 Propagator (org.orekit.propagation.Propagator)92 Frame (org.orekit.frames.Frame)79 OrekitException (org.orekit.errors.OrekitException)74 EquinoctialOrbit (org.orekit.orbits.EquinoctialOrbit)74 NumericalPropagator (org.orekit.propagation.numerical.NumericalPropagator)74 DormandPrince853Integrator (org.hipparchus.ode.nonstiff.DormandPrince853Integrator)70 TimeStampedPVCoordinates (org.orekit.utils.TimeStampedPVCoordinates)64 AbstractLegacyForceModelTest (org.orekit.forces.AbstractLegacyForceModelTest)61 CircularOrbit (org.orekit.orbits.CircularOrbit)58 FieldVector3D (org.hipparchus.geometry.euclidean.threed.FieldVector3D)57