use of org.hipparchus.geometry.euclidean.threed.Vector3D in project Orekit by CS-SI.
the class BatchLSEstimatorTest method testMultiSatWithParameters.
/**
* A modified version of the previous test with a selection of propagation drivers to estimate
* One common (ยต)
* Some specifics for each satellite (Cr and Ca)
*
* @throws OrekitException
*/
@Test
public void testMultiSatWithParameters() throws OrekitException {
// Test: Set the propagator drivers to estimate for each satellite
final boolean muEstimated = true;
final boolean crEstimated1 = true;
final boolean caEstimated1 = true;
final boolean crEstimated2 = true;
final boolean caEstimated2 = false;
// Builder sat 1
final Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
final NumericalPropagatorBuilder propagatorBuilder1 = context.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 1.0, Force.POTENTIAL, Force.SOLAR_RADIATION_PRESSURE);
// Adding selection of parameters
String satName = "sat 1";
for (DelegatingDriver driver : propagatorBuilder1.getPropagationParametersDrivers().getDrivers()) {
if (driver.getName().equals("central attraction coefficient")) {
driver.setSelected(muEstimated);
}
if (driver.getName().equals(RadiationSensitive.REFLECTION_COEFFICIENT)) {
driver.setName(driver.getName() + " " + satName);
driver.setSelected(crEstimated1);
}
if (driver.getName().equals(RadiationSensitive.ABSORPTION_COEFFICIENT)) {
driver.setName(driver.getName() + " " + satName);
driver.setSelected(caEstimated1);
}
}
// Builder for sat 2
final Context context2 = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
final NumericalPropagatorBuilder propagatorBuilder2 = context2.createBuilder(OrbitType.KEPLERIAN, PositionAngle.TRUE, true, 1.0e-6, 60.0, 1.0, Force.POTENTIAL, Force.SOLAR_RADIATION_PRESSURE);
// Adding selection of parameters
satName = "sat 2";
for (ParameterDriver driver : propagatorBuilder2.getPropagationParametersDrivers().getDrivers()) {
if (driver.getName().equals("central attraction coefficient")) {
driver.setSelected(muEstimated);
}
if (driver.getName().equals(RadiationSensitive.REFLECTION_COEFFICIENT)) {
driver.setName(driver.getName() + " " + satName);
driver.setSelected(crEstimated2);
}
if (driver.getName().equals(RadiationSensitive.ABSORPTION_COEFFICIENT)) {
driver.setName(driver.getName() + " " + satName);
driver.setSelected(caEstimated2);
}
}
// 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, propagatorBuilder2);
closePropagator.setEphemerisMode();
closePropagator.propagate(context.initialOrbit.getDate().shiftedBy(3.5 * closeOrbit.getKeplerianPeriod()));
final BoundedPropagator ephemeris = closePropagator.getGeneratedEphemeris();
Propagator propagator1 = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder1);
final List<ObservedMeasurement<?>> r12 = EstimationTestUtils.createMeasurements(propagator1, new InterSatellitesRangeMeasurementCreator(ephemeris), 1.0, 3.0, 300.0);
// create perfect range measurements for first satellite
propagator1 = EstimationTestUtils.createPropagator(context.initialOrbit, propagatorBuilder1);
final List<ObservedMeasurement<?>> r1 = EstimationTestUtils.createMeasurements(propagator1, new RangeMeasurementCreator(context), 1.0, 3.0, 300.0);
// create orbit estimator
final BatchLSEstimator estimator = new BatchLSEstimator(new LevenbergMarquardtOptimizer(), propagatorBuilder1, propagatorBuilder2);
for (final ObservedMeasurement<?> interSat : r12) {
estimator.addMeasurement(interSat);
}
for (final ObservedMeasurement<?> range : r1) {
estimator.addMeasurement(range);
}
estimator.setParametersConvergenceThreshold(1.0e-2);
estimator.setMaxIterations(10);
estimator.setMaxEvaluations(20);
estimator.setObserver(new BatchLSObserver() {
int lastIter = 0;
int lastEval = 0;
/**
* {@inheritDoc}
*/
@Override
public void evaluationPerformed(int iterationsCount, int evaluationscount, Orbit[] orbits, ParameterDriversList estimatedOrbitalParameters, ParameterDriversList estimatedPropagatorParameters, ParameterDriversList estimatedMeasurementsParameters, EstimationsProvider evaluationsProvider, Evaluation lspEvaluation) throws OrekitException {
if (iterationsCount == lastIter) {
Assert.assertEquals(lastEval + 1, evaluationscount);
} else {
Assert.assertEquals(lastIter + 1, iterationsCount);
}
lastIter = iterationsCount;
lastEval = evaluationscount;
AbsoluteDate previous = AbsoluteDate.PAST_INFINITY;
for (int i = 0; i < evaluationsProvider.getNumber(); ++i) {
AbsoluteDate current = evaluationsProvider.getEstimatedMeasurement(i).getDate();
Assert.assertTrue(current.compareTo(previous) >= 0);
previous = current;
}
}
});
List<DelegatingDriver> parameters = estimator.getOrbitalParametersDrivers(true).getDrivers();
ParameterDriver a0Driver = parameters.get(0);
Assert.assertEquals("a[0]", a0Driver.getName());
a0Driver.setValue(a0Driver.getValue() + 1.2);
a0Driver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
ParameterDriver a1Driver = parameters.get(6);
Assert.assertEquals("a[1]", a1Driver.getName());
a1Driver.setValue(a1Driver.getValue() - 5.4);
a1Driver.setReferenceDate(AbsoluteDate.GALILEO_EPOCH);
final Orbit before = new KeplerianOrbit(parameters.get(6).getValue(), parameters.get(7).getValue(), parameters.get(8).getValue(), parameters.get(9).getValue(), parameters.get(10).getValue(), parameters.get(11).getValue(), PositionAngle.TRUE, closeOrbit.getFrame(), closeOrbit.getDate(), closeOrbit.getMu());
Assert.assertEquals(4.7246, Vector3D.distance(closeOrbit.getPVCoordinates().getPosition(), before.getPVCoordinates().getPosition()), 1.0e-3);
Assert.assertEquals(0.0010514, Vector3D.distance(closeOrbit.getPVCoordinates().getVelocity(), before.getPVCoordinates().getVelocity()), 1.0e-6);
EstimationTestUtils.checkFit(context, estimator, 4, 5, 0.0, 6.0e-06, 0.0, 1.7e-05, 0.0, 4.4e-07, 0.0, 1.7e-10);
final Orbit determined = new KeplerianOrbit(parameters.get(6).getValue(), parameters.get(7).getValue(), parameters.get(8).getValue(), parameters.get(9).getValue(), parameters.get(10).getValue(), parameters.get(11).getValue(), PositionAngle.TRUE, closeOrbit.getFrame(), closeOrbit.getDate(), closeOrbit.getMu());
Assert.assertEquals(0.0, Vector3D.distance(closeOrbit.getPVCoordinates().getPosition(), determined.getPVCoordinates().getPosition()), 5.8e-6);
Assert.assertEquals(0.0, Vector3D.distance(closeOrbit.getPVCoordinates().getVelocity(), determined.getPVCoordinates().getVelocity()), 3.5e-9);
// got a default one
for (final ParameterDriver driver : estimator.getOrbitalParametersDrivers(true).getDrivers()) {
if (driver.getName().startsWith("a[")) {
// user-specified reference date
Assert.assertEquals(0, driver.getReferenceDate().durationFrom(AbsoluteDate.GALILEO_EPOCH), 1.0e-15);
} else {
// default reference date
Assert.assertEquals(0, driver.getReferenceDate().durationFrom(propagatorBuilder1.getInitialOrbitDate()), 1.0e-15);
}
}
}
use of org.hipparchus.geometry.euclidean.threed.Vector3D in project Orekit by CS-SI.
the class PVTest method testPVWithCovarianceMatrix.
/**
* Test the PV constructor with one 6x6 covariance matrix as input.
* @throws OrekitException
*/
@Test
public void testPVWithCovarianceMatrix() throws OrekitException {
// Context
Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
// Dummy P, V, T
final Vector3D position = context.initialOrbit.getPVCoordinates().getPosition();
final Vector3D velocity = context.initialOrbit.getPVCoordinates().getVelocity();
final AbsoluteDate date = context.initialOrbit.getDate();
// Initialize standard deviations, weight and corr coeff
final double[] sigmaPV = { 10., 20., 30., 0.1, 0.2, 0.3 };
final double baseWeight = 0.5;
final double[][] corrCoefRef = new double[6][6];
for (int i = 0; i < 6; i++) {
for (int j = 0; j < 6; j++) {
if (i == j) {
corrCoefRef[i][i] = 1.;
} else {
corrCoefRef[i][j] = i + j + 1;
}
}
}
// Reference covariance matrix
final double[][] Pref = new double[6][6];
for (int i = 0; i < 6; i++) {
for (int j = 0; j < 6; j++) {
Pref[i][j] = corrCoefRef[i][j] * sigmaPV[i] * sigmaPV[j];
}
}
// Reference propagator numbers
final int[] propNumRef = { 0, 2 };
// Create PV measurements
final PV[] pvs = new PV[2];
pvs[0] = new PV(date, position, velocity, Pref, baseWeight);
pvs[1] = new PV(date, position, velocity, Pref, baseWeight, propNumRef[1]);
// Tolerance
// tolerance
final double eps = 1.8e-15;
// Check data
for (int k = 0; k < pvs.length; k++) {
final PV pv = pvs[k];
// Propagator numbers
assertEquals(propNumRef[k], pv.getPropagatorsIndices().get(0), eps);
// Weights
for (int i = 0; i < 6; i++) {
assertEquals(baseWeight, pv.getBaseWeight()[i], eps);
}
// Sigmas
for (int i = 0; i < 6; i++) {
assertEquals(sigmaPV[i], pv.getTheoreticalStandardDeviation()[i], eps);
}
// Covariances
final double[][] P = pv.getCovarianceMatrix();
// Substract with ref and get the norm
final double normP = MatrixUtils.createRealMatrix(P).subtract(MatrixUtils.createRealMatrix(Pref)).getNorm();
assertEquals(0., normP, eps);
// Correlation coef
final double[][] corrCoef = pv.getCorrelationCoefficientsMatrix();
// Substract with ref and get the norm
final double normCorrCoef = MatrixUtils.createRealMatrix(corrCoef).subtract(MatrixUtils.createRealMatrix(corrCoefRef)).getNorm();
assertEquals(0., normCorrCoef, eps);
}
}
use of org.hipparchus.geometry.euclidean.threed.Vector3D in project Orekit by CS-SI.
the class PVTest method testExceptions.
/**
* Test exceptions raised if the covariance matrix does not have the proper size.
*/
@Test
public void testExceptions() throws OrekitException {
// Context
Context context = EstimationTestUtils.eccentricContext("regular-data:potential:tides");
// Dummy P, V, T
final Vector3D position = context.initialOrbit.getPVCoordinates().getPosition();
final Vector3D velocity = context.initialOrbit.getPVCoordinates().getVelocity();
final AbsoluteDate date = context.initialOrbit.getDate();
final double weight = 1.;
// Build with two 3-sized vectors
try {
new PV(date, position, velocity, new double[] { 0., 0., 0. }, new double[] { 1. }, weight);
Assert.fail("An OrekitException should have been thrown");
} catch (OrekitException e) {
// An exception should indeed be raised here
}
// Build with one 6-sized vector
try {
new PV(date, position, velocity, new double[] { 0., 0., 0. }, weight);
Assert.fail("An OrekitException should have been thrown");
} catch (OrekitException e) {
// An exception should indeed be raised here
}
// Build with two 3x3 matrices
try {
new PV(date, position, velocity, new double[][] { { 0., 0. }, { 0., 0. } }, new double[][] { { 0., 0. }, { 0., 0. } }, weight);
Assert.fail("An OrekitException should have been thrown");
} catch (OrekitException e) {
// An exception should indeed be raised here
}
// Build with one 6x6 matrix
try {
new PV(date, position, velocity, new double[][] { { 0., 0. }, { 0., 0. } }, weight);
Assert.fail("An OrekitException should have been thrown");
} catch (OrekitException e) {
// An exception should indeed be raised here
}
}
use of org.hipparchus.geometry.euclidean.threed.Vector3D in project Orekit by CS-SI.
the class RangeAnalytic method theoreticalEvaluationAnalytic.
/**
* Analytical version of the theoreticalEvaluation function in Range class
* The derivative structures are not used, an analytical computation is used instead.
* @param iteration current LS estimator iteration
* @param evaluation current LS estimator evaluation
* @param state spacecraft state. At measurement date on first iteration then close to emission date on further iterations
* @param interpolator Orekit step interpolator
* @return
* @throws OrekitException
*/
protected EstimatedMeasurement<Range> theoreticalEvaluationAnalytic(final int iteration, final int evaluation, final SpacecraftState state) throws OrekitException {
// Station attribute from parent Range class
final GroundStation groundStation = this.getStation();
// Station position at signal arrival
final AbsoluteDate downlinkDate = getDate();
final Transform topoToInertDownlink = groundStation.getOffsetToInertial(state.getFrame(), downlinkDate);
final TimeStampedPVCoordinates stationDownlink = topoToInertDownlink.transformPVCoordinates(new TimeStampedPVCoordinates(downlinkDate, PVCoordinates.ZERO));
// Take propagation time into account
// (if state has already been set up to pre-compensate propagation delay,
// we will have offset == downlinkDelay and transitState will be
// the same as state)
// Downlink time of flight
final double tauD = signalTimeOfFlight(state.getPVCoordinates(), stationDownlink.getPosition(), downlinkDate);
final double delta = downlinkDate.durationFrom(state.getDate());
final double dt = delta - tauD;
// Transit state position
final SpacecraftState transitState = state.shiftedBy(dt);
final AbsoluteDate transitDate = transitState.getDate();
final Vector3D transitP = transitState.getPVCoordinates().getPosition();
// Station position at transit state date
final Transform topoToInertAtTransitDate = groundStation.getOffsetToInertial(state.getFrame(), transitDate);
final TimeStampedPVCoordinates stationAtTransitDate = topoToInertAtTransitDate.transformPVCoordinates(new TimeStampedPVCoordinates(transitDate, PVCoordinates.ZERO));
// Uplink time of flight
final double tauU = signalTimeOfFlight(stationAtTransitDate, transitP, transitDate);
final double tau = tauD + tauU;
// Real date and position of station at signal departure
final AbsoluteDate uplinkDate = downlinkDate.shiftedBy(-tau);
final TimeStampedPVCoordinates stationUplink = topoToInertDownlink.shiftedBy(-tau).transformPVCoordinates(new TimeStampedPVCoordinates(uplinkDate, PVCoordinates.ZERO));
// Prepare the evaluation
final EstimatedMeasurement<Range> estimated = new EstimatedMeasurement<Range>(this, iteration, evaluation, new SpacecraftState[] { transitState }, new TimeStampedPVCoordinates[] { stationUplink, transitState.getPVCoordinates(), stationDownlink });
// Set range value
final double cOver2 = 0.5 * Constants.SPEED_OF_LIGHT;
estimated.setEstimatedValue(tau * cOver2);
// Partial derivatives with respect to state
// The formulas below take into account the fact the measurement is at fixed reception date.
// When spacecraft position is changed, the downlink delay is changed, and in order
// to still have the measurement arrive at exactly the same date on ground, we must
// take the spacecraft-station relative velocity into account.
final Vector3D v = state.getPVCoordinates().getVelocity();
final Vector3D qv = stationDownlink.getVelocity();
final Vector3D downInert = stationDownlink.getPosition().subtract(transitP);
final double dDown = Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT * tauD - Vector3D.dotProduct(downInert, v);
final Vector3D upInert = transitP.subtract(stationUplink.getPosition());
// test
// final double dUp = Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT * tauU -
// Vector3D.dotProduct(upInert, qv);
// test
final double dUp = Constants.SPEED_OF_LIGHT * Constants.SPEED_OF_LIGHT * tauU - Vector3D.dotProduct(upInert, stationUplink.getVelocity());
// derivatives of the downlink time of flight
final double dTauDdPx = -downInert.getX() / dDown;
final double dTauDdPy = -downInert.getY() / dDown;
final double dTauDdPz = -downInert.getZ() / dDown;
// Derivatives of the uplink time of flight
final double dTauUdPx = 1. / dUp * upInert.dotProduct(Vector3D.PLUS_I.add((qv.subtract(v)).scalarMultiply(dTauDdPx)));
final double dTauUdPy = 1. / dUp * upInert.dotProduct(Vector3D.PLUS_J.add((qv.subtract(v)).scalarMultiply(dTauDdPy)));
final double dTauUdPz = 1. / dUp * upInert.dotProduct(Vector3D.PLUS_K.add((qv.subtract(v)).scalarMultiply(dTauDdPz)));
// derivatives of the range measurement
final double dRdPx = (dTauDdPx + dTauUdPx) * cOver2;
final double dRdPy = (dTauDdPy + dTauUdPy) * cOver2;
final double dRdPz = (dTauDdPz + dTauUdPz) * cOver2;
estimated.setStateDerivatives(0, new double[] { dRdPx, dRdPy, dRdPz, dRdPx * dt, dRdPy * dt, dRdPz * dt });
if (groundStation.getEastOffsetDriver().isSelected() || groundStation.getNorthOffsetDriver().isSelected() || groundStation.getZenithOffsetDriver().isSelected()) {
// Downlink tme of flight derivatives / station position in topocentric frame
final AngularCoordinates ac = topoToInertDownlink.getAngular().revert();
// final Rotation rotTopoToInert = ac.getRotation();
final Vector3D omega = ac.getRotationRate();
// Inertial frame
final double dTauDdQIx = downInert.getX() / dDown;
final double dTauDdQIy = downInert.getY() / dDown;
final double dTauDdQIz = downInert.getZ() / dDown;
// Uplink tme of flight derivatives / station position in topocentric frame
// Inertial frame
final double dTauUdQIx = 1 / dUp * upInert.dotProduct(Vector3D.MINUS_I.add((qv.subtract(v)).scalarMultiply(dTauDdQIx)).subtract(Vector3D.PLUS_I.crossProduct(omega).scalarMultiply(tau)));
final double dTauUdQIy = 1 / dUp * upInert.dotProduct(Vector3D.MINUS_J.add((qv.subtract(v)).scalarMultiply(dTauDdQIy)).subtract(Vector3D.PLUS_J.crossProduct(omega).scalarMultiply(tau)));
final double dTauUdQIz = 1 / dUp * upInert.dotProduct(Vector3D.MINUS_K.add((qv.subtract(v)).scalarMultiply(dTauDdQIz)).subtract(Vector3D.PLUS_K.crossProduct(omega).scalarMultiply(tau)));
// Range partial derivatives
// with respect to station position in inertial frame
final Vector3D dRdQI = new Vector3D((dTauDdQIx + dTauUdQIx) * cOver2, (dTauDdQIy + dTauUdQIy) * cOver2, (dTauDdQIz + dTauUdQIz) * cOver2);
// convert to topocentric frame, as the station position
// offset parameter is expressed in this frame
final Vector3D dRdQT = ac.getRotation().applyTo(dRdQI);
if (groundStation.getEastOffsetDriver().isSelected()) {
estimated.setParameterDerivatives(groundStation.getEastOffsetDriver(), dRdQT.getX());
}
if (groundStation.getNorthOffsetDriver().isSelected()) {
estimated.setParameterDerivatives(groundStation.getNorthOffsetDriver(), dRdQT.getY());
}
if (groundStation.getZenithOffsetDriver().isSelected()) {
estimated.setParameterDerivatives(groundStation.getZenithOffsetDriver(), dRdQT.getZ());
}
}
return estimated;
}
use of org.hipparchus.geometry.euclidean.threed.Vector3D 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);
}
Aggregations