Search in sources :

Example 46 with Vector3D

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

the class AbstractLegacyForceModelTest method checkStateJacobianVs80Implementation.

protected void checkStateJacobianVs80Implementation(final SpacecraftState state, final ForceModel forceModel, final AttitudeProvider attitudeProvider, final double checkTolerance, final boolean print) throws OrekitException {
    FieldSpacecraftState<DerivativeStructure> fState = toDS(state, attitudeProvider);
    FieldVector3D<DerivativeStructure> dsNew = forceModel.acceleration(fState, forceModel.getParameters(fState.getDate().getField()));
    FieldVector3D<DerivativeStructure> dsOld = accelerationDerivatives(forceModel, fState.getDate().toAbsoluteDate(), fState.getFrame(), fState.getPVCoordinates().getPosition(), fState.getPVCoordinates().getVelocity(), fState.getAttitude().getRotation(), fState.getMass());
    Vector3D dFdPXRef = new Vector3D(dsOld.getX().getPartialDerivative(1, 0, 0, 0, 0, 0), dsOld.getY().getPartialDerivative(1, 0, 0, 0, 0, 0), dsOld.getZ().getPartialDerivative(1, 0, 0, 0, 0, 0));
    Vector3D dFdPXRes = new Vector3D(dsNew.getX().getPartialDerivative(1, 0, 0, 0, 0, 0), dsNew.getY().getPartialDerivative(1, 0, 0, 0, 0, 0), dsNew.getZ().getPartialDerivative(1, 0, 0, 0, 0, 0));
    Vector3D dFdPYRef = new Vector3D(dsOld.getX().getPartialDerivative(0, 1, 0, 0, 0, 0), dsOld.getY().getPartialDerivative(0, 1, 0, 0, 0, 0), dsOld.getZ().getPartialDerivative(0, 1, 0, 0, 0, 0));
    Vector3D dFdPYRes = new Vector3D(dsNew.getX().getPartialDerivative(0, 1, 0, 0, 0, 0), dsNew.getY().getPartialDerivative(0, 1, 0, 0, 0, 0), dsNew.getZ().getPartialDerivative(0, 1, 0, 0, 0, 0));
    Vector3D dFdPZRef = new Vector3D(dsOld.getX().getPartialDerivative(0, 0, 1, 0, 0, 0), dsOld.getY().getPartialDerivative(0, 0, 1, 0, 0, 0), dsOld.getZ().getPartialDerivative(0, 0, 1, 0, 0, 0));
    Vector3D dFdPZRes = new Vector3D(dsNew.getX().getPartialDerivative(0, 0, 1, 0, 0, 0), dsNew.getY().getPartialDerivative(0, 0, 1, 0, 0, 0), dsNew.getZ().getPartialDerivative(0, 0, 1, 0, 0, 0));
    Vector3D dFdVXRef = new Vector3D(dsOld.getX().getPartialDerivative(0, 0, 0, 1, 0, 0), dsOld.getY().getPartialDerivative(0, 0, 0, 1, 0, 0), dsOld.getZ().getPartialDerivative(0, 0, 0, 1, 0, 0));
    Vector3D dFdVXRes = new Vector3D(dsNew.getX().getPartialDerivative(0, 0, 0, 1, 0, 0), dsNew.getY().getPartialDerivative(0, 0, 0, 1, 0, 0), dsNew.getZ().getPartialDerivative(0, 0, 0, 1, 0, 0));
    Vector3D dFdVYRef = new Vector3D(dsOld.getX().getPartialDerivative(0, 0, 0, 0, 1, 0), dsOld.getY().getPartialDerivative(0, 0, 0, 0, 1, 0), dsOld.getZ().getPartialDerivative(0, 0, 0, 0, 1, 0));
    Vector3D dFdVYRes = new Vector3D(dsNew.getX().getPartialDerivative(0, 0, 0, 0, 1, 0), dsNew.getY().getPartialDerivative(0, 0, 0, 0, 1, 0), dsNew.getZ().getPartialDerivative(0, 0, 0, 0, 1, 0));
    Vector3D dFdVZRef = new Vector3D(dsOld.getX().getPartialDerivative(0, 0, 0, 0, 0, 1), dsOld.getY().getPartialDerivative(0, 0, 0, 0, 0, 1), dsOld.getZ().getPartialDerivative(0, 0, 0, 0, 0, 1));
    Vector3D dFdVZRes = new Vector3D(dsNew.getX().getPartialDerivative(0, 0, 0, 0, 0, 1), dsNew.getY().getPartialDerivative(0, 0, 0, 0, 0, 1), dsNew.getZ().getPartialDerivative(0, 0, 0, 0, 0, 1));
    if (print) {
        System.out.println("dF/dPX ref norm: " + dFdPXRef.getNorm() + ", abs error: " + Vector3D.distance(dFdPXRef, dFdPXRes) + ", rel error: " + (Vector3D.distance(dFdPXRef, dFdPXRes) / dFdPXRef.getNorm()));
        System.out.println("dF/dPY ref norm: " + dFdPYRef.getNorm() + ", abs error: " + Vector3D.distance(dFdPYRef, dFdPYRes) + ", rel error: " + (Vector3D.distance(dFdPYRef, dFdPYRes) / dFdPYRef.getNorm()));
        System.out.println("dF/dPZ ref norm: " + dFdPZRef.getNorm() + ", abs error: " + Vector3D.distance(dFdPZRef, dFdPZRes) + ", rel error: " + (Vector3D.distance(dFdPZRef, dFdPZRes) / dFdPZRef.getNorm()));
        System.out.println("dF/dVX ref norm: " + dFdVXRef.getNorm() + ", abs error: " + Vector3D.distance(dFdVXRef, dFdVXRes) + ", rel error: " + (Vector3D.distance(dFdVXRef, dFdVXRes) / dFdVXRef.getNorm()));
        System.out.println("dF/dVY ref norm: " + dFdVYRef.getNorm() + ", abs error: " + Vector3D.distance(dFdVYRef, dFdVYRes) + ", rel error: " + (Vector3D.distance(dFdVYRef, dFdVYRes) / dFdVYRef.getNorm()));
        System.out.println("dF/dVZ ref norm: " + dFdVZRef.getNorm() + ", abs error: " + Vector3D.distance(dFdVZRef, dFdVZRes) + ", rel error: " + (Vector3D.distance(dFdVZRef, dFdVZRes) / dFdVZRef.getNorm()));
    }
    checkdFdP(dFdPXRef, dFdPXRes, checkTolerance);
    checkdFdP(dFdPYRef, dFdPYRes, checkTolerance);
    checkdFdP(dFdPZRef, dFdPZRes, checkTolerance);
    checkdFdP(dFdVXRef, dFdVXRes, checkTolerance);
    checkdFdP(dFdVYRef, dFdVYRes, checkTolerance);
    checkdFdP(dFdVZRef, dFdVZRes, checkTolerance);
}
Also used : Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) DerivativeStructure(org.hipparchus.analysis.differentiation.DerivativeStructure)

Example 47 with Vector3D

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

the class BoxAndSolarArraySpacecraftTest method testNormalSunAlignedField.

@Test
public void testNormalSunAlignedField() throws OrekitException {
    BoxAndSolarArraySpacecraft s = new BoxAndSolarArraySpacecraft(0, 0, 0, (date, frame) -> new TimeStampedPVCoordinates(date, new Vector3D(0, 1e6, 0), Vector3D.ZERO), 20.0, Vector3D.PLUS_J, 0.0, 1.0, 0.0);
    Field<Decimal64> field = Decimal64Field.getInstance();
    FieldVector3D<Decimal64> normal = s.getNormal(FieldAbsoluteDate.getJ2000Epoch(field), FramesFactory.getEME2000(), FieldVector3D.getZero(field), FieldRotation.getIdentity(field));
    Assert.assertEquals(0, FieldVector3D.dotProduct(normal, Vector3D.PLUS_J).getReal(), 1.0e-16);
}
Also used : FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) Decimal64(org.hipparchus.util.Decimal64) TimeStampedPVCoordinates(org.orekit.utils.TimeStampedPVCoordinates) Test(org.junit.Test)

Example 48 with Vector3D

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

the class BoxAndSolarArraySpacecraftTest method testBackwardIllumination.

@Test
public void testBackwardIllumination() throws OrekitException {
    SpacecraftState state = propagator.getInitialState();
    CelestialBody sun = CelestialBodyFactory.getSun();
    BoxAndSolarArraySpacecraft s = new BoxAndSolarArraySpacecraft(0, 0, 0, sun, 20.0, Vector3D.PLUS_J, 0.0, 1.0, 0.0);
    Vector3D n = s.getNormal(state.getDate(), state.getFrame(), state.getPVCoordinates().getPosition(), state.getAttitude().getRotation());
    FieldVector3D<DerivativeStructure> aPlus = s.radiationPressureAcceleration(state.getDate(), state.getFrame(), state.getPVCoordinates().getPosition(), state.getAttitude().getRotation(), state.getMass(), n, getRadiationParameters(s), RadiationSensitive.ABSORPTION_COEFFICIENT);
    FieldVector3D<DerivativeStructure> aMinus = s.radiationPressureAcceleration(state.getDate(), state.getFrame(), state.getPVCoordinates().getPosition(), state.getAttitude().getRotation(), state.getMass(), n.negate(), getRadiationParameters(s), RadiationSensitive.ABSORPTION_COEFFICIENT);
    Assert.assertEquals(0.0, aPlus.add(aMinus).getNorm().getReal(), Double.MIN_VALUE);
}
Also used : SpacecraftState(org.orekit.propagation.SpacecraftState) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) CelestialBody(org.orekit.bodies.CelestialBody) DerivativeStructure(org.hipparchus.analysis.differentiation.DerivativeStructure) Test(org.junit.Test)

