From a72c8c585121a7cab8bb9e97a74954e6a8eee633 Mon Sep 17 00:00:00 2001
From: Luc Maisonobe <luc@orekit.org>
Date: Thu, 24 Mar 2016 15:02:21 +0100
Subject: [PATCH] Fixed Newcomb operators generation for high degree gravity
fields.
Fixes issue #237
---
.../dsst/utilities/NewcombOperators.java | 65 +++++++++++++++++--
src/site/xdoc/changes.xml | 4 ++
.../dsst/DSSTPropagatorTest.java | 31 +++++++++
.../dsst/utilities/GammaMnsFunctionTest.java | 1 +
4 files changed, 97 insertions(+), 4 deletions(-)
diff --git a/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/NewcombOperators.java b/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/NewcombOperators.java
index e32146fa66..9d616486a4 100644
--- a/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/NewcombOperators.java
+++ b/src/main/java/org/orekit/propagation/semianalytical/dsst/utilities/NewcombOperators.java
@@ -23,7 +23,6 @@ import java.util.SortedMap;
import java.util.TreeMap;
import org.apache.commons.math3.analysis.polynomials.PolynomialFunction;
-import org.apache.commons.math3.analysis.polynomials.PolynomialsUtils;
import org.apache.commons.math3.util.FastMath;
/** Implementation of the Modified Newcomb Operators.
@@ -312,8 +311,7 @@ public class NewcombOperators {
}
}
- /** Shift a list of {@link PolynomialFunction}, from the
- * {@link PolynomialsUtils#shift(double[], double)} method.
+ /** Shift a list of {@link PolynomialFunction}.
*
* @param polynomialList list of {@link PolynomialFunction}
* @param shift shift value
@@ -323,11 +321,70 @@ public class NewcombOperators {
final int shift) {
final List<PolynomialFunction> shiftedList = new ArrayList<PolynomialFunction>();
for (PolynomialFunction function : polynomialList) {
- shiftedList.add(new PolynomialFunction(PolynomialsUtils.shift(function.getCoefficients(), shift)));
+ shiftedList.add(new PolynomialFunction(shift(function.getCoefficients(), shift)));
}
return shiftedList;
}
+ /**
+ * Compute the coefficients of the polynomial \(P_s(x)\)
+ * whose values at point {@code x} will be the same as the those from the
+ * original polynomial \(P(x)\) when computed at {@code x + shift}.
+ * <p>
+ * More precisely, let \(\Delta = \) {@code shift} and let
+ * \(P_s(x) = P(x + \Delta)\). The returned array
+ * consists of the coefficients of \(P_s\). So if \(a_0, ..., a_{n-1}\)
+ * are the coefficients of \(P\), then the returned array
+ * \(b_0, ..., b_{n-1}\) satisfies the identity
+ * \(\sum_{i=0}^{n-1} b_i x^i = \sum_{i=0}^{n-1} a_i (x + \Delta)^i\) for all \(x\).
+ * </p>
+ * <p>
+ * This method is a modified version of the method with the same name
+ * in Apache Commons Math {@code PolynomialsUtils} class, simply changing
+ * computation of binomial coefficients so degrees higher than 66 can be used.
+ * </p>
+ *
+ * @param coefficients Coefficients of the original polynomial.
+ * @param shift Shift value.
+ * @return the coefficients \(b_i\) of the shifted
+ * polynomial.
+ */
+ public static double[] shift(final double[] coefficients,
+ final double shift) {
+ final int dp1 = coefficients.length;
+ final double[] newCoefficients = new double[dp1];
+
+ // Pascal triangle.
+ final double[][] coeff = new double[dp1][dp1];
+ coeff[0][0] = 1;
+ for (int i = 1; i < dp1; i++) {
+ coeff[i][0] = 1;
+ for (int j = 1; j < i; j++) {
+ coeff[i][j] = coeff[i - 1][j - 1] + coeff[i - 1][j];
+ }
+ coeff[i][i] = 1;
+ }
+
+ // First polynomial coefficient.
+ double shiftI = 1;
+ for (int i = 0; i < dp1; i++) {
+ newCoefficients[0] += coefficients[i] * shiftI;
+ shiftI *= shift;
+ }
+
+ // Superior order.
+ final int d = dp1 - 1;
+ for (int i = 0; i < d; i++) {
+ double shiftJmI = 1;
+ for (int j = i; j < d; j++) {
+ newCoefficients[i + 1] += coeff[j + 1][j - i] * coefficients[j + 1] * shiftJmI;
+ shiftJmI *= shift;
+ }
+ }
+
+ return newCoefficients;
+ }
+
/** Generate recurrence coefficients.
*
* @param rho Ï value
diff --git a/src/site/xdoc/changes.xml b/src/site/xdoc/changes.xml
index d0090d1677..04bf2d16e6 100644
--- a/src/site/xdoc/changes.xml
+++ b/src/site/xdoc/changes.xml
@@ -22,6 +22,10 @@
<body>
<release version="7.2" date="TBC"
description="">
+ <action dev="luc" type="fix">
+ Fixed Newcomb operators generation in DSST for high degree gravity fields.
+ Fixes issue #237
+ </action>
<action dev="luc" type="update">
Improved tuning of DSST zonal force models. Users can now tune max degree,
max eccentricity power and max frequency in true longitude for short
diff --git a/src/test/java/org/orekit/propagation/semianalytical/dsst/DSSTPropagatorTest.java b/src/test/java/org/orekit/propagation/semianalytical/dsst/DSSTPropagatorTest.java
index f836b64b66..aedc05a03a 100644
--- a/src/test/java/org/orekit/propagation/semianalytical/dsst/DSSTPropagatorTest.java
+++ b/src/test/java/org/orekit/propagation/semianalytical/dsst/DSSTPropagatorTest.java
@@ -44,6 +44,7 @@ import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.errors.OrekitException;
import org.orekit.forces.drag.Atmosphere;
import org.orekit.forces.drag.HarrisPriester;
+import org.orekit.forces.gravity.potential.GRGSFormatReader;
import org.orekit.forces.gravity.potential.GravityFieldFactory;
import org.orekit.forces.gravity.potential.ICGEMFormatReader;
import org.orekit.forces.gravity.potential.TideSystem;
@@ -79,6 +80,7 @@ import org.orekit.time.TimeComponents;
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;
import org.orekit.utils.TimeStampedPVCoordinates;
@@ -86,6 +88,35 @@ public class DSSTPropagatorTest {
private DSSTPropagator dsstProp;
+ @Test
+ public void testHighDegreesSetting() throws OrekitException {
+
+ Utils.setDataRoot("regular-data:potential/grgs-format");
+ GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
+ int earthDegree = 36;
+ int earthOrder = 36;
+ final UnnormalizedSphericalHarmonicsProvider provider =
+ GravityFieldFactory.getUnnormalizedProvider(earthDegree,
+ earthOrder);
+ final org.orekit.frames.Frame earthFrame
+ =FramesFactory.getITRF(IERSConventions.IERS_2010,true); // terrestrial frame
+ final DSSTForceModel force = new
+ DSSTTesseral(earthFrame,Constants.WGS84_EARTH_ANGULAR_VELOCITY, provider,
+ earthDegree, earthOrder, earthDegree, earthOrder);
+ final Collection<DSSTForceModel> forces = new ArrayList<DSSTForceModel>();
+ forces.add(force);
+ TimeScale tai = TimeScalesFactory.getTAI();
+ AbsoluteDate initialDate = new AbsoluteDate("2015-07-01", tai);
+ Frame eci = FramesFactory.getGCRF();
+ KeplerianOrbit orbit = new KeplerianOrbit(
+ 7120000.0, 1.0e-3, FastMath.toRadians(60.0),
+ FastMath.toRadians(120.0), FastMath.toRadians(47.0),
+ FastMath.toRadians(12.0),
+ PositionAngle.TRUE, eci, initialDate, Constants.EIGEN5C_EARTH_MU);
+ SpacecraftState oscuState = DSSTPropagator.computeOsculatingState(new SpacecraftState(orbit), forces);
+ Assert.assertEquals(7119927.097122, oscuState.getA(), 0.001);
+ }
+
@Test
public void testEphemerisDates() throws OrekitException {
//setup
diff --git a/src/test/java/org/orekit/propagation/semianalytical/dsst/utilities/GammaMnsFunctionTest.java b/src/test/java/org/orekit/propagation/semianalytical/dsst/utilities/GammaMnsFunctionTest.java
index e90f52b0d6..9a7bdd3113 100644
--- a/src/test/java/org/orekit/propagation/semianalytical/dsst/utilities/GammaMnsFunctionTest.java
+++ b/src/test/java/org/orekit/propagation/semianalytical/dsst/utilities/GammaMnsFunctionTest.java
@@ -72,6 +72,7 @@ public class GammaMnsFunctionTest {
throws NoSuchFieldException, SecurityException, IllegalArgumentException, IllegalAccessException {
Field precomputedF = GammaMnsFunction.class.getDeclaredField("PRECOMPUTED_RATIOS");
precomputedF.setAccessible(true);
+ precomputedF.set(null, new double[0]);
new GammaMnsFunction(nMax, 0.5, +1);
double[] orginalPrecomputed = (double[]) precomputedF.get(null);
Assert.assertEquals((nMax + 1) * (nMax + 2) * (4 * nMax + 3) / 6, orginalPrecomputed.length);
--
GitLab