Modifications on getMeanElementRate implementation

src/main/java/org/orekit/propagation/semianalytical/dsst/DSSTPartialDerivativeEquations.java

@@ -3,10 +3,7 @@ package org.orekit.propagation.semianalytical.dsst;
 import java.util.Map;
 
 import org.hipparchus.RealFieldElement;
-import org.hipparchus.analysis.differentiation.DerivativeStructure;
-import org.hipparchus.geometry.euclidean.threed.FieldVector3D;
 import org.orekit.errors.OrekitException;
-import org.orekit.forces.ForceModel;
 import org.orekit.propagation.FieldSpacecraftState;
 import org.orekit.utils.ParameterDriver;
 import org.orekit.utils.ParameterDriversList;

src/main/java/org/orekit/propagation/semianalytical/dsst/forces/DSSTThirdBody.java

@@ -18,6 +18,7 @@ package org.orekit.propagation.semianalytical.dsst.forces;
 
 import java.io.NotSerializableException;
 import java.io.Serializable;
+import java.lang.reflect.Field;
 import java.util.ArrayList;
 import java.util.HashMap;
 import java.util.List;
@@ -32,6 +33,7 @@ import org.hipparchus.analysis.differentiation.DerivativeStructure;
 import org.hipparchus.exception.NullArgumentException;
 import org.hipparchus.geometry.euclidean.threed.Vector3D;
 import org.hipparchus.util.FastMath;
+import org.hipparchus.util.MathArrays;
 import org.orekit.attitudes.AttitudeProvider;
 import org.orekit.bodies.CelestialBody;
 import org.orekit.errors.OrekitException;
@@ -438,7 +440,8 @@ public class DSSTThirdBody<T> implements DSSTForceModel {
     }
     
     /** {@inheritDoc} */
-    public <T extends RealFieldElement<T>> T[] getMeanElementRate(final FieldSpacecraftState<T> currentState) {
+    @SuppressWarnings("hiding")
+	public <T extends RealFieldElement<T>> T[] getMeanElementRate(final FieldSpacecraftState<T> currentState) {
 
         // Qns coefficients
         Qns = CoefficientsFactory.computeQns(gamma, maxAR3Pow, maxEccPow);
@@ -450,7 +453,7 @@ public class DSSTThirdBody<T> implements DSSTForceModel {
         }
 
         // Compute potential U derivatives
-        final <T extends RealFieldElement<T>> T[] dU  = computeUDerivatives();
+        final T[] dU  = MathArrays.buildArray(computeUDerivatives(), computeUDerivatives().length);
         final T dUda  = dU[0];
         final T dUdk  = dU[1];
         final T dUdh  = dU[2];
@@ -476,8 +479,7 @@ public class DSSTThirdBody<T> implements DSSTForceModel {
 
         return new T[] {da, dk, dh, dq, dp, dM};
 
-    }
-    
+    }    
 
     /** {@inheritDoc} */
     @Override
@@ -670,6 +672,97 @@ public class DSSTThirdBody<T> implements DSSTForceModel {
         };
 
     }
+    
+    /**
+    private <T extends RealFieldElement<T>> T[] computeUDerivativesDSST() {
+
+        // Gs and Hs coefficients
+        final <T extends RealFieldElement<T>> T[][] GsHs = CoefficientsFactory.computeGsHs(k, h, alpha, beta, maxEccPow);
+
+        // Initialise U.
+        U = 0.;
+
+        // Potential derivatives
+        T dUda  = 0.;
+        T dUdk  = 0.;
+        T dUdh  = 0.;
+        T dUdAl = 0.;
+        T dUdBe = 0.;
+        T dUdGa = 0.;
+
+        for (int s = 0; s <= maxEccPow; s++) {
+            // initialise the Hansen roots
+            this.hansenObjects[s].computeInitValues(B, BB, BBB);
+
+            // Get the current Gs coefficient
+            final T gs = GsHs[0][s];
+
+            // Compute Gs partial derivatives from 3.1-(9)
+            T dGsdh  = 0.;
+            T dGsdk  = 0.;
+            T dGsdAl = 0.;
+            T dGsdBe = 0.;
+            if (s > 0) {
+                // First get the G(s-1) and the H(s-1) coefficients
+                final T sxGsm1 = s * GsHs[0][s - 1];
+                final T sxHsm1 = s * GsHs[1][s - 1];
+                // Then compute derivatives
+                dGsdh  = beta.multiply(sxGsm1).subtract(alpha.multiply(sxHsm1));
+                dGsdk  = alpha.multiply(sxGsm1).add(beta.multiply(sxHsm1));
+                dGsdAl = k.multiply(sxGsm1).subtract(h.multiply(sxHsm1));
+                dGsdBe = h.multiply(sxGsm1).add(k.multiply(sxHsm1));
+            }
+
+            // Kronecker symbol (2 - delta(0,s))
+            final T delta0s = (s == 0) ? 1. : 2.;
+
+            for (int n = FastMath.max(2, s); n <= maxAR3Pow; n++) {
+                // (n - s) must be even
+                if ((n - s) % 2 == 0) {
+                    // Extract data from previous computation :
+                    final T kns   = this.hansenObjects[s].getValue(n, B);
+                    final T dkns  = this.hansenObjects[s].getDerivative(n, B);
+
+                    final T vns   = Vns.get(new NSKey(n, s));
+                    final T coef0 = delta0s.multiply(aoR3Pow[n]).multiply(vns);
+                    final T coef1 = coef0.multiply(Qns[n][s]);
+                    final T coef2 = coef1.multiply(kns);
+                    // dQns/dGamma = Q(n, s + 1) from Equation 3.1-(8)
+                    // for n = s, Q(n, n + 1) = 0. (Cefola & Broucke, 1975)
+                    final T dqns = (n == s) ? 0. : Qns[n][s + 1];
+
+                    //Compute U:
+                    U = U.add(coef2.multiply(gs));
+
+                    // Compute dU / da :
+                    dUda = dUda.add(coef2.multiply(n).multiply(gs));
+                    // Compute dU / dh
+                    dUdh = dUdh.add(coef1.multiply(kns.multiply(dGsdh).add(hXXX.multiply(gs.multiply(dkns)))));
+                    // Compute dU / dk
+                    dUdk = dUdk.add(coef1.multiply(kns.multiply(dGsdk).add(kXXX.multiply(gs.multiply(dkns)))));
+                    // Compute dU / dAlpha
+                    dUdAl = dUdAl.add(coef2.multiply(dGsdAl));
+                    // Compute dU / dBeta
+                    dUdBe = dUdBe.add(coef2.multiply(dGsdBe));
+                    // Compute dU / dGamma
+                    dUdGa = dUdGa.add(coef0.multiply(kns.multiply(dqns.multiply(gs))));
+                }
+            }
+        }
+
+        // multiply by mu3 / R3
+        U = U.multiply(muoR3);
+
+        return new T[] {
+            dUda.multiply(muoR3.divide(a)),
+            dUdk.multiply(muoR3),
+            dUdh.multiply(muoR3),
+            dUdAl.multiply(muoR3),
+            dUdBe.multiply(muoR3),
+            dUdGa.multiply(muoR3)
+        };
+
+    } */
    
     /** {@inheritDoc} */
     @Override