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

Allowed Brouwer-Lyddane model to work for the critical inclination.

parent 30a77101
Pipeline #1577 canceled with stages
in 3 minutes and 2 seconds
Hide whitespace changes
Inline Side-by-side
src/changes/changes.xml
View file @ ab4b46de
......@@ -21,6 +21,9 @@
</properties>
<body>
<release version="11.1" date="TBD" description="TBD">
<action dev="bryan" type="add" issue="869">
Allowed Brouwer-Lyddane model to work for the critical inclination.
</action>
<action dev="bryan" type="fix" issue="868">
Fixed writing of whitespace characters in CPF writer.
</action>
......
src/main/java/org/orekit/propagation/analytical/BrouwerLyddanePropagator.java
View file @ ab4b46de
......@@ -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;
* <p>
* 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.
* </p>
* @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°).
* <p>
* 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.
* </p>
* @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
......
src/main/java/org/orekit/propagation/analytical/FieldBrouwerLyddanePropagator.java
View file @ ab4b46de
......@@ -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).
* <p>
* 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.
* </p>
* @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<T extends CalculusFieldElement<T>> extends FieldAbstractAnalyticalPropagator<T> {
/** Beta constant used by T2 function. */
private static final double BETA = FastMath.scalb(100, -11);
/** Initial Brouwer-Lyddane model. */
private FieldBLModel<T> initialModel;
......@@ -663,7 +668,7 @@ public class FieldBrouwerLyddanePropagator<T extends CalculusFieldElement<T>> 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<T extends CalculusFieldElement<T>> 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<T extends CalculusFieldElement<T>> ex
return ek;
}
/**
* This method is used in Brouwer-Lyddane model to avoid singularity at the
* critical inclination (i = 63.4°).
* <p>
* 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.
* </p>
* @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
......
src/site/markdown/index.md
View file @ ab4b46de
......@@ -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
......
src/test/java/org/orekit/propagation/analytical/BrouwerLyddanePropagatorTest.java
View file @ ab4b46de
......@@ -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);
}
......
src/test/java/org/orekit/propagation/analytical/FieldBrouwerLyddanePropagatorTest.java
View file @ ab4b46de
......@@ -549,28 +549,45 @@ public class FieldBrouwerLyddanePropagatorTest {
doCriticalInclination(Decimal64Field.getInstance());
}
private <T extends CalculusFieldElement<T>> void doCriticalInclination(Field<T> field) {
// for an eccentricity too big for the model
T zero = field.getZero();
FieldAbsoluteDate<T> date = new FieldAbsoluteDate<>(field);
FieldAbsoluteDate<T> initDate = date.shiftedBy(584.);
FieldOrbit<T> 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<T> initDate = new FieldAbsoluteDate<>(field);
final FieldOrbit<T> 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<T> 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<T> extrapDate = initDate.shiftedBy(delta_t);
extrapolator.propagate(extrapDate);
} catch (OrekitException oe) {
Assert.assertEquals(OrekitMessages.ALMOST_CRITICALLY_INCLINED_ORBIT, oe.getSpecifier());
}
final FieldSpacecraftState<T> 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)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment