WIP - RESET THIS

WIP RESET THIS

scenes/orbits/Barycenter.tscn

@@ -1,18 +1,19 @@
 [gd_scene load_steps=3 format=2]
 
-[ext_resource path="res://scripts/orbits/Barycenter.cs" type="Script" id=1]
+[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="Barycenter" type="Spatial"]
 script = ExtResource( 1 )
-SemiMajorAxis = 7.0
-Eccentricity = 0.752
-_pointMassPaths = [ NodePath("Planet A"), NodePath("Planet B") ]
-_speed = 1.0
+SemiMajorAxis = 6.881
+Eccentricity = 0.399
+_a = NodePath("Planet A")
+_b = NodePath("Planet B")
+_barycenter = NodePath("")
 
 [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, -2.55711, 0, -4.25519 )
+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, -11.9998, 0, 1.25565 )
+transform = Transform( 0.5, 0, 0, 0, 0.5, 0, 0, 0, 0.5, -9.62645, 0, -0.0287273 )
 Mass = 0.125

scripts/orbits/IOrbit.cs

@@ -1,7 +1,6 @@
 using Vim.Math3d;
 
-public interface IOrbit
+public interface IOrbit : IEllipse
 {
-    DVector3 GetPosition(double t);
-    IEllipse Ellipse { get; set; }
+    DVector2 GetPosition(float t);
 }

scripts/orbits/IOrbitExtensions.cs

@@ -1,29 +0,0 @@
-using Godot;
-
-public static class IOrbitExtensions
-{
-    public static void Draw(this IOrbit orbit, ImmediateGeometry geo)
-    {
-        geo.Clear();
-        geo.Begin(Mesh.PrimitiveType.LineLoop);
-        geo.SetColor(new Color(1, 0, 0));
-
-        int steps = 100;
-        var ellipse = orbit.Ellipse;
-
-        for (int i = 0; i < steps; i++)
-        {
-            float t = i / (float)steps;
-            float a = t * Mathf.Tau;
-
-            var v2 = ellipse.Foci[0] + ellipse.GetPosition(a);
-
-            geo.AddVertex(new Godot.Vector3
-            {
-                x = (float)v2.X,
-                z = (float)v2.Y
-            });
-        }
-        geo.End();
-    }
-}

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

@@ -4,7 +4,15 @@ using Vim.Math3d;
 
 public class Orbit : IOrbit
 {
+    #region IEllipse
+    public double a => Ellipse.a;
+    public double b => Ellipse.b;
+    public double e => Ellipse.e;
+    public DVector2[] Foci => Ellipse.Foci;
+    #endregion
+
     public IEllipse Ellipse { get; set; }
+    private double U => Kepler.U(_massive.Mass);
 
     private readonly IMassive _massive;
     private readonly ILocation _satellite;
@@ -17,9 +25,9 @@ public class Orbit : IOrbit
     /// <param name="satellite">The orbiting body</param>
     public Orbit(IMassive massive, ILocation satellite)
     {
-        // TODO: determine ellipse better
         // 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);
 
@@ -27,12 +35,16 @@ public class Orbit : IOrbit
         _satellite = satellite;
     }
 
-    public DVector3 GetPosition(double t)
+    public DVector3 GetPosition(float t)
     {
         var a = Ellipse.a;
         var e = Ellipse.e;
 
-        var P = Kepler.GetPeriod(a, _massive.U());
+        Kepler.GetPeriod(a, U);
+
+        // TODO: determine period from mass of orbiting bodies
+        var P = 10f;
+
 
         var M = Kepler.GetMeanAnomaly(t, P);
         var E = Kepler.GetEccentricAnomaly(e, M) + Math.PI;
@@ -40,4 +52,33 @@ public class Orbit : IOrbit
         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.Foci[0] + Ellipse.GetPosition(a);
+
+            geo.AddVertex(new Godot.Vector3
+            {
+                x = (float)v2.X,
+                z = (float)v2.Y
+            });
+        }
+    }
 }

scripts/orbits/OrbitSystem.cs

@@ -4,7 +4,7 @@ using System.Collections.Generic;
 using Vim.Math3d;
 
 [Tool]
-public class Barycenter : Node, IPointMass
+public class OrbitSystem : Node, IPointMass
 {
     [Export]
     private NodePath[] _pointMassPaths;
@@ -40,20 +40,51 @@ public class Barycenter : Node, IPointMass
         }
     }
 
+    private IPointMass[] _pointMasses;
+    public IPointMass[] Orbiters
+    {
+        get
+        {
+            if (_pointMasses == null)
+            {
+                _pointMasses = new IPointMass[_pointMassPaths.Length];
+                // give each body an orbit around ourselves
+                for (int i = 0; i < _pointMassPaths.Length; i++)
+                {
+                    _pointMasses[i] = GetNode<IPointMass>(_pointMassPaths[i]);
+                }
+            }
+            return _pointMasses;
+        }
+    }
+
     public double Mass
     {
         get
         {
             double mass = 0;
-            foreach (var orbiter in Orbiters)
+            foreach (var pointMass in _pointMasses)
             {
-                mass += orbiter.PointMass.Mass;
+                mass += pointMass.Mass;
             }
             return mass;
         }
     }
     public DVector3 Position { get; set; }
 
+    private Orbit _orbit = null;
+    private Orbit Orbit
+    {
+        get
+        {
+            if (_orbit == null)
+            {
+                _orbit = new Orbit();
+            }
+            return _orbit;
+        }
+    }
+
     private ImmediateGeometry _orbitGeometry = null;
     public ImmediateGeometry OrbitGeometry
     {
@@ -72,50 +103,6 @@ public class Barycenter : Node, IPointMass
         }
     }
 
-    private struct Orbiter
-    {
-        public IPointMass PointMass { get; set; }
-        public IOrbit Orbit { get; set; }
-    }
-    private Orbiter[] _orbiters = null;
-    private Orbiter[] Orbiters
-    {
-        get
-        {
-            if (_orbiters == null)
-            {
-                var length = _pointMassPaths.Length;
-                _orbiters = new Orbiter[length];
-                for (int i = 0; i < length; i++)
-                {
-                    var pointMass = GetNode<IPointMass>(_pointMassPaths[i]);
-                    var orbit = new Orbit(this, pointMass);
-
-                    var orbiter = new Orbiter
-                    {
-                        PointMass = pointMass,
-                        Orbit = orbit
-                    };
-                    _orbiters[i] = orbiter;
-                }
-            }
-            return _orbiters;
-        }
-    }
-
-    private IOrbit[] _orbits;
-    private IOrbit[] Orbits
-    {
-        get
-        {
-            if (_orbits == null)
-            {
-
-            }
-            return _orbits;
-        }
-    }
-
     private float _time = 0;
 
     [Export]
@@ -129,15 +116,9 @@ public class Barycenter : Node, IPointMass
 
     private void InvalidateGeometry()
     {
-        for (int i = 0; i < Orbiters.Length; i++)
-        {
-            var orbiter = Orbiters[i];
-            var orbit = orbiter.Orbit;
-            orbit.Ellipse = new Ellipse(SemiMajorAxis, Eccentricity);
-            orbit.Draw(OrbitGeometry);
+        Orbit.Ellipse = new Ellipse(SemiMajorAxis, Eccentricity);
+        Orbit.Draw(OrbitGeometry);
 
-            var t = _time + i * Math.PI;
-            orbiter.PointMass.Position = orbit.GetPosition(t);
-        }
+        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;
+    }
+}