Commit 958336d5 authored by Bryan Cazabonne's avatar Bryan Cazabonne
Browse files

Fixed wrong equations in Field Brouwer-Lyddane orbit propagator.

parent 74a5b9a2
......@@ -698,7 +698,6 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator {
0.375 * yp2 * yp2 * ((-5.0 + 12.0 * n + 9.0 * n2) * cosI1 + (-35.0 - 36.0 * n - 5.0 * n2) * cosI3) +
1.25 * yp4 * (5.0 - 3.0 * n2) * cosI1 * (3.0 - 7.0 * cosI2);
final double cA = 1.0 - 11.0 * cosI2 - 40.0 * cosI4 / C5c2;
final double cB = 1.0 - 3.0 * cosI2 - 8.0 * cosI4 / C5c2;
final double cC = 1.0 - 9.0 * cosI2 - 24.0 * cosI4 / C5c2;
......@@ -716,7 +715,6 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator {
(4.0 + 3.0 * epp2) *
(3.0 + 16.0 * cosI2 / C5c2 + 40.0 * cosI4 / C5c2 / C5c2);
// long periodic multiplicative
dei3sg = 35.0 / 96.0 * yp5 / yp2 * epp2 * n2 * cD * sinI1;
de2sg = -1.0 / 12.0 * epp * n2 / yp2 * qyp2_4;
......@@ -724,7 +722,6 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator {
1.0 / 4.0 * n2 / yp2 * (yp3 + 5.0 / 16.0 * yp5 * (4.0 + 3.0 * epp2) * cC)) * sinI1;
de = epp2 * n2 / 24.0 / yp2 * qyp2_4;
final double qyp52quotient = epp * (-32.0 + 81.0 * epp4) / (4.0 + 3.0 * epp2 + n * (4.0 + 9.0 * epp2));
dlgs2g = 1.0 / 48.0 / yp2 * (-3.0 * yp2 * yp2 * qyp22 +
10.0 * yp4 * qyp42 ) +
......@@ -746,14 +743,11 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator {
epp * cosI1 / (4.0 * yp2 * sinI1) * (yp3 + 0.3125 * yp5 * (4.0 + 3.0 * epp2) * cC) +
1.875 / (4.0 * yp2) * yp5 * qyp52bis;
// short periodic multiplicative
aC = -yp2 * C3c2 * app / n3;
aCbis = y2 * app * C3c2;
ac2g2f = y2 * app * 3.0 * sinI2;
double qe = 0.5 * n2 * y2 * C3c2 / n6;
eC = qe * epp / (1.0 + n3) * (3.0 - epp2 * (3.0 - epp2));
ecf = 3.0 * qe;
......@@ -775,7 +769,6 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator {
qi = yp2 * cosI1 * sinI1;
ic2f2g = 1.5 * qi;
double qgl1 = 0.25 * yp2;
double qgl2 = 0.25 * yp2 * epp * n2 / (1.0 + n);
glf = qgl1 * -6.0 * C5c2;
......@@ -793,7 +786,6 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator {
1.0 / 3.0 * qgl2;
glos2g3f = qgl2;
final double qh = 3.0 * yp2 * cosI1;
hf = -qh;
hl = qh;
......@@ -802,7 +794,6 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator {
hs2g2f = 1.5 * yp2 * cosI1;
hsfc2g2f = -epp * yp2 * cosI1;
final double qedl = -0.25 * yp2 * n3;
edls2g = 1.0 / 24.0 * epp * n3 / yp2 * qyp2_4;
edlcg = -0.25 * yp3 / yp2 * n3 * sinI1 -
......
......@@ -639,14 +639,13 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
this.mean = mean;
this.mass = mass;
final T zero = mass.getField().getZero();
final T one = mass.getField().getOne();
final T app = mean.getA();
xnotDot = mu.divide(mean.getA()).sqrt().divide(mean.getA());
xnotDot = mu.divide(app).sqrt().divide(app);
// preliminary processing
final T q = zero.add(referenceRadius).divide(app);
final T q = app.divide(referenceRadius).reciprocal();
T ql = q.multiply(q);
final T y2 = ql.multiply(-0.5 * ck0[2]);
......@@ -710,7 +709,6 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
add(cosI2.multiply(n2.multiply(126.0).add(-270.0))).
add(cosI4.multiply(n2.multiply(-189.0).add(385.0)))));
ht = yp2.multiply(-3.0).multiply(cosI1).
add(yp2.multiply(0.375).multiply(yp2).multiply(cosI1.multiply(n2.multiply(9.0).add(n.multiply(12.0)).add(-5.0)).
add(cosI3.multiply(n2.multiply(-5.0).add(n.multiply(-36.0)).add(-35.0))))).
......@@ -721,16 +719,15 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
final T cC = one.subtract(cosI2.multiply(9.0)).subtract(cosI4.multiply(24.0).divide(C5c2));
final T cD = one.subtract(cosI2.multiply(5.0)).subtract(cosI4.multiply(16.0).divide(C5c2));
final T qyp2_4 = yp2.multiply(3.0).multiply(yp2).multiply(cA).
subtract(yp4.multiply(10.0).multiply(cB));
final T qyp52 = cosI1.multiply(epp3).multiply(cD.multiply(0.5).divide(sinI1).
add(sinI1.multiply(cosI4.divide(C5c2).divide(C5c2).multiply(80.0).
add(cosI2.divide(C5c2).multiply(32.0).
add(5.0)))));
final T qyp22 = epp2.add(2.0).subtract(cosI2.multiply(11.0).multiply(epp2.multiply(3.0).add(2.0))).
subtract(cosI4.divide(C5c2).multiply(40.0).multiply(epp2.multiply(3.0).add(2.0))).
subtract(cosI6.divide(C5c2).divide(C5c2).multiply(400.0).multiply(epp2));
final T qyp22 = epp2.add(2.0).subtract(epp2.multiply(3.0).add(2.0).multiply(11.0).multiply(cosI2)).
subtract(epp2.multiply(5.0).add(2.0).multiply(40.0).multiply(cosI4.divide(C5c2))).
subtract(epp2.multiply(400.0).multiply(cosI6).divide(C5c2).divide(C5c2));
final T qyp42 = one.divide(5.0).multiply(qyp22.
add(one.multiply(4.0).multiply(epp2.
add(2.0).
......@@ -739,21 +736,19 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
multiply(cosI4.divide(C5c2).divide(C5c2).multiply(40.0).
add(cosI2.divide(C5c2).multiply(16.0)).
add(3.0));
// long periodic multiplicative
dei3sg = yp5.divide(yp2).multiply(35.0 / 96.0).multiply(epp2).multiply(n2).multiply(cD).multiply(sinI1);
de2sg = qyp2_4.divide(yp2).multiply(epp).multiply(n2).multiply(-1.0 / 12.0);
deisg = sinI1.multiply(yp5.divide(yp2).multiply(-35.0 / 128.0).multiply(epp2).multiply(n2).multiply(cD).
add(yp2.multiply(n2).multiply(1.0 / 4.0).multiply(yp3.
add(yp5.multiply(5.0 / 16.0).multiply(cC).multiply(epp2.multiply(3.0).add(4.0))))));
deisg = sinI1.multiply(yp5.multiply(-35.0 / 128.0).divide(yp2).multiply(epp2).multiply(n2).multiply(cD).
add(n2.multiply(0.25).divide(yp2).multiply(yp3.add(yp5.multiply(5.0 / 16.0).multiply(epp2.multiply(3.0).add(4.0)).multiply(cC)))));
de = epp2.multiply(n2).multiply(qyp2_4).divide(24.0).divide(yp2);
final T qyp52quotient = epp.multiply(epp4.multiply(81.0).add(-32.0)).divide(n.multiply(epp2.multiply(9.0).add(4.0)).add(epp2.multiply(3.0)).add(4.0));
dlgs2g = yp2.multiply(1.0 / 48.0).multiply(qyp22.multiply(yp2).multiply(yp2).multiply(-3.0).add(qyp42.multiply(yp4).multiply(10.0))).
add(qyp2_4.divide(yp2).multiply(n3).divide(24.0));
dlgc3g = yp5.divide(yp2).multiply(35.0 / 384.0).multiply(n3).multiply(epp).multiply(cD).multiply(sinI1).
add(yp5.divide(yp2).multiply(35.0 / 1152.0)).multiply(qyp52.multiply(2.0).multiply(cosI1).
subtract(epp.multiply(cD).multiply(sinI1).multiply(epp2.multiply(2.0).add(3.0))));
dlgs2g = yp2.multiply(48.0).reciprocal().multiply(yp2.multiply(-3.0).multiply(yp2).multiply(qyp22).
add(yp4.multiply(10.0).multiply(qyp42))).
add(n3.divide(yp2).multiply(qyp2_4).divide(24.0));
dlgc3g = yp5.multiply(35.0 / 384.0).divide(yp2).multiply(n3).multiply(epp).multiply(cD).multiply(sinI1).
add(yp5.multiply(35.0 / 1152.0).divide(yp2).multiply(qyp52.multiply(2.0).multiply(cosI1).subtract(epp.multiply(cD).multiply(sinI1).multiply(epp2.multiply(2.0).add(3.0)))));
dlgcg = yp3.negate().multiply(cosI2).multiply(epp).divide(yp2.multiply(sinI1).multiply(4.0)).
add(yp5.divide(yp2).multiply(0.078125).multiply(cC).multiply(cosI2.divide(sinI1).multiply(epp.negate()).multiply(epp2.multiply(3.0).add(4.0)).
add(sinI1.multiply(epp2).multiply(epp2.multiply(9.0).add(26.0)))).
......@@ -768,13 +763,11 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
add(cosI1.multiply(epp).divide(yp2.multiply(sinI1).multiply(4.0)).multiply(yp3.add(yp5.multiply(0.3125).multiply(cC).multiply(epp2.multiply(3.0).add(4.0))))).
add(yp5.multiply(qyp52bis).multiply(1.875).divide(yp2.multiply(4.0)));
// short periodic multiplicative
aC = yp2.negate().multiply(C3c2).multiply(app).divide(n3);
aCbis = y2.multiply(app).multiply(C3c2);
ac2g2f = y2.multiply(app).multiply(sinI2).multiply(3.0);
T qe = y2.multiply(C3c2).multiply(0.5).multiply(n2).divide(n6);
eC = qe.multiply(epp).divide(n3.add(1.0)).multiply(epp2.multiply(epp2.subtract(3.0)).add(3.0));
ecf = qe.multiply(3.0);
......@@ -792,10 +785,9 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
T qi = yp2.multiply(epp).multiply(cosI1).multiply(sinI1);
ide = cosI1.multiply(epp.negate()).divide(sinI1.multiply(n2));
isfs2f2g = qi;
icfc2f2g = qe.multiply(2.0);
icfc2f2g = qi.multiply(2.0);
qi = yp2.multiply(cosI1).multiply(sinI1);
ic2f2g = qe.multiply(1.5);
ic2f2g = qi.multiply(1.5);
T qgl1 = yp2.multiply(0.25);
T qgl2 = yp2.multiply(epp).multiply(n2).multiply(0.25).divide(n.add(1.0));
......@@ -814,7 +806,6 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
add(qgl2.multiply(1.0 / 3.0));
glos2g3f = qgl2;
final T qh = yp2.multiply(cosI1).multiply(3.0);
hf = qh.negate();
hl = qh;
......@@ -823,14 +814,13 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
hs2g2f = yp2.multiply(cosI1).multiply(1.5);
hsfc2g2f = yp2.multiply(cosI1).multiply(epp).negate();
final T qedl = yp2.multiply(n3).multiply(-0.25);
edls2g = qyp2_4.multiply(1.0 / 24.0).multiply(epp).multiply(n3).divide(yp2);
edlcg = yp3.divide(yp2).multiply(-0.25).multiply(n3).multiply(sinI1).
subtract(yp5.divide(yp2).multiply(0.078125).multiply(n3).multiply(sinI1).multiply(cC).multiply(epp2.multiply(9.0).add(4.0)));
edlc3g = yp5.divide(yp2).multiply(n3).multiply(epp2).multiply(cD).multiply(sinI1).multiply(35.0 / 384.0);
edlsf = qedl.multiply(C3c2).multiply(2.0);
edls2gf = qedl.multiply(cosI1.negate().add(1.0)).multiply(3.0);
edls2gf = qedl.multiply(3.0).multiply(cosI2.negate().add(1.0));
edls2g3f = qedl.multiply(1.0 / 3.0);
}
......
......@@ -230,7 +230,7 @@ public class BrouwerLyddanePropagatorTest {
SpacecraftState BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift));
final KeplerianOrbit BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit());
Assert.assertEquals(NumOrbit.getA(), BLOrbit.getA(), 0.073);
Assert.assertEquals(NumOrbit.getA(), BLOrbit.getA(), 0.072);
Assert.assertEquals(NumOrbit.getE(), BLOrbit.getE(), 0.00000028);
Assert.assertEquals(NumOrbit.getI(), BLOrbit.getI(), 0.000004);
Assert.assertEquals(MathUtils.normalizeAngle(NumOrbit.getPerigeeArgument(), FastMath.PI),
......
......@@ -40,7 +40,7 @@ import org.orekit.utils.FieldPVCoordinates;
import org.orekit.utils.IERSConventions;
public class FieldBrouwerLyddanePropagatorTest {
private static final AttitudeProvider DEFAULT_LAW = Utils.defaultLaw();
private static final AttitudeProvider DEFAULT_LAW = Utils.defaultLaw();
@Test
public void sameDateCartesian() {
......@@ -106,12 +106,12 @@ public class FieldBrouwerLyddanePropagatorTest {
Assert.assertEquals(0.0,
FieldVector3D.distance(initialOrbit.getPVCoordinates().getPosition(),
finalOrbit.getPVCoordinates().getPosition()).getReal(),
1.0e-8);
1.1e-8);
Assert.assertEquals(0.0,
FieldVector3D.distance(initialOrbit.getPVCoordinates().getVelocity(),
finalOrbit.getPVCoordinates().getVelocity()).getReal(),
1.0e-11);
1.4e-11);
Assert.assertEquals(0.0, finalOrbit.getA().getReal() - initialOrbit.getA().getReal(), 0.0);
}
......@@ -259,9 +259,9 @@ public class FieldBrouwerLyddanePropagatorTest {
final KeplerianOrbit BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit().toOrbit());
Assert.assertEquals(NumOrbit.getA(), BLOrbit.getA(), 0.073);
Assert.assertEquals(NumOrbit.getA(), BLOrbit.getA(), 0.072);
Assert.assertEquals(NumOrbit.getE(), BLOrbit.getE(), 0.00000028);
Assert.assertEquals(NumOrbit.getI(), BLOrbit.getI(), 0.000005);
Assert.assertEquals(NumOrbit.getI(), BLOrbit.getI(), 0.000004);
Assert.assertEquals(MathUtils.normalizeAngle(NumOrbit.getPerigeeArgument(), FastMath.PI),
MathUtils.normalizeAngle(BLOrbit.getPerigeeArgument(), FastMath.PI), 0.119);
Assert.assertEquals(MathUtils.normalizeAngle(NumOrbit.getRightAscensionOfAscendingNode(), FastMath.PI),
......@@ -340,7 +340,7 @@ public class FieldBrouwerLyddanePropagatorTest {
Assert.assertEquals(NumOrbit.getA(), BLOrbit.getA(), 0.17);
Assert.assertEquals(NumOrbit.getE(), BLOrbit.getE(), 0.00000028);
Assert.assertEquals(NumOrbit.getI(), BLOrbit.getI(), 0.000005);
Assert.assertEquals(NumOrbit.getI(), BLOrbit.getI(), 0.000004);
Assert.assertEquals(MathUtils.normalizeAngle(NumOrbit.getPerigeeArgument(), FastMath.PI),
MathUtils.normalizeAngle(BLOrbit.getPerigeeArgument(), FastMath.PI), 0.197);
Assert.assertEquals(MathUtils.normalizeAngle(NumOrbit.getRightAscensionOfAscendingNode(), FastMath.PI),
......@@ -573,6 +573,58 @@ public class FieldBrouwerLyddanePropagatorTest {
}
}
@Test
public void testMeanComparisonWithNonField() {
doTestMeanComparisonWithNonField(Decimal64Field.getInstance());
}
private <T extends CalculusFieldElement<T>> void doTestMeanComparisonWithNonField(Field<T> field) {
T zero = field.getZero();
FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
final Frame inertialFrame = FramesFactory.getEME2000();
FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
double timeshift = -60000.0;
// Initial orbit
final double a = 24396159; // semi major axis in meters
final double e = 0.9; // eccentricity
final double i = FastMath.toRadians(7); // inclination
final double omega = FastMath.toRadians(180); // perigee argument
final double raan = FastMath.toRadians(261); // right ascention of ascending node
final double lM = FastMath.toRadians(0); // mean anomaly
final FieldOrbit<T> initialOrbit = new FieldKeplerianOrbit<>(zero.add(a), zero.add(e), zero.add(i), zero.add(omega),
zero.add(raan), zero.add(lM), PositionAngle.TRUE,
inertialFrame, initDate, zero.add(provider.getMu()));
// Initial state definition
final FieldSpacecraftState<T> initialStateField = new FieldSpacecraftState<>(initialOrbit);
final SpacecraftState initialState = initialStateField.toSpacecraftState();
// Field propagation
final FieldBrouwerLyddanePropagator<T> blField = new FieldBrouwerLyddanePropagator<>(initialStateField.getOrbit(),
GravityFieldFactory.getUnnormalizedProvider(provider),
PropagationType.MEAN);
final FieldSpacecraftState<T> finalStateField = blField.propagate(initialStateField.getDate().shiftedBy(timeshift));
final FieldKeplerianOrbit<T> finalOrbitField = (FieldKeplerianOrbit<T>) OrbitType.KEPLERIAN.convertType(finalStateField.getOrbit());
final KeplerianOrbit finalOrbitFieldReal = finalOrbitField.toOrbit();
// Classical propagation
final BrouwerLyddanePropagator bl = new BrouwerLyddanePropagator(initialState.getOrbit(),
GravityFieldFactory.getUnnormalizedProvider(provider),
PropagationType.MEAN);
final SpacecraftState finalState = bl.propagate(initialState.getDate().shiftedBy(timeshift));
final KeplerianOrbit finalOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(finalState.getOrbit());
Assert.assertEquals(finalOrbitFieldReal.getA(), finalOrbit.getA(), Double.MIN_VALUE);
Assert.assertEquals(finalOrbitFieldReal.getE(), finalOrbit.getE(), Double.MIN_VALUE);
Assert.assertEquals(finalOrbitFieldReal.getI(), finalOrbit.getI(), Double.MIN_VALUE);
Assert.assertEquals(finalOrbitFieldReal.getRightAscensionOfAscendingNode(), + finalOrbit.getRightAscensionOfAscendingNode(), Double.MIN_VALUE);
Assert.assertEquals(finalOrbitFieldReal.getPerigeeArgument(), finalOrbit.getPerigeeArgument(), Double.MIN_VALUE);
Assert.assertEquals(finalOrbitFieldReal.getMeanAnomaly(), finalOrbit.getMeanAnomaly(), Double.MIN_VALUE);
Assert.assertEquals(0.0, finalOrbitFieldReal.getPVCoordinates().getPosition().distance(finalOrbit.getPVCoordinates().getPosition()), Double.MIN_VALUE);
Assert.assertEquals(0.0, finalOrbitFieldReal.getPVCoordinates().getVelocity().distance(finalOrbit.getPVCoordinates().getVelocity()), Double.MIN_VALUE);
}
@Before
public void setUp() {
......
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