Search in sources :

Example 26 with Vector3D

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

the class HolmesFeatherstoneAttractionModelTest method testParameterDerivative.

@Test
public void testParameterDerivative() throws OrekitException {
    Utils.setDataRoot("regular-data:potential/grgs-format");
    GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
    // pos-vel (from a ZOOM ephemeris reference)
    final Vector3D pos = new Vector3D(6.46885878304673824e+06, -1.88050918456274318e+06, -1.32931592294715829e+04);
    final Vector3D vel = new Vector3D(2.14718074509906819e+03, 7.38239351251748485e+03, -1.14097953925384523e+01);
    final SpacecraftState state = new SpacecraftState(new CartesianOrbit(new PVCoordinates(pos, vel), FramesFactory.getGCRF(), new AbsoluteDate(2005, 3, 5, 0, 24, 0.0, TimeScalesFactory.getTAI()), GravityFieldFactory.getUnnormalizedProvider(1, 1).getMu()));
    final HolmesFeatherstoneAttractionModel holmesFeatherstoneModel = new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), GravityFieldFactory.getNormalizedProvider(20, 20));
    final String name = NewtonianAttraction.CENTRAL_ATTRACTION_COEFFICIENT;
    checkParameterDerivative(state, holmesFeatherstoneModel, name, 1.0e-5, 5.0e-12);
}
Also used : SpacecraftState(org.orekit.propagation.SpacecraftState) FieldSpacecraftState(org.orekit.propagation.FieldSpacecraftState) CartesianOrbit(org.orekit.orbits.CartesianOrbit) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) PVCoordinates(org.orekit.utils.PVCoordinates) FieldPVCoordinates(org.orekit.utils.FieldPVCoordinates) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) AbsoluteDate(org.orekit.time.AbsoluteDate) GRGSFormatReader(org.orekit.forces.gravity.potential.GRGSFormatReader) AbstractLegacyForceModelTest(org.orekit.forces.AbstractLegacyForceModelTest) Test(org.junit.Test)

Example 27 with Vector3D

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

the class HolmesFeatherstoneAttractionModelTest method testCompleteWithCunninghamReference.

@Test
@Deprecated
public void testCompleteWithCunninghamReference() throws OrekitException {
    Utils.setDataRoot("regular-data:potential/grgs-format");
    GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
    // initialization
    AbsoluteDate date = new AbsoluteDate(new DateComponents(2000, 07, 01), new TimeComponents(13, 59, 27.816), TimeScalesFactory.getUTC());
    double i = FastMath.toRadians(98.7);
    double omega = FastMath.toRadians(93.0);
    double OMEGA = FastMath.toRadians(15.0 * 22.5);
    Orbit orbit = new KeplerianOrbit(7201009.7124401, 1e-3, i, omega, OMEGA, 0, PositionAngle.MEAN, FramesFactory.getEME2000(), date, mu);
    double[][] tolerances = NumericalPropagator.tolerances(0.01, orbit, OrbitType.CARTESIAN);
    AbsoluteDate targetDate = date.shiftedBy(3 * Constants.JULIAN_DAY);
    propagator = new NumericalPropagator(new DormandPrince853Integrator(1.0e-3, 120, tolerances[0], tolerances[1]));
    propagator.setOrbitType(OrbitType.CARTESIAN);
    propagator.addForceModel(new HolmesFeatherstoneAttractionModel(itrf, GravityFieldFactory.getNormalizedProvider(69, 69)));
    propagator.setInitialState(new SpacecraftState(orbit));
    SpacecraftState hfOrb = propagator.propagate(targetDate);
    propagator.removeForceModels();
    propagator.addForceModel(new CunninghamAttractionModel(itrf, GravityFieldFactory.getUnnormalizedProvider(69, 69)));
    propagator.setInitialState(new SpacecraftState(orbit));
    SpacecraftState cOrb = propagator.propagate(targetDate);
    Vector3D dif = hfOrb.getPVCoordinates().getPosition().subtract(cOrb.getPVCoordinates().getPosition());
    Assert.assertEquals(0, dif.getNorm(), 4e-5);
}
Also used : EquinoctialOrbit(org.orekit.orbits.EquinoctialOrbit) CartesianOrbit(org.orekit.orbits.CartesianOrbit) FieldKeplerianOrbit(org.orekit.orbits.FieldKeplerianOrbit) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) Orbit(org.orekit.orbits.Orbit) DateComponents(org.orekit.time.DateComponents) TimeComponents(org.orekit.time.TimeComponents) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) AbsoluteDate(org.orekit.time.AbsoluteDate) GRGSFormatReader(org.orekit.forces.gravity.potential.GRGSFormatReader) SpacecraftState(org.orekit.propagation.SpacecraftState) FieldSpacecraftState(org.orekit.propagation.FieldSpacecraftState) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) FieldNumericalPropagator(org.orekit.propagation.numerical.FieldNumericalPropagator) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) FieldKeplerianOrbit(org.orekit.orbits.FieldKeplerianOrbit) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) DormandPrince853Integrator(org.hipparchus.ode.nonstiff.DormandPrince853Integrator) AbstractLegacyForceModelTest(org.orekit.forces.AbstractLegacyForceModelTest) Test(org.junit.Test)

Example 28 with Vector3D

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

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

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

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