move according to kepler

scenes/Main.tscn

@@ -14,4 +14,5 @@ transform = Transform( 0.515898, 0.606099, -0.605386, -0.393123, 0.795389, 0.461
 transform = Transform( 0.856869, 0.3292, -0.396741, 0.0949996, 0.655565, 0.749139, 0.506706, -0.679604, 0.53046, -0.599122, 0, 0 )
 SemiMajorAxis = 6.166
 Eccentricity = 0.239
+Period = 1.0
 _speed = 0.877

scenes/orbits/OrbitSystem.tscn

@@ -5,8 +5,9 @@
 
 [node name="Orbit System" type="Spatial"]
 script = ExtResource( 1 )
-SemiMajorAxis = 5.994
-Eccentricity = 0.518
+SemiMajorAxis = 6.881
+Eccentricity = 0.399
+Period = 10.0
 _a = NodePath("Planet A")
 _b = NodePath("Planet B")
 _barycenter = NodePath("Barycenter")
@@ -17,5 +18,5 @@ _barycenter = NodePath("Barycenter")
 transform = Transform( 0.999977, 0.00603502, 0.00319089, -0.00604282, 0.999979, 0.0024369, -0.00317609, -0.00245615, 0.999992, 0, 0, 0 )
 
 [node name="Planet B" parent="." instance=ExtResource( 2 )]
-transform = Transform( 0.5, 0, 0, 0, 0.5, 0, 0, 0, 0.5, -6.57314, 0, 4.18169 )
+transform = Transform( 0.5, 0, 0, 0, 0.5, 0, 0, 0, 0.5, -8.44634, 0, -3.53339 )
 Mass = 0.125

scripts/orbits/Orbit.cs

@@ -4,16 +4,25 @@ using Vim.Math3d;
 
 public class Orbit
 {
-    // the position of the primary body relative to the ellipse - the 
+    // The position of the primary body relative to the ellipse - the 
     // orbit itself is presented as being centred on the primary, but
-    // all our math is ellipse stuff
+    // all our math is elliptical.
     private DVector2 PrimaryPosition => Ellipse.Focus0;
 
     public Ellipse Ellipse { get; set; }
 
     public DVector3 GetPosition(float t)
     {
-        var p = PrimaryPosition + Ellipse.GetPosition(t);
+        // TODO: determine period from mass of orbiting bodies
+        var P = 10f;
+
+        var a = Ellipse.a;
+        var e = Ellipse.e;
+
+        var M = Kepler.GetMeanAnomaly(t, P);
+        var E = Kepler.GetEccentricAnomaly(e, M) + Math.PI;
+
+        var p = Kepler.GetEccentricPosition2D(e, a, t);
         return new DVector3(p.X, 0, p.Y);
     }
 

scripts/orbits/math/Kepler.cs

@@ -0,0 +1,127 @@
+using static System.Math;
+using Vim.Math3d;
+
+public static class Kepler
+{
+    /// <summary>
+    /// 2 * PI
+    /// </summary>
+    private static double TAU => PI * 2.0;
+    /// <summary>
+    /// Universal gravitational constant
+    /// </summary>
+    private static double G => 6.677e-11;
+
+    /// <summary>
+    /// Get the mean anomaly of an orbit given at a given time, given
+    /// the period of the orbit.
+    /// </summary>
+    /// <param name="time">Time of mean anomaly</param>
+    /// <param name="P">Orbital period</param>
+    /// <returns></returns>
+    public static double GetMeanAnomaly(double time, double P)
+    {
+        var M = ((time % P / P) * TAU) % TAU;
+        return M % TAU;
+    }
+
+    /// <summary>
+    /// 
+    /// </summary>
+    /// <param name="P">Orbital period</param>
+    /// <returns></returns>
+    public static double GetMeanAngularVelocity(double P) => TAU / P;
+
+    /// <summary>
+    /// 
+    /// </summary>
+    /// <param name="a">Semimajor axis in AU</param>
+    /// <param name="U">Standard gravitational parameter</param>
+    /// <returns></returns>
+    public static double GetPeriod(double a, double U)
+    {
+        var a3 = a * a * a;
+        var T = TAU * Sqrt(a3 / U);
+
+        return T;
+    }
+
+    /// <summary>
+    /// Get the standard gravitational parameter of a mass
+    /// </summary>
+    /// <param name="M">Mass in kilograms</param>
+    /// <returns>Standard gravitational parameter</returns>
+    public static double U(double M) => M * G;
+    /// <summary>
+    /// Get the standard gravitational parameter of a mass
+    /// </summary>
+    /// <param name="M">Mass in kilograms</param>
+    /// <returns>Standard gravitational parameter</returns>
+    public static double GetStandardGravitationalParameter(double M) => M * G;
+
+    /// <summary>
+    /// Get the eccentric anomaly for a given mean anomaly.
+    /// </summary>
+    /// <param name="e">Eccentricity</param>
+    /// <param name="M">Mean anomaly</param>
+    /// <returns>Eccentric anomaly</returns>
+    public static double GetEccentricAnomaly(double e, double M) =>
+        GetEccentricAnomalyBS(e, M);
+
+    /// <summary>
+    /// Get the eccentric anomaly for a given mean anomaly 
+    /// using a binary search.
+    /// </summary>
+    /// <param name="e">Eccentricity</param>
+    /// <param name="M">Mean anomaly</param>
+    /// <returns>Eccentric anomaly</returns>
+    private static double GetEccentricAnomalyBS(double e, double M)
+    {
+        // from Astronomical Algorithms, p. 206
+        const int iterations = 33;
+        double E0, D;
+
+        double F = Sign(M);
+        M = Abs(M) / TAU;
+        M = (M - (int)M) * TAU * F;
+        if (M < 0)
+        {
+            M = M + TAU;
+        }
+        F = 1.0;
+        if (M > PI)
+        {
+            F = -1.0;
+            M = TAU - M;
+        }
+        E0 = PI / 2;
+        D = PI / 4;
+        for (int J = 0; J < iterations; J++)
+        {
+            double M1 = E0 - e * Sin(E0);
+            E0 = E0 + D * Sign(M - M1);
+            D *= 0.5;
+        }
+        E0 *= F;
+
+        return E0;
+    }
+
+    /// <summary>
+    /// Get the Cartesian position on an orbit given its eccentricity, 
+    /// its semi-major axis and a given mean anomaly.
+    /// </summary>
+    /// <param name="e">Eccentricity</param>
+    /// <param name="a">Semi-major axis</param>
+    /// <param name="M">Mean anomaly</param>
+    /// <returns>Position on the orbit</returns>
+    public static DVector2 GetEccentricPosition2D(double e, double a, double M)
+    {
+        var E = Kepler.GetEccentricAnomaly(e, M);
+
+        var P = a * (Cos(E) - e);
+        var Q = a * Sin(E) * Sqrt(1.0 - Pow(e, 2));
+
+        return new DVector2(P, Q);
+    }
+}