use of org.orekit.utils.FieldPVCoordinates in project Orekit by CS-SI.
the class FieldTransformTest method doTestJacobianP.
private <T extends RealFieldElement<T>> void doTestJacobianP(Field<T> field) {
// base directions for finite differences
@SuppressWarnings("unchecked") FieldPVCoordinates<T>[] directions = (FieldPVCoordinates<T>[]) Array.newInstance(FieldPVCoordinates.class, 3);
directions[0] = new FieldPVCoordinates<>(FieldVector3D.getPlusI(field), FieldVector3D.getZero(field), FieldVector3D.getZero(field));
directions[1] = new FieldPVCoordinates<>(FieldVector3D.getPlusJ(field), FieldVector3D.getZero(field), FieldVector3D.getZero(field));
directions[2] = new FieldPVCoordinates<>(FieldVector3D.getPlusK(field), FieldVector3D.getZero(field), FieldVector3D.getZero(field));
double h = 0.01;
RandomGenerator random = new Well19937a(0x47fd0d6809f4b173l);
for (int i = 0; i < 20; ++i) {
// generate a random transform
FieldTransform<T> combined = randomTransform(field, random);
// compute Jacobian
T[][] jacobian = MathArrays.buildArray(field, 9, 9);
for (int l = 0; l < jacobian.length; ++l) {
for (int c = 0; c < jacobian[l].length; ++c) {
jacobian[l][c] = field.getZero().add(l + 0.1 * c);
}
}
combined.getJacobian(CartesianDerivativesFilter.USE_P, jacobian);
for (int j = 0; j < 100; ++j) {
FieldPVCoordinates<T> pv0 = new FieldPVCoordinates<>(randomVector(field, 1e3, random), randomVector(field, 1.0, random), randomVector(field, 1.0e-3, random));
double epsilonP = 2.0e-12 * pv0.getPosition().getNorm().getReal();
for (int c = 0; c < directions.length; ++c) {
// eight points finite differences estimation of a Jacobian column
FieldPVCoordinates<T> pvm4h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, -4 * h, directions[c]));
FieldPVCoordinates<T> pvm3h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, -3 * h, directions[c]));
FieldPVCoordinates<T> pvm2h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, -2 * h, directions[c]));
FieldPVCoordinates<T> pvm1h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, -1 * h, directions[c]));
FieldPVCoordinates<T> pvp1h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, +1 * h, directions[c]));
FieldPVCoordinates<T> pvp2h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, +2 * h, directions[c]));
FieldPVCoordinates<T> pvp3h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, +3 * h, directions[c]));
FieldPVCoordinates<T> pvp4h = combined.transformPVCoordinates(new FieldPVCoordinates<>(1.0, pv0, +4 * h, directions[c]));
FieldPVCoordinates<T> d4 = new FieldPVCoordinates<>(pvm4h, pvp4h);
FieldPVCoordinates<T> d3 = new FieldPVCoordinates<>(pvm3h, pvp3h);
FieldPVCoordinates<T> d2 = new FieldPVCoordinates<>(pvm2h, pvp2h);
FieldPVCoordinates<T> d1 = new FieldPVCoordinates<>(pvm1h, pvp1h);
double d = 1.0 / (840 * h);
FieldPVCoordinates<T> estimatedColumn = new FieldPVCoordinates<>(-3 * d, d4, 32 * d, d3, -168 * d, d2, 672 * d, d1);
// check analytical Jacobian against finite difference reference
Assert.assertEquals(estimatedColumn.getPosition().getX().getReal(), jacobian[0][c].getReal(), epsilonP);
Assert.assertEquals(estimatedColumn.getPosition().getY().getReal(), jacobian[1][c].getReal(), epsilonP);
Assert.assertEquals(estimatedColumn.getPosition().getZ().getReal(), jacobian[2][c].getReal(), epsilonP);
// check the rest of the matrix remains untouched
for (int l = 3; l < jacobian.length; ++l) {
Assert.assertEquals(l + 0.1 * c, jacobian[l][c].getReal(), 1.0e-15);
}
}
// check the rest of the matrix remains untouched
for (int c = directions.length; c < jacobian[0].length; ++c) {
for (int l = 0; l < jacobian.length; ++l) {
Assert.assertEquals(l + 0.1 * c, jacobian[l][c].getReal(), 1.0e-15);
}
}
}
}
}
use of org.orekit.utils.FieldPVCoordinates in project Orekit by CS-SI.
the class FieldEquinoctialOrbitTest method doTestNonInertialFrame.
private <T extends RealFieldElement<T>> void doTestNonInertialFrame(Field<T> field) throws IllegalArgumentException {
T zero = field.getZero();
FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
FieldVector3D<T> position = new FieldVector3D<>(zero.add(4512.9), zero.add(18260.), zero.add(-5127.));
FieldVector3D<T> velocity = new FieldVector3D<>(zero.add(134664.6), zero.add(90066.8), zero.add(72047.6));
FieldPVCoordinates<T> FieldPVCoordinates = new FieldPVCoordinates<>(position, velocity);
new FieldEquinoctialOrbit<>(FieldPVCoordinates, new Frame(FramesFactory.getEME2000(), Transform.IDENTITY, "non-inertial", false), date, mu);
}
use of org.orekit.utils.FieldPVCoordinates in project Orekit by CS-SI.
the class FieldKeplerianOrbitTest method doTestVeryLargeEccentricity.
private <T extends RealFieldElement<T>> void doTestVeryLargeEccentricity(final Field<T> field) {
FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
final Frame eme2000 = FramesFactory.getEME2000();
final double meanAnomaly = 1.;
final FieldKeplerianOrbit<T> orb0 = new FieldKeplerianOrbit<>(field.getZero().add(42600e3), field.getZero().add(0.9), field.getZero().add(0.00001), field.getZero().add(0), field.getZero().add(0), field.getZero().add(FastMath.toRadians(meanAnomaly)), PositionAngle.MEAN, eme2000, date, mu);
// big dV along Y
final FieldVector3D<T> deltaV = new FieldVector3D<>(field.getZero().add(0.0), field.getZero().add(110000.0), field.getZero());
final FieldPVCoordinates<T> pv1 = new FieldPVCoordinates<>(orb0.getPVCoordinates().getPosition(), orb0.getPVCoordinates().getVelocity().add(deltaV));
final FieldKeplerianOrbit<T> orb1 = new FieldKeplerianOrbit<>(pv1, eme2000, date, mu);
// Despite large eccentricity, the conversion of mean anomaly to hyperbolic eccentric anomaly
// converges in less than 50 iterations (issue #114)
final FieldPVCoordinates<T> pvTested = orb1.shiftedBy(field.getZero()).getPVCoordinates();
final FieldVector3D<T> pTested = pvTested.getPosition();
final FieldVector3D<T> vTested = pvTested.getVelocity();
final FieldPVCoordinates<T> pvReference = orb1.getPVCoordinates();
final FieldVector3D<T> pReference = pvReference.getPosition();
final FieldVector3D<T> vReference = pvReference.getVelocity();
final double threshold = 1.e-15;
Assert.assertEquals(0, pTested.subtract(pReference).getNorm().getReal(), pReference.getNorm().multiply(threshold).getReal());
Assert.assertEquals(0, vTested.subtract(vReference).getNorm().getReal(), vReference.getNorm().multiply(threshold).getReal());
}
use of org.orekit.utils.FieldPVCoordinates in project Orekit by CS-SI.
the class FieldKeplerianOrbitTest method doTestCopyNonKeplerianAcceleration.
private <T extends RealFieldElement<T>> void doTestCopyNonKeplerianAcceleration(Field<T> field) throws OrekitException {
final Frame eme2000 = FramesFactory.getEME2000();
// Define GEO satellite position
final FieldVector3D<T> position = new FieldVector3D<>(field.getZero().add(42164140), field.getZero(), field.getZero());
// Build PVCoodrinates starting from its position and computing the corresponding circular velocity
final FieldPVCoordinates<T> pv = new FieldPVCoordinates<>(position, new FieldVector3D<>(field.getZero(), position.getNorm().reciprocal().multiply(mu).sqrt(), field.getZero()));
// Build a KeplerianOrbit in eme2000
final FieldOrbit<T> orbit = new FieldKeplerianOrbit<>(pv, eme2000, FieldAbsoluteDate.getJ2000Epoch(field), mu);
// Build another KeplerianOrbit as a copy of the first one
final FieldOrbit<T> orbitCopy = new FieldKeplerianOrbit<>(orbit);
// Shift the orbit of a time-interval
// This works good
final FieldOrbit<T> shiftedOrbit = orbit.shiftedBy(10);
// This does not work
final FieldOrbit<T> shiftedOrbitCopy = orbitCopy.shiftedBy(10);
Assert.assertEquals(0.0, FieldVector3D.distance(shiftedOrbit.getPVCoordinates().getPosition(), shiftedOrbitCopy.getPVCoordinates().getPosition()).getReal(), 1.0e-10);
Assert.assertEquals(0.0, FieldVector3D.distance(shiftedOrbit.getPVCoordinates().getVelocity(), shiftedOrbitCopy.getPVCoordinates().getVelocity()).getReal(), 1.0e-10);
}
use of org.orekit.utils.FieldPVCoordinates in project Orekit by CS-SI.
the class FieldSpacecraftStateTest method doTestShiftVsEcksteinHechlerError.
private <T extends RealFieldElement<T>> void doTestShiftVsEcksteinHechlerError(final Field<T> field) throws OrekitException {
T zero = field.getZero();
T mass = zero.add(2500.);
T a = zero.add(rOrbit.getA());
T e = zero.add(rOrbit.getE());
T i = zero.add(rOrbit.getI());
T pa = zero.add(1.9674147913622104);
T raan = zero.add(FastMath.toRadians(261));
T lv = zero.add(0);
final double ae = 6.378137e6;
final double c20 = -1.08263e-3;
final double c30 = 2.54e-6;
final double c40 = 1.62e-6;
final double c50 = 2.3e-7;
final double c60 = -5.5e-7;
// polynomial models for interpolation error in position, velocity, acceleration and attitude
// these models grow as follows
// interpolation time (s) position error (m) velocity error (m/s) acceleration error (m/s²) attitude error (°)
// 60 2 0.07 0.002 0.00002
// 120 12 0.3 0.005 0.00009
// 300 170 1.6 0.012 0.0009
// 600 1200 5.7 0.024 0.006
// 900 3600 10.6 0.034 0.02
// the expected maximum residuals with respect to these models are about 0.4m, 0.5mm/s, 8μm/s² and 3e-6°
PolynomialFunction pModel = new PolynomialFunction(new double[] { 1.5664070631933846e-01, 7.5504722733047560e-03, -8.2460562451009510e-05, 6.9546332080305580e-06, -1.7045365367533077e-09, -4.2187860791066264e-13 });
PolynomialFunction vModel = new PolynomialFunction(new double[] { -3.5472364019908720e-04, 1.6568103861124980e-05, 1.9637913327830596e-05, -3.4248792843039766e-09, -5.6565135131014254e-12, 1.4730170946808630e-15 });
PolynomialFunction aModel = new PolynomialFunction(new double[] { 3.0731707577766896e-06, 3.9770746399850350e-05, 1.9779039254538660e-09, 8.0263328220724900e-12, -1.5600835252366078e-14, 1.1785257001549687e-18 });
PolynomialFunction rModel = new PolynomialFunction(new double[] { -2.7689062063188115e-06, 1.7406542538258334e-07, 2.5109795349592287e-09, 2.0399322661074575e-11, 9.9126348912426750e-15, -3.5015638905729510e-18 });
FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field, new DateComponents(2004, 01, 01), TimeComponents.H00, TimeScalesFactory.getUTC());
FieldKeplerianOrbit<T> orbit = new FieldKeplerianOrbit<>(a, e, i, pa, raan, lv, PositionAngle.TRUE, FramesFactory.getEME2000(), date, mu);
BodyCenterPointing attitudeLaw = new BodyCenterPointing(orbit.getFrame(), earth);
FieldPropagator<T> propagator = new FieldEcksteinHechlerPropagator<>(orbit, attitudeLaw, mass, ae, mu, c20, c30, c40, c50, c60);
FieldAbsoluteDate<T> centerDate = orbit.getDate().shiftedBy(100.0);
FieldSpacecraftState<T> centerState = propagator.propagate(centerDate);
double maxResidualP = 0;
double maxResidualV = 0;
double maxResidualA = 0;
double maxResidualR = 0;
for (T dt = field.getZero(); dt.getReal() < 900.0; dt = dt.add(5)) {
FieldSpacecraftState<T> shifted = centerState.shiftedBy(dt);
FieldSpacecraftState<T> propagated = propagator.propagate(centerDate.shiftedBy(dt));
FieldPVCoordinates<T> dpv = new FieldPVCoordinates<>(propagated.getPVCoordinates(), shifted.getPVCoordinates());
double residualP = pModel.value(dt.getReal()) - dpv.getPosition().getNorm().getReal();
double residualV = vModel.value(dt.getReal()) - dpv.getVelocity().getNorm().getReal();
double residualA = aModel.value(dt.getReal()) - dpv.getAcceleration().getNorm().getReal();
double residualR = rModel.value(dt.getReal()) - FastMath.toDegrees(FieldRotation.distance(shifted.getAttitude().getRotation(), propagated.getAttitude().getRotation()).getReal());
maxResidualP = FastMath.max(maxResidualP, FastMath.abs(residualP));
maxResidualV = FastMath.max(maxResidualV, FastMath.abs(residualV));
maxResidualA = FastMath.max(maxResidualA, FastMath.abs(residualA));
maxResidualR = FastMath.max(maxResidualR, FastMath.abs(residualR));
}
Assert.assertEquals(0.40, maxResidualP, 0.01);
Assert.assertEquals(4.9e-4, maxResidualV, 1.0e-5);
Assert.assertEquals(2.8e-6, maxResidualR, 1.0e-1);
}
Aggregations