Skip to content

GitLab

Commit 11de4863 authored by Bryan Cazabonne's avatar Bryan Cazabonne
Browse files

Fixed eD and eY equation in ECOM2 model.

parent a86752a5
Hide whitespace changes
Inline Side-by-side
src/changes/changes.xml
View file @ 11de4863
......@@ -21,6 +21,9 @@
</properties>
<body>
<release version="11.2" date="TBD" description="TBD">
<action dev="bryan" type="fix" issue="910">
Fixed eD and eY equation in ECOM2 model.
</action>
<action dev="pascal" type="fix" issue="908">
Fixed unmanaged comment in OMM.
</action>
......
src/main/java/org/orekit/forces/radiation/ECOM2.java
View file @ 11de4863
......@@ -148,19 +148,23 @@ public class ECOM2 extends AbstractRadiationForceModel {
/** {@inheritDoc} */
@Override
public Vector3D acceleration(final SpacecraftState s, final double[] parameters) {
// Spacecraft and Sun position vectors (expressed in the spacecraft's frame)
final Vector3D satPos = s.getPVCoordinates().getPosition();
final Vector3D sunPos = sun.getPVCoordinates(s.getDate(), s.getFrame()).getPosition();
// Build the coordinate system
final Vector3D Z = s.getPVCoordinates().getMomentum().normalize();
final Vector3D sunPos = sun.getPVCoordinates(s.getDate(), s.getFrame()).getPosition().normalize();
final Vector3D Y = Z.crossProduct(sunPos);
final Vector3D X = Y.crossProduct(Z);
final Vector3D Z = s.getPVCoordinates().getMomentum();
final Vector3D Y = Z.crossProduct(sunPos).normalize();
final Vector3D X = Y.crossProduct(Z).normalize();
// Build eD, eY, eB vectors
final Vector3D position = s.getPVCoordinates().getPosition().normalize();
final Vector3D eD = sunPos.add(-1.0, position);
final Vector3D eY = position.crossProduct(eD);
final Vector3D eD = sunPos.subtract(satPos).normalize();
final Vector3D eY = eD.crossProduct(satPos).normalize();
final Vector3D eB = eD.crossProduct(eY);
// Angular argument difference u_s - u
final double delta_u = FastMath.atan2(position.dotProduct(Y), position.dotProduct(X));
final double delta_u = FastMath.atan2(satPos.dotProduct(Y), satPos.dotProduct(X));
// Compute B(u)
double b_u = parameters[0];
......@@ -181,20 +185,23 @@ public class ECOM2 extends AbstractRadiationForceModel {
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> FieldVector3D<T> acceleration(final FieldSpacecraftState<T> s, final T[] parameters) {
// Spacecraft and Sun position vectors (expressed in the spacecraft's frame)
final FieldVector3D<T> satPos = s.getPVCoordinates().getPosition();
final FieldVector3D<T> sunPos = sun.getPVCoordinates(s.getDate(), s.getFrame()).getPosition();
// Build the coordinate system
final FieldVector3D<T> Z = s.getPVCoordinates().getMomentum().normalize();
final FieldVector3D<T> sunPos = sun.getPVCoordinates(s.getDate(), s.getFrame()).getPosition().normalize();
final FieldVector3D<T> Y = Z.crossProduct(sunPos);
final FieldVector3D<T> X = Y.crossProduct(Z);
final FieldVector3D<T> Z = s.getPVCoordinates().getMomentum();
final FieldVector3D<T> Y = Z.crossProduct(sunPos).normalize();
final FieldVector3D<T> X = Y.crossProduct(Z).normalize();
// Build eD, eY, eB vectors
final FieldVector3D<T> position = s.getPVCoordinates().getPosition().normalize();
final FieldVector3D<T> eD = sunPos.add(-1.0, position);
final FieldVector3D<T> eY = position.crossProduct(eD);
final FieldVector3D<T> eD = sunPos.subtract(satPos).normalize();
final FieldVector3D<T> eY = eD.crossProduct(satPos).normalize();
final FieldVector3D<T> eB = eD.crossProduct(eY);
// Angular argument difference u_s - u
final T delta_u = FastMath.atan2(position.dotProduct(Y), position.dotProduct(X));
final T delta_u = FastMath.atan2(satPos.dotProduct(Y), satPos.dotProduct(X));
// Compute B(u)
T b_u = parameters[0];
......
src/test/java/org/orekit/forces/radiation/ECOM2Test.java
View file @ 11de4863
......@@ -496,6 +496,45 @@ public class ECOM2Test extends AbstractForceModelTest {
Assert.assertTrue(calc.getCalls() < 7100);
}
@Test
public void testRealAndFieldComparison() {
// Orbital parameters from GNSS almanac
final int freeParameters = 6;
final Gradient sma = Gradient.variable(freeParameters, 0, 26559614.1);
final Gradient ecc = Gradient.variable(freeParameters, 1, 0.00522136);
final Gradient inc = Gradient.variable(freeParameters, 2, 0.963785748);
final Gradient aop = Gradient.variable(freeParameters, 3, 0.451712027);
final Gradient raan = Gradient.variable(freeParameters, 4, -1.159458779);
final Gradient lm = Gradient.variable(freeParameters, 4, -2.105941778);
// Field and zero
final Field<Gradient> field = sma.getField();
// Epoch
final FieldAbsoluteDate<Gradient> epoch = FieldAbsoluteDate.getJ2000Epoch(field);
// Create a Keplerian orbit
FieldKeplerianOrbit<Gradient> orbit = new FieldKeplerianOrbit<>(sma, ecc, inc, aop, raan, lm,
PositionAngle.MEAN,
FramesFactory.getEME2000(),
epoch,
field.getZero().add(Constants.EIGEN5C_EARTH_MU));
// Model
final ECOM2 forceModel = new ECOM2(2, 2, 1e-7, CelestialBodyFactory.getSun(), Constants.EGM96_EARTH_EQUATORIAL_RADIUS);
// Field acceleration
final FieldVector3D<Gradient> accField = forceModel.acceleration(new FieldSpacecraftState<>(orbit), forceModel.getParameters(field));
// Real acceleration
final Vector3D accReal = forceModel.acceleration(new SpacecraftState(orbit.toOrbit()), forceModel.getParameters());
// Verify
Assert.assertEquals(0.0, accReal.distance(accField.toVector3D()), 1.0e-20);
}
private static class SolarStepHandler implements OrekitFixedStepHandler {
public void handleStep(SpacecraftState currentState) {
......@@ -504,7 +543,7 @@ public class ECOM2Test extends AbstractForceModelTest {
final double alpha = FastMath.toDegrees(FastMath.atan2(dey, dex));
Assert.assertTrue(alpha > 100.0);
Assert.assertTrue(alpha < 112.0);
checkRadius(FastMath.sqrt(dex * dex + dey * dey), 0.003482, 0.003525);
checkRadius(FastMath.sqrt(dex * dex + dey * dey), 0.003482, 0.003529);
}
}
......
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