From 87add4f59302b11c5d9fc4e72d21ebaba316d79d Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Thu, 2 Jun 2022 10:49:45 +0200 Subject: [PATCH 1/2] =?UTF-8?q?Use=20consistent=20values=20for=20=C2=B5=20?= =?UTF-8?q?in=20tests.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../EcksteinHechlerPropagatorTest.java | 16 ++++++++-------- .../FieldEcksteinHechlerPropagatorTest.java | 16 ++++++++-------- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/src/test/java/org/orekit/propagation/analytical/EcksteinHechlerPropagatorTest.java b/src/test/java/org/orekit/propagation/analytical/EcksteinHechlerPropagatorTest.java index a635178ed..424df7f63 100644 --- a/src/test/java/org/orekit/propagation/analytical/EcksteinHechlerPropagatorTest.java +++ b/src/test/java/org/orekit/propagation/analytical/EcksteinHechlerPropagatorTest.java @@ -504,7 +504,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 +514,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 +614,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 +644,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 +658,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 +679,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 +778,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 +814,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<>(); diff --git a/src/test/java/org/orekit/propagation/analytical/FieldEcksteinHechlerPropagatorTest.java b/src/test/java/org/orekit/propagation/analytical/FieldEcksteinHechlerPropagatorTest.java index e734d3705..2d9d23808 100644 --- a/src/test/java/org/orekit/propagation/analytical/FieldEcksteinHechlerPropagatorTest.java +++ b/src/test/java/org/orekit/propagation/analytical/FieldEcksteinHechlerPropagatorTest.java @@ -598,7 +598,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 +619,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 +754,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 +798,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 +821,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 +851,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 +890,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 +931,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<>(); -- GitLab From b7c6b72f256f9634f08c526efc84861db233d831 Mon Sep 17 00:00:00 2001 From: Luc Maisonobe Date: Thu, 2 Jun 2022 11:50:11 +0200 Subject: [PATCH 2/2] Added a way to compute mean parameters in Eckstein-Hechler model. Fixes #933 --- .../analytical/EcksteinHechlerPropagator.java | 70 ++++++++++++++- .../FieldEcksteinHechlerPropagator.java | 85 +++++++++++++++++-- .../EcksteinHechlerPropagatorTest.java | 43 ++++++++++ .../FieldEcksteinHechlerPropagatorTest.java | 57 ++++++++++++- 4 files changed, 246 insertions(+), 9 deletions(-) diff --git a/src/main/java/org/orekit/propagation/analytical/EcksteinHechlerPropagator.java b/src/main/java/org/orekit/propagation/analytical/EcksteinHechlerPropagator.java index 4ce46c82a..b51696425 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 836d26afe..5cb1e9b0d 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 424df7f63..c27eafec2 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; @@ -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 2d9d23808..673878dba 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; @@ -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( -- GitLab