update main scene

scenes/orbits/Barycenter.tscn

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

scripts/orbits/Barycenter.cs

@@ -4,7 +4,7 @@ using System.Collections.Generic;
 using Vim.Math3d;
 
 [Tool]
-public class OrbitSystem : Node, IPointMass
+public class Barycenter : Node, IPointMass
 {
     [Export]
     private NodePath[] _pointMassPaths;
@@ -40,51 +40,20 @@ public class OrbitSystem : 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 pointMass in _pointMasses)
+            foreach (var orbiter in Orbiters)
             {
-                mass += pointMass.Mass;
+                mass += orbiter.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
     {
@@ -103,6 +72,50 @@ public class OrbitSystem : 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]
@@ -116,9 +129,15 @@ public class OrbitSystem : Node, IPointMass
 
     private void InvalidateGeometry()
     {
-        Orbit.Ellipse = new Ellipse(SemiMajorAxis, Eccentricity);
-        Orbit.Draw(OrbitGeometry);
+        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);
 
-        Orbiter.Position = Orbit.GetPosition(_time);
+            var t = _time + i * Math.PI;
+            orbiter.PointMass.Position = orbit.GetPosition(t);
+        }
     }
 }

scripts/orbits/IOrbit.cs

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

scripts/orbits/IOrbitExtensions.cs

@@ -0,0 +1,29 @@
+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

@@ -1,13 +0,0 @@
-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,15 +4,7 @@ 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;
@@ -25,9 +17,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);
 
@@ -35,16 +27,12 @@ public class Orbit : IOrbit
         _satellite = satellite;
     }
 
-    public DVector3 GetPosition(float t)
+    public DVector3 GetPosition(double 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 P = Kepler.GetPeriod(a, _massive.U());
 
         var M = Kepler.GetMeanAnomaly(t, P);
         var E = Kepler.GetEccentricAnomaly(e, M) + Math.PI;
@@ -52,33 +40,4 @@ 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/OrbitalElements.cs

@@ -1,46 +0,0 @@
-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;
-    }
-}