+ * 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