Search in sources :

Example 6 with FieldVector3D

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

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

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

the class HolmesFeatherstoneAttractionModel method accelerationWrtState.

/**
 * Compute acceleration derivatives with respect to state parameters.
 * <p>
 * From a theoretical point of view, this method computes the same values
 * as {@link #acceleration(FieldSpacecraftState, RealFieldElement[])} in the
 * specific case of {@link DerivativeStructure} with respect to state, so
 * it is less general. However, it is *much* faster in this important case.
 * <p>
 * <p>
 * The derivatives should be computed with respect to position. The input
 * parameters already take into account the free parameters (6 or 7 depending
 * on derivation with respect to mass being considered or not) and order
 * (always 1). Free parameters at indices 0, 1 and 2 correspond to derivatives
 * with respect to position. Free parameters at indices 3, 4 and 5 correspond
 * to derivatives with respect to velocity (these derivatives will remain zero
 * as acceleration due to gravity does not depend on velocity). Free parameter
 * at index 6 (if present) corresponds to to derivatives with respect to mass
 * (this derivative will remain zero as acceleration due to gravity does not
 * depend on mass).
 * </p>
 * @param date current date
 * @param frame inertial reference frame for state (both orbit and attitude)
 * @param position position of spacecraft in inertial frame
 * @param mu central attraction coefficient to use
 * @return acceleration with all derivatives specified by the input parameters
 * own derivatives
 * @exception OrekitException if derivatives cannot be computed
 * @since 6.0
 */
private FieldVector3D<DerivativeStructure> accelerationWrtState(final AbsoluteDate date, final Frame frame, final FieldVector3D<DerivativeStructure> position, final DerivativeStructure mu) throws OrekitException {
    // get the position in body frame
    final Transform fromBodyFrame = bodyFrame.getTransformTo(frame, date);
    final Transform toBodyFrame = fromBodyFrame.getInverse();
    final Vector3D positionBody = toBodyFrame.transformPosition(position.toVector3D());
    // compute gradient and Hessian
    final GradientHessian gh = gradientHessian(date, positionBody, mu.getReal());
    // gradient of the non-central part of the gravity field
    final double[] gInertial = fromBodyFrame.transformVector(new Vector3D(gh.getGradient())).toArray();
    // Hessian of the non-central part of the gravity field
    final RealMatrix hBody = new Array2DRowRealMatrix(gh.getHessian(), false);
    final RealMatrix rot = new Array2DRowRealMatrix(toBodyFrame.getRotation().getMatrix());
    final RealMatrix hInertial = rot.transpose().multiply(hBody).multiply(rot);
    // distribute all partial derivatives in a compact acceleration vector
    final double[] derivatives = new double[1 + position.getX().getFreeParameters()];
    final DerivativeStructure[] accDer = new DerivativeStructure[3];
    for (int i = 0; i < 3; ++i) {
        // first element is value of acceleration (i.e. gradient of field)
        derivatives[0] = gInertial[i];
        // next three elements are one row of the Jacobian of acceleration (i.e. Hessian of field)
        derivatives[1] = hInertial.getEntry(i, 0);
        derivatives[2] = hInertial.getEntry(i, 1);
        derivatives[3] = hInertial.getEntry(i, 2);
        // next element is derivative with respect to parameter mu
        if (derivatives.length > 4 && isVariable(mu, 3)) {
            derivatives[4] = gInertial[i] / mu.getReal();
        }
        accDer[i] = position.getX().getFactory().build(derivatives);
    }
    return new FieldVector3D<>(accDer);
}
Also used : RealMatrix(org.hipparchus.linear.RealMatrix) Array2DRowRealMatrix(org.hipparchus.linear.Array2DRowRealMatrix) Array2DRowRealMatrix(org.hipparchus.linear.Array2DRowRealMatrix) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) DerivativeStructure(org.hipparchus.analysis.differentiation.DerivativeStructure) Transform(org.orekit.frames.Transform) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D)

Example 9 with FieldVector3D

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

the class JB2008 method getDensity.

/**
 * {@inheritDoc}
 */
@Override
public <T extends RealFieldElement<T>> T getDensity(final FieldAbsoluteDate<T> date, final FieldVector3D<T> position, final Frame frame) throws OrekitException {
    // check if data are available :
    final AbsoluteDate dateD = date.toAbsoluteDate();
    if ((dateD.compareTo(inputParams.getMaxDate()) > 0) || (dateD.compareTo(inputParams.getMinDate()) < 0)) {
        throw new OrekitException(OrekitMessages.NO_SOLAR_ACTIVITY_AT_DATE, dateD, inputParams.getMinDate(), inputParams.getMaxDate());
    }
    // compute MJD date
    final T dateMJD = date.durationFrom(AbsoluteDate.MODIFIED_JULIAN_EPOCH).divide(Constants.JULIAN_DAY);
    // compute geodetic position (km and °)
    final FieldGeodeticPoint<T> inBody = earth.transform(position, frame, date);
    // compute sun position
    final Frame ecef = earth.getBodyFrame();
    final FieldVector3D<T> sunPos = new FieldVector3D<>(date.getField(), sun.getPVCoordinates(dateD, ecef).getPosition());
    final FieldGeodeticPoint<T> sunInBody = earth.transform(sunPos, ecef, date);
    return getDensity(dateMJD, sunInBody.getLongitude(), sunInBody.getLatitude(), inBody.getLongitude(), inBody.getLatitude(), inBody.getAltitude(), inputParams.getF10(dateD), inputParams.getF10B(dateD), inputParams.getS10(dateD), inputParams.getS10B(dateD), inputParams.getXM10(dateD), inputParams.getXM10B(dateD), inputParams.getY10(dateD), inputParams.getY10B(dateD), inputParams.getDSTDTC(dateD));
}
Also used : Frame(org.orekit.frames.Frame) OrekitException(org.orekit.errors.OrekitException) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) FieldAbsoluteDate(org.orekit.time.FieldAbsoluteDate) AbsoluteDate(org.orekit.time.AbsoluteDate)

Example 10 with FieldVector3D

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

the class DragForce method getDensityWrtStateUsingFiniteDifferences.

/**
 * Compute density and its derivatives.
 * Using finite differences for the derivatives.
 * And doing the actual computation only for the derivatives with respect to position (others are set to 0.).
 * <p>
 * From a theoretical point of view, this method computes the same values
 * as {@link Atmosphere#getDensity(FieldAbsoluteDate, FieldVector3D, Frame)} in the
 * specific case of {@link DerivativeStructure} with respect to state, so
 * it is less general. However, it is *much* faster in this important case.
 * <p>
 * <p>
 * The derivatives should be computed with respect to position. The input
 * parameters already take into account the free parameters (6, 7 or 8 depending
 * on derivation with respect to drag coefficient and lift ratio being considered or not)
 * and order (always 1). Free parameters at indices 0, 1 and 2 correspond to derivatives
 * with respect to position. Free parameters at indices 3, 4 and 5 correspond
 * to derivatives with respect to velocity (these derivatives will remain zero
 * as the atmospheric density does not depend on velocity). Free parameter
 * at indexes 6 and 7 (if present) corresponds to derivatives with respect to drag coefficient
 * and/or lift ratio (one of these or both).
 * This 2 last derivatives will remain zero as atmospheric density does not depend on them.
 * </p>
 * @param date current date
 * @param frame inertial reference frame for state (both orbit and attitude)
 * @param position position of spacecraft in inertial frame
 * @param <T> type of the elements
 * @return the density and its derivatives
 * @exception OrekitException if derivatives cannot be computed
 * @since 9.0
 */
