diff --git a/src/changes/changes.xml b/src/changes/changes.xml index 2d484900ebd80eb31ea5fb2fb31d93e477c0a5bb..3a3a102fdef1b62b0405d22f5b367bf9d028b09e 100644 --- a/src/changes/changes.xml +++ b/src/changes/changes.xml @@ -21,6 +21,9 @@ + + Added atmospheric drag effect for Brouwer-Lyddane model. + Allowed Brouwer-Lyddane model to work for the critical inclination. diff --git a/src/main/java/org/orekit/propagation/analytical/BrouwerLyddanePropagator.java b/src/main/java/org/orekit/propagation/analytical/BrouwerLyddanePropagator.java index d54ebfbb3f8584721dff5c05b8775531332e0253..82d901be4fdde2d51f224ba716dc6cec5e908961 100644 --- a/src/main/java/org/orekit/propagation/analytical/BrouwerLyddanePropagator.java +++ b/src/main/java/org/orekit/propagation/analytical/BrouwerLyddanePropagator.java @@ -16,7 +16,8 @@ */ package org.orekit.propagation.analytical; -import java.io.Serializable; +import java.util.Collections; +import java.util.List; import org.hipparchus.analysis.differentiation.UnivariateDerivative2; import org.hipparchus.util.CombinatoricsUtils; @@ -36,11 +37,11 @@ import org.orekit.orbits.OrbitType; import org.orekit.orbits.PositionAngle; import org.orekit.propagation.PropagationType; import org.orekit.propagation.SpacecraftState; +import org.orekit.propagation.analytical.tle.TLE; import org.orekit.time.AbsoluteDate; +import org.orekit.utils.ParameterDriver; import org.orekit.utils.TimeSpanMap; - - /** * This class propagates a {@link org.orekit.propagation.SpacecraftState} * using the analytical Brouwer-Lyddane model (from J2 to J5 zonal harmonics). @@ -49,7 +50,24 @@ import org.orekit.utils.TimeSpanMap; * suited for elliptical orbits, there is no problem having a rather small eccentricity or inclination * (Lyddane helped to solve this issue with the Brouwer model). Singularity for the critical * inclination i = 63.4° is avoided using the method developed in Warren Phipps' 1992 thesis. - *

+ *

+ * By default, Brouwer-Lyddane model considers only the perturbations due to zonal harmonics. + * However, for low Earth orbits, the magnitude of the perturbative acceleration due to + * atmospheric drag can be significant. Warren Phipps' 1992 thesis considered the atmospheric + * drag by time derivatives of the mean mean anomaly using the catch-all coefficient + * {@link #M2Driver}. + * + * Usually, M2 is adjusted during an orbit determination process and it represents the + * combination of all unmodeled secular along-track effects (i.e. not just the atmospheric drag). + * The behavior of M2 is closed to the {@link TLE#getBStar()} parameter for the TLE. + * + * If the value of M2 is equal to {@link #M2 0.0}, the along-track secular effects are not + * considered in the dynamical model. Typical values for M2 are not known. It depends on the + * orbit type. However, the value of M2 must be very small (e.g. between 1.0e-14 and 1.0e-15). + * + * The along-track effects, represented by the secular rates of the mean semi-major axis + * and eccentricity, are computed following Eq. 2.38, 2.41, and 2.45 of Warren Phipps' thesis. + * * @see "Brouwer, Dirk. Solution of the problem of artificial satellite theory without drag. * YALE UNIV NEW HAVEN CT NEW HAVEN United States, 1959." * @@ -65,6 +83,20 @@ import org.orekit.utils.TimeSpanMap; */ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { + /** Parameter name for M2 coefficient. */ + public static final String M2_NAME = "M2"; + + /** Default value for M2 coefficient. */ + public static final double M2 = 0.0; + + /** Parameters scaling factor. + *

+ * We use a power of 2 to avoid numeric noise introduction + * in the multiplications/divisions sequences. + *

+ */ + private static final double SCALE = FastMath.scalb(1.0, -20); + /** Beta constant used by T2 function. */ private static final double BETA = FastMath.scalb(100, -11); @@ -83,6 +115,9 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { /** Un-normalized zonal coefficients. */ private double[] ck0; + /** Empirical coefficient used in the drag modeling. */ + private final ParameterDriver M2Driver; + /** Build a propagator from orbit and potential provider. *

Mass and attitude provider are set to unspecified non-null arbitrary values.

* @@ -90,15 +125,16 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { * * @param initialOrbit initial orbit * @param provider for un-normalized zonal coefficients - * @see #BrouwerLyddanePropagator(Orbit, AttitudeProvider, - * UnnormalizedSphericalHarmonicsProvider) - * @see #BrouwerLyddanePropagator(Orbit, UnnormalizedSphericalHarmonicsProvider, - * PropagationType) + * @param M2 value of empirical drag coefficient. + * If equal to {@link #M2} drag is not computed + * @see #BrouwerLyddanePropagator(Orbit, AttitudeProvider, UnnormalizedSphericalHarmonicsProvider, double) + * @see #BrouwerLyddanePropagator(Orbit, UnnormalizedSphericalHarmonicsProvider, PropagationType, double) */ public BrouwerLyddanePropagator(final Orbit initialOrbit, - final UnnormalizedSphericalHarmonicsProvider provider) { + final UnnormalizedSphericalHarmonicsProvider provider, + final double M2) { this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()), - DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate())); + DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate()), M2); } /** @@ -109,21 +145,25 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { * @param mass spacecraft mass * @param provider for un-normalized zonal coefficients * @param harmonics {@code provider.onDate(initialOrbit.getDate())} + * @param M2 value of empirical drag coefficient. + * If equal to {@link #M2} drag is not computed * @see #BrouwerLyddanePropagator(Orbit, AttitudeProvider, double, * UnnormalizedSphericalHarmonicsProvider, * UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics, - * PropagationType) + * PropagationType, double) */ public BrouwerLyddanePropagator(final Orbit initialOrbit, final AttitudeProvider attitude, final double mass, final UnnormalizedSphericalHarmonicsProvider provider, - final UnnormalizedSphericalHarmonics harmonics) { + final UnnormalizedSphericalHarmonics harmonics, + final double M2) { this(initialOrbit, attitude, mass, provider.getAe(), provider.getMu(), harmonics.getUnnormalizedCnm(2, 0), harmonics.getUnnormalizedCnm(3, 0), harmonics.getUnnormalizedCnm(4, 0), - harmonics.getUnnormalizedCnm(5, 0)); + harmonics.getUnnormalizedCnm(5, 0), + M2); } /** Build a propagator from orbit and potential. @@ -147,16 +187,18 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { * @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 M2 value of empirical drag coefficient. + * If equal to {@link #M2} drag is not computed * @see org.orekit.utils.Constants * @see #BrouwerLyddanePropagator(Orbit, AttitudeProvider, double, double, double, - * double, double, double, double) + * double, double, double, double, double) */ public BrouwerLyddanePropagator(final Orbit initialOrbit, final double referenceRadius, final double mu, final double c20, final double c30, final double c40, - final double c50) { + final double c50, final double M2) { this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()), - DEFAULT_MASS, referenceRadius, mu, c20, c30, c40, c50); + DEFAULT_MASS, referenceRadius, mu, c20, c30, c40, c50, M2); } /** Build a propagator from orbit, mass and potential provider. @@ -167,13 +209,15 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { * @param initialOrbit initial orbit * @param mass spacecraft mass * @param provider for un-normalized zonal coefficients - * @see #BrouwerLyddanePropagator(Orbit, AttitudeProvider, double, - * UnnormalizedSphericalHarmonicsProvider) + * @param M2 value of empirical drag coefficient. + * If equal to {@link #M2} drag is not computed + * @see #BrouwerLyddanePropagator(Orbit, AttitudeProvider, double, UnnormalizedSphericalHarmonicsProvider, double) */ public BrouwerLyddanePropagator(final Orbit initialOrbit, final double mass, - final UnnormalizedSphericalHarmonicsProvider provider) { + final UnnormalizedSphericalHarmonicsProvider provider, + final double M2) { this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()), - mass, provider, provider.onDate(initialOrbit.getDate())); + mass, provider, provider.onDate(initialOrbit.getDate()), M2); } /** Build a propagator from orbit, mass and potential. @@ -198,15 +242,17 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { * @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 M2 value of empirical drag coefficient. + * If equal to {@link #M2} drag is not computed * @see #BrouwerLyddanePropagator(Orbit, AttitudeProvider, double, double, double, - * double, double, double, double) + * double, double, double, double, double) */ public BrouwerLyddanePropagator(final Orbit initialOrbit, final double mass, final double referenceRadius, final double mu, final double c20, final double c30, final double c40, - final double c50) { + final double c50, final double M2) { this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()), - mass, referenceRadius, mu, c20, c30, c40, c50); + mass, referenceRadius, mu, c20, c30, c40, c50, M2); } /** Build a propagator from orbit, attitude provider and potential provider. @@ -215,11 +261,14 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { * @param initialOrbit initial orbit * @param attitudeProv attitude provider * @param provider for un-normalized zonal coefficients + * @param M2 value of empirical drag coefficient. + * If equal to {@link #M2} drag is not computed */ public BrouwerLyddanePropagator(final Orbit initialOrbit, final AttitudeProvider attitudeProv, - final UnnormalizedSphericalHarmonicsProvider provider) { - this(initialOrbit, attitudeProv, DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate())); + final UnnormalizedSphericalHarmonicsProvider provider, + final double M2) { + this(initialOrbit, attitudeProv, DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate()), M2); } /** Build a propagator from orbit, attitude provider and potential. @@ -244,13 +293,15 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { * @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 M2 value of empirical drag coefficient. + * If equal to {@link #M2} drag is not computed */ public BrouwerLyddanePropagator(final Orbit initialOrbit, final AttitudeProvider attitudeProv, final double referenceRadius, final double mu, final double c20, final double c30, final double c40, - final double c50) { - this(initialOrbit, attitudeProv, DEFAULT_MASS, referenceRadius, mu, c20, c30, c40, c50); + final double c50, final double M2) { + this(initialOrbit, attitudeProv, DEFAULT_MASS, referenceRadius, mu, c20, c30, c40, c50, M2); } /** Build a propagator from orbit, attitude provider, mass and potential provider. @@ -259,14 +310,17 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { * @param attitudeProv attitude provider * @param mass spacecraft mass * @param provider for un-normalized zonal coefficients + * @param M2 value of empirical drag coefficient. + * If equal to {@link #M2} drag is not computed * @see #BrouwerLyddanePropagator(Orbit, AttitudeProvider, double, - * UnnormalizedSphericalHarmonicsProvider, PropagationType) + * UnnormalizedSphericalHarmonicsProvider, PropagationType, double) */ public BrouwerLyddanePropagator(final Orbit initialOrbit, final AttitudeProvider attitudeProv, final double mass, - final UnnormalizedSphericalHarmonicsProvider provider) { - this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate())); + final UnnormalizedSphericalHarmonicsProvider provider, + final double M2) { + this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate()), M2); } /** Build a propagator from orbit, attitude provider, mass and potential. @@ -291,17 +345,19 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { * @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 M2 value of empirical drag coefficient. + * If equal to {@link #M2} drag is not computed * @see #BrouwerLyddanePropagator(Orbit, AttitudeProvider, double, double, double, - * double, double, double, double, PropagationType) + * double, double, double, double, PropagationType, double) */ public BrouwerLyddanePropagator(final Orbit initialOrbit, final AttitudeProvider attitudeProv, final double mass, final double referenceRadius, final double mu, final double c20, final double c30, final double c40, - final double c50) { + final double c50, final double M2) { this(initialOrbit, attitudeProv, mass, referenceRadius, mu, c20, c30, c40, c50, - PropagationType.OSCULATING); + PropagationType.OSCULATING, M2); } @@ -314,12 +370,14 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { * @param initialOrbit initial orbit * @param provider for un-normalized zonal coefficients * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit) + * @param M2 value of empirical drag coefficient. + * If equal to {@link #M2} drag is not computed */ public BrouwerLyddanePropagator(final Orbit initialOrbit, final UnnormalizedSphericalHarmonicsProvider provider, - final PropagationType initialType) { + final PropagationType initialType, final double M2) { this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()), - DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate()), initialType); + DEFAULT_MASS, provider, provider.onDate(initialOrbit.getDate()), initialType, M2); } /** Build a propagator from orbit, attitude provider, mass and potential provider. @@ -330,13 +388,15 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { * @param mass spacecraft mass * @param provider for un-normalized zonal coefficients * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit) + * @param M2 value of empirical drag coefficient. + * If equal to {@link #M2} drag is not computed */ public BrouwerLyddanePropagator(final Orbit initialOrbit, final AttitudeProvider attitudeProv, final double mass, final UnnormalizedSphericalHarmonicsProvider provider, - final PropagationType initialType) { - this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate()), initialType); + final PropagationType initialType, final double M2) { + this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate()), initialType, M2); } /** @@ -349,20 +409,21 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { * @param provider for un-normalized zonal coefficients * @param harmonics {@code provider.onDate(initialOrbit.getDate())} * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit) - + * @param M2 value of empirical drag coefficient. + * If equal to {@link #M2} drag is not computed */ public BrouwerLyddanePropagator(final Orbit initialOrbit, final AttitudeProvider attitude, final double mass, final UnnormalizedSphericalHarmonicsProvider provider, final UnnormalizedSphericalHarmonics harmonics, - final PropagationType initialType) { + final PropagationType initialType, final double M2) { this(initialOrbit, attitude, mass, provider.getAe(), provider.getMu(), harmonics.getUnnormalizedCnm(2, 0), harmonics.getUnnormalizedCnm(3, 0), harmonics.getUnnormalizedCnm(4, 0), harmonics.getUnnormalizedCnm(5, 0), - initialType); + initialType, M2); } /** Build a propagator from orbit, attitude provider, mass and potential. @@ -389,6 +450,8 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { * @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 initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit) + * @param M2 value of empirical drag coefficient. + * If equal to {@link #M2} drag is not computed */ public BrouwerLyddanePropagator(final Orbit initialOrbit, final AttitudeProvider attitudeProv, @@ -396,7 +459,7 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { final double referenceRadius, final double mu, final double c20, final double c30, final double c40, final double c50, - final PropagationType initialType) { + final PropagationType initialType, final double M2) { super(attitudeProv); @@ -405,14 +468,19 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { this.mu = mu; this.ck0 = new double[] {0.0, 0.0, c20, c30, c40, c50}; + // initialize M2 driver + this.M2Driver = new ParameterDriver(M2_NAME, M2, SCALE, + Double.NEGATIVE_INFINITY, + Double.POSITIVE_INFINITY); + // compute mean parameters if needed - // transform into circular adapted parameters used by the Brouwer-Lyddane model resetInitialState(new SpacecraftState(initialOrbit, attitudeProv.getAttitude(initialOrbit, initialOrbit.getDate(), initialOrbit.getFrame()), mass), initialType); + } /** {@inheritDoc} @@ -527,13 +595,27 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { return new KeplerianOrbit(propOrb_parameters[0].getValue(), propOrb_parameters[1].getValue(), propOrb_parameters[2].getValue(), propOrb_parameters[3].getValue(), propOrb_parameters[4].getValue(), propOrb_parameters[5].getValue(), - PositionAngle.MEAN, current.mean.getFrame(), date, mu); } + PositionAngle.MEAN, current.mean.getFrame(), date, mu); + } - /** Local class for Brouwer-Lyddane model. */ - private static class BLModel implements Serializable { + /** + * Get the value of the M2 drag parameter. + * @return the value of the M2 drag parameter + */ + public double getM2() { + return M2Driver.getValue(); + } - /** Serializable UID. */ - private static final long serialVersionUID = 20211123L; + /** + * Get the parameters driver for propagation model. + * @return drivers for propagation model + */ + public List getParametersDrivers() { + return Collections.singletonList(M2Driver); + } + + /** Local class for Brouwer-Lyddane model. */ + private class BLModel { /** Mean orbit. */ private final KeplerianOrbit mean; @@ -621,6 +703,10 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { private final double edls2gf; private final double edls2g3f; + // Drag terms + private final double aRate; + private final double eRate; + // CHECKSTYLE: resume JavadocVariable check /** Create a model for specified mean orbit. @@ -802,6 +888,11 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { edls2gf = 3.0 * qedl * (1.0 - cosI2); edls2g3f = 1.0 / 3.0 * qedl; + // secular rates of the mean semi-major axis and eccentricity + // Eq. 2.41 and Eq. 2.45 of Phipps' 1992 thesis + aRate = -4.0 * app / (3.0 * xnotDot); + eRate = -4.0 * epp * n * n / (3.0 * xnotDot); + } /** @@ -941,6 +1032,9 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { */ public UnivariateDerivative2[] propagateParameters(final AbsoluteDate date) { + // Empirical drag coefficient M2 + final double m2 = M2Driver.getValue(); + // Keplerian evolution final UnivariateDerivative2 dt = new UnivariateDerivative2(date.durationFrom(mean.getDate()), 1.0, 0.0); final UnivariateDerivative2 xnot = dt.multiply(xnotDot); @@ -948,11 +1042,13 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { //____________________________________ // secular effects - // mean mean anomaly - final UnivariateDerivative2 lpp = new UnivariateDerivative2(MathUtils.normalizeAngle(mean.getMeanAnomaly() + lt * xnot.getValue(), + // mean mean anomaly (with drag Eq. 2.38 of Phipps' 1992 thesis) + final UnivariateDerivative2 dtM2 = dt.multiply(m2); + final UnivariateDerivative2 dt2M2 = dt.multiply(dtM2); + final UnivariateDerivative2 lpp = new UnivariateDerivative2(MathUtils.normalizeAngle(mean.getMeanAnomaly() + lt * xnot.getValue() + dt2M2.getValue(), FastMath.PI), - lt * xnotDot, - 0.0); + lt * xnotDot + 2.0 * dtM2.getValue(), + 2.0 * m2); // mean argument of perigee final UnivariateDerivative2 gpp = new UnivariateDerivative2(MathUtils.normalizeAngle(mean.getPerigeeArgument() + gt * xnot.getValue(), FastMath.PI), @@ -964,6 +1060,15 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { ht * xnotDot, 0.0); + // ________________________________________________ + // secular rates of the mean semi-major axis and eccentricity + + // semi-major axis + final UnivariateDerivative2 appDrag = dt.multiply(aRate * m2); + + // eccentricity + final UnivariateDerivative2 eppDrag = dt.multiply(eRate * m2); + //____________________________________ // Long periodical terms final UnivariateDerivative2 cg1 = gpp.cos(); @@ -1019,13 +1124,13 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { final UnivariateDerivative2 eE3 = eE.multiply(eE).multiply(eE); final UnivariateDerivative2 sigma = eE.multiply(n * n).multiply(eE).add(eE); - // Semi-major axis - final UnivariateDerivative2 a = eE3.multiply(aCbis).add(mean.getA()). + // Semi-major axis (with drag Eq. 2.41 of Phipps' 1992 thesis) + final UnivariateDerivative2 a = eE3.multiply(aCbis).add(appDrag.add(mean.getA())). add(aC). add(eE3.multiply(c2g2f).multiply(ac2g2f)); - // Eccentricity - final UnivariateDerivative2 e = d1e.add(mean.getE()). + // Eccentricity (with drag Eq. 2.45 of Phipps' 1992 thesis) + final UnivariateDerivative2 e = d1e.add(eppDrag.add(mean.getE())). add(eC). add(cf1.multiply(ecf)). add(cf2.multiply(e2cf)). diff --git a/src/main/java/org/orekit/propagation/analytical/FieldBrouwerLyddanePropagator.java b/src/main/java/org/orekit/propagation/analytical/FieldBrouwerLyddanePropagator.java index 87bca775f61d92070a5617b067a84cd3eafc3142..c633b8749426d6f82ed4ae18167e17d772225f7b 100644 --- a/src/main/java/org/orekit/propagation/analytical/FieldBrouwerLyddanePropagator.java +++ b/src/main/java/org/orekit/propagation/analytical/FieldBrouwerLyddanePropagator.java @@ -40,6 +40,7 @@ import org.orekit.orbits.OrbitType; import org.orekit.orbits.PositionAngle; import org.orekit.propagation.FieldSpacecraftState; import org.orekit.propagation.PropagationType; +import org.orekit.propagation.analytical.tle.FieldTLE; import org.orekit.time.FieldAbsoluteDate; import org.orekit.utils.FieldTimeSpanMap; import org.orekit.utils.ParameterDriver; @@ -51,7 +52,24 @@ import org.orekit.utils.ParameterDriver; * suited for elliptical orbits, there is no problem having a rather small eccentricity or inclination * (Lyddane helped to solve this issue with the Brouwer model). Singularity for the critical * inclination i = 63.4° is avoided using the method developed in Warren Phipps' 1992 thesis. - *

+ *

+ * By default, Brouwer-Lyddane model considers only the perturbations due to zonal harmonics. + * However, for low Earth orbits, the magnitude of the perturbative acceleration due to + * atmospheric drag can be significant. Warren Phipps' 1992 thesis considered the atmospheric + * drag by time derivatives of the mean mean anomaly using the catch-all coefficient + * {@link #M2Driver}. + * + * Usually, M2 is adjusted during an orbit determination process and it represents the + * combination of all unmodeled secular along-track effects (i.e. not just the atmospheric drag). + * The behavior of M2 is closed to the {@link FieldTLE#getBStar()} parameter for the TLE. + * + * If the value of M2 is equal to {@link BrouwerLyddanePropagator#M2 0.0}, the along-track secular + * effects are not considered in the dynamical model. Typical values for M2 are not known. + * It depends on the orbit type. However, the value of M2 must be very small (e.g. between 1.0e-14 and 1.0e-15). + * + * The along-track effects, represented by the secular rates of the mean semi-major axis + * and eccentricity, are computed following Eq. 2.38, 2.41, and 2.45 of Warren Phipps' thesis. + * * @see "Brouwer, Dirk. Solution of the problem of artificial satellite theory without drag. * YALE UNIV NEW HAVEN CT NEW HAVEN United States, 1959." * @@ -67,6 +85,14 @@ import org.orekit.utils.ParameterDriver; */ public class FieldBrouwerLyddanePropagator> extends FieldAbstractAnalyticalPropagator { + /** Parameters scaling factor. + *

+ * We use a power of 2 to avoid numeric noise introduction + * in the multiplications/divisions sequences. + *

+ */ + private static final double SCALE = FastMath.scalb(1.0, -20); + /** Beta constant used by T2 function. */ private static final double BETA = FastMath.scalb(100, -11); @@ -85,6 +111,9 @@ public class FieldBrouwerLyddanePropagator> ex /** Un-normalized zonal coefficients. */ private double[] ck0; + /** Empirical coefficient used in the drag modeling. */ + private final ParameterDriver M2Driver; + /** Build a propagator from orbit and potential provider. *

Mass and attitude provider are set to unspecified non-null arbitrary values.

* @@ -92,13 +121,16 @@ public class FieldBrouwerLyddanePropagator> ex * * @param initialOrbit initial orbit * @param provider for un-normalized zonal coefficients - * @see #FieldBrouwerLyddanePropagator(FieldOrbit, UnnormalizedSphericalHarmonicsProvider, PropagationType) + * @param M2 value of empirical drag coefficient. + * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed + * @see #FieldBrouwerLyddanePropagator(FieldOrbit, UnnormalizedSphericalHarmonicsProvider, PropagationType, double) */ public FieldBrouwerLyddanePropagator(final FieldOrbit initialOrbit, - final UnnormalizedSphericalHarmonicsProvider provider) { + final UnnormalizedSphericalHarmonicsProvider provider, + final double M2) { this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()), initialOrbit.getA().getField().getZero().add(DEFAULT_MASS), provider, - provider.onDate(initialOrbit.getDate().toAbsoluteDate())); + provider.onDate(initialOrbit.getDate().toAbsoluteDate()), M2); } /** @@ -109,19 +141,23 @@ public class FieldBrouwerLyddanePropagator> ex * @param mass spacecraft mass * @param provider for un-normalized zonal coefficients * @param harmonics {@code provider.onDate(initialOrbit.getDate())} + * @param M2 value of empirical drag coefficient. + * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement, - * UnnormalizedSphericalHarmonicsProvider, UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics, PropagationType) + * UnnormalizedSphericalHarmonicsProvider, UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics, PropagationType, double) */ public FieldBrouwerLyddanePropagator(final FieldOrbit initialOrbit, - final AttitudeProvider attitude, + final AttitudeProvider attitude, final T mass, final UnnormalizedSphericalHarmonicsProvider provider, - final UnnormalizedSphericalHarmonics harmonics) { + final UnnormalizedSphericalHarmonics harmonics, + final double M2) { this(initialOrbit, attitude, mass, provider.getAe(), initialOrbit.getA().getField().getZero().add(provider.getMu()), harmonics.getUnnormalizedCnm(2, 0), harmonics.getUnnormalizedCnm(3, 0), harmonics.getUnnormalizedCnm(4, 0), - harmonics.getUnnormalizedCnm(5, 0)); + harmonics.getUnnormalizedCnm(5, 0), + M2); } /** Build a propagator from orbit and potential. @@ -145,16 +181,18 @@ public class FieldBrouwerLyddanePropagator> ex * @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 M2 value of empirical drag coefficient. + * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed * @see org.orekit.utils.Constants - * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, double, CalculusFieldElement, double, double, double, double) + * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, double, CalculusFieldElement, double, double, double, double, double) */ public FieldBrouwerLyddanePropagator(final FieldOrbit initialOrbit, final double referenceRadius, final T mu, final double c20, final double c30, final double c40, - final double c50) { + final double c50, final double M2) { this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()), initialOrbit.getDate().getField().getZero().add(DEFAULT_MASS), - referenceRadius, mu, c20, c30, c40, c50); + referenceRadius, mu, c20, c30, c40, c50, M2); } /** Build a propagator from orbit, mass and potential provider. @@ -165,12 +203,15 @@ public class FieldBrouwerLyddanePropagator> ex * @param initialOrbit initial orbit * @param mass spacecraft mass * @param provider for un-normalized zonal coefficients - * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement, UnnormalizedSphericalHarmonicsProvider) + * @param M2 value of empirical drag coefficient. + * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed + * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement, UnnormalizedSphericalHarmonicsProvider, double) */ public FieldBrouwerLyddanePropagator(final FieldOrbit initialOrbit, final T mass, - final UnnormalizedSphericalHarmonicsProvider provider) { + final UnnormalizedSphericalHarmonicsProvider provider, + final double M2) { this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()), - mass, provider, provider.onDate(initialOrbit.getDate().toAbsoluteDate())); + mass, provider, provider.onDate(initialOrbit.getDate().toAbsoluteDate()), M2); } /** Build a propagator from orbit, mass and potential. @@ -195,14 +236,16 @@ public class FieldBrouwerLyddanePropagator> ex * @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) - * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement, double, CalculusFieldElement, double, double, double, double) + * @param M2 value of empirical drag coefficient. + * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed + * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement, double, CalculusFieldElement, double, double, double, double, double) */ public FieldBrouwerLyddanePropagator(final FieldOrbit initialOrbit, final T mass, final double referenceRadius, final T mu, final double c20, final double c30, final double c40, - final double c50) { + final double c50, final double M2) { this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()), - mass, referenceRadius, mu, c20, c30, c40, c50); + mass, referenceRadius, mu, c20, c30, c40, c50, M2); } /** Build a propagator from orbit, attitude provider and potential provider. @@ -211,12 +254,15 @@ public class FieldBrouwerLyddanePropagator> ex * @param initialOrbit initial orbit * @param attitudeProv attitude provider * @param provider for un-normalized zonal coefficients + * @param M2 value of empirical drag coefficient. + * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed */ public FieldBrouwerLyddanePropagator(final FieldOrbit initialOrbit, final AttitudeProvider attitudeProv, - final UnnormalizedSphericalHarmonicsProvider provider) { + final UnnormalizedSphericalHarmonicsProvider provider, + final double M2) { this(initialOrbit, attitudeProv, initialOrbit.getA().getField().getZero().add(DEFAULT_MASS), provider, - provider.onDate(initialOrbit.getDate().toAbsoluteDate())); + provider.onDate(initialOrbit.getDate().toAbsoluteDate()), M2); } /** Build a propagator from orbit, attitude provider and potential. @@ -241,14 +287,16 @@ public class FieldBrouwerLyddanePropagator> ex * @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 M2 value of empirical drag coefficient. + * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed */ public FieldBrouwerLyddanePropagator(final FieldOrbit initialOrbit, final AttitudeProvider attitudeProv, final double referenceRadius, final T mu, final double c20, final double c30, final double c40, - final double c50) { + final double c50, final double M2) { this(initialOrbit, attitudeProv, initialOrbit.getDate().getField().getZero().add(DEFAULT_MASS), - referenceRadius, mu, c20, c30, c40, c50); + referenceRadius, mu, c20, c30, c40, c50, M2); } /** Build a propagator from orbit, attitude provider, mass and potential provider. @@ -257,13 +305,16 @@ public class FieldBrouwerLyddanePropagator> ex * @param attitudeProv attitude provider * @param mass spacecraft mass * @param provider for un-normalized zonal coefficients - * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement, UnnormalizedSphericalHarmonicsProvider, PropagationType) + * @param M2 value of empirical drag coefficient. + * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed + * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement, UnnormalizedSphericalHarmonicsProvider, PropagationType, double) */ public FieldBrouwerLyddanePropagator(final FieldOrbit initialOrbit, final AttitudeProvider attitudeProv, final T mass, - final UnnormalizedSphericalHarmonicsProvider provider) { - this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate().toAbsoluteDate())); + final UnnormalizedSphericalHarmonicsProvider provider, + final double M2) { + this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate().toAbsoluteDate()), M2); } /** Build a propagator from orbit, attitude provider, mass and potential. @@ -288,15 +339,17 @@ public class FieldBrouwerLyddanePropagator> ex * @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) - * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement, double, CalculusFieldElement, double, double, double, double, PropagationType) + * @param M2 value of empirical drag coefficient. + * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed + * @see #FieldBrouwerLyddanePropagator(FieldOrbit, AttitudeProvider, CalculusFieldElement, double, CalculusFieldElement, double, double, double, double, PropagationType, double) */ public FieldBrouwerLyddanePropagator(final FieldOrbit initialOrbit, final AttitudeProvider attitudeProv, final T mass, final double referenceRadius, final T mu, final double c20, final double c30, final double c40, - final double c50) { - this(initialOrbit, attitudeProv, mass, referenceRadius, mu, c20, c30, c40, c50, PropagationType.OSCULATING); + final double c50, final double M2) { + this(initialOrbit, attitudeProv, mass, referenceRadius, mu, c20, c30, c40, c50, PropagationType.OSCULATING, M2); } @@ -309,13 +362,16 @@ public class FieldBrouwerLyddanePropagator> ex * @param initialOrbit initial orbit * @param provider for un-normalized zonal coefficients * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit) + * @param M2 value of empirical drag coefficient. + * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed */ public FieldBrouwerLyddanePropagator(final FieldOrbit initialOrbit, final UnnormalizedSphericalHarmonicsProvider provider, - final PropagationType initialType) { + final PropagationType initialType, + final double M2) { this(initialOrbit, InertialProvider.of(initialOrbit.getFrame()), initialOrbit.getA().getField().getZero().add(DEFAULT_MASS), provider, - provider.onDate(initialOrbit.getDate().toAbsoluteDate()), initialType); + provider.onDate(initialOrbit.getDate().toAbsoluteDate()), initialType, M2); } /** Build a propagator from orbit, attitude provider, mass and potential provider. @@ -326,13 +382,16 @@ public class FieldBrouwerLyddanePropagator> ex * @param mass spacecraft mass * @param provider for un-normalized zonal coefficients * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit) + * @param M2 value of empirical drag coefficient. + * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed */ public FieldBrouwerLyddanePropagator(final FieldOrbit initialOrbit, final AttitudeProvider attitudeProv, final T mass, final UnnormalizedSphericalHarmonicsProvider provider, - final PropagationType initialType) { - this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate().toAbsoluteDate()), initialType); + final PropagationType initialType, + final double M2) { + this(initialOrbit, attitudeProv, mass, provider, provider.onDate(initialOrbit.getDate().toAbsoluteDate()), initialType, M2); } /** @@ -345,20 +404,22 @@ public class FieldBrouwerLyddanePropagator> ex * @param provider for un-normalized zonal coefficients * @param harmonics {@code provider.onDate(initialOrbit.getDate())} * @param initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit) - + * @param M2 value of empirical drag coefficient. + * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed */ public FieldBrouwerLyddanePropagator(final FieldOrbit initialOrbit, final AttitudeProvider attitude, final T mass, final UnnormalizedSphericalHarmonicsProvider provider, final UnnormalizedSphericalHarmonics harmonics, - final PropagationType initialType) { + final PropagationType initialType, + final double M2) { this(initialOrbit, attitude, mass, provider.getAe(), initialOrbit.getA().getField().getZero().add(provider.getMu()), harmonics.getUnnormalizedCnm(2, 0), harmonics.getUnnormalizedCnm(3, 0), harmonics.getUnnormalizedCnm(4, 0), harmonics.getUnnormalizedCnm(5, 0), - initialType); + initialType, M2); } /** Build a propagator from orbit, attitude provider, mass and potential. @@ -385,6 +446,8 @@ public class FieldBrouwerLyddanePropagator> ex * @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 initialType initial orbit type (mean Brouwer-Lyddane orbit or osculating orbit) + * @param M2 value of empirical drag coefficient. + * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed */ public FieldBrouwerLyddanePropagator(final FieldOrbit initialOrbit, final AttitudeProvider attitudeProv, @@ -392,27 +455,29 @@ public class FieldBrouwerLyddanePropagator> ex final double referenceRadius, final T mu, final double c20, final double c30, final double c40, final double c50, - final PropagationType initialType) { + final PropagationType initialType, + final double M2) { super(mass.getField(), attitudeProv); - try { // store model coefficients - this.referenceRadius = referenceRadius; - this.mu = mu; - this.ck0 = new double[] {0.0, 0.0, c20, c30, c40, c50}; - - // compute mean parameters if needed - // transform into circular adapted parameters used by the Brouwer-Lyddane model - resetInitialState(new FieldSpacecraftState<>(initialOrbit, - attitudeProv.getAttitude(initialOrbit, - initialOrbit.getDate(), - initialOrbit.getFrame()), - mass), - initialType); - } catch (OrekitException oe) { - throw new OrekitException(oe); - } + this.referenceRadius = referenceRadius; + this.mu = mu; + this.ck0 = new double[] {0.0, 0.0, c20, c30, c40, c50}; + + // initialize M2 driver + this.M2Driver = new ParameterDriver(BrouwerLyddanePropagator.M2_NAME, M2, SCALE, + Double.NEGATIVE_INFINITY, + Double.POSITIVE_INFINITY); + + // compute mean parameters if needed + resetInitialState(new FieldSpacecraftState<>(initialOrbit, + attitudeProv.getAttitude(initialOrbit, + initialOrbit.getDate(), + initialOrbit.getFrame()), + mass), + initialType); + } @@ -483,7 +548,7 @@ public class FieldBrouwerLyddanePropagator> ex while (i++ < 200) { // recompute the osculating parameters from the current mean parameters - final FieldUnivariateDerivative2[] parameters = current.propagateParameters(current.mean.getDate()); + final FieldUnivariateDerivative2[] parameters = current.propagateParameters(current.mean.getDate(), getParameters(mass.getField())); // adapted parameters residuals final T deltaA = osculating.getA() .subtract(parameters[0].getValue()); @@ -524,11 +589,20 @@ public class FieldBrouwerLyddanePropagator> ex // compute Cartesian parameters, taking derivatives into account // to make sure velocity and acceleration are consistent final FieldBLModel current = models.get(date); - final FieldUnivariateDerivative2[] propOrb_parameters = current.propagateParameters(date); + final FieldUnivariateDerivative2[] propOrb_parameters = current.propagateParameters(date, parameters); return new FieldKeplerianOrbit(propOrb_parameters[0].getValue(), propOrb_parameters[1].getValue(), propOrb_parameters[2].getValue(), propOrb_parameters[3].getValue(), propOrb_parameters[4].getValue(), propOrb_parameters[5].getValue(), - PositionAngle.MEAN, current.mean.getFrame(), date, mu); } + PositionAngle.MEAN, current.mean.getFrame(), date, mu); + } + + /** + * Get the value of the M2 drag parameter. + * @return the value of the M2 drag parameter + */ + public double getM2() { + return M2Driver.getValue(); + } /** Local class for Brouwer-Lyddane model. */ private static class FieldBLModel> { @@ -619,6 +693,10 @@ public class FieldBrouwerLyddanePropagator> ex private final T edls2gf; private final T edls2g3f; + // Drag terms + private final T aRate; + private final T eRate; + // CHECKSTYLE: resume JavadocVariable check /** Create a model for specified mean orbit. @@ -811,6 +889,11 @@ public class FieldBrouwerLyddanePropagator> ex edls2gf = qedl.multiply(3.0).multiply(cosI2.negate().add(1.0)); edls2g3f = qedl.multiply(1.0 / 3.0); + // secular rates of the mean semi-major axis and eccentricity + // Eq. 2.41 and Eq. 2.45 of Phipps' 1992 thesis + aRate = app.multiply(-4.0).divide(xnotDot.multiply(3.0)); + eRate = epp.multiply(n).multiply(n).multiply(-4.0).divide(xnotDot.multiply(3.0)); + } /** @@ -949,12 +1032,19 @@ public class FieldBrouwerLyddanePropagator> ex /** Extrapolate an orbit up to a specific target date. * @param date target date for the orbit + * @param parameters model parameters * @return propagated parameters */ - public FieldUnivariateDerivative2[] propagateParameters(final FieldAbsoluteDate date) { - final Field field = date.durationFrom(mean.getDate()).getField(); - final T one = field.getOne(); + public FieldUnivariateDerivative2[] propagateParameters(final FieldAbsoluteDate date, final T[] parameters) { + + // Field + final Field field = date.getField(); + final T one = field.getOne(); final T zero = field.getZero(); + + // Empirical drag coefficient M2 + final T m2 = parameters[0]; + // Keplerian evolution final FieldUnivariateDerivative2 dt = new FieldUnivariateDerivative2<>(date.durationFrom(mean.getDate()), one, zero); final FieldUnivariateDerivative2 xnot = dt.multiply(xnotDot); @@ -963,10 +1053,12 @@ public class FieldBrouwerLyddanePropagator> ex // secular effects // mean mean anomaly - final FieldUnivariateDerivative2 lpp = new FieldUnivariateDerivative2(MathUtils.normalizeAngle(mean.getMeanAnomaly().add(lt.multiply(xnot.getValue())), + final FieldUnivariateDerivative2 dtM2 = dt.multiply(m2); + final FieldUnivariateDerivative2 dt2M2 = dt.multiply(dtM2); + final FieldUnivariateDerivative2 lpp = new FieldUnivariateDerivative2(MathUtils.normalizeAngle(mean.getMeanAnomaly().add(lt.multiply(xnot.getValue())).add(dt2M2.getValue()), one.getPi()), - lt.multiply(xnotDot), - zero); + lt.multiply(xnotDot).add(dtM2.multiply(2.0).getValue()), + m2.multiply(2.0)); // mean argument of perigee final FieldUnivariateDerivative2 gpp = new FieldUnivariateDerivative2(MathUtils.normalizeAngle(mean.getPerigeeArgument().add(gt.multiply(xnot.getValue())), one.getPi()), @@ -978,6 +1070,15 @@ public class FieldBrouwerLyddanePropagator> ex ht.multiply(xnotDot), zero); + // ________________________________________________ + // secular rates of the mean semi-major axis and eccentricity + + // semi-major axis + final FieldUnivariateDerivative2 appDrag = dt.multiply(aRate.multiply(m2)); + + // eccentricity + final FieldUnivariateDerivative2 eppDrag = dt.multiply(eRate.multiply(m2)); + //____________________________________ // Long periodical terms final FieldUnivariateDerivative2 cg1 = gpp.cos(); @@ -1034,12 +1135,12 @@ public class FieldBrouwerLyddanePropagator> ex final FieldUnivariateDerivative2 sigma = eE.multiply(n.multiply(n)).multiply(eE).add(eE); // Semi-major axis - final FieldUnivariateDerivative2 a = eE3.multiply(aCbis).add(mean.getA()). + final FieldUnivariateDerivative2 a = eE3.multiply(aCbis).add(appDrag.add(mean.getA())). add(aC). add(eE3.multiply(c2g2f).multiply(ac2g2f)); // Eccentricity - final FieldUnivariateDerivative2 e = d1e.add(mean.getE()). + final FieldUnivariateDerivative2 e = d1e.add(eppDrag.add(mean.getE())). add(eC). add(cf1.multiply(ecf)). add(cf2.multiply(e2cf)). @@ -1117,8 +1218,7 @@ public class FieldBrouwerLyddanePropagator> ex /** {@inheritDoc} */ @Override protected List getParametersDrivers() { - // Brouwer Lyddane propagation model does not have parameter drivers. - return Collections.emptyList(); + return Collections.singletonList(M2Driver); } } diff --git a/src/main/java/org/orekit/propagation/conversion/BrouwerLyddanePropagatorBuilder.java b/src/main/java/org/orekit/propagation/conversion/BrouwerLyddanePropagatorBuilder.java index f65c2056b726696b07cffa3405911032fab59825..13ae760576ced8353ea481c489d143504eae4288 100644 --- a/src/main/java/org/orekit/propagation/conversion/BrouwerLyddanePropagatorBuilder.java +++ b/src/main/java/org/orekit/propagation/conversion/BrouwerLyddanePropagatorBuilder.java @@ -17,6 +17,7 @@ package org.orekit.propagation.conversion; +import org.hipparchus.util.FastMath; import org.orekit.attitudes.AttitudeProvider; import org.orekit.attitudes.InertialProvider; import org.orekit.forces.gravity.potential.GravityFieldFactory; @@ -25,15 +26,48 @@ import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvide import org.orekit.orbits.Orbit; import org.orekit.orbits.OrbitType; import org.orekit.orbits.PositionAngle; -import org.orekit.propagation.Propagator; import org.orekit.propagation.analytical.BrouwerLyddanePropagator; +import org.orekit.propagation.analytical.tle.TLE; +import org.orekit.utils.ParameterDriver; /** Builder for Brouwer-Lyddane propagator. + *

+ * By default, Brouwer-Lyddane model considers only the perturbations due to zonal harmonics. + * However, for low Earth orbits, the magnitude of the perturbative acceleration due to + * atmospheric drag can be significant. Warren Phipps' 1992 thesis considered the atmospheric + * drag by time derivatives of the mean mean anomaly using the catch-all coefficient M2. + * + * Usually, M2 is adjusted during an orbit determination process and it represents the + * combination of all unmodeled secular along-track effects (i.e. not just the atmospheric drag). + * The behavior of M2 is closed to the {@link TLE#getBStar()} parameter for the TLE. + * + * If the value of M2 is equal to {@link BrouwerLyddanePropagator#M2 0.0}, the along-track + * secular effects are not considered in the dynamical model. Typical values for M2 are not known. + * It depends on the orbit type. However, the value of M2 must be very small (e.g. between 1.0e-14 and 1.0e-15). + *

+ * To estimate the M2 parameter, it is necessary to call the {@link #getPropagationParametersDrivers()} method + * as follow: + * 

+ *  for (ParameterDriver driver : builder.getPropagationParametersDrivers().getDrivers()) {
+ *     if (BrouwerLyddanePropagator.M2_NAME.equals(driver.getName())) {
+ *        driver.setSelected(true);
+ *     }
+ *  }
+ *
* @author Melina Vanel + * @author Bryan Cazabonne * @since 11.1 */ public class BrouwerLyddanePropagatorBuilder extends AbstractPropagatorBuilder { + /** Parameters scaling factor. + *

+ * We use a power of 2 to avoid numeric noise introduction + * in the multiplications/divisions sequences. + *

+ */ + private static final double SCALE = FastMath.scalb(1.0, -20); + /** Provider for un-normalized coefficients. */ private final UnnormalizedSphericalHarmonicsProvider provider; @@ -54,43 +88,17 @@ public class BrouwerLyddanePropagatorBuilder extends AbstractPropagatorBuilder { * @param positionAngle position angle type to use * @param positionScale scaling factor used for orbital parameters normalization * (typically set to the expected standard deviation of the position) + * @param M2 value of empirical drag coefficient. + * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed * @see #BrouwerLyddanePropagatorBuilder(Orbit, - * UnnormalizedSphericalHarmonicsProvider, PositionAngle, double, AttitudeProvider) - */ - public BrouwerLyddanePropagatorBuilder(final Orbit templateOrbit, - final UnnormalizedSphericalHarmonicsProvider provider, - final PositionAngle positionAngle, - final double positionScale) { - this(templateOrbit, provider, positionAngle, positionScale, - InertialProvider.of(templateOrbit.getFrame())); - } - - /** Build a new instance. - *

- * The template orbit is used as a model to {@link - * #createInitialOrbit() create initial orbit}. It defines the - * inertial frame, the central attraction coefficient, the orbit type, and is also - * used together with the {@code positionScale} to convert from the {@link - * org.orekit.utils.ParameterDriver#setNormalizedValue(double) normalized} parameters used by the - * callers of this builder to the real orbital parameters. - *

- * @param templateOrbit reference orbit from which real orbits will be built - * (note that the mu from this orbit will be overridden with the mu from the - * {@code provider}) - * @param provider for un-normalized zonal coefficients - * @param positionAngle position angle type to use - * @param positionScale scaling factor used for orbital parameters normalization - * (typically set to the expected standard deviation of the position) - * @param attitudeProvider attitude law to use. + * UnnormalizedSphericalHarmonicsProvider, PositionAngle, double, AttitudeProvider, double) */ public BrouwerLyddanePropagatorBuilder(final Orbit templateOrbit, final UnnormalizedSphericalHarmonicsProvider provider, final PositionAngle positionAngle, final double positionScale, - final AttitudeProvider attitudeProvider) { - super(overrideMu(templateOrbit, provider, positionAngle), positionAngle, - positionScale, true, attitudeProvider); - this.provider = provider; + final double M2) { + this(templateOrbit, provider, positionAngle, positionScale, InertialProvider.of(templateOrbit.getFrame()), M2); } /** Build a new instance. @@ -117,8 +125,10 @@ public class BrouwerLyddanePropagatorBuilder extends AbstractPropagatorBuilder { * @param positionAngle position angle type to use * @param positionScale scaling factor used for orbital parameters normalization * (typically set to the expected standard deviation of the position) + * @param M2 value of empirical drag coefficient. + * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed * @see #BrouwerLyddanePropagatorBuilder(Orbit, - * UnnormalizedSphericalHarmonicsProvider, PositionAngle, double, AttitudeProvider) + * UnnormalizedSphericalHarmonicsProvider, PositionAngle, double, AttitudeProvider, double) */ public BrouwerLyddanePropagatorBuilder(final Orbit templateOrbit, final double referenceRadius, @@ -130,7 +140,8 @@ public class BrouwerLyddanePropagatorBuilder extends AbstractPropagatorBuilder { final double c50, final OrbitType orbitType, final PositionAngle positionAngle, - final double positionScale) { + final double positionScale, + final double M2) { this(templateOrbit, GravityFieldFactory.getUnnormalizedProvider(referenceRadius, mu, tideSystem, new double[][] { @@ -162,7 +173,42 @@ public class BrouwerLyddanePropagatorBuilder extends AbstractPropagatorBuilder { 0 } }), - positionAngle, positionScale); + positionAngle, positionScale, M2); + } + + /** Build a new instance. + *

+ * The template orbit is used as a model to {@link + * #createInitialOrbit() create initial orbit}. It defines the + * inertial frame, the central attraction coefficient, the orbit type, and is also + * used together with the {@code positionScale} to convert from the {@link + * org.orekit.utils.ParameterDriver#setNormalizedValue(double) normalized} parameters used by the + * callers of this builder to the real orbital parameters. + *

+ * @param templateOrbit reference orbit from which real orbits will be built + * (note that the mu from this orbit will be overridden with the mu from the + * {@code provider}) + * @param provider for un-normalized zonal coefficients + * @param positionAngle position angle type to use + * @param positionScale scaling factor used for orbital parameters normalization + * (typically set to the expected standard deviation of the position) + * @param M2 value of empirical drag coefficient. + * If equal to {@link BrouwerLyddanePropagator#M2} drag is not computed + * @param attitudeProvider attitude law to use + */ + public BrouwerLyddanePropagatorBuilder(final Orbit templateOrbit, + final UnnormalizedSphericalHarmonicsProvider provider, + final PositionAngle positionAngle, + final double positionScale, + final AttitudeProvider attitudeProvider, + final double M2) { + super(overrideMu(templateOrbit, provider, positionAngle), positionAngle, positionScale, true, attitudeProvider); + this.provider = provider; + // initialize M2 driver + final ParameterDriver M2Driver = new ParameterDriver(BrouwerLyddanePropagator.M2_NAME, M2, SCALE, + Double.NEGATIVE_INFINITY, + Double.POSITIVE_INFINITY); + addSupportedParameter(M2Driver); } /** Override central attraction coefficient. @@ -184,9 +230,26 @@ public class BrouwerLyddanePropagatorBuilder extends AbstractPropagatorBuilder { } /** {@inheritDoc} */ - public Propagator buildPropagator(final double[] normalizedParameters) { + public BrouwerLyddanePropagator buildPropagator(final double[] normalizedParameters) { setParameters(normalizedParameters); - return new BrouwerLyddanePropagator(createInitialOrbit(), getAttitudeProvider(), - provider); + + // Update M2 value and selection + double newM2 = 0.0; + boolean isSelected = false; + for (final ParameterDriver driver : getPropagationParametersDrivers().getDrivers()) { + if (BrouwerLyddanePropagator.M2_NAME.equals(driver.getName())) { + newM2 = driver.getValue(); + isSelected = driver.isSelected(); + } + } + + // Initialize propagator + final BrouwerLyddanePropagator propagator = new BrouwerLyddanePropagator(createInitialOrbit(), getAttitudeProvider(), provider, newM2); + propagator.getParametersDrivers().get(0).setSelected(isSelected); + + // Return + return propagator; + } + } diff --git a/src/site/markdown/index.md b/src/site/markdown/index.md index e8e90a4819444c961e83ec52900eb70debad779a..6a561625214895a1b3882cc6718fb13a729ae055 100644 --- a/src/site/markdown/index.md +++ b/src/site/markdown/index.md @@ -78,6 +78,7 @@ * Kepler * Eckstein-Heschler * Brouwer-Lyddane with Warren Phipps' correction for the critical inclination of 63.4° + and the perturbative acceleration due to atmospheric drag * SDP4/SGP4 with 2006 corrections * GNSS: GPS, QZSS, Galileo, GLONASS, Beidou, IRNSS and SBAS * numerical propagators diff --git a/src/test/java/org/orekit/propagation/analytical/BrouwerLyddanePropagatorTest.java b/src/test/java/org/orekit/propagation/analytical/BrouwerLyddanePropagatorTest.java index 9db7aee303a055b94d7fdc5fb4423c6c192789e8..97125da99ef663e07b4829f4f2d2e019467847be 100644 --- a/src/test/java/org/orekit/propagation/analytical/BrouwerLyddanePropagatorTest.java +++ b/src/test/java/org/orekit/propagation/analytical/BrouwerLyddanePropagatorTest.java @@ -1,5 +1,7 @@ package org.orekit.propagation.analytical; +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; @@ -11,8 +13,13 @@ import org.junit.Before; import org.junit.Test; import org.orekit.Utils; import org.orekit.attitudes.AttitudeProvider; +import org.orekit.bodies.CelestialBodyFactory; +import org.orekit.bodies.OneAxisEllipsoid; +import org.orekit.data.DataContext; import org.orekit.errors.OrekitException; import org.orekit.forces.ForceModel; +import org.orekit.forces.drag.DragForce; +import org.orekit.forces.drag.IsotropicDrag; import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel; import org.orekit.forces.gravity.potential.GravityFieldFactory; import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider; @@ -20,6 +27,8 @@ import org.orekit.forces.gravity.potential.TideSystem; import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider; import org.orekit.frames.Frame; import org.orekit.frames.FramesFactory; +import org.orekit.models.earth.atmosphere.DTM2000; +import org.orekit.models.earth.atmosphere.data.MarshallSolarActivityFutureEstimation; import org.orekit.orbits.EquinoctialOrbit; import org.orekit.orbits.KeplerianOrbit; import org.orekit.orbits.Orbit; @@ -30,6 +39,9 @@ import org.orekit.propagation.Propagator; import org.orekit.propagation.SpacecraftState; import org.orekit.propagation.numerical.NumericalPropagator; import org.orekit.time.AbsoluteDate; +import org.orekit.time.TimeScale; +import org.orekit.time.TimeScalesFactory; +import org.orekit.utils.Constants; import org.orekit.utils.IERSConventions; import org.orekit.utils.PVCoordinates; @@ -54,7 +66,7 @@ public class BrouwerLyddanePropagatorTest { // Extrapolation at the initial date // --------------------------------- BrouwerLyddanePropagator extrapolator = - new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider)); + new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2); SpacecraftState finalOrbit = extrapolator.propagate(initDate); // positions velocity and semi major axis match perfectly @@ -82,7 +94,7 @@ public class BrouwerLyddanePropagatorTest { FramesFactory.getEME2000(), initDate, provider.getMu()); BrouwerLyddanePropagator extrapolator = - new BrouwerLyddanePropagator(initialOrbit, DEFAULT_LAW, GravityFieldFactory.getUnnormalizedProvider(provider)); + new BrouwerLyddanePropagator(initialOrbit, DEFAULT_LAW, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2); SpacecraftState finalOrbit = extrapolator.propagate(initDate); @@ -131,7 +143,7 @@ public class BrouwerLyddanePropagatorTest { // Extrapolators definitions // ------------------------- BrouwerLyddanePropagator extrapolatorAna = - new BrouwerLyddanePropagator(initialOrbit, 1000.0, kepProvider); + new BrouwerLyddanePropagator(initialOrbit, 1000.0, kepProvider, BrouwerLyddanePropagator.M2); KeplerianPropagator extrapolatorKep = new KeplerianPropagator(initialOrbit); // Extrapolation at a final date different from initial date @@ -225,7 +237,7 @@ public class BrouwerLyddanePropagatorTest { //_______________________________________________________________________________________________ BrouwerLyddanePropagator BLextrapolator = - new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider)); + new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2); SpacecraftState BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift)); final KeplerianOrbit BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit()); @@ -241,6 +253,103 @@ public class BrouwerLyddanePropagatorTest { MathUtils.normalizeAngle(BLOrbit.getTrueAnomaly(), FastMath.PI), 0.12); } + @Test + public void compareToNumericalPropagationWithDrag() { + + final Frame inertialFrame = FramesFactory.getEME2000(); + final TimeScale utc = TimeScalesFactory.getUTC(); + final AbsoluteDate initDate = new AbsoluteDate(2003, 1, 1, 00, 00, 00.000, utc); + double timeshift = 60000. ; + + // Initial orbit + final double a = Constants.WGS84_EARTH_EQUATORIAL_RADIUS + 400e3; // semi major axis in meters + final double e = 0.01; // eccentricity + final double i = FastMath.toRadians(7); // inclination + final double omega = FastMath.toRadians(180); // perigee argument + final double raan = FastMath.toRadians(261); // right ascention of ascending node + final double lM = 0; // mean anomaly + final Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.TRUE, + inertialFrame, initDate, provider.getMu()); + // Initial state definition + final SpacecraftState initialState = new SpacecraftState(initialOrbit); + + //_______________________________________________________________________________________________ + // SET UP A REFERENCE NUMERICAL PROPAGATION + //_______________________________________________________________________________________________ + + // Adaptive step integrator with a minimum step of 0.001 and a maximum step of 1000 + final double minStep = 0.001; + final double maxstep = 1000.0; + final double positionTolerance = 10.0; + final OrbitType propagationType = OrbitType.KEPLERIAN; + final double[][] tolerances = + NumericalPropagator.tolerances(positionTolerance, initialOrbit, propagationType); + final AdaptiveStepsizeIntegrator integrator = + new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]); + + // Numerical Propagator + final NumericalPropagator NumPropagator = new NumericalPropagator(integrator); + NumPropagator.setOrbitType(propagationType); + + // Atmosphere + final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, + FramesFactory.getITRF(IERSConventions.IERS_2010, true)); + MarshallSolarActivityFutureEstimation msafe = + new MarshallSolarActivityFutureEstimation("Jan2000F10-edited-data\\.txt", + MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE); + DataContext.getDefault().getDataProvidersManager().feed(msafe.getSupportedNames(), msafe); + DTM2000 atmosphere = new DTM2000(msafe, CelestialBodyFactory.getSun(), earth); + + // Force model + final ForceModel holmesFeatherstone = + new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider); + final ForceModel drag = + new DragForce(atmosphere, new IsotropicDrag(1.0, 1.0)); + NumPropagator.addForceModel(holmesFeatherstone); + NumPropagator.addForceModel(drag); + + // Set up initial state in the propagator + NumPropagator.setInitialState(initialState); + + // Extrapolate from the initial to the final date + final SpacecraftState NumFinalState = NumPropagator.propagate(initDate.shiftedBy(timeshift)); + final KeplerianOrbit NumOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(NumFinalState.getOrbit()); + + //_______________________________________________________________________________________________ + // SET UP A BROUWER LYDDANE PROPAGATION WITHOUT DRAG + //_______________________________________________________________________________________________ + + BrouwerLyddanePropagator BLextrapolator = + new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2); + + SpacecraftState BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift)); + KeplerianOrbit BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit()); + + // Verify a and e differences without the drag effect on Brouwer-Lyddane + final double deltaSmaBefore = 20.44; + final double deltaEccBefore = 1.0125e-4; + Assert.assertEquals(NumOrbit.getA(), BLOrbit.getA(), deltaSmaBefore); + Assert.assertEquals(NumOrbit.getE(), BLOrbit.getE(), deltaEccBefore); + + //_______________________________________________________________________________________________ + // SET UP A BROUWER LYDDANE PROPAGATION WITH DRAG + //_______________________________________________________________________________________________ + + double M2 = 1.0e-14; + BLextrapolator = new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), M2); + BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift)); + BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit()); + + // Verify a and e differences without the drag effect on Brouwer-Lyddane + final double deltaSmaAfter = 15.66; + final double deltaEccAfter = 1.0121e-4; + Assert.assertEquals(NumOrbit.getA(), BLOrbit.getA(), deltaSmaAfter); + Assert.assertEquals(NumOrbit.getE(), BLOrbit.getE(), deltaEccAfter); + assertTrue(deltaSmaAfter < deltaSmaBefore); + assertTrue(deltaEccAfter < deltaEccBefore); + + } + @Test public void compareToNumericalPropagationMeanInitialOrbit() { @@ -262,7 +371,7 @@ public class BrouwerLyddanePropagatorTest { BrouwerLyddanePropagator BLextrapolator = new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), - PropagationType.MEAN); + PropagationType.MEAN, BrouwerLyddanePropagator.M2); SpacecraftState initialOsculatingState = BLextrapolator.propagate(initDate); final KeplerianOrbit InitOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(initialOsculatingState.getOrbit()); @@ -341,14 +450,14 @@ public class BrouwerLyddanePropagatorTest { BrouwerLyddanePropagator BLextrapolator1 = new BrouwerLyddanePropagator(initialOrbit, DEFAULT_LAW, Propagator.DEFAULT_MASS, GravityFieldFactory.getUnnormalizedProvider(provider), - PropagationType.OSCULATING); + PropagationType.OSCULATING, BrouwerLyddanePropagator.M2); //_______________________________________________________________________________________________ // SET UP ANOTHER BROUWER LYDDANE PROPAGATOR //_______________________________________________________________________________________________ BrouwerLyddanePropagator BLextrapolator2 = new BrouwerLyddanePropagator( new KeplerianOrbit(a + 3000, e + 0.001, i - FastMath.toRadians(12.0), omega, raan, lM, PositionAngle.TRUE, - inertialFrame, initDate, provider.getMu()),DEFAULT_LAW, Propagator.DEFAULT_MASS, GravityFieldFactory.getUnnormalizedProvider(provider)); + inertialFrame, initDate, provider.getMu()),DEFAULT_LAW, Propagator.DEFAULT_MASS, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2); // Reset BL2 with BL1 initial state BLextrapolator2.resetInitialState(initialState); @@ -389,11 +498,11 @@ public class BrouwerLyddanePropagatorTest { BrouwerLyddanePropagator BLPropagator1 = new BrouwerLyddanePropagator(initialOrbit, DEFAULT_LAW, - provider.getAe(), provider.getMu(), -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7); + provider.getAe(), provider.getMu(), -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2); BrouwerLyddanePropagator BLPropagator2 = new BrouwerLyddanePropagator(initialOrbit, - provider.getAe(), provider.getMu(), -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7); + provider.getAe(), provider.getMu(), -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2); BrouwerLyddanePropagator BLPropagator3 = new BrouwerLyddanePropagator(initialOrbit, Propagator.DEFAULT_MASS, - provider.getAe(), provider.getMu(), -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7); + provider.getAe(), provider.getMu(), -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2); SpacecraftState BLFinalState1 = BLPropagator1.propagate(initDate.shiftedBy(timeshift)); final KeplerianOrbit BLOrbit1 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState1.getOrbit()); @@ -438,7 +547,7 @@ public class BrouwerLyddanePropagatorTest { // ----------------------- BrouwerLyddanePropagator extrapolator = new BrouwerLyddanePropagator(initialOrbit, DEFAULT_LAW, provider.getAe(), provider.getMu(), - -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7); + -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2); // Extrapolation at the initial date // --------------------------------- @@ -459,7 +568,7 @@ public class BrouwerLyddanePropagatorTest { // ----------------------- BrouwerLyddanePropagator extrapolator = new BrouwerLyddanePropagator(initialOrbit, provider.getAe(), provider.getMu(), - -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7); + -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2); // Extrapolation at the initial date // --------------------------------- @@ -488,7 +597,7 @@ public class BrouwerLyddanePropagatorTest { // Extrapolator definition // ----------------------- BrouwerLyddanePropagator extrapolator = - new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider)); + new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2); // Extrapolation at the initial date // --------------------------------- @@ -547,7 +656,7 @@ public class BrouwerLyddanePropagatorTest { // Extrapolator definition // ----------------------- BrouwerLyddanePropagator extrapolator = - new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider)); + new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2); // Extrapolation at the initial date // --------------------------------- @@ -559,7 +668,7 @@ public class BrouwerLyddanePropagatorTest { @Before public void setUp() { - Utils.setDataRoot("regular-data:potential/icgem-format"); + Utils.setDataRoot("regular-data:atmosphere:potential/icgem-format"); provider = GravityFieldFactory.getNormalizedProvider(5, 0); } diff --git a/src/test/java/org/orekit/propagation/analytical/FieldBrouwerLyddanePropagatorTest.java b/src/test/java/org/orekit/propagation/analytical/FieldBrouwerLyddanePropagatorTest.java index 0470ecc83179795300cc35d92493033cdd42f5b5..3ce417eb9c70e9edfbfa805f742d79b3450cedfa 100644 --- a/src/test/java/org/orekit/propagation/analytical/FieldBrouwerLyddanePropagatorTest.java +++ b/src/test/java/org/orekit/propagation/analytical/FieldBrouwerLyddanePropagatorTest.java @@ -14,9 +14,14 @@ import org.junit.Before; import org.junit.Test; import org.orekit.Utils; import org.orekit.attitudes.AttitudeProvider; +import org.orekit.bodies.CelestialBodyFactory; +import org.orekit.bodies.OneAxisEllipsoid; +import org.orekit.data.DataContext; import org.orekit.errors.OrekitException; import org.orekit.errors.OrekitMessages; import org.orekit.forces.ForceModel; +import org.orekit.forces.drag.DragForce; +import org.orekit.forces.drag.IsotropicDrag; import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel; import org.orekit.forces.gravity.potential.GravityFieldFactory; import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider; @@ -24,6 +29,8 @@ import org.orekit.forces.gravity.potential.TideSystem; import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider; import org.orekit.frames.Frame; import org.orekit.frames.FramesFactory; +import org.orekit.models.earth.atmosphere.DTM2000; +import org.orekit.models.earth.atmosphere.data.MarshallSolarActivityFutureEstimation; import org.orekit.orbits.FieldEquinoctialOrbit; import org.orekit.orbits.FieldKeplerianOrbit; import org.orekit.orbits.FieldOrbit; @@ -35,7 +42,11 @@ import org.orekit.propagation.PropagationType; import org.orekit.propagation.Propagator; import org.orekit.propagation.SpacecraftState; import org.orekit.propagation.numerical.NumericalPropagator; +import org.orekit.time.AbsoluteDate; import org.orekit.time.FieldAbsoluteDate; +import org.orekit.time.TimeScale; +import org.orekit.time.TimeScalesFactory; +import org.orekit.utils.Constants; import org.orekit.utils.FieldPVCoordinates; import org.orekit.utils.IERSConventions; @@ -64,7 +75,7 @@ public class FieldBrouwerLyddanePropagatorTest { // Extrapolation at the initial date // --------------------------------- FieldBrouwerLyddanePropagator extrapolator = - new FieldBrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider)); + new FieldBrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2); FieldSpacecraftState finalOrbit = extrapolator.propagate(initDate); // positions velocity and semi major axis match perfectly @@ -98,7 +109,7 @@ public class FieldBrouwerLyddanePropagatorTest { FramesFactory.getEME2000(), initDate, zero.add(provider.getMu())); FieldBrouwerLyddanePropagator extrapolator = - new FieldBrouwerLyddanePropagator(initialOrbit, DEFAULT_LAW, GravityFieldFactory.getUnnormalizedProvider(provider)); + new FieldBrouwerLyddanePropagator(initialOrbit, DEFAULT_LAW, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2); FieldSpacecraftState finalOrbit = extrapolator.propagate(initDate); @@ -154,7 +165,7 @@ public class FieldBrouwerLyddanePropagatorTest { // Extrapolators definitions // ------------------------- FieldBrouwerLyddanePropagator extrapolatorAna = - new FieldBrouwerLyddanePropagator<>(initialOrbit, zero.add(1000.0), kepProvider); + new FieldBrouwerLyddanePropagator<>(initialOrbit, zero.add(1000.0), kepProvider, BrouwerLyddanePropagator.M2); FieldKeplerianPropagator extrapolatorKep = new FieldKeplerianPropagator<>(initialOrbit); // Extrapolation at a final date different from initial date @@ -253,7 +264,7 @@ public class FieldBrouwerLyddanePropagatorTest { //_______________________________________________________________________________________________ FieldBrouwerLyddanePropagator BLextrapolator = - new FieldBrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider)); + new FieldBrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2); FieldSpacecraftState BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift)); final KeplerianOrbit BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit().toOrbit()); @@ -270,6 +281,110 @@ public class FieldBrouwerLyddanePropagatorTest { MathUtils.normalizeAngle(BLOrbit.getTrueAnomaly(), FastMath.PI), 0.12); } + @Test + public void compareToNumericalPropagationWithDrag() { + doCompareToNumericalPropagationWithDrag(Decimal64Field.getInstance()); + } + + private > void doCompareToNumericalPropagationWithDrag(Field field) { + + T zero = field.getZero(); + final Frame inertialFrame = FramesFactory.getEME2000(); + final TimeScale utc = TimeScalesFactory.getUTC(); + final AbsoluteDate date = new AbsoluteDate(2003, 1, 1, 00, 00, 00.000, utc); + final FieldAbsoluteDate initDate = new FieldAbsoluteDate<>(date, zero); + double timeshift = 60000. ; + + // Initial orbit + final double a = Constants.WGS84_EARTH_EQUATORIAL_RADIUS + 400e3; // semi major axis in meters + final double e = 0.01; // eccentricity + final double i = FastMath.toRadians(7); // inclination + final double omega = FastMath.toRadians(180); // perigee argument + final double raan = FastMath.toRadians(261); // right ascention of ascending node + final double lM = 0; // mean anomaly + final FieldOrbit initialOrbit = new FieldKeplerianOrbit<>(zero.add(a), zero.add(e), zero.add(i), zero.add(omega), + zero.add(raan), zero.add(lM), PositionAngle.TRUE, + inertialFrame, initDate, zero.add(provider.getMu())); + // Initial state definition + final FieldSpacecraftState initialState = new FieldSpacecraftState<>(initialOrbit); + + //_______________________________________________________________________________________________ + // SET UP A REFERENCE NUMERICAL PROPAGATION + //_______________________________________________________________________________________________ + + // Adaptive step integrator with a minimum step of 0.001 and a maximum step of 1000 + final double minStep = 0.001; + final double maxstep = 1000.0; + final double positionTolerance = 10.0; + final OrbitType propagationType = OrbitType.KEPLERIAN; + final double[][] tolerances = + NumericalPropagator.tolerances(positionTolerance, initialOrbit.toOrbit(), propagationType); + final AdaptiveStepsizeIntegrator integrator = + new DormandPrince853Integrator(minStep, maxstep, tolerances[0], tolerances[1]); + + // Numerical Propagator + final NumericalPropagator NumPropagator = new NumericalPropagator(integrator); + NumPropagator.setOrbitType(propagationType); + + // Atmosphere + final OneAxisEllipsoid earth = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, + FramesFactory.getITRF(IERSConventions.IERS_2010, true)); + MarshallSolarActivityFutureEstimation msafe = + new MarshallSolarActivityFutureEstimation("Jan2000F10-edited-data\\.txt", + MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE); + DataContext.getDefault().getDataProvidersManager().feed(msafe.getSupportedNames(), msafe); + DTM2000 atmosphere = new DTM2000(msafe, CelestialBodyFactory.getSun(), earth); + + // Force model + final ForceModel holmesFeatherstone = + new HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, true), provider); + final ForceModel drag = + new DragForce(atmosphere, new IsotropicDrag(1.0, 1.0)); + NumPropagator.addForceModel(holmesFeatherstone); + NumPropagator.addForceModel(drag); + + // Set up initial state in the propagator + NumPropagator.setInitialState(initialState.toSpacecraftState()); + + // Extrapolate from the initial to the final date + final SpacecraftState NumFinalState = NumPropagator.propagate(initDate.toAbsoluteDate().shiftedBy(timeshift)); + final KeplerianOrbit NumOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(NumFinalState.getOrbit()); + + //_______________________________________________________________________________________________ + // SET UP A BROUWER LYDDANE PROPAGATION WITHOUT DRAG + //_______________________________________________________________________________________________ + + FieldBrouwerLyddanePropagator BLextrapolator = + new FieldBrouwerLyddanePropagator<>(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2); + + FieldSpacecraftState BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift)); + KeplerianOrbit BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit().toOrbit()); + + // Verify a and e differences without the drag effect on Brouwer-Lyddane + final double deltaSmaBefore = 20.44; + final double deltaEccBefore = 1.0125e-4; + Assert.assertEquals(NumOrbit.getA(), BLOrbit.getA(), deltaSmaBefore); + Assert.assertEquals(NumOrbit.getE(), BLOrbit.getE(), deltaEccBefore); + + //_______________________________________________________________________________________________ + // SET UP A BROUWER LYDDANE PROPAGATION WITH DRAG + //_______________________________________________________________________________________________ + + double M2 = 1.0e-14; + BLextrapolator = new FieldBrouwerLyddanePropagator<>(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), M2); + BLFinalState = BLextrapolator.propagate(initDate.shiftedBy(timeshift)); + BLOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState.getOrbit().toOrbit()); + + // Verify a and e differences without the drag effect on Brouwer-Lyddane + final double deltaSmaAfter = 15.66; + final double deltaEccAfter = 1.0121e-4; + Assert.assertEquals(NumOrbit.getA(), BLOrbit.getA(), deltaSmaAfter); + Assert.assertEquals(NumOrbit.getE(), BLOrbit.getE(), deltaEccAfter); + Assert.assertTrue(deltaSmaAfter < deltaSmaBefore); + Assert.assertTrue(deltaEccAfter < deltaEccBefore); + Assert.assertEquals(M2, BLextrapolator.getM2(), Double.MIN_VALUE); + + } @Test public void compareToNumericalPropagationMeanInitialOrbit() { @@ -297,7 +412,7 @@ public class FieldBrouwerLyddanePropagatorTest { FieldBrouwerLyddanePropagator BLextrapolator = new FieldBrouwerLyddanePropagator<>(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), - PropagationType.MEAN); + PropagationType.MEAN, BrouwerLyddanePropagator.M2); FieldSpacecraftState initialOsculatingState = BLextrapolator.propagate(initDate); final KeplerianOrbit InitOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(initialOsculatingState.getOrbit().toOrbit()); @@ -381,7 +496,7 @@ public class FieldBrouwerLyddanePropagatorTest { FieldBrouwerLyddanePropagator BLextrapolator1 = new FieldBrouwerLyddanePropagator<>(initialOrbit, DEFAULT_LAW, zero.add(Propagator.DEFAULT_MASS), GravityFieldFactory.getUnnormalizedProvider(provider), - PropagationType.OSCULATING); + PropagationType.OSCULATING, BrouwerLyddanePropagator.M2); //_______________________________________________________________________________________________ // SET UP ANOTHER BROUWER LYDDANE PROPAGATOR //_______________________________________________________________________________________________ @@ -392,7 +507,7 @@ public class FieldBrouwerLyddanePropagatorTest { zero.add(raan), zero.add(lM), PositionAngle.TRUE, inertialFrame, initDate, zero.add(provider.getMu())),DEFAULT_LAW, zero.add(Propagator.DEFAULT_MASS), - GravityFieldFactory.getUnnormalizedProvider(provider)); + GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2); // Reset BL2 with BL1 initial state BLextrapolator2.resetInitialState(initialState); @@ -440,12 +555,12 @@ public class FieldBrouwerLyddanePropagatorTest { FieldBrouwerLyddanePropagator BLPropagator1 = new FieldBrouwerLyddanePropagator(initialOrbit, DEFAULT_LAW, - provider.getAe(), zero.add(provider.getMu()), -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7); + provider.getAe(), zero.add(provider.getMu()), -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2); FieldBrouwerLyddanePropagator BLPropagator2 = new FieldBrouwerLyddanePropagator<>(initialOrbit, - provider.getAe(), zero.add(provider.getMu()), -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7); + provider.getAe(), zero.add(provider.getMu()), -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2); FieldBrouwerLyddanePropagator BLPropagator3 = new FieldBrouwerLyddanePropagator<>(initialOrbit, zero.add(Propagator.DEFAULT_MASS), provider.getAe(), zero.add(provider.getMu()), -1.08263e-3, - 2.54e-6, 1.62e-6, 2.3e-7); + 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2); FieldSpacecraftState BLFinalState1 = BLPropagator1.propagate(initDate.shiftedBy(timeshift)); final KeplerianOrbit BLOrbit1 = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(BLFinalState1.getOrbit().toOrbit()); @@ -499,7 +614,7 @@ public class FieldBrouwerLyddanePropagatorTest { FieldBrouwerLyddanePropagator extrapolator = new FieldBrouwerLyddanePropagator<>(initialOrbit, DEFAULT_LAW, provider.getAe(), zero.add(provider.getMu()), - -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7); + -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2); // Extrapolation at the initial date // --------------------------------- @@ -530,7 +645,7 @@ public class FieldBrouwerLyddanePropagatorTest { // ----------------------- FieldBrouwerLyddanePropagator extrapolator = new FieldBrouwerLyddanePropagator<>(initialOrbit, provider.getAe(), zero.add(provider.getMu()), - -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7); + -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7, BrouwerLyddanePropagator.M2); // Extrapolation at the initial date // --------------------------------- @@ -569,7 +684,7 @@ public class FieldBrouwerLyddanePropagatorTest { // Extrapolator definition // ----------------------- FieldBrouwerLyddanePropagator extrapolator = - new FieldBrouwerLyddanePropagator<>(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider)); + new FieldBrouwerLyddanePropagator<>(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2); // Extrapolation at the initial date // --------------------------------- @@ -615,7 +730,7 @@ public class FieldBrouwerLyddanePropagatorTest { // Extrapolator definition // ----------------------- - final FieldBrouwerLyddanePropagator blField = new FieldBrouwerLyddanePropagator<>(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider)); + final FieldBrouwerLyddanePropagator blField = new FieldBrouwerLyddanePropagator<>(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2); // Extrapolation at the initial date // --------------------------------- @@ -656,7 +771,7 @@ public class FieldBrouwerLyddanePropagatorTest { // Field propagation final FieldBrouwerLyddanePropagator blField = new FieldBrouwerLyddanePropagator<>(initialStateField.getOrbit(), GravityFieldFactory.getUnnormalizedProvider(provider), - PropagationType.MEAN); + PropagationType.MEAN, BrouwerLyddanePropagator.M2); final FieldSpacecraftState finalStateField = blField.propagate(initialStateField.getDate().shiftedBy(timeshift)); final FieldKeplerianOrbit finalOrbitField = (FieldKeplerianOrbit) OrbitType.KEPLERIAN.convertType(finalStateField.getOrbit()); final KeplerianOrbit finalOrbitFieldReal = finalOrbitField.toOrbit(); @@ -664,7 +779,7 @@ public class FieldBrouwerLyddanePropagatorTest { // Classical propagation final BrouwerLyddanePropagator bl = new BrouwerLyddanePropagator(initialState.getOrbit(), GravityFieldFactory.getUnnormalizedProvider(provider), - PropagationType.MEAN); + PropagationType.MEAN, BrouwerLyddanePropagator.M2); final SpacecraftState finalState = bl.propagate(initialState.getDate().shiftedBy(timeshift)); final KeplerianOrbit finalOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(finalState.getOrbit()); @@ -708,14 +823,15 @@ public class FieldBrouwerLyddanePropagatorTest { // Field propagation final FieldBrouwerLyddanePropagator blField = new FieldBrouwerLyddanePropagator<>(initialStateField.getOrbit(), - GravityFieldFactory.getUnnormalizedProvider(provider)); + GravityFieldFactory.getUnnormalizedProvider(provider), BrouwerLyddanePropagator.M2); final FieldSpacecraftState finalStateField = blField.propagate(initialStateField.getDate().shiftedBy(timeshift)); final FieldKeplerianOrbit finalOrbitField = (FieldKeplerianOrbit) OrbitType.KEPLERIAN.convertType(finalStateField.getOrbit()); final KeplerianOrbit finalOrbitFieldReal = finalOrbitField.toOrbit(); // Classical propagation final BrouwerLyddanePropagator bl = new BrouwerLyddanePropagator(initialState.getOrbit(), - GravityFieldFactory.getUnnormalizedProvider(provider)); + GravityFieldFactory.getUnnormalizedProvider(provider), + BrouwerLyddanePropagator.M2); final SpacecraftState finalState = bl.propagate(initialState.getDate().shiftedBy(timeshift)); final KeplerianOrbit finalOrbit = (KeplerianOrbit) OrbitType.KEPLERIAN.convertType(finalState.getOrbit()); @@ -731,7 +847,7 @@ public class FieldBrouwerLyddanePropagatorTest { @Before public void setUp() { - Utils.setDataRoot("regular-data:potential/icgem-format"); + Utils.setDataRoot("regular-data:atmosphere:potential/icgem-format"); provider = GravityFieldFactory.getNormalizedProvider(5, 0); } diff --git a/src/test/java/org/orekit/propagation/conversion/BrouwerLyddanePropagatorBuilderTest.java b/src/test/java/org/orekit/propagation/conversion/BrouwerLyddanePropagatorBuilderTest.java index 88bae96393c42cc47bd66547e81b6434f2b00004..9b39a9aa60cf886219f542dab9720d675d027538 100644 --- a/src/test/java/org/orekit/propagation/conversion/BrouwerLyddanePropagatorBuilderTest.java +++ b/src/test/java/org/orekit/propagation/conversion/BrouwerLyddanePropagatorBuilderTest.java @@ -22,8 +22,6 @@ import org.junit.Before; import org.junit.Test; import org.orekit.Utils; import org.orekit.forces.gravity.potential.GravityFieldFactory; -import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider; -import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider.NormalizedSphericalHarmonics; import org.orekit.forces.gravity.potential.TideSystem; import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider; import org.orekit.forces.gravity.potential.UnnormalizedSphericalHarmonicsProvider.UnnormalizedSphericalHarmonics; @@ -33,9 +31,9 @@ import org.orekit.orbits.KeplerianOrbit; import org.orekit.orbits.Orbit; import org.orekit.orbits.OrbitType; import org.orekit.orbits.PositionAngle; -import org.orekit.propagation.Propagator; import org.orekit.propagation.analytical.BrouwerLyddanePropagator; import org.orekit.time.AbsoluteDate; +import org.orekit.utils.ParameterDriver; public class BrouwerLyddanePropagatorBuilderTest { @@ -49,15 +47,12 @@ public class BrouwerLyddanePropagatorBuilderTest { // Define initial state and BrouwerLyddane Propagator AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH.shiftedBy(583.); - BrouwerLyddanePropagator propagator = new BrouwerLyddanePropagator(orbit, provider); + BrouwerLyddanePropagator propagator = new BrouwerLyddanePropagator(orbit, provider, BrouwerLyddanePropagator.M2); // We propagate using directly the propagator of the set up final Orbit orbitWithPropagator = propagator.propagate(initDate.shiftedBy(60000)).getOrbit(); - // Convert povider to normalized provider to be able to build a Brouwer Lyddane propagator UnnormalizedSphericalHarmonics harmonics = provider.onDate(orbit.getDate()); - NormalizedSphericalHarmonicsProvider normProv = GravityFieldFactory.getNormalizedProvider(provider); - NormalizedSphericalHarmonics normHarmonics = normProv.onDate(orbit.getDate()); // We propagate using a build version of the propagator // We shall have the same results than before @@ -71,11 +66,10 @@ public class BrouwerLyddanePropagatorBuilderTest { harmonics.getUnnormalizedCnm(5, 0), OrbitType.KEPLERIAN, PositionAngle.TRUE, - 1.0); + 1.0, + BrouwerLyddanePropagator.M2); - double[] normCnm = new double[] - { 0, 0, normHarmonics.getNormalizedCnm(2, 0), normHarmonics.getNormalizedCnm(3, 0), normHarmonics.getNormalizedCnm(4, 0), normHarmonics.getNormalizedCnm(5, 0)}; - final Propagator prop = builder.buildPropagator(normCnm); + final BrouwerLyddanePropagator prop = builder.buildPropagator(builder.getSelectedNormalizedParameters()); final Orbit orbitWithBuilder = prop.propagate(initDate.shiftedBy(60000)).getOrbit(); // Verify @@ -85,8 +79,47 @@ public class BrouwerLyddanePropagatorBuilderTest { Assert.assertEquals(orbitWithPropagator.getHx(), orbitWithBuilder.getHx(), eps); Assert.assertEquals(orbitWithPropagator.getHy(), orbitWithBuilder.getHy(), eps); Assert.assertEquals(orbitWithPropagator.getLM(), orbitWithBuilder.getLM(), 8.0e-10); + } + @Test + public void doTestBuildPropagatorWithDrag() { + + // M2 + final double M2 = 1.0e-15; + + // Convert provider to normalized provider to be able to build a Brouwer Lyddane propagator + UnnormalizedSphericalHarmonics harmonics = provider.onDate(orbit.getDate()); + + // Initialize propagator builder + BrouwerLyddanePropagatorBuilder builder = new BrouwerLyddanePropagatorBuilder(orbit, + provider.getAe(), + provider.getMu(), + provider.getTideSystem(), + harmonics.getUnnormalizedCnm(2, 0), + harmonics.getUnnormalizedCnm(3, 0), + harmonics.getUnnormalizedCnm(4, 0), + harmonics.getUnnormalizedCnm(5, 0), + OrbitType.KEPLERIAN, + PositionAngle.TRUE, + 1.0, + M2); + + // Set the M2 parameter to selected + for (ParameterDriver driver : builder.getPropagationParametersDrivers().getDrivers()) { + if (BrouwerLyddanePropagator.M2_NAME.equals(driver.getName())) { + driver.setSelected(true); + } + } + + // Build the propagator + final BrouwerLyddanePropagator prop = builder.buildPropagator(builder.getSelectedNormalizedParameters()); + + // Verify + Assert.assertEquals(M2, prop.getM2(), Double.MIN_VALUE); + Assert.assertTrue(prop.getParametersDrivers().get(0).isSelected()); + + } @Before public void setUp() {