Search in sources :

Example 31 with Vector3D

use of org.hipparchus.geometry.euclidean.threed.Vector3D in project Orekit by CS-SI.

the class GroundStationTest method testEstimateStationPosition.

@Test
public void testEstimateStationPosition() throws OrekitException, IOException, ClassNotFoundException {
    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);
    // move one station
    final RandomGenerator random = new Well19937a(0x4adbecfc743bda60l);
    final TopocentricFrame base = context.stations.get(0).getBaseFrame();
    final BodyShape parent = base.getParentShape();
    final Vector3D baseOrigin = parent.transform(base.getPoint());
    final Vector3D deltaTopo = new Vector3D(2 * random.nextDouble() - 1, 2 * random.nextDouble() - 1, 2 * random.nextDouble() - 1);
    final Transform topoToParent = base.getTransformTo(parent.getBodyFrame(), (AbsoluteDate) null);
    final Vector3D deltaParent = topoToParent.transformVector(deltaTopo);
    final String movedSuffix = "-moved";
    final GroundStation moved = new GroundStation(new TopocentricFrame(parent, parent.transform(baseOrigin.subtract(deltaParent), parent.getBodyFrame(), null), base.getName() + movedSuffix), context.ut1.getEOPHistory(), context.stations.get(0).getDisplacements());
    // create orbit estimator
    final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(), propagatorBuilder);
    for (final ObservedMeasurement<?> measurement : measurements) {
        final Range range = (Range) measurement;
        final String name = range.getStation().getBaseFrame().getName() + movedSuffix;
        if (moved.getBaseFrame().getName().equals(name)) {
            estimator.addMeasurement(new Range(moved, range.getDate(), range.getObservedValue()[0], range.getTheoreticalStandardDeviation()[0], range.getBaseWeight()[0]));
        } else {
            estimator.addMeasurement(range);
        }
    }
    estimator.setParametersConvergenceThreshold(1.0e-3);
    estimator.setMaxIterations(100);
    estimator.setMaxEvaluations(200);
    // we want to estimate station offsets
    moved.getEastOffsetDriver().setSelected(true);
    moved.getNorthOffsetDriver().setSelected(true);
    moved.getZenithOffsetDriver().setSelected(true);
    EstimationTestUtils.checkFit(context, estimator, 2, 3, 0.0, 5.6e-7, 0.0, 1.4e-6, 0.0, 4.8e-7, 0.0, 2.6e-10);
    Assert.assertEquals(deltaTopo.getX(), moved.getEastOffsetDriver().getValue(), 4.5e-7);
    Assert.assertEquals(deltaTopo.getY(), moved.getNorthOffsetDriver().getValue(), 6.2e-7);
    Assert.assertEquals(deltaTopo.getZ(), moved.getZenithOffsetDriver().getValue(), 2.6e-7);
    GeodeticPoint result = moved.getOffsetGeodeticPoint(null);
    GeodeticPoint reference = context.stations.get(0).getBaseFrame().getPoint();
    Assert.assertEquals(reference.getLatitude(), result.getLatitude(), 1.4e-14);
    Assert.assertEquals(reference.getLongitude(), result.getLongitude(), 2.9e-14);
    Assert.assertEquals(reference.getAltitude(), result.getAltitude(), 2.6e-7);
    RealMatrix normalizedCovariances = estimator.getOptimum().getCovariances(1.0e-10);
    RealMatrix physicalCovariances = estimator.getPhysicalCovariances(1.0e-10);
    Assert.assertEquals(9, normalizedCovariances.getRowDimension());
    Assert.assertEquals(9, normalizedCovariances.getColumnDimension());
    Assert.assertEquals(9, physicalCovariances.getRowDimension());
    Assert.assertEquals(9, physicalCovariances.getColumnDimension());
    Assert.assertEquals(0.55431, physicalCovariances.getEntry(6, 6), 1.0e-5);
    Assert.assertEquals(0.22694, physicalCovariances.getEntry(7, 7), 1.0e-5);
    Assert.assertEquals(0.13106, physicalCovariances.getEntry(8, 8), 1.0e-5);
    ByteArrayOutputStream bos = new ByteArrayOutputStream();
    ObjectOutputStream oos = new ObjectOutputStream(bos);
    oos.writeObject(moved.getEstimatedEarthFrame().getTransformProvider());
    Assert.assertTrue(bos.size() > 155000);
    Assert.assertTrue(bos.size() < 160000);
    ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
    ObjectInputStream ois = new ObjectInputStream(bis);
    EstimatedEarthFrameProvider deserialized = (EstimatedEarthFrameProvider) ois.readObject();
    Assert.assertEquals(moved.getPrimeMeridianOffsetDriver().getValue(), deserialized.getPrimeMeridianOffsetDriver().getValue(), 1.0e-15);
    Assert.assertEquals(moved.getPrimeMeridianDriftDriver().getValue(), deserialized.getPrimeMeridianDriftDriver().getValue(), 1.0e-15);
    Assert.assertEquals(moved.getPolarOffsetXDriver().getValue(), deserialized.getPolarOffsetXDriver().getValue(), 1.0e-15);
    Assert.assertEquals(moved.getPolarDriftXDriver().getValue(), deserialized.getPolarDriftXDriver().getValue(), 1.0e-15);
    Assert.assertEquals(moved.getPolarOffsetYDriver().getValue(), deserialized.getPolarOffsetYDriver().getValue(), 1.0e-15);
    Assert.assertEquals(moved.getPolarDriftYDriver().getValue(), deserialized.getPolarDriftYDriver().getValue(), 1.0e-15);
}
Also used : TopocentricFrame(org.orekit.frames.TopocentricFrame) Well19937a(org.hipparchus.random.Well19937a) ObjectOutputStream(java.io.ObjectOutputStream) BodyShape(org.orekit.bodies.BodyShape) RandomGenerator(org.hipparchus.random.RandomGenerator) BatchLSEstimator(org.orekit.estimation.leastsquares.BatchLSEstimator) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) Propagator(org.orekit.propagation.Propagator) GeodeticPoint(org.orekit.bodies.GeodeticPoint) Context(org.orekit.estimation.Context) ByteArrayOutputStream(java.io.ByteArrayOutputStream) LevenbergMarquardtOptimizer(org.hipparchus.optim.nonlinear.vector.leastsquares.LevenbergMarquardtOptimizer) RealMatrix(org.hipparchus.linear.RealMatrix) ByteArrayInputStream(java.io.ByteArrayInputStream) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) FieldTransform(org.orekit.frames.FieldTransform) Transform(org.orekit.frames.Transform) ObjectInputStream(java.io.ObjectInputStream) Test(org.junit.Test)

Example 32 with Vector3D

use of org.hipparchus.geometry.euclidean.threed.Vector3D 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 33 with Vector3D

use of org.hipparchus.geometry.euclidean.threed.Vector3D 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)

Example 34 with Vector3D

use of org.hipparchus.geometry.euclidean.threed.Vector3D in project Orekit by CS-SI.

the class KalmanEstimatorTest method testKeplerianRangeWithOnBoardAntennaOffset.

/**
 * Perfect range measurements with a biased start and an on-board antenna range offset
 * Keplerian formalism
 * @throws OrekitException
 */
@Test
public void testKeplerianRangeWithOnBoardAntennaOffset() throws OrekitException {
    // Create context
    Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
    // Create initial orbit and propagator builder
    final OrbitType orbitType = OrbitType.KEPLERIAN;
    final PositionAngle positionAngle = PositionAngle.TRUE;
    final boolean perfectStart = true;
    final double minStep = 1.e-6;
    final double maxStep = 60.;
    final double dP = 1.;
    final NumericalPropagatorBuilder propagatorBuilder = context.createBuilder(orbitType, positionAngle, perfectStart, minStep, maxStep, dP);
    propagatorBuilder.setAttitudeProvider(new LofOffset(propagatorBuilder.getFrame(), LOFType.LVLH));
    // Antenna phase center definition
    final Vector3D antennaPhaseCenter = new Vector3D(-1.2, 2.3, -0.7);
    // Create perfect range measurements with antenna offset
    final Propagator propagator = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder);
    final List<ObservedMeasurement<?>> measurements = EstimationTestUtils.createMeasurements(propagator, new RangeMeasurementCreator(context, antennaPhaseCenter), 1.0, 3.0, 300.0);
    // Add antenna offset to the measurements
    final OnBoardAntennaRangeModifier obaModifier = new OnBoardAntennaRangeModifier(antennaPhaseCenter);
    for (final ObservedMeasurement<?> range : measurements) {
        ((Range) range).addModifier(obaModifier);
    }
    // Reference propagator for estimation performances
    final NumericalPropagator referencePropagator = propagatorBuilder.buildPropagator(propagatorBuilder.getSelectedNormalizedParameters());
    // Reference position/velocity at last measurement date
    final Orbit refOrbit = referencePropagator.propagate(measurements.get(measurements.size() - 1).getDate()).getOrbit();
    // Change semi-major axis of 1.2m as in the batch test
    ParameterDriver aDriver = propagatorBuilder.getOrbitalParametersDrivers().getDrivers().get(0);
    aDriver.setValue(aDriver.getValue() + 1.2);
    aDriver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
    // Cartesian covariance matrix initialization
    // 100m on position / 1e-2m/s on velocity
    final RealMatrix cartesianP = MatrixUtils.createRealDiagonalMatrix(new double[] { 10., 10., 10., 1e-3, 1e-3, 1e-3 });
    // Jacobian of the orbital parameters w/r to Cartesian
    final Orbit initialOrbit = OrbitType.KEPLERIAN.convertType(context.initialOrbit);
    final double[][] dYdC = new double[6][6];
    initialOrbit.getJacobianWrtCartesian(PositionAngle.TRUE, dYdC);
    final RealMatrix Jac = MatrixUtils.createRealMatrix(dYdC);
    // Keplerian initial covariance matrix
    final RealMatrix initialP = Jac.multiply(cartesianP.multiply(Jac.transpose()));
    // Process noise matrix is set to 0 here
    RealMatrix Q = MatrixUtils.createRealMatrix(6, 6);
    // Build the Kalman filter
    final KalmanEstimatorBuilder kalmanBuilder = new KalmanEstimatorBuilder();
    kalmanBuilder.builder(propagatorBuilder);
    kalmanBuilder.estimatedMeasurementsParameters(new ParameterDriversList());
    kalmanBuilder.initialCovarianceMatrix(initialP);
    kalmanBuilder.processNoiseMatrixProvider(new ConstantProcessNoise(Q));
    final KalmanEstimator kalman = kalmanBuilder.build();
    // Filter the measurements and check the results
    final double expectedDeltaPos = 0.;
    final double posEps = 4.57e-3;
    final double expectedDeltaVel = 0.;
    final double velEps = 7.29e-6;
    final double[] expectedSigmasPos = { 1.105194, 0.930785, 1.254579 };
    final double sigmaPosEps = 1e-6;
    final double[] expectedSigmasVel = { 6.193718e-4, 4.088774e-4, 3.299135e-4 };
    final double sigmaVelEps = 1e-10;
    EstimationTestUtils.checkKalmanFit(context, kalman, measurements, refOrbit, positionAngle, expectedDeltaPos, posEps, expectedDeltaVel, velEps, expectedSigmasPos, sigmaPosEps, expectedSigmasVel, sigmaVelEps);
}
Also used : PositionAngle(org.orekit.orbits.PositionAngle) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) ParameterDriversList(org.orekit.utils.ParameterDriversList) Propagator(org.orekit.propagation.Propagator) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) LofOffset(org.orekit.attitudes.LofOffset) ObservedMeasurement(org.orekit.estimation.measurements.ObservedMeasurement) Context(org.orekit.estimation.Context) Orbit(org.orekit.orbits.Orbit) Range(org.orekit.estimation.measurements.Range) ParameterDriver(org.orekit.utils.ParameterDriver) OnBoardAntennaRangeModifier(org.orekit.estimation.measurements.modifiers.OnBoardAntennaRangeModifier) RealMatrix(org.hipparchus.linear.RealMatrix) NumericalPropagatorBuilder(org.orekit.propagation.conversion.NumericalPropagatorBuilder) OrbitType(org.orekit.orbits.OrbitType) RangeMeasurementCreator(org.orekit.estimation.measurements.RangeMeasurementCreator) Test(org.junit.Test)

Example 35 with Vector3D

use of org.hipparchus.geometry.euclidean.threed.Vector3D in project Orekit by CS-SI.

the class KalmanOrbitDeterminationTest 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) TimeStampedPVCoordinates(org.orekit.utils.TimeStampedPVCoordinates) PVCoordinates(org.orekit.utils.PVCoordinates) 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)

Aggregations

Vector3D (org.hipparchus.geometry.euclidean.threed.Vector3D)750 Test (org.junit.Test)466 AbsoluteDate (org.orekit.time.AbsoluteDate)323 PVCoordinates (org.orekit.utils.PVCoordinates)280 FieldVector3D (org.hipparchus.geometry.euclidean.threed.FieldVector3D)259 TimeStampedPVCoordinates (org.orekit.utils.TimeStampedPVCoordinates)187 SpacecraftState (org.orekit.propagation.SpacecraftState)152 FieldAbsoluteDate (org.orekit.time.FieldAbsoluteDate)124 Rotation (org.hipparchus.geometry.euclidean.threed.Rotation)119 Frame (org.orekit.frames.Frame)115 KeplerianOrbit (org.orekit.orbits.KeplerianOrbit)105 Orbit (org.orekit.orbits.Orbit)100 GeodeticPoint (org.orekit.bodies.GeodeticPoint)84 OrekitException (org.orekit.errors.OrekitException)83 CartesianOrbit (org.orekit.orbits.CartesianOrbit)75 EquinoctialOrbit (org.orekit.orbits.EquinoctialOrbit)68 DateComponents (org.orekit.time.DateComponents)67 Transform (org.orekit.frames.Transform)61 OneAxisEllipsoid (org.orekit.bodies.OneAxisEllipsoid)59 CircularOrbit (org.orekit.orbits.CircularOrbit)59