WIP - RESET THIS

WIP RESET THIS

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/orbits/OrbitSystem.tscn" type="PackedScene" id=4]
+[ext_resource path="res://scenes/orbits/Barycenter.tscn" type="PackedScene" id=4]
 
 [node name="Main" type="Node2D"]
 
@@ -10,9 +10,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 )]
+[node name="Barycenter" 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/orbits/Barycenter.tscn

@@ -3,20 +3,17 @@
 [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"]
+[node name="Barycenter" 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")
+_barycenter = NodePath("")
-
-[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 )
+transform = Transform( 0.5, 0, 0, 0, 0.5, 0, 0, 0, 0.5, -9.62645, 0, -0.0287273 )
 Mass = 0.125

scripts/orbits/IMassive.cs

@@ -3,5 +3,5 @@ using System;
 
 public interface IMassive
 {
-    float Mass { get; }
+    double Mass { get; }
 }

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,6 @@
+using Vim.Math3d;
+
+public interface IOrbit : IEllipse
+{
+    DVector2 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/Orbit.cs

@@ -2,22 +2,49 @@ using Godot;
 using System;
 using Vim.Math3d;
 
-public class Orbit
+public class Orbit : IOrbit
 {
-    // The position of the primary body relative to the ellipse - the 
+    #region IEllipse
-    // orbit itself is presented as being centred on the primary, but
+    public double a => Ellipse.a;
-    // all our math is elliptical.
+    public double b => Ellipse.b;
-    private DVector2 PrimaryPosition => Ellipse.Focus0;
+    public double e => Ellipse.e;
+    public DVector2[] Foci => Ellipse.Foci;
+    #endregion
 
-    public Ellipse Ellipse { get; set; }
+    public IEllipse Ellipse { get; set; }
+    private double U => Kepler.U(_massive.Mass);
+
+    private readonly IMassive _massive;
+    private readonly ILocation _satellite;
+
+    /// <summary>
+    /// An orbit is treated as one stationary body at one focus of the orbit's ellipse. 
+    /// The satellite is treated as a massless point.
+    /// </summary>
+    /// <param name="massive">Mass of the object being orbited</param>
+    /// <param name="satellite">The orbiting body</param>
+    public Orbit(IMassive massive, ILocation satellite)
+    {
+        // semi-major axis is distance of satellite from origin
+        var a = satellite.Position.Magnitude();
+        // TODO: determine eccentricity better
+        var e = 0.8f;
+        this.Ellipse = new Ellipse(a, e);
+
+        _massive = massive;
+        _satellite = satellite;
+    }
 
     public DVector3 GetPosition(float t)
     {
+        var a = Ellipse.a;
+        var e = Ellipse.e;
+
+        Kepler.GetPeriod(a, U);
+
         // 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;
@@ -45,7 +72,7 @@ public class Orbit
             float t = i / (float)steps;
             float a = t * Mathf.Tau;
 
-            var v2 = Ellipse.Focus0 + Ellipse.GetPosition(a);
+            var v2 = Ellipse.Foci[0] + Ellipse.GetPosition(a);
 
             geo.AddVertex(new Godot.Vector3
             {

scripts/orbits/OrbitSystem.cs

@@ -1,13 +1,13 @@
 using Godot;
 using System;
+using System.Collections.Generic;
 using Vim.Math3d;
 
 [Tool]
-public class OrbitSystem : Node, IMassive, ILocation
+public class OrbitSystem : Node, IPointMass
 {
-    [Export] private NodePath _a;
+    [Export]
-    [Export] private NodePath _b;
+    private NodePath[] _pointMassPaths;
-    [Export] private NodePath _barycenter;
 
     private double _semiMajorAxis;
     [Export]
@@ -40,66 +40,37 @@ public class OrbitSystem : Node, IMassive, ILocation
         }
     }
 
-    #region Point Masses
+    private IPointMass[] _pointMasses;
-    private IPointMass Primary
+    public IPointMass[] Orbiters
     {
         get
         {
-            if (_pointMasses[0] == null)
+            if (_pointMasses == null)
             {
-                InitPointMasses();
+                _pointMasses = new IPointMass[_pointMassPaths.Length];
-            }
+                // give each body an orbit around ourselves
-            return _pointMasses[0];
+                for (int i = 0; i < _pointMassPaths.Length; i++)
+                {
+                    _pointMasses[i] = GetNode<IPointMass>(_pointMassPaths[i]);
                 }
             }
-    private IPointMass Secondary
+            return _pointMasses;
+        }
+    }
+
+    public double Mass
     {
         get
         {
-            if (_pointMasses[1] == null)
+            double mass = 0;
+            foreach (var pointMass in _pointMasses)
             {
-                InitPointMasses();
+                mass += pointMass.Mass;
             }
-            return _pointMasses[1];
+            return mass;
         }
     }
-    private readonly IPointMass[] _pointMasses = new IPointMass[2];
+    public DVector3 Position { get; set; }
-    private void InitPointMasses()
-    {
-        var a = GetNode<IPointMass>(_a);
-        var b = GetNode<IPointMass>(_b);
-
-        if (a.Mass > b.Mass)
-        {
-            _pointMasses[0] = a;
-            _pointMasses[1] = b;
-        }
-        else
-        {
-            _pointMasses[0] = b;
-            _pointMasses[1] = a;
-        }
-    }
-    #endregion
-
-    public float Mass => Primary.Mass + Secondary.Mass;
-    public DVector3 Position
-    {
-        get => Barycenter;
-        set => Barycenter = value;
-    }
-
-    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
@@ -148,6 +119,6 @@ public class OrbitSystem : Node, IMassive, ILocation
         Orbit.Ellipse = new Ellipse(SemiMajorAxis, Eccentricity);
         Orbit.Draw(OrbitGeometry);
 
-        Secondary.Position = Orbit.GetPosition(_time);
+        Orbiter.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

@@ -6,7 +6,7 @@ using Vim.Math3d;
 public class Planet : Spatial, IPointMass
 {
     [Export]
-    public float Mass { get; set; }
+    public double Mass { get; set; }
 
     private DVector3? _position;
     public DVector3 Position

scripts/orbits/math/Ellipse.cs

@@ -2,7 +2,7 @@ using Godot;
 using System;
 using Vim.Math3d;
 
-public struct Ellipse
+public struct Ellipse : IEllipse
 {
     public Ellipse(double a = 1, double e = 0)
     {
@@ -12,6 +12,16 @@ public struct Ellipse
         // TODO: this is an immutable struct, so initialise everything else 
         // in the constructor to avoid recalculating properties whenever
         // they are accessed
+
+        Apisides = new Apisides(a, e);
+        b = Math.Sqrt(Apisides.max * Apisides.min);
+
+        var d = Math.Sqrt(a * a - b * b);
+        Foci = new DVector2[]
+        {
+            new DVector2(-d, 0),
+            new DVector2(d, 0)
+        };
     }
 
     /// <summary>
@@ -29,8 +39,7 @@ public struct Ellipse
     /// <summary>
     /// Semi-minor axis
     /// </summary>
-    public double b => SemiMinorAxis;
+    public double b { get; }
-    public double SemiMinorAxis => Math.Sqrt(Apisides.max * Apisides.min);
 
     /// <summary>
     /// Get a position on the auxiliary circle
@@ -40,13 +49,10 @@ public struct Ellipse
     public DVector2 GetAuxiliaryPosition2D(double t = 0) =>
         new DVector2(Math.Cos(t), Math.Sin(t)) * a;
 
-    private Apisides Apisides => new Apisides(a, e);
+    private Apisides Apisides { get; }
 
-    public DVector2 GetPosition(double t) =>
+    public DVector2[] Foci { get; }
-        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);
 }
 

scripts/orbits/math/IEllipse.cs

@@ -0,0 +1,20 @@
+using Godot;
+using System;
+using Vim.Math3d;
+
+public interface IEllipse
+{
+    /// <summary>
+    /// Semi-major axis
+    /// </summary>
+    double a { get; }
+    double b { get; }
+    /// <summary>
+    /// Eccentricity
+    /// </summary>
+    double e { get; }
+    /// <summary>
+    /// Positions of the two foci of the ellipse
+    /// </summary>
+    DVector2[] Foci { get; }
+}

scripts/orbits/math/IEllipseExtensions.cs

@@ -0,0 +1,8 @@
+using Vim.Math3d;
+using System;
+
+public static class IEllipseExtensions
+{
+    public static DVector2 GetPosition(this IEllipse e, double t) =>
+        new DVector2(Math.Cos(t) * e.a, Math.Sin(t) * e.b);
+}