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);
}
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;
}
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);
}
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));
}
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;
}
Aggregations