diff --git a/src/changes/changes.xml b/src/changes/changes.xml index 115f408ba86aeb6d57f1eaae5550777fe4ae0281..2d484900ebd80eb31ea5fb2fb31d93e477c0a5bb 100644 --- a/src/changes/changes.xml +++ b/src/changes/changes.xml @@ -21,6 +21,9 @@ + + Allowed Brouwer-Lyddane model to work for the critical inclination. + Fixed handling of multiple historical eccentricities for a same station. diff --git a/src/main/java/org/orekit/propagation/analytical/BrouwerLyddanePropagator.java b/src/main/java/org/orekit/propagation/analytical/BrouwerLyddanePropagator.java index e0530c07cfabef73a32e6ffe4e31667a086220ef..d54ebfbb3f8584721dff5c05b8775531332e0253 100644 --- a/src/main/java/org/orekit/propagation/analytical/BrouwerLyddanePropagator.java +++ b/src/main/java/org/orekit/propagation/analytical/BrouwerLyddanePropagator.java @@ -19,6 +19,7 @@ package org.orekit.propagation.analytical; import java.io.Serializable; import org.hipparchus.analysis.differentiation.UnivariateDerivative2; +import org.hipparchus.util.CombinatoricsUtils; import org.hipparchus.util.FastMath; import org.hipparchus.util.MathUtils; import org.hipparchus.util.Precision; @@ -46,8 +47,8 @@ import org.orekit.utils.TimeSpanMap; *

* At the opposite of the {@link EcksteinHechlerPropagator}, the Brouwer-Lyddane model is * 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). One needs still to be careful with eccentricities - * lower than 5e-4. The computation should not be done for the critical inclination : 63.4°. + * (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. *

* @see "Brouwer, Dirk. Solution of the problem of artificial satellite theory without drag. * YALE UNIV NEW HAVEN CT NEW HAVEN United States, 1959." @@ -55,14 +56,18 @@ import org.orekit.utils.TimeSpanMap; * @see "Lyddane, R. H. Small eccentricities or inclinations in the Brouwer theory of the * artificial satellite. The Astronomical Journal 68 (1963): 555." * - * @see "Parks, A. D. A drag-augmented Brouwer-Lyddane artificial satellite theory and its - * application to long-term station alert predictions. NAVAL SURFACE WEAPONS CENTER DAHLGREN VA, 1983." + * @see "Phipps Jr, Warren E. Parallelization of the Navy Space Surveillance Center + * (NAVSPASUR) Satellite Model. NAVAL POSTGRADUATE SCHOOL MONTEREY CA, 1992." * * @author Melina Vanel + * @author Bryan Cazabonne * @since 11.1 */ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { + /** Beta constant used by T2 function. */ + private static final double BETA = FastMath.scalb(100, -11); + /** Initial Brouwer-Lyddane model. */ private BLModel initialModel; @@ -664,7 +669,7 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { final double cosI3 = cosI2 * cosI1; final double cosI4 = cosI2 * cosI2; final double cosI6 = cosI4 * cosI2; - final double C5c2 = 1.0 - 5.0 * cosI2; + final double C5c2 = 1.0 / T2(cosI1); final double C3c2 = 3.0 * cosI2 - 1.0; final double epp = mean.getE(); @@ -678,12 +683,6 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { mean.getE()); } - if (FastMath.abs(cosI2 - 1.0 / 5.0) < 1.0e-3) { - throw new OrekitException(OrekitMessages.ALMOST_CRITICALLY_INCLINED_ORBIT, - FastMath.toDegrees(mean.getI())); - } - - // secular multiplicative lt = 1 + 1.5 * yp2 * n * C3c2 + @@ -901,6 +900,41 @@ public class BrouwerLyddanePropagator extends AbstractAnalyticalPropagator { return ek; } + /** + * This method is used in Brouwer-Lyddane model to avoid singularity at the + * critical inclination (i = 63.4°). + *

+ * This method, based on Warren Phipps's 1992 thesis (Eq. 2.47 and 2.48), + * approximate the factor (1.0 - 5.0 * cos²(inc))^-1 (causing the singularity) + * by a function, named T2 in the thesis. + *

+ * @param cosInc cosine of the mean inclination + * @return an approximation of (1.0 - 5.0 * cos²(inc))^-1 term + */ + private double T2(final double cosInc) { + + // X = (1.0 - 5.0 * cos²(inc)) + final double x = 1.0 - 5.0 * cosInc * cosInc; + final double x2 = x * x; + + // Eq. 2.48 + double sum = 0.0; + for (int i = 0; i <= 12; i++) { + final double sign = i % 2 == 0 ? +1.0 : -1.0; + sum += sign * FastMath.pow(BETA, i) * FastMath.pow(x2, i) / CombinatoricsUtils.factorialDouble(i + 1); + } + + // Right term of equation 2.47 + double product = 1.0; + for (int i = 0; i <= 10; i++) { + product *= 1 + FastMath.exp(FastMath.scalb(-1.0, i) * BETA * x2); + } + + // Return (Eq. 2.47) + return BETA * x * sum * product; + + } + /** Extrapolate an orbit up to a specific target date. * @param date target date for the orbit * @return propagated parameters diff --git a/src/main/java/org/orekit/propagation/analytical/FieldBrouwerLyddanePropagator.java b/src/main/java/org/orekit/propagation/analytical/FieldBrouwerLyddanePropagator.java index 4604681c4af93f564ca3c400843db502242770b0..87bca775f61d92070a5617b067a84cd3eafc3142 100644 --- a/src/main/java/org/orekit/propagation/analytical/FieldBrouwerLyddanePropagator.java +++ b/src/main/java/org/orekit/propagation/analytical/FieldBrouwerLyddanePropagator.java @@ -22,6 +22,7 @@ import java.util.List; import org.hipparchus.CalculusFieldElement; import org.hipparchus.Field; import org.hipparchus.analysis.differentiation.FieldUnivariateDerivative2; +import org.hipparchus.util.CombinatoricsUtils; import org.hipparchus.util.FastMath; import org.hipparchus.util.FieldSinCos; import org.hipparchus.util.MathArrays; @@ -46,10 +47,10 @@ import org.orekit.utils.ParameterDriver; /** This class propagates a {@link org.orekit.propagation.FieldSpacecraftState} * using the analytical Brouwer-Lyddane model (from J2 to J5 zonal harmonics). *

- * At the opposite of the {@link EcksteinHechlerPropagator}, the Brouwer-Lyddane model is + * At the opposite of the {@link FieldEcksteinHechlerPropagator}, the Brouwer-Lyddane model is * 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). One needs still to be careful with eccentricities - * lower than 5e-4. The computation should not be done for the critical inclination : 63.4°. + * (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. *

* @see "Brouwer, Dirk. Solution of the problem of artificial satellite theory without drag. * YALE UNIV NEW HAVEN CT NEW HAVEN United States, 1959." @@ -57,14 +58,18 @@ import org.orekit.utils.ParameterDriver; * @see "Lyddane, R. H. Small eccentricities or inclinations in the Brouwer theory of the * artificial satellite. The Astronomical Journal 68 (1963): 555." * - * @see "Parks, A. D. A drag-augmented Brouwer-Lyddane artificial satellite theory and its - * application to long-term station alert predictions. NAVAL SURFACE WEAPONS CENTER DAHLGREN VA, 1983." + * @see "Phipps Jr, Warren E. Parallelization of the Navy Space Surveillance Center + * (NAVSPASUR) Satellite Model. NAVAL POSTGRADUATE SCHOOL MONTEREY CA, 1992." * * @author Melina Vanel + * @author Bryan Cazabonne * @since 11.1 */ public class FieldBrouwerLyddanePropagator> extends FieldAbstractAnalyticalPropagator { + /** Beta constant used by T2 function. */ + private static final double BETA = FastMath.scalb(100, -11); + /** Initial Brouwer-Lyddane model. */ private FieldBLModel initialModel; @@ -663,7 +668,7 @@ public class FieldBrouwerLyddanePropagator> ex final T cosI3 = cosI2.multiply(cosI1); final T cosI4 = cosI2.multiply(cosI2); final T cosI6 = cosI4.multiply(cosI2); - final T C5c2 = cosI2.multiply(-5.0).add(1.0); + final T C5c2 = T2(cosI1).reciprocal(); final T C3c2 = cosI2.multiply(3.0).subtract(1.0); final T epp = mean.getE(); @@ -677,12 +682,6 @@ public class FieldBrouwerLyddanePropagator> ex mean.getE().getReal()); } - if (cosI2.subtract(-0.2).abs().getReal() < 1.0e-3) { - throw new OrekitException(OrekitMessages.ALMOST_CRITICALLY_INCLINED_ORBIT, - FastMath.toDegrees(mean.getI().getReal())); - } - - // secular multiplicative lt = one.add(yp2.multiply(n).multiply(C3c2).multiply(1.5)). add(yp2.multiply(0.09375).multiply(yp2).multiply(n).multiply(n2.multiply(25.0).add(n.multiply(16.0)).add(-15.0). @@ -913,6 +912,41 @@ public class FieldBrouwerLyddanePropagator> ex return ek; } + /** + * This method is used in Brouwer-Lyddane model to avoid singularity at the + * critical inclination (i = 63.4°). + *

+ * This method, based on Warren Phipps's 1992 thesis (Eq. 2.47 and 2.48), + * approximate the factor (1.0 - 5.0 * cos²(inc))^-1 (causing the singularity) + * by a function, named T2 in the thesis. + *

+ * @param cosInc cosine of the mean inclination + * @return an approximation of (1.0 - 5.0 * cos²(inc))^-1 term + */ + private T T2(final T cosInc) { + + // X = (1.0 - 5.0 * cos²(inc)) + final T x = cosInc.multiply(cosInc).multiply(-5.0).add(1.0); + final T x2 = x.multiply(x); + + // Eq. 2.48 + T sum = x.getField().getZero(); + for (int i = 0; i <= 12; i++) { + final double sign = i % 2 == 0 ? +1.0 : -1.0; + sum = sum.add(FastMath.pow(x2, i).multiply(FastMath.pow(BETA, i)).multiply(sign).divide(CombinatoricsUtils.factorialDouble(i + 1))); + } + + // Right term of equation 2.47 + T product = x.getField().getOne(); + for (int i = 0; i <= 10; i++) { + product = product.multiply(FastMath.exp(x2.multiply(BETA).multiply(FastMath.scalb(-1.0, i))).add(1.0)); + } + + // Return (Eq. 2.47) + return x.multiply(BETA).multiply(sum).multiply(product); + + } + /** Extrapolate an orbit up to a specific target date. * @param date target date for the orbit * @return propagated parameters diff --git a/src/site/markdown/index.md b/src/site/markdown/index.md index 097ce86061a1549c6b718718a27f5b5e106fe781..e8e90a4819444c961e83ec52900eb70debad779a 100644 --- a/src/site/markdown/index.md +++ b/src/site/markdown/index.md @@ -77,7 +77,7 @@ * analytical propagation models * Kepler * Eckstein-Heschler - * Brouwer-Lyddane + * Brouwer-Lyddane with Warren Phipps' correction for the critical inclination of 63.4° * 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 aca1e4d352c39c27676662bca604a304276f89ed..9db7aee303a055b94d7fdc5fb4423c6c192789e8 100644 --- a/src/test/java/org/orekit/propagation/analytical/BrouwerLyddanePropagatorTest.java +++ b/src/test/java/org/orekit/propagation/analytical/BrouwerLyddanePropagatorTest.java @@ -469,24 +469,43 @@ public class BrouwerLyddanePropagatorTest { } - @Test(expected = OrekitException.class) + @Test public void criticalInclination() { - // for an eccentricity too big for the model + + final Frame inertialFrame = FramesFactory.getEME2000(); + + final double a = 24396159; // semi major axis in meters + final double e = 0.01; // eccentricity + final double i = FastMath.acos(1.0 / FastMath.sqrt(5.0)); // critical 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 + AbsoluteDate initDate = AbsoluteDate.J2000_EPOCH; - Orbit initialOrbit = new KeplerianOrbit(67679244.0, 0.01, FastMath.toRadians(63.44), 2.1, 2.9, - 6.2, PositionAngle.TRUE, FramesFactory.getEME2000(), - initDate, provider.getMu()); + final Orbit initialOrbit = new KeplerianOrbit(a, e, i, omega, raan, lM, PositionAngle.TRUE, + inertialFrame, initDate, provider.getMu()); + // Extrapolator definition // ----------------------- BrouwerLyddanePropagator extrapolator = - new BrouwerLyddanePropagator(initialOrbit, DEFAULT_LAW, provider.getAe(), provider.getMu(), - -1.08263e-3, 2.54e-6, 1.62e-6, 2.3e-7); + new BrouwerLyddanePropagator(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider)); // Extrapolation at the initial date // --------------------------------- - double delta_t = 0.0; - AbsoluteDate extrapDate = initDate.shiftedBy(delta_t); - extrapolator.propagate(extrapDate); + SpacecraftState finalOrbit = extrapolator.propagate(initDate); + + // Verify + Assert.assertEquals(0.0, + Vector3D.distance(initialOrbit.getPVCoordinates().getPosition(), + finalOrbit.getPVCoordinates().getPosition()), + 1.9e-8); + + Assert.assertEquals(0.0, + Vector3D.distance(initialOrbit.getPVCoordinates().getVelocity(), + finalOrbit.getPVCoordinates().getVelocity()), + 3.0e-12); + Assert.assertEquals(0.0, finalOrbit.getA() - initialOrbit.getA(), 0.0); + } diff --git a/src/test/java/org/orekit/propagation/analytical/FieldBrouwerLyddanePropagatorTest.java b/src/test/java/org/orekit/propagation/analytical/FieldBrouwerLyddanePropagatorTest.java index e20349004eb036283177ba39a91650092fd5c9b7..0470ecc83179795300cc35d92493033cdd42f5b5 100644 --- a/src/test/java/org/orekit/propagation/analytical/FieldBrouwerLyddanePropagatorTest.java +++ b/src/test/java/org/orekit/propagation/analytical/FieldBrouwerLyddanePropagatorTest.java @@ -549,28 +549,45 @@ public class FieldBrouwerLyddanePropagatorTest { doCriticalInclination(Decimal64Field.getInstance()); } private > void doCriticalInclination(Field field) { - // for an eccentricity too big for the model - T zero = field.getZero(); - FieldAbsoluteDate date = new FieldAbsoluteDate<>(field); - FieldAbsoluteDate initDate = date.shiftedBy(584.); - FieldOrbit initialOrbit = new FieldKeplerianOrbit<>(zero.add(67679244.0), zero.add(0.01), zero.add(FastMath.toRadians(63.44)), - zero.add(2.1), zero.add(2.9), zero.add(6.2), PositionAngle.TRUE, - FramesFactory.getEME2000(), initDate, zero.add(provider.getMu())); - try { + + final Frame inertialFrame = FramesFactory.getEME2000(); + + // Initial orbit + final double a = 24396159; // 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 + + T zero = field.getZero(); + FieldAbsoluteDate initDate = new FieldAbsoluteDate<>(field); + 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())); + // Extrapolator definition // ----------------------- 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); + new FieldBrouwerLyddanePropagator<>(initialOrbit, GravityFieldFactory.getUnnormalizedProvider(provider)); // Extrapolation at the initial date // --------------------------------- - double delta_t = 0.0; - FieldAbsoluteDate extrapDate = initDate.shiftedBy(delta_t); - extrapolator.propagate(extrapDate); - } catch (OrekitException oe) { - Assert.assertEquals(OrekitMessages.ALMOST_CRITICALLY_INCLINED_ORBIT, oe.getSpecifier()); - } + final FieldSpacecraftState finalOrbit = extrapolator.propagate(initDate); + + // Verify + Assert.assertEquals(0.0, + FieldVector3D.distance(initialOrbit.getPVCoordinates().getPosition(), + finalOrbit.getPVCoordinates().getPosition()).getReal(), + 7.0e-8); + + Assert.assertEquals(0.0, + FieldVector3D.distance(initialOrbit.getPVCoordinates().getVelocity(), + finalOrbit.getPVCoordinates().getVelocity()).getReal(), + 1.2e-11); + + Assert.assertEquals(0.0, finalOrbit.getA().getReal() - initialOrbit.getA().getReal(), 0.0); + } @Test(expected = OrekitException.class)