Example 49 with Vector3D

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

the class EstimationTestUtils method geoStationnaryContext.

public static Context geoStationnaryContext(final String dataRoot) throws OrekitException {
    Utils.setDataRoot(dataRoot);
    Context context = new Context();
    context.conventions = IERSConventions.IERS_2010;
    context.utc = TimeScalesFactory.getUTC();
    context.ut1 = TimeScalesFactory.getUT1(context.conventions, true);
    context.displacements = new StationDisplacement[0];
    String Myframename = "MyEarthFrame";
    final AbsoluteDate datedef = new AbsoluteDate(2000, 1, 1, 12, 0, 0.0, context.utc);
    final double omega = Constants.WGS84_EARTH_ANGULAR_VELOCITY;
    final Vector3D rotationRate = new Vector3D(0.0, 0.0, omega);
    TransformProvider MyEarthFrame = new TransformProvider() {

        private static final long serialVersionUID = 1L;

        public Transform getTransform(final AbsoluteDate date) {
            final double rotationduration = date.durationFrom(datedef);
            final Vector3D alpharot = new Vector3D(rotationduration, rotationRate);
            final Rotation rotation = new Rotation(Vector3D.PLUS_K, -alpharot.getZ(), RotationConvention.VECTOR_OPERATOR);
            return new Transform(date, rotation, rotationRate);
        }

        public <T extends RealFieldElement<T>> FieldTransform<T> getTransform(final FieldAbsoluteDate<T> date) {
            final T rotationduration = date.durationFrom(datedef);
            final FieldVector3D<T> alpharot = new FieldVector3D<>(rotationduration, rotationRate);
            final FieldRotation<T> rotation = new FieldRotation<>(FieldVector3D.getPlusK(date.getField()), alpharot.getZ().negate(), RotationConvention.VECTOR_OPERATOR);
            return new FieldTransform<>(date, rotation, new FieldVector3D<>(date.getField(), rotationRate));
        }
    };
    Frame FrameTest = new Frame(FramesFactory.getEME2000(), MyEarthFrame, Myframename, true);
    // Earth is spherical, rotating in one sidereal day
    context.earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 0.0, FrameTest);
    context.sun = CelestialBodyFactory.getSun();
    context.moon = CelestialBodyFactory.getMoon();
    context.radiationSensitive = new IsotropicRadiationClassicalConvention(2.0, 0.2, 0.8);
    context.dragSensitive = new IsotropicDrag(2.0, 1.2);
    GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
    AstronomicalAmplitudeReader aaReader = new AstronomicalAmplitudeReader("hf-fes2004.dat", 5, 2, 3, 1.0);
    DataProvidersManager.getInstance().feed(aaReader.getSupportedNames(), aaReader);
    Map<Integer, Double> map = aaReader.getAstronomicalAmplitudesMap();
    GravityFieldFactory.addOceanTidesReader(new FESCHatEpsilonReader("fes2004-7x7.dat", 0.01, FastMath.toRadians(1.0), OceanLoadDeformationCoefficients.IERS_2010, map));
    context.gravity = GravityFieldFactory.getNormalizedProvider(20, 20);
    // semimajor axis for a geostationnary satellite
    double da = FastMath.cbrt(context.gravity.getMu() / (omega * omega));
    // context.stations = Arrays.asList(context.createStation(  0.0,  0.0, 0.0, "Lat0_Long0"),
    // context.createStation( 62.29639,   -7.01250,  880.0, "Slættaratindur")
    // );
    context.stations = Arrays.asList(context.createStation(0.0, 0.0, 0.0, "Lat0_Long0"));
    // Station position & velocity in EME2000
    final Vector3D geovelocity = new Vector3D(0., 0., 0.);
    // Compute the frames transformation from station frame to EME2000
    Transform topoToEME = context.stations.get(0).getBaseFrame().getTransformTo(FramesFactory.getEME2000(), new AbsoluteDate(2000, 1, 1, 12, 0, 0.0, context.utc));
    // Station position in EME2000 at reference date
    Vector3D stationPositionEME = topoToEME.transformPosition(Vector3D.ZERO);
    // Satellite position and velocity in Station Frame
    final Vector3D sat_pos = new Vector3D(0., 0., da - stationPositionEME.getNorm());
    final Vector3D acceleration = new Vector3D(-context.gravity.getMu(), sat_pos);
    final PVCoordinates pv_sat_topo = new PVCoordinates(sat_pos, geovelocity, acceleration);
    // satellite position in EME2000
    final PVCoordinates pv_sat_iner = topoToEME.transformPVCoordinates(pv_sat_topo);
    // Geo-stationary Satellite Orbit, tightly above the station (l0-L0)
    context.initialOrbit = new KeplerianOrbit(pv_sat_iner, FramesFactory.getEME2000(), new AbsoluteDate(2000, 1, 1, 12, 0, 0.0, context.utc), context.gravity.getMu());
    context.stations = Arrays.asList(context.createStation(10.0, 45.0, 0.0, "Lat10_Long45"));
    // Turn-around range stations
    // Map entry = master station
    // Map value = slave station associated
    context.TARstations = new HashMap<GroundStation, GroundStation>();
    context.TARstations.put(context.createStation(41.977, 13.600, 671.354, "Fucino"), context.createStation(43.604, 1.444, 263.0, "Toulouse"));
    context.TARstations.put(context.createStation(49.867, 8.65, 144.0, "Darmstadt"), context.createStation(-25.885, 27.707, 1566.633, "Pretoria"));
    return context;
}
Also used : Frame(org.orekit.frames.Frame) OneAxisEllipsoid(org.orekit.bodies.OneAxisEllipsoid) IsotropicDrag(org.orekit.forces.drag.IsotropicDrag) PVCoordinates(org.orekit.utils.PVCoordinates) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) AbsoluteDate(org.orekit.time.AbsoluteDate) GRGSFormatReader(org.orekit.forces.gravity.potential.GRGSFormatReader) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) TransformProvider(org.orekit.frames.TransformProvider) KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) AstronomicalAmplitudeReader(org.orekit.forces.gravity.potential.AstronomicalAmplitudeReader) GroundStation(org.orekit.estimation.measurements.GroundStation) RealFieldElement(org.hipparchus.RealFieldElement) FieldTransform(org.orekit.frames.FieldTransform) Rotation(org.hipparchus.geometry.euclidean.threed.Rotation) FieldRotation(org.hipparchus.geometry.euclidean.threed.FieldRotation) FieldRotation(org.hipparchus.geometry.euclidean.threed.FieldRotation) FESCHatEpsilonReader(org.orekit.forces.gravity.potential.FESCHatEpsilonReader) IsotropicRadiationClassicalConvention(org.orekit.forces.radiation.IsotropicRadiationClassicalConvention) FieldTransform(org.orekit.frames.FieldTransform) Transform(org.orekit.frames.Transform) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate)