private <T extends RealFieldElement<T>> T getDensityWrtStateUsingFiniteDifferences(final AbsoluteDate date, final Frame frame, final FieldVector3D<T> position) throws OrekitException {
    // Retrieve derivation properties for parameter T
    // It is implied here that T is a DerivativeStructure
    // With order 1 and 6, 7 or 8 free parameters
    // This is all checked before in method isStateDerivatives
    final DSFactory factory = ((DerivativeStructure) position.getX()).getFactory();
    // Build a DerivativeStructure using only derivatives with respect to position
    final DSFactory factory3 = new DSFactory(3, 1);
    final FieldVector3D<DerivativeStructure> position3 = new FieldVector3D<>(factory3.variable(0, position.getX().getReal()), factory3.variable(1, position.getY().getReal()), factory3.variable(2, position.getZ().getReal()));
    // Get atmosphere properties in atmosphere own frame
    final Frame atmFrame = atmosphere.getFrame();
    final Transform toBody = frame.getTransformTo(atmFrame, date);
    final FieldVector3D<DerivativeStructure> posBodyDS = toBody.transformPosition(position3);
    final Vector3D posBody = posBodyDS.toVector3D();
    // Estimate density model by finite differences and composition
    // Using a delta of 1m
    final double delta = 1.0;
    final double x = posBody.getX();
    final double y = posBody.getY();
    final double z = posBody.getZ();
    final double rho0 = atmosphere.getDensity(date, posBody, atmFrame);
    final double dRhodX = (atmosphere.getDensity(date, new Vector3D(x + delta, y, z), atmFrame) - rho0) / delta;
    final double dRhodY = (atmosphere.getDensity(date, new Vector3D(x, y + delta, z), atmFrame) - rho0) / delta;
    final double dRhodZ = (atmosphere.getDensity(date, new Vector3D(x, y, z + delta), atmFrame) - rho0) / delta;
    final double[] dXdQ = posBodyDS.getX().getAllDerivatives();
    final double[] dYdQ = posBodyDS.getY().getAllDerivatives();
    final double[] dZdQ = posBodyDS.getZ().getAllDerivatives();
    // Density with derivatives:
    // - The value and only the 3 first derivatives (those with respect to spacecraft position) are computed
    // - Others are set to 0.
    final int p = factory.getCompiler().getFreeParameters();
    final double[] rhoAll = new double[p + 1];
    rhoAll[0] = rho0;
    for (int i = 1; i < 4; ++i) {
        rhoAll[i] = dRhodX * dXdQ[i] + dRhodY * dYdQ[i] + dRhodZ * dZdQ[i];
    }
    @SuppressWarnings("unchecked") final T rho = (T) (factory.build(rhoAll));
    return rho;
}
Also used : Frame(org.orekit.frames.Frame) Vector3D(org.hipparchus.geometry.euclidean.threed.Vector3D) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D) DerivativeStructure(org.hipparchus.analysis.differentiation.DerivativeStructure) DSFactory(org.hipparchus.analysis.differentiation.DSFactory) Transform(org.orekit.frames.Transform) FieldVector3D(org.hipparchus.geometry.euclidean.threed.FieldVector3D)

Aggregations

FieldVector3D (org.hipparchus.geometry.euclidean.threed.FieldVector3D)124 Vector3D (org.hipparchus.geometry.euclidean.threed.Vector3D)53 FieldAbsoluteDate (org.orekit.time.FieldAbsoluteDate)49 DerivativeStructure (org.hipparchus.analysis.differentiation.DerivativeStructure)38 Test (org.junit.Test)38 Frame (org.orekit.frames.Frame)36 TimeStampedFieldPVCoordinates (org.orekit.utils.TimeStampedFieldPVCoordinates)31 OrekitException (org.orekit.errors.OrekitException)23 DSFactory (org.hipparchus.analysis.differentiation.DSFactory)20 FieldPVCoordinates (org.orekit.utils.FieldPVCoordinates)20 Decimal64 (org.hipparchus.util.Decimal64)18 FieldEquinoctialOrbit (org.orekit.orbits.FieldEquinoctialOrbit)15 OrbitType (org.orekit.orbits.OrbitType)15 AbsoluteDate (org.orekit.time.AbsoluteDate)15 DormandPrince853FieldIntegrator (org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator)14 Transform (org.orekit.frames.Transform)14 FieldDerivativeStructure (org.hipparchus.analysis.differentiation.FieldDerivativeStructure)12 FieldEcksteinHechlerPropagator (org.orekit.propagation.analytical.FieldEcksteinHechlerPropagator)10 TimeStampedPVCoordinates (org.orekit.utils.TimeStampedPVCoordinates)9 ArrayList (java.util.ArrayList)8