diff --git a/src/main/java/org/orekit/propagation/analytical/EcksteinHechlerPropagator.java b/src/main/java/org/orekit/propagation/analytical/EcksteinHechlerPropagator.java
index 4ce46c82a4ea8254042eb84c5e4db17882e59519..b5169642509b263414fc7c3fb242a20ee8704fff 100644
--- a/src/main/java/org/orekit/propagation/analytical/EcksteinHechlerPropagator.java
+++ b/src/main/java/org/orekit/propagation/analytical/EcksteinHechlerPropagator.java
@@ -442,6 +442,74 @@ public class EcksteinHechlerPropagator extends AbstractAnalyticalPropagator {
}
+ /** Conversion from osculating to mean orbit.
+ *
+ * Compute mean orbit in a Eckstein-Hechler sense, corresponding to the
+ * osculating SpacecraftState in input.
+ *
+ *
+ * Since the osculating orbit is obtained with the computation of
+ * short-periodic variation, the resulting output will depend on
+ * the gravity field parameterized in input.
+ *
+ *
+ * The computation is done through a fixed-point iteration process.
+ *
+ * @param osculating osculating orbit to convert
+ * @param provider for un-normalized zonal coefficients
+ * @param harmonics {@code provider.onDate(osculating.getDate())}
+ * @return mean orbit in a Eckstein-Hechler sense
+ * @since 11.2
+ */
+ public static CircularOrbit computeMeanOrbit(final Orbit osculating,
+ final UnnormalizedSphericalHarmonicsProvider provider,
+ final UnnormalizedSphericalHarmonics harmonics) {
+ return computeMeanOrbit(osculating,
+ provider.getAe(), provider.getMu(),
+ harmonics.getUnnormalizedCnm(2, 0),
+ harmonics.getUnnormalizedCnm(3, 0),
+ harmonics.getUnnormalizedCnm(4, 0),
+ harmonics.getUnnormalizedCnm(5, 0),
+ harmonics.getUnnormalizedCnm(6, 0));
+ }
+
+ /** Conversion from osculating to mean orbit.
+ *
+ * Compute mean orbit in a Eckstein-Hechler sense, corresponding to the
+ * osculating SpacecraftState in input.
+ *
+ *
+ * Since the osculating orbit is obtained with the computation of
+ * short-periodic variation, the resulting output will depend on
+ * the gravity field parameterized in input.
+ *
+ *
+ * The computation is done through a fixed-point iteration process.
+ *
+ * @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 c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
+ * @return mean orbit in a Eckstein-Hechler sense
+ * @since 11.2
+ */
+ public static CircularOrbit computeMeanOrbit(final Orbit osculating,
+ final double referenceRadius, final double mu,
+ final double c20, final double c30, final double c40,
+ final double c50, final double c60) {
+ final EcksteinHechlerPropagator propagator =
+ new EcksteinHechlerPropagator(osculating,
+ InertialProvider.of(osculating.getFrame()),
+ DEFAULT_MASS,
+ referenceRadius, mu, c20, c30, c40, c50, c60,
+ PropagationType.OSCULATING);
+ return propagator.initialModel.mean;
+ }
+
/** {@inheritDoc}
* The new initial state to consider
* must be defined with an osculating orbit.
@@ -487,7 +555,7 @@ public class EcksteinHechlerPropagator extends AbstractAnalyticalPropagator {
// sanity check
if (osculating.getA() < referenceRadius) {
throw new OrekitException(OrekitMessages.TRAJECTORY_INSIDE_BRILLOUIN_SPHERE,
- osculating.getA());
+ osculating.getA());
}
// rough initialization of the mean parameters
diff --git a/src/main/java/org/orekit/propagation/analytical/FieldEcksteinHechlerPropagator.java b/src/main/java/org/orekit/propagation/analytical/FieldEcksteinHechlerPropagator.java
index 836d26afe0c58d0bc6dbc12e8427b59882d19272..5cb1e9b0d258ad17c1161c904f825620b14132f2 100644
--- a/src/main/java/org/orekit/propagation/analytical/FieldEcksteinHechlerPropagator.java
+++ b/src/main/java/org/orekit/propagation/analytical/FieldEcksteinHechlerPropagator.java
@@ -86,7 +86,7 @@ public class FieldEcksteinHechlerPropagator> e
public FieldEcksteinHechlerPropagator(final FieldOrbit initialOrbit,
final UnnormalizedSphericalHarmonicsProvider provider) {
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()));
}
@@ -106,7 +106,7 @@ public class FieldEcksteinHechlerPropagator> e
final T mass,
final UnnormalizedSphericalHarmonicsProvider provider,
final UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics harmonics) {
- 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),
@@ -145,7 +145,7 @@ public class FieldEcksteinHechlerPropagator> e
final double c20, final double c30, final double c40,
final double c50, final double c60) {
this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()),
- initialOrbit.getDate().getField().getZero().add(DEFAULT_MASS),
+ initialOrbit.getMu().newInstance(DEFAULT_MASS),
referenceRadius, mu, c20, c30, c40, c50, c60);
}
@@ -210,7 +210,7 @@ public class FieldEcksteinHechlerPropagator> e
public FieldEcksteinHechlerPropagator(final FieldOrbit initialOrbit,
final AttitudeProvider attitudeProv,
final UnnormalizedSphericalHarmonicsProvider provider) {
- this(initialOrbit, attitudeProv, initialOrbit.getA().getField().getZero().add(DEFAULT_MASS), provider, provider.onDate(initialOrbit.getDate().toAbsoluteDate()));
+ this(initialOrbit, attitudeProv, initialOrbit.getMu().newInstance(DEFAULT_MASS), provider, provider.onDate(initialOrbit.getDate().toAbsoluteDate()));
}
/** Build a propagator from FieldOrbit, attitude provider and potential.
@@ -242,7 +242,7 @@ public class FieldEcksteinHechlerPropagator> e
final double referenceRadius, final T mu,
final double c20, final double c30, final double c40,
final double c50, final double c60) {
- this(initialOrbit, attitudeProv, initialOrbit.getDate().getField().getZero().add(DEFAULT_MASS),
+ this(initialOrbit, attitudeProv, initialOrbit.getMu().newInstance(DEFAULT_MASS),
referenceRadius, mu, c20, c30, c40, c50, c60);
}
@@ -312,7 +312,7 @@ public class FieldEcksteinHechlerPropagator> e
final UnnormalizedSphericalHarmonicsProvider provider,
final PropagationType initialType) {
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);
}
@@ -352,7 +352,7 @@ public class FieldEcksteinHechlerPropagator> e
final UnnormalizedSphericalHarmonicsProvider provider,
final UnnormalizedSphericalHarmonics harmonics,
final PropagationType initialType) {
- 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),
@@ -420,6 +420,77 @@ public class FieldEcksteinHechlerPropagator> e
}
}
+ /** Conversion from osculating to mean orbit.
+ *
+ * Compute mean orbit in a Eckstein-Hechler sense, corresponding to the
+ * osculating SpacecraftState in input.
+ *
+ *
+ * Since the osculating orbit is obtained with the computation of
+ * short-periodic variation, the resulting output will depend on
+ * the gravity field parameterized in input.
+ *
+ *
+ * The computation is done through a fixed-point iteration process.
+ *
+ * @param 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())}
+ * @return mean orbit in a Eckstein-Hechler sense
+ * @since 11.2
+ */
+ public static > FieldCircularOrbit computeMeanOrbit(final FieldOrbit osculating,
+ final UnnormalizedSphericalHarmonicsProvider provider,
+ final UnnormalizedSphericalHarmonics harmonics) {
+ return computeMeanOrbit(osculating,
+ provider.getAe(), provider.getMu(),
+ harmonics.getUnnormalizedCnm(2, 0),
+ harmonics.getUnnormalizedCnm(3, 0),
+ harmonics.getUnnormalizedCnm(4, 0),
+ harmonics.getUnnormalizedCnm(5, 0),
+ harmonics.getUnnormalizedCnm(6, 0));
+ }
+
+ /** Conversion from osculating to mean orbit.
+ *
+ * Compute mean orbit in a Eckstein-Hechler sense, corresponding to the
+ * osculating SpacecraftState in input.
+ *
+ *
+ * Since the osculating orbit is obtained with the computation of
+ * short-periodic variation, the resulting output will depend on
+ * the gravity field parameterized in input.
+ *
+ *
+ * The computation is done through a fixed-point iteration process.
+ *
+ * @param 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 c60 un-normalized zonal coefficient (about -5.41e-7 for Earth)
+ * @return mean orbit in a Eckstein-Hechler sense
+ * @since 11.2
+ */
+ public static > FieldCircularOrbit computeMeanOrbit(final FieldOrbit osculating,
+ final double referenceRadius, final double mu,
+ final double c20, final double c30, final double c40,
+ final double c50, final double c60) {
+ final FieldEcksteinHechlerPropagator propagator =
+ new FieldEcksteinHechlerPropagator<>(osculating,
+ InertialProvider.of(osculating.getFrame()),
+ osculating.getMu().newInstance(DEFAULT_MASS),
+ referenceRadius, osculating.getMu().newInstance(mu),
+ c20, c30, c40, c50, c60,
+ PropagationType.OSCULATING);
+ return propagator.initialModel.mean;
+ }
+
/** {@inheritDoc}
* The new initial state to consider
* must be defined with an osculating orbit.
diff --git a/src/test/java/org/orekit/propagation/analytical/EcksteinHechlerPropagatorTest.java b/src/test/java/org/orekit/propagation/analytical/EcksteinHechlerPropagatorTest.java
index a635178ed28a9ff3c81557eb978ce0130425001e..c27eafec26e21d703eb681b93b54be12f5b96782 100644
--- a/src/test/java/org/orekit/propagation/analytical/EcksteinHechlerPropagatorTest.java
+++ b/src/test/java/org/orekit/propagation/analytical/EcksteinHechlerPropagatorTest.java
@@ -30,6 +30,9 @@ import org.hipparchus.geometry.euclidean.threed.RotationOrder;
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;
@@ -48,6 +51,7 @@ import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
import org.orekit.forces.gravity.potential.GravityFieldFactory;
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.frames.LOFType;
@@ -82,6 +86,7 @@ import org.orekit.time.FieldAbsoluteDate;
import org.orekit.time.TimeComponents;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.CartesianDerivativesFilter;
+import org.orekit.utils.Constants;
import org.orekit.utils.FieldPVCoordinatesProvider;
import org.orekit.utils.IERSConventions;
import org.orekit.utils.PVCoordinates;
@@ -504,7 +509,7 @@ public class EcksteinHechlerPropagatorTest {
public void hyperbolic() {
KeplerianOrbit hyperbolic =
new KeplerianOrbit(-1.0e10, 2, 0, 0, 0, 0, PositionAngle.TRUE,
- FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, 3.986004415e14);
+ FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, provider.getMu());
EcksteinHechlerPropagator propagator =
new EcksteinHechlerPropagator(hyperbolic, provider);
propagator.propagate(AbsoluteDate.J2000_EPOCH.shiftedBy(10.0));
@@ -514,7 +519,7 @@ public class EcksteinHechlerPropagatorTest {
public void wrongAttitude() {
KeplerianOrbit orbit =
new KeplerianOrbit(1.0e10, 1.0e-4, 1.0e-2, 0, 0, 0, PositionAngle.TRUE,
- FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, 3.986004415e14);
+ FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, provider.getMu());
AttitudeProvider wrongLaw = new AttitudeProvider() {
public Attitude getAttitude(PVCoordinatesProvider pvProv, AbsoluteDate date, Frame frame) {
throw new OrekitException(new DummyLocalizable("gasp"), new RuntimeException());
@@ -614,7 +619,7 @@ public class EcksteinHechlerPropagatorTest {
public void stopAtTargetDate() {
final KeplerianOrbit orbit =
new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngle.TRUE,
- FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, 3.986004415e14);
+ FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, provider.getMu());
EcksteinHechlerPropagator propagator =
new EcksteinHechlerPropagator(orbit, provider);
Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
@@ -644,7 +649,7 @@ public class EcksteinHechlerPropagatorTest {
public void date() {
final KeplerianOrbit orbit =
new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngle.TRUE,
- FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, 3.986004415e14);
+ FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, provider.getMu());
EcksteinHechlerPropagator propagator =
new EcksteinHechlerPropagator(orbit, provider);
final AbsoluteDate stopDate = AbsoluteDate.J2000_EPOCH.shiftedBy(500.0);
@@ -658,7 +663,7 @@ public class EcksteinHechlerPropagatorTest {
public void fixedStep() {
final KeplerianOrbit orbit =
new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngle.TRUE,
- FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, 3.986004415e14);
+ FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, provider.getMu());
EcksteinHechlerPropagator propagator =
new EcksteinHechlerPropagator(orbit, provider);
final double step = 100.0;
@@ -679,7 +684,7 @@ public class EcksteinHechlerPropagatorTest {
public void setting() {
final KeplerianOrbit orbit =
new KeplerianOrbit(7.8e6, 0.032, 0.4, 0.1, 0.2, 0.3, PositionAngle.TRUE,
- FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, 3.986004415e14);
+ FramesFactory.getEME2000(), AbsoluteDate.J2000_EPOCH, provider.getMu());
EcksteinHechlerPropagator propagator =
new EcksteinHechlerPropagator(orbit, provider);
final OneAxisEllipsoid earthShape =
@@ -778,7 +783,7 @@ public class EcksteinHechlerPropagatorTest {
final SpacecraftState initialState = new SpacecraftState(new EquinoctialOrbit(new PVCoordinates(position, velocity),
FramesFactory.getEME2000(),
initDate,
- 3.986004415E14));
+ provider.getMu()));
// Mean state computation
final List models = new ArrayList<>();
@@ -814,7 +819,7 @@ public class EcksteinHechlerPropagatorTest {
final SpacecraftState initialState = new SpacecraftState(new EquinoctialOrbit(new PVCoordinates(position, velocity),
FramesFactory.getEME2000(),
initDate,
- 3.986004415E14));
+ provider.getMu()));
// Mean state computation
final List models = new ArrayList<>();
@@ -841,6 +846,44 @@ public class EcksteinHechlerPropagatorTest {
4.2e-2);
}
+ @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 UnnormalizedSphericalHarmonics ush = provider.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.CIRCULAR);
+ 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, GravityFieldFactory.getNormalizedProvider(provider)));
+ num.setInitialState(new SpacecraftState(initialOsculating));
+ num.setOrbitType(OrbitType.CIRCULAR);
+ 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 = EcksteinHechlerPropagator.computeMeanOrbit(state.getOrbit(), provider, ush);
+ meanMin.increment(mean.getA());
+ meanMax.increment(mean.getA());
+ });
+ num.propagate(initialOsculating.getDate().shiftedBy(Constants.JULIAN_DAY));
+
+ Assert.assertEquals(3190.029, oscMax.getResult() - oscMin.getResult(), 1.0e-3);
+ Assert.assertEquals( 49.638, meanMax.getResult() - meanMin.getResult(), 1.0e-3);
+
+ }
+
private static double tangLEmLv(double Lv, double ex, double ey) {
// tan ((LE - Lv) /2)) =
return (ey * FastMath.cos(Lv) - ex * FastMath.sin(Lv))
diff --git a/src/test/java/org/orekit/propagation/analytical/FieldEcksteinHechlerPropagatorTest.java b/src/test/java/org/orekit/propagation/analytical/FieldEcksteinHechlerPropagatorTest.java
index e734d370509a25e42ab9e286b5eba316aac9f0a9..673878dba0218f4f43d2bae4ec9dd3b8e21db961 100644
--- a/src/test/java/org/orekit/propagation/analytical/FieldEcksteinHechlerPropagatorTest.java
+++ b/src/test/java/org/orekit/propagation/analytical/FieldEcksteinHechlerPropagatorTest.java
@@ -17,16 +17,22 @@
package org.orekit.propagation.analytical;
+import java.io.IOException;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.Collection;
import java.util.List;
-import org.hipparchus.Field;
import org.hipparchus.CalculusFieldElement;
+import org.hipparchus.Field;
import org.hipparchus.exception.DummyLocalizable;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
import org.hipparchus.geometry.euclidean.threed.RotationOrder;
+import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
+import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
+import org.hipparchus.stat.descriptive.StorelessUnivariateStatistic;
+import org.hipparchus.stat.descriptive.rank.Max;
+import org.hipparchus.stat.descriptive.rank.Min;
import org.hipparchus.util.Decimal64Field;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathUtils;
@@ -43,9 +49,11 @@ import org.orekit.bodies.GeodeticPoint;
import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;
+import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
import org.orekit.forces.gravity.potential.GravityFieldFactory;
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.frames.LOFType;
@@ -65,6 +73,7 @@ import org.orekit.propagation.events.FieldElevationDetector;
import org.orekit.propagation.events.FieldEventDetector;
import org.orekit.propagation.events.FieldNodeDetector;
import org.orekit.propagation.events.handlers.FieldContinueOnEvent;
+import org.orekit.propagation.numerical.FieldNumericalPropagator;
import org.orekit.propagation.sampling.FieldOrekitFixedStepHandler;
import org.orekit.propagation.semianalytical.dsst.FieldDSSTPropagator;
import org.orekit.propagation.semianalytical.dsst.forces.DSSTForceModel;
@@ -75,6 +84,7 @@ import org.orekit.time.FieldAbsoluteDate;
import org.orekit.time.TimeComponents;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.CartesianDerivativesFilter;
+import org.orekit.utils.Constants;
import org.orekit.utils.FieldPVCoordinates;
import org.orekit.utils.FieldPVCoordinatesProvider;
import org.orekit.utils.IERSConventions;
@@ -598,7 +608,7 @@ public class FieldEcksteinHechlerPropagatorTest {
FieldAbsoluteDate date = new FieldAbsoluteDate<>(field);
FieldKeplerianOrbit hyperbolic =
new FieldKeplerianOrbit<>(zero.add(-1.0e10), zero.add(2), zero, zero, zero, zero, PositionAngle.TRUE,
- FramesFactory.getEME2000(), date, zero.add(3.986004415e14));
+ FramesFactory.getEME2000(), date, zero.add(provider.getMu()));
try {
FieldEcksteinHechlerPropagator propagator =
new FieldEcksteinHechlerPropagator<>(hyperbolic, provider);
@@ -619,7 +629,7 @@ public class FieldEcksteinHechlerPropagatorTest {
FieldAbsoluteDate date = new FieldAbsoluteDate<>(field);
FieldKeplerianOrbit orbit =
new FieldKeplerianOrbit<>(zero.add(1.0e10), zero.add(1.0e-4), zero.add(1.0e-2), zero, zero, zero, PositionAngle.TRUE,
- FramesFactory.getEME2000(), date, zero.add(3.986004415e14));
+ FramesFactory.getEME2000(), date, zero.add(provider.getMu()));
final DummyLocalizable gasp = new DummyLocalizable("gasp");
AttitudeProvider wrongLaw = new AttitudeProvider() {
@@ -754,7 +764,7 @@ public class FieldEcksteinHechlerPropagatorTest {
FieldAbsoluteDate date = new FieldAbsoluteDate<>(field);
final FieldKeplerianOrbit orbit =
new FieldKeplerianOrbit<>(zero.add(7.8e6), zero.add(0.032), zero.add(0.4), zero.add(0.1), zero.add(0.2), zero.add(0.3), PositionAngle.TRUE,
- FramesFactory.getEME2000(), date, zero.add(3.986004415e14));
+ FramesFactory.getEME2000(), date, zero.add(provider.getMu()));
FieldEcksteinHechlerPropagator propagator =
new FieldEcksteinHechlerPropagator<>(orbit, provider);
Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
@@ -798,7 +808,7 @@ public class FieldEcksteinHechlerPropagatorTest {
FieldAbsoluteDate date = new FieldAbsoluteDate<>(field);
final FieldKeplerianOrbit orbit =
new FieldKeplerianOrbit<>(zero.add(7.8e6), zero.add(0.032), zero.add(0.4), zero.add(0.1), zero.add(0.2), zero.add(0.3), PositionAngle.TRUE,
- FramesFactory.getEME2000(), date, zero.add(3.986004415e14));
+ FramesFactory.getEME2000(), date, zero.add(provider.getMu()));
FieldEcksteinHechlerPropagator propagator =
new FieldEcksteinHechlerPropagator<>(orbit, provider);
final FieldAbsoluteDate stopDate = date.shiftedBy(500.0);
@@ -821,7 +831,7 @@ public class FieldEcksteinHechlerPropagatorTest {
FieldAbsoluteDate date = new FieldAbsoluteDate<>(field);
final FieldKeplerianOrbit orbit =
new FieldKeplerianOrbit<>(zero.add(7.8e6), zero.add(0.032), zero.add(0.4), zero.add(0.1), zero.add(0.2), zero.add(0.3), PositionAngle.TRUE,
- FramesFactory.getEME2000(), date, zero.add(3.986004415e14));
+ FramesFactory.getEME2000(), date, zero.add(provider.getMu()));
FieldEcksteinHechlerPropagator propagator =
new FieldEcksteinHechlerPropagator<>(orbit, provider);
final T step = zero.add(100.0);
@@ -851,7 +861,7 @@ public class FieldEcksteinHechlerPropagatorTest {
FieldAbsoluteDate date = new FieldAbsoluteDate<>(field);
final FieldKeplerianOrbit orbit =
new FieldKeplerianOrbit<>(zero.add(7.8e6), zero.add(0.032), zero.add(0.4), zero.add(0.1), zero.add(0.2), zero.add(0.3), PositionAngle.TRUE,
- FramesFactory.getEME2000(), date, zero.add(3.986004415e14));
+ FramesFactory.getEME2000(), date, zero.add(provider.getMu()));
FieldEcksteinHechlerPropagator propagator =
new FieldEcksteinHechlerPropagator<>(orbit, provider);
final OneAxisEllipsoid earthShape =
@@ -890,7 +900,7 @@ public class FieldEcksteinHechlerPropagatorTest {
final FieldSpacecraftState initialState = new FieldSpacecraftState<>(new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
FramesFactory.getEME2000(),
initDate,
- zero.add(3.986004415E14)));
+ zero.add(provider.getMu())));
// Mean state computation
final List models = new ArrayList<>();
@@ -931,7 +941,7 @@ public class FieldEcksteinHechlerPropagatorTest {
final FieldSpacecraftState initialState = new FieldSpacecraftState<>(new FieldEquinoctialOrbit<>(new FieldPVCoordinates<>(position, velocity),
FramesFactory.getEME2000(),
initDate,
- zero.add(3.986004415E14)));
+ zero.add(provider.getMu())));
// Mean state computation
final List models = new ArrayList<>();
@@ -958,6 +968,51 @@ public class FieldEcksteinHechlerPropagatorTest {
4.2e-2);
}
+ @Test
+ public void testMeanOrbit() throws IOException {
+ doTestMeanOrbit(Decimal64Field.getInstance());
+ }
+
+ private > void doTestMeanOrbit(Field field) {
+ T zero = field.getZero();
+ FieldAbsoluteDate date = new FieldAbsoluteDate<>(field);
+ final FieldKeplerianOrbit initialOsculating =
+ new FieldKeplerianOrbit<>(zero.newInstance(7.8e6), zero.newInstance(0.032), zero.newInstance(0.4),
+ zero.newInstance(0.1), zero.newInstance(0.2), zero.newInstance(0.3),
+ PositionAngle.TRUE,
+ FramesFactory.getEME2000(), date, zero.add(provider.getMu()));
+ final UnnormalizedSphericalHarmonics ush = provider.onDate(initialOsculating.getDate().toAbsoluteDate());
+
+ // 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 = FieldNumericalPropagator.tolerances(zero.newInstance(0.1), initialOsculating, OrbitType.CIRCULAR);
+ AdaptiveStepsizeFieldIntegrator integrator = new DormandPrince853FieldIntegrator<>(field, 0.001, 1000, tol[0], tol[1]);
+ integrator.setInitialStepSize(60);
+ FieldNumericalPropagator num = new FieldNumericalPropagator<>(field, integrator);
+ Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
+ num.addForceModel(new HolmesFeatherstoneAttractionModel(itrf, GravityFieldFactory.getNormalizedProvider(provider)));
+ num.setInitialState(new FieldSpacecraftState<>(initialOsculating));
+ num.setOrbitType(OrbitType.CIRCULAR);
+ final StorelessUnivariateStatistic oscMin = new Min();
+ final StorelessUnivariateStatistic oscMax = new Max();
+ final StorelessUnivariateStatistic meanMin = new Min();
+ final StorelessUnivariateStatistic meanMax = new Max();
+ num.getMultiplexer().add(zero.newInstance(60), state -> {
+ final FieldOrbit osc = state.getOrbit();
+ oscMin.increment(osc.getA().getReal());
+ oscMax.increment(osc.getA().getReal());
+ // compute mean orbit at current date (this is what we test)
+ final FieldOrbit mean = FieldEcksteinHechlerPropagator.computeMeanOrbit(state.getOrbit(), provider, ush);
+ meanMin.increment(mean.getA().getReal());
+ meanMax.increment(mean.getA().getReal());
+ });
+ num.propagate(initialOsculating.getDate().shiftedBy(Constants.JULIAN_DAY));
+
+ Assert.assertEquals(3190.029, oscMax.getResult() - oscMin.getResult(), 1.0e-3);
+ Assert.assertEquals( 49.638, meanMax.getResult() - meanMin.getResult(), 1.0e-3);
+
+ }
+
private > T tangLEmLv(T Lv, T ex, T ey) {
// tan ((LE - Lv) /2)) =
return ey.multiply(Lv.cos()).subtract(ex.multiply(Lv.sin())).divide(