Fixed Newcomb operators generation for high degree gravity fields.

Fixes issue #237

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

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

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

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);