Commit 7e736694 authored by Luc Maisonobe's avatar Luc Maisonobe
Browse files

Fixed an issue with non-Keplerian accelerations.

Thanks to roberto Alacevich for finding the issue.
parent 485b4746
......@@ -429,8 +429,8 @@ public class CircularOrbit
final double h2 = hx * hx + hy * hy;
final double h = FastMath.sqrt(h2);
raan = FastMath.atan2(hy, hx);
final double cosRaan = hx / h;
final double sinRaan = hy / h;
final double cosRaan = h == 0 ? FastMath.cos(raan) : hx / h;
final double sinRaan = h == 0 ? FastMath.sin(raan) : hy / h;
final double equiEx = op.getEquinoctialEx();
final double equiEy = op.getEquinoctialEy();
ex = equiEx * cosRaan + equiEy * sinRaan;
......
......@@ -407,8 +407,8 @@ public class FieldCircularOrbit<T extends RealFieldElement<T>>
final T h2 = hx.multiply(hx).add(hy.multiply(hy));
final T h = h2.sqrt();
raan = hy.atan2(hx);
final T cosRaan = hx.divide(h);
final T sinRaan = hy.divide(h);
final T cosRaan = h.getReal() == 0 ? raan.cos() : hx.divide(h);
final T sinRaan = h.getReal() == 0 ? raan.sin() : hy.divide(h);
final T equiEx = op.getEquinoctialEx();
final T equiEy = op.getEquinoctialEy();
ex = equiEx.multiply(cosRaan).add(equiEy.multiply(sinRaan));
......
......@@ -312,6 +312,28 @@ public class FieldKeplerianOrbit<T extends RealFieldElement<T>> extends FieldOrb
public FieldKeplerianOrbit(final TimeStampedFieldPVCoordinates<T> pvCoordinates,
final Frame frame, final double mu)
throws IllegalArgumentException {
this(pvCoordinates, frame, mu, hasNonKeplerianAcceleration(pvCoordinates, mu));
}
/** Constructor from Cartesian parameters.
*
* <p> The acceleration provided in {@code FieldPVCoordinates} is accessible using
* {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
* use {@code mu} and the position to compute the acceleration, including
* {@link #shiftedBy(RealFieldElement)} and {@link #getPVCoordinates(FieldAbsoluteDate, Frame)}.
*
* @param pvCoordinates the PVCoordinates of the satellite
* @param frame the frame in which are defined the {@link FieldPVCoordinates}
* (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
* @param mu central attraction coefficient (m³/s²)
* @param reliableAcceleration if true, the acceleration is considered to be reliable
* @exception IllegalArgumentException if frame is not a {@link
* Frame#isPseudoInertial pseudo-inertial frame}
*/
private FieldKeplerianOrbit(final TimeStampedFieldPVCoordinates<T> pvCoordinates,
final Frame frame, final double mu,
final boolean reliableAcceleration)
throws IllegalArgumentException {
super(pvCoordinates, frame, mu);
......@@ -372,7 +394,7 @@ public class FieldKeplerianOrbit<T extends RealFieldElement<T>> extends FieldOrb
FACTORIES.put(a.getField(), new FDSFactory<>(a.getField(), 1, 1));
}
if (hasNonKeplerianAcceleration(pvCoordinates, mu)) {
if (reliableAcceleration) {
// we have a relevant acceleration, we can compute derivatives
final T[][] jacobian = MathArrays.buildArray(a.getField(), 6, 6);
......@@ -432,7 +454,7 @@ public class FieldKeplerianOrbit<T extends RealFieldElement<T>> extends FieldOrb
* Frame#isPseudoInertial pseudo-inertial frame}
*/
public FieldKeplerianOrbit(final FieldPVCoordinates<T> FieldPVCoordinates,
final Frame frame, final FieldAbsoluteDate<T> date, final double mu)
final Frame frame, final FieldAbsoluteDate<T> date, final double mu)
throws IllegalArgumentException {
this(new TimeStampedFieldPVCoordinates<>(date, FieldPVCoordinates), frame, mu);
}
......@@ -441,7 +463,7 @@ public class FieldKeplerianOrbit<T extends RealFieldElement<T>> extends FieldOrb
* @param op orbital parameters to copy
*/
public FieldKeplerianOrbit(final FieldOrbit<T> op) {
this(op.getPVCoordinates(), op.getFrame(), op.getDate(), op.getMu());
this(op.getPVCoordinates(), op.getFrame(), op.getMu(), op.hasDerivatives());
}
/** {@inheritDoc} */
......
......@@ -282,6 +282,28 @@ public class KeplerianOrbit extends Orbit {
public KeplerianOrbit(final TimeStampedPVCoordinates pvCoordinates,
final Frame frame, final double mu)
throws IllegalArgumentException {
this(pvCoordinates, frame, mu, hasNonKeplerianAcceleration(pvCoordinates, mu));
}
/** Constructor from Cartesian parameters.
*
* <p> The acceleration provided in {@code pvCoordinates} is accessible using
* {@link #getPVCoordinates()} and {@link #getPVCoordinates(Frame)}. All other methods
* use {@code mu} and the position to compute the acceleration, including
* {@link #shiftedBy(double)} and {@link #getPVCoordinates(AbsoluteDate, Frame)}.
*
* @param pvCoordinates the PVCoordinates of the satellite
* @param frame the frame in which are defined the {@link PVCoordinates}
* (<em>must</em> be a {@link Frame#isPseudoInertial pseudo-inertial frame})
* @param mu central attraction coefficient (m³/s²)
* @param reliableAcceleration if true, the acceleration is considered to be reliable
* @exception IllegalArgumentException if frame is not a {@link
* Frame#isPseudoInertial pseudo-inertial frame}
*/
private KeplerianOrbit(final TimeStampedPVCoordinates pvCoordinates,
final Frame frame, final double mu,
final boolean reliableAcceleration)
throws IllegalArgumentException {
super(pvCoordinates, frame, mu);
// compute inclination
......@@ -328,7 +350,7 @@ public class KeplerianOrbit extends Orbit {
partialPV = pvCoordinates;
if (hasNonKeplerianAcceleration(pvCoordinates, mu)) {
if (reliableAcceleration) {
// we have a relevant acceleration, we can compute derivatives
final double[][] jacobian = new double[6][6];
......@@ -395,7 +417,7 @@ public class KeplerianOrbit extends Orbit {
* @param op orbital parameters to copy
*/
public KeplerianOrbit(final Orbit op) {
this(op.getPVCoordinates(), op.getFrame(), op.getDate(), op.getMu());
this(op.getPVCoordinates(), op.getFrame(), op.getMu(), op.hasDerivatives());
}
/** {@inheritDoc} */
......
......@@ -711,6 +711,37 @@ public class CartesianOrbitTest {
orbit.toString());
}
@Test
public void testCopyNonKeplerianAcceleration() throws OrekitException {
final Frame eme2000 = FramesFactory.getEME2000();
// Define GEO satellite position
final Vector3D position = new Vector3D(42164140, 0, 0);
// Build PVCoodrinates starting from its position and computing the corresponding circular velocity
final PVCoordinates pv = new PVCoordinates(position,
new Vector3D(0, FastMath.sqrt(mu / position.getNorm()), 0));
// Build a KeplerianOrbit in eme2000
final Orbit orbit = new CartesianOrbit(pv, eme2000, date, mu);
// Build another KeplerianOrbit as a copy of the first one
final Orbit orbitCopy = new CartesianOrbit(orbit);
// Shift the orbit of a time-interval
final Orbit shiftedOrbit = orbit.shiftedBy(10); // This works good
final Orbit shiftedOrbitCopy = orbitCopy.shiftedBy(10); // This does not work
Assert.assertEquals(0.0,
Vector3D.distance(shiftedOrbit.getPVCoordinates().getPosition(),
shiftedOrbitCopy.getPVCoordinates().getPosition()),
1.0e-10);
Assert.assertEquals(0.0,
Vector3D.distance(shiftedOrbit.getPVCoordinates().getVelocity(),
shiftedOrbitCopy.getPVCoordinates().getVelocity()),
1.0e-10);
}
@Before
public void setUp() {
......
......@@ -1168,6 +1168,37 @@ public class CircularOrbitTest {
orbit.toString());
}
@Test
public void testCopyNonKeplerianAcceleration() throws OrekitException {
final Frame eme2000 = FramesFactory.getEME2000();
// Define GEO satellite position
final Vector3D position = new Vector3D(42164140, 0, 0);
// Build PVCoodrinates starting from its position and computing the corresponding circular velocity
final PVCoordinates pv = new PVCoordinates(position,
new Vector3D(0, FastMath.sqrt(mu / position.getNorm()), 0));
// Build a KeplerianOrbit in eme2000
final Orbit orbit = new CircularOrbit(pv, eme2000, date, mu);
// Build another KeplerianOrbit as a copy of the first one
final Orbit orbitCopy = new CircularOrbit(orbit);
// Shift the orbit of a time-interval
final Orbit shiftedOrbit = orbit.shiftedBy(10); // This works good
final Orbit shiftedOrbitCopy = orbitCopy.shiftedBy(10); // This does not work
Assert.assertEquals(0.0,
Vector3D.distance(shiftedOrbit.getPVCoordinates().getPosition(),
shiftedOrbitCopy.getPVCoordinates().getPosition()),
1.0e-10);
Assert.assertEquals(0.0,
Vector3D.distance(shiftedOrbit.getPVCoordinates().getVelocity(),
shiftedOrbitCopy.getPVCoordinates().getVelocity()),
1.0e-10);
}
@Before
public void setUp() {
......
/* Copyright 2002-2017 CS Systèmes d'Information
/* Copyright 2002-2017 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
......@@ -1103,6 +1103,37 @@ public class EquinoctialOrbitTest {
orbit.toString());
}
@Test
public void testCopyNonKeplerianAcceleration() throws OrekitException {
final Frame eme2000 = FramesFactory.getEME2000();
// Define GEO satellite position
final Vector3D position = new Vector3D(42164140, 0, 0);
// Build PVCoodrinates starting from its position and computing the corresponding circular velocity
final PVCoordinates pv = new PVCoordinates(position,
new Vector3D(0, FastMath.sqrt(mu / position.getNorm()), 0));
// Build a KeplerianOrbit in eme2000
final Orbit orbit = new EquinoctialOrbit(pv, eme2000, date, mu);
// Build another KeplerianOrbit as a copy of the first one
final Orbit orbitCopy = new EquinoctialOrbit(orbit);
// Shift the orbit of a time-interval
final Orbit shiftedOrbit = orbit.shiftedBy(10); // This works good
final Orbit shiftedOrbitCopy = orbitCopy.shiftedBy(10); // This does not work
Assert.assertEquals(0.0,
Vector3D.distance(shiftedOrbit.getPVCoordinates().getPosition(),
shiftedOrbitCopy.getPVCoordinates().getPosition()),
1.0e-10);
Assert.assertEquals(0.0,
Vector3D.distance(shiftedOrbit.getPVCoordinates().getVelocity(),
shiftedOrbitCopy.getPVCoordinates().getVelocity()),
1.0e-10);
}
@Before
public void setUp() {
......
......@@ -191,6 +191,11 @@ public class FieldCartesianOrbitTest {
doTestEquatorialRetrograde(Decimal64Field.getInstance());
}
@Test
public void testCopyNonKeplerianAcceleration() throws OrekitException {
doTestCopyNonKeplerianAcceleration(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>> void doTestCartesianToCartesian(Field<T> field)
throws NoSuchMethodException, SecurityException, IllegalAccessException, IllegalArgumentException, InvocationTargetException {
......@@ -832,5 +837,41 @@ public class FieldCartesianOrbitTest {
orbit.toString());
}
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 FieldCartesianOrbit<>(pv, eme2000, FieldAbsoluteDate.getJ2000Epoch(field), mu);
// Build another KeplerianOrbit as a copy of the first one
final FieldOrbit<T> orbitCopy = new FieldCartesianOrbit<>(orbit);
// Shift the orbit of a time-interval
final FieldOrbit<T> shiftedOrbit = orbit.shiftedBy(10); // This works good
final FieldOrbit<T> shiftedOrbitCopy = orbitCopy.shiftedBy(10); // This does not work
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);
}
}
......@@ -215,6 +215,11 @@ public class FieldCircularOrbitTest {
doTestEquatorialRetrograde(Decimal64Field.getInstance());
}
@Test
public void testCopyNonKeplerianAcceleration() throws OrekitException {
doTestCopyNonKeplerianAcceleration(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>> void doTestCircularToEquinoctialEll(Field<T> field) {
T zero = field.getZero();
......@@ -1331,4 +1336,40 @@ public class FieldCircularOrbitTest {
orbit.toString());
}
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 FieldCircularOrbit<>(pv, eme2000, FieldAbsoluteDate.getJ2000Epoch(field), mu);
// Build another KeplerianOrbit as a copy of the first one
final FieldOrbit<T> orbitCopy = new FieldCircularOrbit<>(orbit);
// Shift the orbit of a time-interval
final FieldOrbit<T> shiftedOrbit = orbit.shiftedBy(10); // This works good
final FieldOrbit<T> shiftedOrbitCopy = orbitCopy.shiftedBy(10); // This does not work
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);
}
}
......@@ -186,6 +186,11 @@ public class FieldEquinoctialOrbitTest {
doTestEquatorialRetrograde(Decimal64Field.getInstance());
}
@Test
public void testCopyNonKeplerianAcceleration() throws OrekitException {
doTestCopyNonKeplerianAcceleration(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>> void doTestEquinoctialToEquinoctialEll(Field<T> field) {
T zero = field.getZero();
FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
......@@ -1245,4 +1250,40 @@ public class FieldEquinoctialOrbitTest {
orbit.toString());
}
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
final FieldOrbit<T> shiftedOrbit = orbit.shiftedBy(10); // This works good
final FieldOrbit<T> shiftedOrbitCopy = orbitCopy.shiftedBy(10); // This does not work
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);
}
}
......@@ -236,6 +236,11 @@ public class FieldKeplerianOrbitTest {
doTestPerfectlyEquatorialConversion(Decimal64Field.getInstance());
}
@Test
public void testCopyNonKeplerianAcceleration() throws OrekitException {
doTestCopyNonKeplerianAcceleration(Decimal64Field.getInstance());
}
private <T extends RealFieldElement<T>> void doTestKeplerianToKeplerian(final Field<T> field) {
FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
......@@ -250,7 +255,7 @@ public class FieldKeplerianOrbitTest {
// elliptic orbit
FieldKeplerianOrbit<T> kep =
new FieldKeplerianOrbit<>(a, e, i, pa, raan, anomaly, PositionAngle.MEAN,
new FieldKeplerianOrbit<>(a, e, i, pa, raan, anomaly, PositionAngle.MEAN,
FramesFactory.getEME2000(), date, mu);
FieldVector3D<T> pos = kep.getPVCoordinates().getPosition();
......@@ -1015,7 +1020,7 @@ public class FieldKeplerianOrbitTest {
FieldAbsoluteDate<T> dateTca = new FieldAbsoluteDate<>(field, 2000, 04, 01, 0, 0, 0.000, TimeScalesFactory.getUTC());
double mu = 3.986004415e+14;
FieldKeplerianOrbit<T> orbKep = new FieldKeplerianOrbit<>(field.getZero().add(7000000.0), field.getZero().add(0.01), field.getZero().add(FastMath.toRadians(80.)), field.getZero().add(FastMath.toRadians(80.)), field.getZero().add(FastMath.toRadians(20.)),
FieldKeplerianOrbit<T> orbKep = new FieldKeplerianOrbit<>(field.getZero().add(7000000.0), field.getZero().add(0.01), field.getZero().add(FastMath.toRadians(80.)), field.getZero().add(FastMath.toRadians(80.)), field.getZero().add(FastMath.toRadians(20.)),
field.getZero().add(FastMath.toRadians(40.)), PositionAngle.MEAN,
FramesFactory.getEME2000(), dateTca, mu);
......@@ -1823,6 +1828,42 @@ public class FieldKeplerianOrbitTest {
orbit.toString());
}
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
final FieldOrbit<T> shiftedOrbit = orbit.shiftedBy(10); // This works good
final FieldOrbit<T> shiftedOrbitCopy = orbitCopy.shiftedBy(10); // This does not work
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);
}
@Before
public void setUp() {
......
......@@ -1551,6 +1551,37 @@ public class KeplerianOrbitTest {
orbit.toString());
}
@Test
public void testCopyNonKeplerianAcceleration() throws OrekitException {
final Frame eme2000 = FramesFactory.getEME2000();
// Define GEO satellite position
final Vector3D position = new Vector3D(42164140, 0, 0);
// Build PVCoodrinates starting from its position and computing the corresponding circular velocity
final PVCoordinates pv = new PVCoordinates(position,
new Vector3D(0, FastMath.sqrt(mu / position.getNorm()), 0));
// Build a KeplerianOrbit in eme2000
final Orbit orbit = new KeplerianOrbit(pv, eme2000, date, mu);
// Build another KeplerianOrbit as a copy of the first one
final Orbit orbitCopy = new KeplerianOrbit(orbit);
// Shift the orbit of a time-interval
final Orbit shiftedOrbit = orbit.shiftedBy(10); // This works good
final Orbit shiftedOrbitCopy = orbitCopy.shiftedBy(10); // This does not work
Assert.assertEquals(0.0,
Vector3D.distance(shiftedOrbit.getPVCoordinates().getPosition(),
shiftedOrbitCopy.getPVCoordinates().getPosition()),
1.0e-10);
Assert.assertEquals(0.0,
Vector3D.distance(shiftedOrbit.getPVCoordinates().getVelocity(),
shiftedOrbitCopy.getPVCoordinates().getVelocity()),
1.0e-10);
}
@Before
public void setUp() {
......
......@@ -810,8 +810,8 @@ public class KeplerianPropagatorTest {
ObjectOutputStream oos = new ObjectOutputStream(bos);
oos.writeObject(propagator.getGeneratedEphemeris());
Assert.assertTrue(bos.size() > 2400);
Assert.assertTrue(bos.size() < 2500);
Assert.assertTrue(bos.size() > 2300);
Assert.assertTrue(bos.size() < 2400);
ByteArrayInputStream bis = new ByteArrayInputStream(bos.toByteArray());
ObjectInputStream ois = new ObjectInputStream(bis);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment