Skip to content
GitLab
Commit 88e2e9fd authored by Bryan Cazabonne's avatar Bryan Cazabonne
Browse files

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

parent de540de7
Pipeline #2087 failed with stages
in 22 minutes and 46 seconds
Hide whitespace changes
Inline Side-by-side
src/changes/changes.xml
View file @ 88e2e9fd
......@@ -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 @ 88e2e9fd
......@@ -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>
......
src/main/java/org/orekit/propagation/analytical/FieldBrouwerLyddanePropagator.java
View file @ 88e2e9fd
......@@ -129,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);
}
......@@ -152,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),
......@@ -191,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);
}
......@@ -261,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);
}
......@@ -295,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);
}
......@@ -370,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);
}
......@@ -414,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),
......@@ -480,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
......
src/test/java/org/orekit/propagation/analytical/BrouwerLyddanePropagatorTest.java
View file @ 88e2e9fd
......@@ -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;
......@@ -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/FieldBrouwerLyddanePropagatorTest.java
View file @ 88e2e9fd
package org.orekit.propagation.analytical;
import java.io.IOException;
import org.hipparchus.CalculusFieldElement;
import org.hipparchus.Field;
import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
import org.hipparchus.ode.nonstiff.AdaptiveStepsizeFieldIntegrator;
import org.hipparchus.ode.nonstiff.AdaptiveStepsizeIntegrator;
import org.hipparchus.ode.nonstiff.DormandPrince853FieldIntegrator;
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.Decimal64Field;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathUtils;
......@@ -27,6 +34,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;
......@@ -41,6 +49,7 @@ import org.orekit.propagation.FieldSpacecraftState;
import org.orekit.propagation.PropagationType;
import org.orekit.propagation.Propagator;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.numerical.FieldNumericalPropagator;
import org.orekit.propagation.numerical.NumericalPropagator;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.FieldAbsoluteDate;
......@@ -845,6 +854,53 @@ public class FieldBrouwerLyddanePropagatorTest {
Assert.assertEquals(0.0, finalOrbitFieldReal.getPVCoordinates().getVelocity().distance(finalOrbit.getPVCoordinates().getVelocity()), Double.MIN_VALUE);
}
@Test
public void testMeanOrbit() throws IOException {
doTestMeanOrbit(Decimal64Field.getInstance());
}
private <T extends CalculusFieldElement<T>> void doTestMeanOrbit(Field<T> field) {
T zero = field.getZero();
FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
final UnnormalizedSphericalHarmonicsProvider ushp = GravityFieldFactory.getUnnormalizedProvider(provider);
final FieldKeplerianOrbit<T> 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 = ushp.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.KEPLERIAN);
AdaptiveStepsizeFieldIntegrator<T> integrator = new DormandPrince853FieldIntegrator<>(field, 0.001, 1000, tol[0], tol[1]);
integrator.setInitialStepSize(60);
FieldNumericalPropagator<T> num = new FieldNumericalPropagator<>(field, integrator);
Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
num.addForceModel(new HolmesFeatherstoneAttractionModel(itrf, provider));
num.setInitialState(new FieldSpacecraftState<>(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(zero.newInstance(60), state -> {
final FieldOrbit<T> 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<T> mean = FieldBrouwerLyddanePropagator.computeMeanOrbit(state.getOrbit(), ushp, ush, BrouwerLyddanePropagator.M2);
meanMin.increment(mean.getA().getReal());
meanMax.increment(mean.getA().getReal());
});
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");
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment