WIP - RESET THIS

move planet around ellipse

adjust orbit speed

scenes/Main.tscn

@@ -1,7 +1,7 @@
 [gd_scene load_steps=3 format=2]
 
 [ext_resource path="res://scenes/OrbitCamera.tscn" type="PackedScene" id=3]
-[ext_resource path="res://scenes/OrbitSystem.tscn" type="PackedScene" id=4]
+[ext_resource path="res://scenes/orbits/OrbitSystem.tscn" type="PackedScene" id=4]
 
 [node name="Main" type="Node2D"]
 
@@ -11,3 +11,8 @@ transform = Transform( 0.515898, 0.606099, -0.605386, -0.393123, 0.795389, 0.461
 [node name="OrbitCamera" parent="." instance=ExtResource( 3 )]
 
 [node name="Orbit System" parent="." instance=ExtResource( 4 )]
+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/OrbitSystem.tscn

@@ -1,15 +0,0 @@
-[gd_scene load_steps=2 format=2]
-
-[ext_resource path="res://scripts/orbits/OrbitSystem.cs" type="Script" id=1]
-
-[node name="Orbit System" type="Spatial"]
-script = ExtResource( 1 )
-_a = NodePath("A")
-_b = NodePath("B")
-_barycenter = NodePath("Barycenter")
-
-[node name="A" type="Spatial" parent="."]
-
-[node name="B" type="Spatial" parent="."]
-
-[node name="Barycenter" type="Spatial" parent="."]

scenes/orbits/OrbitSystem.tscn

@@ -0,0 +1,22 @@
+[gd_scene load_steps=3 format=2]
+
+[ext_resource path="res://scripts/orbits/OrbitSystem.cs" type="Script" id=1]
+[ext_resource path="res://scenes/orbits/Planet.tscn" type="PackedScene" id=2]
+
+[node name="Orbit System" type="Spatial"]
+script = ExtResource( 1 )
+SemiMajorAxis = 6.881
+Eccentricity = 0.399
+Period = 10.0
+_a = NodePath("Planet A")
+_b = NodePath("Planet B")
+_barycenter = NodePath("Barycenter")
+
+[node name="Barycenter" type="Spatial" parent="."]
+
+[node name="Planet A" parent="." instance=ExtResource( 2 )]
+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, -8.44634, 0, -3.53339 )
+Mass = 0.125

scenes/orbits/Planet.tscn

@@ -0,0 +1,11 @@
+[gd_scene load_steps=2 format=2]
+
+[ext_resource path="res://scripts/orbits/Planet.cs" type="Script" id=1]
+
+[node name="Planet" type="Spatial"]
+script = ExtResource( 1 )
+Mass = 1.0
+
+[node name="CSGSphere" type="CSGSphere" parent="."]
+radial_segments = 20
+rings = 20

scripts/orbits/ILocation.cs

@@ -1,7 +1,6 @@
-using Godot;
+using Vim.Math3d;
-using System;
 
 public interface ILocation
 {
-    Vector3 Position { get; }
+    DVector3 Position { get; set; }
 }

scripts/orbits/IMassiveExtensions.cs

@@ -0,0 +1,4 @@
+public static class IMassiveExtensions
+{
+    public static double U(this IMassive m) => Kepler.U(m.Mass);
+}

scripts/orbits/IOrbit.cs

@@ -0,0 +1,7 @@
+using Vim.Math3d;
+
+public interface IOrbit
+{
+    Ellipse Ellipse { get; }
+    DVector3 GetPosition(float t);
+}

scripts/orbits/IOrbitalElementsExtensions.cs

@@ -0,0 +1,13 @@
+using Vim.Math3d;
+using static System.Math;
+
+public static class IOrbitalElementsExtensions
+{
+    /// <summary>
+    /// Get the Cartesian position on the orbit for a given mean anomaly.
+    /// </summary>
+    /// <param name="M">Mean anomaly</param>
+    /// <returns>Position</returns>
+    public static DVector2 GetEccentricPosition2D(this OrbitalElements elements, double M) =>
+        Kepler.GetEccentricPosition2D(elements.e, elements.a, M);
+}

scripts/orbits/IPointMass.cs

@@ -0,0 +1,3 @@
+public interface IPointMass : IMassive, ILocation
+{
+}

scripts/orbits/Orbit.cs

@@ -0,0 +1,57 @@
+using Godot;
+using System;
+using Vim.Math3d;
+
+public class Orbit : IOrbit
+{
+    // 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 elliptical.
+    private DVector2 PrimaryPosition => Ellipse.Focus0;
+
+    public Ellipse Ellipse { get; set; }
+
+    public DVector3 GetPosition(float 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);
+    }
+
+    public void Draw(ImmediateGeometry geo)
+    {
+        geo.Clear();
+        geo.Begin(Mesh.PrimitiveType.LineLoop);
+        DrawEllipse(geo);
+        geo.End();
+    }
+
+    private void DrawEllipse(ImmediateGeometry geo)
+    {
+        geo.SetColor(new Color(1, 0, 0));
+
+        int steps = 100;
+
+        for (int i = 0; i < steps; i++)
+        {
+            float t = i / (float)steps;
+            float a = t * Mathf.Tau;
+
+            var v2 = Ellipse.Focus0 + Ellipse.GetPosition(a);
+
+            geo.AddVertex(new Godot.Vector3
+            {
+                x = (float)v2.X,
+                z = (float)v2.Y
+            });
+        }
+    }
+}

scripts/orbits/OrbitSystem.cs

@@ -1,79 +1,73 @@
 using Godot;
 using System;
+using Vim.Math3d;
 
+[Tool]
 public class OrbitSystem : Node, IMassive, ILocation
 {
     [Export] private NodePath _a;
     [Export] private NodePath _b;
     [Export] private NodePath _barycenter;
 
-    private readonly PointMass[] _pointMasses = new PointMass[2];
+    private double _semiMajorAxis;
-
-    public Vector3 Position => Barycenter;
-
     [Export]
-    public Vector3 Barycenter
+    private double SemiMajorAxis
+    {
+        get => _semiMajorAxis;
+        set
+        {
+            _semiMajorAxis = value;
+
+            if (!Engine.EditorHint) return;
+
+            InvalidateGeometry();
+        }
+    }
+
+
+    private double _eccentricity;
+    [Export(PropertyHint.Range, "0,1")]
+    private double Eccentricity
+    {
+        get => _eccentricity;
+        set
+        {
+            _eccentricity = value;
+
+            if (!Engine.EditorHint) return;
+
+            InvalidateGeometry();
+        }
+    }
+
+    #region Point Masses
+    private IPointMass Primary
     {
         get
         {
-            var p0 = _pointMasses[0].Position;
+            if (_pointMasses[0] == null)
-            var p1 = _pointMasses[1].Position;
-
-            return p0.LinearInterpolate(p1, .5f);
-        }
-        set => _ = value;
-    }
-
-    public float Mass { get; } = 1;
-
-    private ImmediateGeometry _orbit;
-
-    public override void _Ready()
             {
                 InitPointMasses();
-        InitGeometry();
             }
-
+            return _pointMasses[0];
-    public override void _Process(float delta)
-    {
-        DrawOrbit();
         }
-
-    private void DrawOrbit()
-    {
-        int steps = 100;
-
-        _orbit.Clear();
-        _orbit.Begin(Mesh.PrimitiveType.LineLoop);
-        _orbit.SetColor(new Color(1, 0, 0));
-        for (int i = 0; i < steps; i++)
-        {
-            float t = i / (float)steps;
-            float a = t * Mathf.Tau;
-
-            _orbit.AddVertex(new Vector3
-            {
-                x = Mathf.Sin(a),
-                z = Mathf.Cos(a)
-            });
     }
-        _orbit.End();
+    private IPointMass Secondary
-    }
-
-    private void InitGeometry()
     {
-        _orbit = new ImmediateGeometry();
+        get
-        var m = new SpatialMaterial();
+        {
-        m.VertexColorUseAsAlbedo = true;
+            if (_pointMasses[1] == null)
-        m.FlagsUnshaded = true;
+            {
-        _orbit.MaterialOverride = m;
+                InitPointMasses();
-        AddChild(_orbit);
             }
-
+            return _pointMasses[1];
+        }
+    }
+    private readonly IPointMass[] _pointMasses = new IPointMass[2];
     private void InitPointMasses()
     {
-        var a = GetPointMass(_a);
+        var a = GetNode<IPointMass>(_a);
-        var b = GetPointMass(_b);
+        var b = GetNode<IPointMass>(_b);
 
         if (a.Mass > b.Mass)
         {
@@ -86,12 +80,74 @@ public class OrbitSystem : Node, IMassive, ILocation
             _pointMasses[1] = a;
         }
     }
+    #endregion
 
-    private PointMass GetPointMass(NodePath path)
+    public float Mass => Primary.Mass + Secondary.Mass;
+    public DVector3 Position
     {
-        var spatial = GetNode<Spatial>(path);
+        get => Barycenter;
-        //var massive = node.GetChild<IMassive>(0);
+        set => Barycenter = value;
+    }
 
-        return new PointMass(spatial, 1f);
+    public DVector3 Barycenter
+    {
+        get
+        {
+            var p0 = Primary.Position;
+            var p1 = Secondary.Position;
+            return p0.Lerp(p1, .5f);
+        }
+        // TODO - make setting the berycenter do something sensible?
+        set => _ = value;
+    }
+
+    private Orbit _orbit = null;
+    private Orbit Orbit
+    {
+        get
+        {
+            if (_orbit == null)
+            {
+                _orbit = new Orbit();
+            }
+            return _orbit;
+        }
+    }
+
+    private ImmediateGeometry _orbitGeometry = null;
+    public ImmediateGeometry OrbitGeometry
+    {
+        get
+        {
+            if (_orbitGeometry == null)
+            {
+                _orbitGeometry = new ImmediateGeometry();
+                var m = new SpatialMaterial();
+                m.VertexColorUseAsAlbedo = true;
+                m.FlagsUnshaded = true;
+                _orbitGeometry.MaterialOverride = m;
+                AddChild(_orbitGeometry);
+            }
+            return _orbitGeometry;
+        }
+    }
+
+    private float _time = 0;
+
+    [Export]
+    private float _speed = 3f;
+
+    public override void _Process(float delta)
+    {
+        _time += delta * _speed;
+        InvalidateGeometry();
+    }
+
+    private void InvalidateGeometry()
+    {
+        Orbit.Ellipse = new Ellipse(SemiMajorAxis, Eccentricity);
+        Orbit.Draw(OrbitGeometry);
+
+        Secondary.Position = Orbit.GetPosition(_time);
     }
 }

scripts/orbits/OrbitalElements.cs

@@ -0,0 +1,46 @@
+public struct OrbitalElements
+{
+    private Ellipse _ellipse;
+    /// <summary>
+    /// Semi-major axis in AU
+    /// </summary>
+    public double a => _ellipse.a;
+    /// <summary>
+    /// Eccentricity
+    /// </summary>
+    public double e => _ellipse.e;
+    /// <summary>
+    /// Inclination in degrees
+    /// </summary>
+    public double i { get; private set; }
+    /// <summary>
+    /// Argument of ascending node in degrees
+    /// </summary>
+    public double W { get; private set; }
+    /// <summary>
+    /// Argument of periapsis in degrees
+    /// </summary>
+    public double p { get; private set; }
+    /// <summary>
+    /// Mean anomaly at epoch
+    /// </summary>
+    public double M0 { get; private set; }
+
+    /// <summary>
+    ///
+    /// </summary>
+    /// <param name="a">Semi-major axis in AU</param>
+    /// <param name="e">eccentricity</param>
+    /// <param name="i">Inclination in degrees</param>
+    /// <param name="W">Argument of ascending node in degrees</param>
+    /// <param name="p">Argument of periapsis in degrees</param>
+    /// <param name="M0">Mean anomaly at epoch</param>
+    public OrbitalElements(double a, double e, double i, double W, double p, double M0)
+    {
+        _ellipse = new Ellipse(a, e);
+        this.i = i;
+        this.W = W;
+        this.p = p;
+        this.M0 = M0;
+    }
+}

scripts/orbits/Planet.cs

@@ -1,18 +1,32 @@
 using Godot;
 using System;
+using Vim.Math3d;
 
-public class Planet : Node
+[Tool]
+public class Planet : Spatial, IPointMass
 {
     [Export]
     public float Mass { get; set; }
 
-    public override void _Ready()
+    private DVector3? _position;
+    public DVector3 Position
     {
+        get
+        {
+            if (_position.HasValue) return _position.Value;
 
+            var t = Translation;
+            return new DVector3(t.x, t.y, t.z);
+        }
+        set
+        {
+            _position = value;
+            Translation = new Godot.Vector3
+            {
+                x = (float)value.X,
+                y = (float)value.Y,
+                z = (float)value.Z
+            };
+        }
     }
-
-    //  public override void _Process(float delta)
-    //  {
-    //
-    //  }
 }

scripts/orbits/PointMass.cs

@@ -1,25 +0,0 @@
-using Godot;
-
-struct PointMass : IMassive, ILocation
-{
-    public float Mass => _massive == null ? _mass : _massive.Mass;
-    public Vector3 Position => _spatial.Translation;
-
-    private IMassive _massive;
-    private Spatial _spatial;
-    private float _mass;
-
-    public PointMass(Spatial spatial, IMassive massive)
-    {
-        _spatial = spatial;
-        _massive = massive;
-        _mass = 0;
-    }
-
-    public PointMass(Spatial spatial, float mass)
-    {
-        _spatial = spatial;
-        _massive = null;
-        _mass = mass;
-    }
-}

scripts/orbits/math/Ellipse.cs

@@ -0,0 +1,63 @@
+using Godot;
+using System;
+using Vim.Math3d;
+
+public struct Ellipse
+{
+    public Ellipse(double a = 1, double e = 0)
+    {
+        SemiMajorAxis = a;
+        Eccentricity = e;
+
+        // TODO: this is an immutable struct, so initialise everything else 
+        // in the constructor to avoid recalculating properties whenever
+        // they are accessed
+    }
+
+    /// <summary>
+    /// Semi-major axis
+    /// </summary>
+    public double a => SemiMajorAxis;
+    public double SemiMajorAxis { get; }
+
+    /// <summary>
+    /// Eccentricity
+    /// </summary>
+    public double e => Eccentricity;
+    public double Eccentricity { get; }
+
+    /// <summary>
+    /// Semi-minor axis
+    /// </summary>
+    public double b => SemiMinorAxis;
+    public double SemiMinorAxis => Math.Sqrt(Apisides.max * Apisides.min);
+
+    /// <summary>
+    /// Get a position on the auxiliary circle
+    /// </summary>
+    /// <param name="t">Angle in radians around the circle</param>
+    /// <returns>2D position on the circle scaled by the semi-major axis</returns>
+    public DVector2 GetAuxiliaryPosition2D(double t = 0) =>
+        new DVector2(Math.Cos(t), Math.Sin(t)) * a;
+
+    private Apisides Apisides => new Apisides(a, e);
+
+    public DVector2 GetPosition(double t) =>
+        new DVector2(Math.Cos(t) * a, Math.Sin(t) * b);
+
+    public DVector2 Focus0 => new DVector2(-GetFocusDistance(), 0);
+    public DVector2 Focus1 => new DVector2(GetFocusDistance(), 0);
+    private double GetFocusDistance() => Math.Sqrt(a * a - b * b);
+}
+
+public struct Apisides
+{
+    public readonly double min;
+    public readonly double max;
+
+    public Apisides(double a, double e)
+    {
+        min = a * (1 - e);
+        max = a * (1 + e);
+    }
+}

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