Skip to content
GitLab
Commit 75a6e52f authored by Bryan Cazabonne's avatar Bryan Cazabonne
Browse files

Merge branch 'issue-932' into 'develop' 

Added a way to compute mean parameters in Brouwer-Lyddane model.

Closes #932

See merge request !276
parents 448bd1d0 6ac2b8db
Pipeline #2089 passed with stages
in 14 minutes and 38 seconds
Hide whitespace changes
Inline Side-by-side
src/changes/changes.xml
View file @ 75a6e52f
......@@ -21,6 +21,9 @@
</properties>
<body>
<release version="11.2" date="TBD" description="TBD">
<action dev="bryan" type="add" issue="932">
Added a way to compute mean parameters in Brouwer-Lyddane model.
</action>
<action dev="markrutten" type="add" issue="922">
Added bistatic range measurement.
</action>
......
src/main/java/org/orekit/propagation/analytical/BrouwerLyddanePropagator.java
View file @ 75a6e52f
......@@ -489,6 +489,80 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator {
}
/** Conversion from osculating to mean orbit.
* <p>
* Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
* osculating SpacecraftState in input.
* </p>
* <p>
* Since the osculating orbit is obtained with the computation of
* short-periodic variation, the resulting output will depend on
* both the gravity field parameterized in input and the
* atmospheric drag represented by the {@code m2} parameter.
* </p>
* <p>
* The computation is done through a fixed-point iteration process.
* </p>
* @param osculating osculating orbit to convert
* @param provider for un-normalized zonal coefficients
* @param harmonics {@code provider.onDate(osculating.getDate())}
* @param M2Value value of empirical drag coefficient in rad/s².
* If equal to {@code BrouwerLyddanePropagator.M2} drag is not considered
* @return mean orbit in a Brouwer-Lyddane sense
* @since 11.2
*/
public static KeplerianOrbit computeMeanOrbit(final Orbit osculating,
final UnnormalizedSphericalHarmonicsProvider provider,
final UnnormalizedSphericalHarmonics harmonics,
final double M2Value) {
return computeMeanOrbit(osculating,
provider.getAe(), provider.getMu(),
harmonics.getUnnormalizedCnm(2, 0),
harmonics.getUnnormalizedCnm(3, 0),
harmonics.getUnnormalizedCnm(4, 0),
harmonics.getUnnormalizedCnm(5, 0),
M2Value);
}
/** Conversion from osculating to mean orbit.
* <p>
* Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
* osculating SpacecraftState in input.
* </p>
* <p>
* Since the osculating orbit is obtained with the computation of
* short-periodic variation, the resulting output will depend on
* both the gravity field parameterized in input and the
* atmospheric drag represented by the {@code m2} parameter.
* </p>
* <p>
* The computation is done through a fixed-point iteration process.
* </p>
* @param osculating osculating orbit to convert
* @param referenceRadius reference radius of the Earth for the potential model (m)
* @param mu central attraction coefficient (m³/s²)
* @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
* @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
* @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
* @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
* @param M2Value value of empirical drag coefficient in rad/s².
* If equal to {@code BrouwerLyddanePropagator.M2} drag is not considered
* @return mean orbit in a Brouwer-Lyddane sense
* @since 11.2
*/
public static KeplerianOrbit computeMeanOrbit(final Orbit osculating,
final double referenceRadius, final double mu,
final double c20, final double c30, final double c40,
final double c50, final double M2Value) {
final BrouwerLyddanePropagator propagator =
new BrouwerLyddanePropagator(osculating,
InertialProvider.of(osculating.getFrame()),
DEFAULT_MASS,
referenceRadius, mu, c20, c30, c40, c50,
PropagationType.OSCULATING, M2Value);
return propagator.initialModel.mean;
}
/** {@inheritDoc}
* <p>The new initial state to consider
* must be defined with an osculating orbit.</p>
......@@ -550,21 +624,21 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator {
while (i++ < 200) {
// recompute the osculating parameters from the current mean parameters
final UnivariateDerivative2[] parameters = current.propagateParameters(current.mean.getDate());
final KeplerianOrbit parameters = current.propagateParameters(current.mean.getDate());
// adapted parameters residuals
final double deltaA = osculating.getA() - parameters[0].getValue();
final double deltaE = osculating.getE() - parameters[1].getValue();
final double deltaI = osculating.getI() - parameters[2].getValue();
final double deltaA = osculating.getA() - parameters.getA();
final double deltaE = osculating.getE() - parameters.getE();
final double deltaI = osculating.getI() - parameters.getI();
final double deltaOmega = MathUtils.normalizeAngle(osculating.getPerigeeArgument() -
parameters[3].getValue(),
0.0);
parameters.getPerigeeArgument(),
0.0);
final double deltaRAAN = MathUtils.normalizeAngle(osculating.getRightAscensionOfAscendingNode() -
parameters[4].getValue(),
0.0);
parameters.getRightAscensionOfAscendingNode(),
0.0);
final double deltaAnom = MathUtils.normalizeAngle(osculating.getMeanAnomaly() -
parameters[5].getValue(),
0.0);
parameters.getMeanAnomaly(),
0.0);
// update mean parameters
......@@ -594,14 +668,9 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator {
/** {@inheritDoc} */
public KeplerianOrbit propagateOrbit(final AbsoluteDate date) {
// compute Cartesian parameters, taking derivatives into account
// to make sure velocity and acceleration are consistent
// compute keplerian parameters, taking derivatives into account
final BLModel current = models.get(date);
final UnivariateDerivative2[] propOrb_parameters = current.propagateParameters(date);
return new KeplerianOrbit(propOrb_parameters[0].getValue(), propOrb_parameters[1].getValue(),
propOrb_parameters[2].getValue(), propOrb_parameters[3].getValue(),
propOrb_parameters[4].getValue(), propOrb_parameters[5].getValue(),
PositionAngle.MEAN, current.mean.getFrame(), date, mu);
return current.propagateParameters(date);
}
/**
......@@ -1082,7 +1151,7 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator {
* @param date target date for the orbit
* @return propagated parameters
*/
public UnivariateDerivative2[] propagateParameters(final AbsoluteDate date) {
public KeplerianOrbit propagateParameters(final AbsoluteDate date) {
// Empirical drag coefficient M2
final double m2 = M2Driver.getValue();
......@@ -1240,7 +1309,12 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator {
// Argument of perigee
final UnivariateDerivative2 g = g_p_l.subtract(l);
return new UnivariateDerivative2[] { a, e, i, g, h, l };
// Return a Keplerian orbit
return new KeplerianOrbit(a.getValue(), e.getValue(), i.getValue(),
g.getValue(), h.getValue(), l.getValue(),
a.getFirstDerivative(), e.getFirstDerivative(), i.getFirstDerivative(),
g.getFirstDerivative(), h.getFirstDerivative(), l.getFirstDerivative(),
PositionAngle.MEAN, mean.getFrame(), date, mu);
}
......
src/main/java/org/orekit/propagation/analytical/FieldBrouwerLyddanePropagator.java
View file @ 75a6e52f
......@@ -25,7 +25,6 @@ import org.hipparchus.analysis.differentiation.FieldUnivariateDerivative2;
import org.hipparchus.util.CombinatoricsUtils;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.FieldSinCos;
import org.hipparchus.util.MathArrays;
import org.hipparchus.util.MathUtils;
import org.hipparchus.util.Precision;
import org.orekit.attitudes.AttitudeProvider;
......@@ -130,7 +129,7 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
final UnnormalizedSphericalHarmonicsProvider provider,
final double M2) {
this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()),
initialOrbit.getA().getField().getZero().add(DEFAULT_MASS), provider,
initialOrbit.getMu().newInstance(DEFAULT_MASS), provider,
provider.onDate(initialOrbit.getDate().toAbsoluteDate()), M2);
}
......@@ -153,7 +152,7 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
final UnnormalizedSphericalHarmonicsProvider provider,
final UnnormalizedSphericalHarmonics harmonics,
final double M2) {
this(initialOrbit, attitude, mass, provider.getAe(), initialOrbit.getA().getField().getZero().add(provider.getMu()),
this(initialOrbit, attitude, mass, provider.getAe(), initialOrbit.getMu().newInstance(provider.getMu()),
harmonics.getUnnormalizedCnm(2, 0),
harmonics.getUnnormalizedCnm(3, 0),
harmonics.getUnnormalizedCnm(4, 0),
......@@ -192,7 +191,7 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
final double c20, final double c30, final double c40,
final double c50, final double M2) {
this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()),
initialOrbit.getDate().getField().getZero().add(DEFAULT_MASS),
initialOrbit.getMu().newInstance(DEFAULT_MASS),
referenceRadius, mu, c20, c30, c40, c50, M2);
}
......@@ -262,7 +261,7 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
final AttitudeProvider attitudeProv,
final UnnormalizedSphericalHarmonicsProvider provider,
final double M2) {
this(initialOrbit, attitudeProv, initialOrbit.getA().getField().getZero().add(DEFAULT_MASS), provider,
this(initialOrbit, attitudeProv, initialOrbit.getMu().newInstance(DEFAULT_MASS), provider,
provider.onDate(initialOrbit.getDate().toAbsoluteDate()), M2);
}
......@@ -296,7 +295,7 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
final double referenceRadius, final T mu,
final double c20, final double c30, final double c40,
final double c50, final double M2) {
this(initialOrbit, attitudeProv, initialOrbit.getDate().getField().getZero().add(DEFAULT_MASS),
this(initialOrbit, attitudeProv, initialOrbit.getMu().newInstance(DEFAULT_MASS),
referenceRadius, mu, c20, c30, c40, c50, M2);
}
......@@ -371,7 +370,7 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
final PropagationType initialType,
final double M2) {
this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()),
initialOrbit.getA().getField().getZero().add(DEFAULT_MASS), provider,
initialOrbit.getMu().newInstance(DEFAULT_MASS), provider,
provider.onDate(initialOrbit.getDate().toAbsoluteDate()), initialType, M2);
}
......@@ -415,7 +414,7 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
final UnnormalizedSphericalHarmonics harmonics,
final PropagationType initialType,
final double M2) {
this(initialOrbit, attitude, mass, provider.getAe(), initialOrbit.getA().getField().getZero().add(provider.getMu()),
this(initialOrbit, attitude, mass, provider.getAe(), initialOrbit.getMu().newInstance(provider.getMu()),
harmonics.getUnnormalizedCnm(2, 0),
harmonics.getUnnormalizedCnm(3, 0),
harmonics.getUnnormalizedCnm(4, 0),
......@@ -481,6 +480,82 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
}
/** Conversion from osculating to mean orbit.
* <p>
* Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
* osculating SpacecraftState in input.
* </p>
* <p>
* Since the osculating orbit is obtained with the computation of
* short-periodic variation, the resulting output will depend on
* both the gravity field parameterized in input and the
* atmospheric drag represented by the {@code m2} parameter.
* </p>
* <p>
* The computation is done through a fixed-point iteration process.
* </p>
* @param <T> type of the filed elements
* @param osculating osculating orbit to convert
* @param provider for un-normalized zonal coefficients
* @param harmonics {@code provider.onDate(osculating.getDate())}
* @param M2Value value of empirical drag coefficient in rad/s².
* If equal to {@code BrouwerLyddanePropagator.M2} drag is not considered
* @return mean orbit in a Brouwer-Lyddane sense
* @since 11.2
*/
public static <T extends CalculusFieldElement<T>> FieldKeplerianOrbit<T> computeMeanOrbit(final FieldOrbit<T> osculating,
final UnnormalizedSphericalHarmonicsProvider provider,
final UnnormalizedSphericalHarmonics harmonics,
final double M2Value) {
return computeMeanOrbit(osculating,
provider.getAe(), provider.getMu(),
harmonics.getUnnormalizedCnm(2, 0),
harmonics.getUnnormalizedCnm(3, 0),
harmonics.getUnnormalizedCnm(4, 0),
harmonics.getUnnormalizedCnm(5, 0),
M2Value);
}
/** Conversion from osculating to mean orbit.
* <p>
* Compute mean orbit <b>in a Brouwer-Lyddane sense</b>, corresponding to the
* osculating SpacecraftState in input.
* </p>
* <p>
* Since the osculating orbit is obtained with the computation of
* short-periodic variation, the resulting output will depend on
* both the gravity field parameterized in input and the
* atmospheric drag represented by the {@code m2} parameter.
* </p>
* <p>
* The computation is done through a fixed-point iteration process.
* </p>
* @param <T> type of the filed elements
* @param osculating osculating orbit to convert
* @param referenceRadius reference radius of the Earth for the potential model (m)
* @param mu central attraction coefficient (m³/s²)
* @param c20 un-normalized zonal coefficient (about -1.08e-3 for Earth)
* @param c30 un-normalized zonal coefficient (about +2.53e-6 for Earth)
* @param c40 un-normalized zonal coefficient (about +1.62e-6 for Earth)
* @param c50 un-normalized zonal coefficient (about +2.28e-7 for Earth)
* @param M2Value value of empirical drag coefficient in rad/s².
* If equal to {@code BrouwerLyddanePropagator.M2} drag is not considered
* @return mean orbit in a Brouwer-Lyddane sense
* @since 11.2
*/
public static <T extends CalculusFieldElement<T>> FieldKeplerianOrbit<T> computeMeanOrbit(final FieldOrbit<T> osculating,
final double referenceRadius, final double mu,
final double c20, final double c30, final double c40,
final double c50, final double M2Value) {
final FieldBrouwerLyddanePropagator<T> propagator =
new FieldBrouwerLyddanePropagator<>(osculating,
InertialProvider.of(osculating.getFrame()),
osculating.getMu().newInstance(DEFAULT_MASS),
referenceRadius, osculating.getMu().newInstance(mu),
c20, c30, c40, c50,
PropagationType.OSCULATING, M2Value);
return propagator.initialModel.mean;
}
/** {@inheritDoc}
* <p>The new initial state to consider
......@@ -549,15 +624,15 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
while (i++ < 200) {
// recompute the osculating parameters from the current mean parameters
final FieldUnivariateDerivative2<T>[] parameters = current.propagateParameters(current.mean.getDate(), getParameters(mass.getField()));
final FieldKeplerianOrbit<T> parameters = current.propagateParameters(current.mean.getDate(), getParameters(mass.getField()));
// adapted parameters residuals
final T deltaA = osculating.getA() .subtract(parameters[0].getValue());
final T deltaE = osculating.getE() .subtract(parameters[1].getValue());
final T deltaI = osculating.getI() .subtract(parameters[2].getValue());
final T deltaOmega = MathUtils.normalizeAngle(osculating.getPerigeeArgument().subtract(parameters[3].getValue()), zero);
final T deltaRAAN = MathUtils.normalizeAngle(osculating.getRightAscensionOfAscendingNode().subtract(parameters[4].getValue()), zero);
final T deltaAnom = MathUtils.normalizeAngle(osculating.getMeanAnomaly().subtract(parameters[5].getValue()), zero);
final T deltaA = osculating.getA() .subtract(parameters.getA());
final T deltaE = osculating.getE() .subtract(parameters.getE());
final T deltaI = osculating.getI() .subtract(parameters.getI());
final T deltaOmega = MathUtils.normalizeAngle(osculating.getPerigeeArgument().subtract(parameters.getPerigeeArgument()), zero);
final T deltaRAAN = MathUtils.normalizeAngle(osculating.getRightAscensionOfAscendingNode().subtract(parameters.getRightAscensionOfAscendingNode()), zero);
final T deltaAnom = MathUtils.normalizeAngle(osculating.getMeanAnomaly().subtract(parameters.getMeanAnomaly()), zero);
// update mean parameters
......@@ -588,13 +663,8 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
/** {@inheritDoc} */
public FieldKeplerianOrbit<T> propagateOrbit(final FieldAbsoluteDate<T> date, final T[] parameters) {
// compute Cartesian parameters, taking derivatives into account
// to make sure velocity and acceleration are consistent
final FieldBLModel<T> current = models.get(date);
final FieldUnivariateDerivative2<T>[] propOrb_parameters = current.propagateParameters(date, parameters);
return new FieldKeplerianOrbit<T>(propOrb_parameters[0].getValue(), propOrb_parameters[1].getValue(),
propOrb_parameters[2].getValue(), propOrb_parameters[3].getValue(),
propOrb_parameters[4].getValue(), propOrb_parameters[5].getValue(),
PositionAngle.MEAN, current.mean.getFrame(), date, mu);
return current.propagateParameters(date, parameters);
}
/**
......@@ -614,6 +684,9 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
/** Constant mass. */
private final T mass;
/** Central attraction coefficient. */
private final T mu;
// CHECKSTYLE: stop JavadocVariable check
// preprocessed values
......@@ -712,6 +785,7 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
this.mean = mean;
this.mass = mass;
this.mu = mu;
final T one = mass.getField().getOne();
final T app = mean.getA();
......@@ -1036,7 +1110,7 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
* @param parameters model parameters
* @return propagated parameters
*/
public FieldUnivariateDerivative2<T>[] propagateParameters(final FieldAbsoluteDate<T> date, final T[] parameters) {
public FieldKeplerianOrbit<T> propagateParameters(final FieldAbsoluteDate<T> date, final T[] parameters) {
// Field
final Field<T> field = date.getField();
......@@ -1199,14 +1273,13 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> ex
// Argument of perigee
final FieldUnivariateDerivative2<T> g = g_p_l.subtract(l);
final FieldUnivariateDerivative2<T>[] FTD = MathArrays.buildArray(g.getField(), 6);
FTD[0] = a;
FTD[1] = e;
FTD[2] = i;
FTD[3] = g;
FTD[4] = h;
FTD[5] = l;
return FTD;
// Return a keplerian orbit
return new FieldKeplerianOrbit<>(a.getValue(), e.getValue(), i.getValue(),
g.getValue(), h.getValue(), l.getValue(),
a.getFirstDerivative(), e.getFirstDerivative(), i.getFirstDerivative(),
g.getFirstDerivative(), h.getFirstDerivative(), l.getFirstDerivative(),
PositionAngle.MEAN, mean.getFrame(), date, this.mu);
}
}
......
src/test/java/org/orekit/estimation/leastsquares/BrouwerLyddaneBatchLSEstimatorTest.java
View file @ 75a6e52f
......@@ -252,7 +252,7 @@ public class BrouwerLyddaneBatchLSEstimatorTest {
BrouwerLyddaneEstimationTestUtils.checkFit(context, estimator, 1, 2,
0.0, 7.9e-8,
0.0, 1.1e-7,
0.0, 1.4e-5,
0.0, 1.5e-5,
0.0, 1.4e-8);
}
......
src/test/java/org/orekit/estimation/sequential/BrouwerLyddaneKalmanEstimatorTest.java
View file @ 75a6e52f
......@@ -83,7 +83,7 @@ public class BrouwerLyddaneKalmanEstimatorTest {
// Filter the measurements and check the results
final double expectedDeltaPos = 0.;
final double posEps = 2.33e-8;
final double posEps = 2.70e-8;
final double expectedDeltaVel = 0.;
final double velEps = 6.59e-11;
final double[] expectedsigmasPos = {0.998894, 0.933798, 0.997346};
......
src/test/java/org/orekit/propagation/analytical/BrouwerLyddaneParametersDerivativesTest.java
View file @ 75a6e52f
......@@ -120,7 +120,7 @@ public class BrouwerLyddaneParametersDerivativesTest {
sM4h, sM3h, sM2h, sM1h, sP1h, sP2h, sP3h, sP4h);
for (int i = 0; i < 6; ++i) {
Assert.assertEquals(0.0, (dYdPRef[i][0] - dYdP.getEntry(i, 0)) / dYdPRef[i][0], 8.49e-13);
Assert.assertEquals(0.0, (dYdPRef[i][0] - dYdP.getEntry(i, 0)) / dYdPRef[i][0], 8.78e-13);
}
}
......
src/test/java/org/orekit/propagation/analytical/BrouwerLyddanePropagatorTest.java
View file @ 75a6e52f
......@@ -5,6 +5,9 @@ import static org.junit.Assert.assertTrue;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
import org.hipparchus.stat.descriptive.StorelessUnivariateStatistic;
import org.hipparchus.stat.descriptive.rank.Max;
import org.hipparchus.stat.descriptive.rank.Min;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathUtils;
import org.junit.After;
......@@ -25,6 +28,7 @@ import org.orekit.forces.gravity.potential.GravityFieldFactory;
import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
import org.orekit.forces.gravity.potential.TideSystem;
import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider;
import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.models.earth.atmosphere.DTM2000;
......@@ -73,12 +77,12 @@ public class BrouwerLyddanePropagatorTest {
Assert.assertEquals(0.0,
Vector3D.distance(initialOrbit.getPVCoordinates().getPosition(),
finalOrbit.getPVCoordinates().getPosition()),
1.0e-8);
3.6e-9);
Assert.assertEquals(0.0,
Vector3D.distance(initialOrbit.getPVCoordinates().getVelocity(),
finalOrbit.getPVCoordinates().getVelocity()),
1.0e-11);
3.8e-12);
Assert.assertEquals(0.0, finalOrbit.getA() - initialOrbit.getA(), 0.0);
}
......@@ -102,12 +106,12 @@ public class BrouwerLyddanePropagatorTest {
Assert.assertEquals(0.0,
Vector3D.distance(initialOrbit.getPVCoordinates().getPosition(),
finalOrbit.getPVCoordinates().getPosition()),
1.0e-8);
7.0e-9);
Assert.assertEquals(0.0,
Vector3D.distance(initialOrbit.getPVCoordinates().getVelocity(),
finalOrbit.getPVCoordinates().getVelocity()),
1.0e-11);
7.0e-12);
Assert.assertEquals(0.0, finalOrbit.getA() - initialOrbit.getA(), 0.0);
}
......@@ -242,15 +246,15 @@ 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.072);
Assert.assertEquals(NumOrbit.getA(), BLOrbit.getA(), 0.070);
Assert.assertEquals(NumOrbit.getE(), BLOrbit.getE(), 0.00000028);
Assert.assertEquals(NumOrbit.getI(), BLOrbit.getI(), 0.000004);
Assert.assertEquals(NumOrbit.getI(), BLOrbit.getI(), 0.000000053);
Assert.assertEquals(MathUtils.normalizeAngle(NumOrbit.getPerigeeArgument(), FastMath.PI),
MathUtils.normalizeAngle(BLOrbit.getPerigeeArgument(), FastMath.PI), 0.119);
MathUtils.normalizeAngle(BLOrbit.getPerigeeArgument(), FastMath.PI), 0.0021);
Assert.assertEquals(MathUtils.normalizeAngle(NumOrbit.getRightAscensionOfAscendingNode(), FastMath.PI),
MathUtils.normalizeAngle(BLOrbit.getRightAscensionOfAscendingNode(), FastMath.PI), 0.000072);
MathUtils.normalizeAngle(BLOrbit.getRightAscensionOfAscendingNode(), FastMath.PI), 0.0000013);
Assert.assertEquals(MathUtils.normalizeAngle(NumOrbit.getTrueAnomaly(), FastMath.PI),
MathUtils.normalizeAngle(BLOrbit.getTrueAnomaly(), FastMath.PI), 0.12);
MathUtils.normalizeAngle(BLOrbit.getTrueAnomaly(), FastMath.PI), 0.0021);
}
@Test
......@@ -666,6 +670,45 @@ public class BrouwerLyddanePropagatorTest {
}
@Test
public void testMeanOrbit() {
final KeplerianOrbit initialOsculating =
new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngle.TRUE,
FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH,
provider.getMu());
final UnnormalizedSphericalHarmonicsProvider ushp = GravityFieldFactory.getUnnormalizedProvider(provider);
final UnnormalizedSphericalHarmonics ush = ushp.onDate(initialOsculating.getDate());
// set up a reference numerical propagator starting for the specified start orbit
// using the same force models (i.e. the first few zonal terms)
double[][] tol = NumericalPropagator.tolerances(0.1, initialOsculating, OrbitType.KEPLERIAN);
AdaptiveStepsizeIntegrator integrator = new DormandPrince853Integrator(0.001, 1000, tol[0], tol[1]);
integrator.setInitialStepSize(60);
NumericalPropagator num = new NumericalPropagator(integrator);
Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
num.addForceModel(new HolmesFeatherstoneAttractionModel(itrf, provider));
num.setInitialState(new SpacecraftState(initialOsculating));
num.setOrbitType(OrbitType.KEPLERIAN);
final StorelessUnivariateStatistic oscMin = new Min();
final StorelessUnivariateStatistic oscMax = new Max();
final StorelessUnivariateStatistic meanMin = new Min();
final StorelessUnivariateStatistic meanMax = new Max();
num.getMultiplexer().add(60, state -> {
final Orbit osc = state.getOrbit();
oscMin.increment(osc.getA());
oscMax.increment(osc.getA());
// compute mean orbit at current date (this is what we test)
final Orbit mean = BrouwerLyddanePropagator.computeMeanOrbit(state.getOrbit(), ushp, ush, BrouwerLyddanePropagator.M2);
meanMin.increment(mean.getA());
meanMax.increment(mean.getA());
});
num.propagate(initialOsculating.getDate().shiftedBy(Constants.JULIAN_DAY));
Assert.assertEquals(3188.347, oscMax.getResult() - oscMin.getResult(), 1.0e-3);
Assert.assertEquals( 55.540, meanMax.getResult() - meanMin.getResult(), 1.0e-3);
}
@Before
public void setUp() {
Utils.setDataRoot("regular-data:atmosphere:potential/icgem-format");
......
src/test/java/org/orekit/propagation/analytical/BrouwerLyddaneStateTransitionMatrixTest.java
View file @ 75a6e52f
......@@ -109,7 +109,7 @@ public class BrouwerLyddaneStateTransitionMatrixTest {
for (int j = 0; j < 6; ++j) {
if (stateVector[i] != 0) {
double error = FastMath.abs((dYdY0.getEntry(i, j) - dYdY0Ref[i][j]) / stateVector[i]) * steps[j];
Assert.assertEquals(0, error, 5.15e-14);
Assert.assertEquals(0, error, 1.42e-13);
}
}
}
......