Example 50 with Vector3D

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

the class EstimationTestUtils method checkKalmanFit.

/**
 * Checker for Kalman estimator validation
 * @param context context used for the test
 * @param kalman Kalman filter
 * @param measurements List of observed measurements to be processed by the Kalman
 * @param refOrbit Reference orbit at last measurement date
 * @param expectedDeltaPos Expected position difference between estimated orbit and reference orbit
 * @param posEps Tolerance on expected position difference
 * @param expectedDeltaVel Expected velocity difference between estimated orbit and reference orbit
 * @param velEps Tolerance on expected velocity difference
 * @param expectedSigmasPos Expected values for covariance matrix on position
 * @param sigmaPosEps Tolerance on expected covariance matrix on position
 * @param expectedSigmasVel Expected values for covariance matrix on velocity
 * @param sigmaVelEps Tolerance on expected covariance matrix on velocity
 * @throws OrekitException
 */
public static void checkKalmanFit(final Context context, final KalmanEstimator kalman, final List<ObservedMeasurement<?>> measurements, final Orbit refOrbit, final PositionAngle positionAngle, final double expectedDeltaPos, final double posEps, final double expectedDeltaVel, final double velEps, final double[] expectedSigmasPos, final double sigmaPosEps, final double[] expectedSigmasVel, final double sigmaVelEps) throws OrekitException {
    // Add the measurements to the Kalman filter
    NumericalPropagator estimated = kalman.processMeasurements(measurements);
    // Check the number of measurements processed by the filter
    Assert.assertEquals(measurements.size(), kalman.getCurrentMeasurementNumber());
    // Get the last estimation
    final Orbit estimatedOrbit = estimated.getInitialState().getOrbit();
    final Vector3D estimatedPosition = estimatedOrbit.getPVCoordinates().getPosition();
    final Vector3D estimatedVelocity = estimatedOrbit.getPVCoordinates().getVelocity();
    // Get the last covariance matrix estimation
    final RealMatrix estimatedP = kalman.getPhysicalEstimatedCovarianceMatrix();
    // Convert the orbital part to Cartesian formalism
    // Assuming all 6 orbital parameters are estimated by the filter
    final double[][] dCdY = new double[6][6];
    estimatedOrbit.getJacobianWrtParameters(positionAngle, dCdY);
    final RealMatrix Jacobian = MatrixUtils.createRealMatrix(dCdY);
    final RealMatrix estimatedCartesianP = Jacobian.multiply(estimatedP.getSubMatrix(0, 5, 0, 5)).multiply(Jacobian.transpose());
    // Get the final sigmas (ie.sqrt of the diagonal of the Cartesian orbital covariance matrix)
    final double[] sigmas = new double[6];
    for (int i = 0; i < 6; i++) {
        sigmas[i] = FastMath.sqrt(estimatedCartesianP.getEntry(i, i));
    }
    // // FIXME: debug
    // final double dPos = Vector3D.distance(refOrbit.getPVCoordinates().getPosition(), estimatedPosition);
    // final double dVel = Vector3D.distance(refOrbit.getPVCoordinates().getVelocity(), estimatedVelocity);
    // System.out.println("Nb Meas = " + kalman.getCurrentMeasurementNumber());
    // System.out.println("dPos    = " + dPos + " m");
    // System.out.println("dVel    = " + dVel + " m/s");
    // System.out.println("sigmas  = " + sigmas[0] + " "
    // + sigmas[1] + " "
    // + sigmas[2] + " "
    // + sigmas[3] + " "
    // + sigmas[4] + " "
    // + sigmas[5]);
    // //debug
    // Check the final orbit estimation & PV sigmas
    Assert.assertEquals(expectedDeltaPos, Vector3D.distance(refOrbit.getPVCoordinates().getPosition(), estimatedPosition), posEps);
    Assert.assertEquals(expectedDeltaVel, Vector3D.distance(refOrbit.getPVCoordinates().getVelocity(), estimatedVelocity), velEps);
    for (int i = 0; i < 3; i++) {
        Assert.assertEquals(expectedSigmasPos[i], sigmas[i], sigmaPosEps);
        Assert.assertEquals(expectedSigmasVel[i], sigmas[i + 3], sigmaVelEps);
    }
}
Also used : KeplerianOrbit(org.orekit.orbits.KeplerianOrbit) Orbit(org.orekit.orbits.Orbit) RealMatrix(org.hipparchus.linear.RealMatrix) NumericalPropagator(org.orekit.propagation.numerical.NumericalPropagator) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D)

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