Fixed issue #996 : HolmesFeatherstoneAttractionModel error with a degree 0 gravity field

src/changes/changes.xml

@@ -24,6 +24,9 @@
       <action dev="luc" type="update" issue="1017">
         Removed generics in EventHandler.
       </action>
+      <action dev="pascal" type="fix" issue="996">
+        Fixed HolmesFeatherstoneAttractionModel error with a degree 0 gravity field.
+      </action>
       <action dev="luc" type="update" issue="1013">
         Use SI units (i.e. seconds) in GNSSDate.
       </action>

src/main/java/org/orekit/forces/gravity/HolmesFeatherstoneAttractionModel.java

@@ -21,8 +21,8 @@ import java.util.Collections;
 import java.util.List;
 import java.util.stream.Stream;
 
-import org.hipparchus.Field;
 import org.hipparchus.CalculusFieldElement;
+import org.hipparchus.Field;
 import org.hipparchus.analysis.differentiation.DerivativeStructure;
 import org.hipparchus.analysis.differentiation.Gradient;
 import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
@@ -153,7 +153,9 @@ public class HolmesFeatherstoneAttractionModel extends AbstractForceModel implem
         // scaled sectorial terms corresponding to equation 28 in Holmes and Featherstone paper
         sectorial    = new double[degree + 1];
         sectorial[0] = FastMath.scalb(1.0, -SCALING);
-        sectorial[1] = FastMath.sqrt(3) * sectorial[0];
+        if (degree > 0) {
+            sectorial[1] = FastMath.sqrt(3) * sectorial[0];
+        }
         for (int m = 2; m < sectorial.length; ++m) {
             sectorial[m] = FastMath.sqrt((2 * m + 1) / (2.0 * m)) * sectorial[m - 1];
         }
@@ -720,7 +722,9 @@ public class HolmesFeatherstoneAttractionModel extends AbstractForceModel implem
         // initialize array
         final double[] aOrN = new double[provider.getMaxDegree() + 1];
         aOrN[0] = 1;
-        aOrN[1] = aOr;
+        if (provider.getMaxDegree() > 0) {
+            aOrN[1] = aOr;
+        }
 
         // fill up array
         for (int n = 2; n < aOrN.length; ++n) {
@@ -742,7 +746,9 @@ public class HolmesFeatherstoneAttractionModel extends AbstractForceModel implem
         // initialize array
         final T[] aOrN = MathArrays.buildArray(aOr.getField(), provider.getMaxDegree() + 1);
         aOrN[0] = aOr.getField().getOne();
-        aOrN[1] = aOr;
+        if (provider.getMaxDegree() > 0) {
+            aOrN[1] = aOr;
+        }
 
         // fill up array
         for (int n = 2; n < aOrN.length; ++n) {

src/test/java/org/orekit/forces/gravity/HolmesFeatherstoneAttractionModelTest.java

@@ -16,6 +16,9 @@
  */
 package org.orekit.forces.gravity;
 
+import java.io.Serializable;
+import java.lang.reflect.InvocationTargetException;
+
 import org.hipparchus.Field;
 import org.hipparchus.analysis.differentiation.DSFactory;
 import org.hipparchus.analysis.differentiation.DerivativeStructure;
@@ -80,9 +83,6 @@ import org.orekit.utils.IERSConventions;
 import org.orekit.utils.PVCoordinates;
 import org.orekit.utils.PVCoordinatesProvider;
 
-import java.io.Serializable;
-import java.lang.reflect.InvocationTargetException;
-
 
 public class HolmesFeatherstoneAttractionModelTest extends AbstractLegacyForceModelTest {
 
@@ -198,7 +198,7 @@ public class HolmesFeatherstoneAttractionModelTest extends AbstractLegacyForceMo
     }
 
     private GradientHessian gradientHessian(final HolmesFeatherstoneAttractionModel hfModel,
-                                           final AbsoluteDate date, final Vector3D position)
+                                            final AbsoluteDate date, final Vector3D position)
         {
         try {
 
@@ -1187,6 +1187,48 @@ public class HolmesFeatherstoneAttractionModelTest extends AbstractLegacyForceMo
 
     }
 
+    @Test
+    public void testIssue996() {
+
+        Utils.setDataRoot("regular-data:potential/grgs-format");
+        GravityFieldFactory.addPotentialCoefficientsReader(new GRGSFormatReader("grim4s4_gr", true));
+
+        // initialization
+        AbsoluteDate date = new AbsoluteDate(new DateComponents(2000, 07, 01),
+                                             new TimeComponents(13, 59, 27.816),
+                                             TimeScalesFactory.getUTC());
+        double i     = FastMath.toRadians(98.7);
+        double omega = FastMath.toRadians(93.0);
+        double OMEGA = FastMath.toRadians(15.0 * 22.5);
+        Orbit orbit  = new KeplerianOrbit(7201009.7124401, 1e-3, i , omega, OMEGA,
+                                          0, PositionAngle.MEAN, FramesFactory.getEME2000(), date, mu);
+        Vector3D pos = orbit.getPosition(itrf);
+
+        double[] zeroGradient = new double[] {-0., 0., 0.};
+
+        NormalizedSphericalHarmonicsProvider provider00 = GravityFieldFactory.getNormalizedProvider(0,  0);
+        HolmesFeatherstoneAttractionModel hfModel00 = new HolmesFeatherstoneAttractionModel(itrf, provider00);
+
+        Assertions.assertEquals(0., hfModel00.nonCentralPart(date, pos, mu));
+        Assertions.assertEquals(mu / pos.getNorm(), hfModel00.value(date, pos, mu));
+        Assertions.assertArrayEquals(zeroGradient, hfModel00.gradient(date, pos, mu));
+
+        NormalizedSphericalHarmonicsProvider provider10 = GravityFieldFactory.getNormalizedProvider(1,  0);
+        HolmesFeatherstoneAttractionModel hfModel10 = new HolmesFeatherstoneAttractionModel(itrf, provider10);
+
+        Assertions.assertEquals(0., hfModel10.nonCentralPart(date, pos, mu));
+        Assertions.assertEquals(mu / pos.getNorm(), hfModel10.value(date, pos, mu));
+        Assertions.assertArrayEquals(zeroGradient, hfModel10.gradient(date, pos, mu));
+
+        NormalizedSphericalHarmonicsProvider provider11 = GravityFieldFactory.getNormalizedProvider(1,  1);
+        HolmesFeatherstoneAttractionModel hfModel11 = new HolmesFeatherstoneAttractionModel(itrf, provider11);
+
+        Assertions.assertEquals(0., hfModel11.nonCentralPart(date, pos, mu));
+        Assertions.assertEquals(mu / pos.getNorm(), hfModel11.value(date, pos, mu));
+        Assertions.assertArrayEquals(zeroGradient, hfModel11.gradient(date, pos, mu));
+
+    }
+
     @BeforeEach
     public void setUp() {
         itrf   = null;