diff --git a/scenes/Main.tscn b/scenes/Main.tscn
index 7a192cf..13bd289 100644
--- a/scenes/Main.tscn
+++ b/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
diff --git a/scenes/orbits/OrbitSystem.tscn b/scenes/orbits/Barycenter.tscn
similarity index 72%
rename from scenes/orbits/OrbitSystem.tscn
rename to scenes/orbits/Barycenter.tscn
index 3a5c227..140eba6 100644
--- a/scenes/orbits/OrbitSystem.tscn
+++ b/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")
-
-[node name="Barycenter" type="Spatial" parent="."]
+_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, 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
diff --git a/scripts/orbits/IMassive.cs b/scripts/orbits/IMassive.cs
index a7d5f80..f5fac4b 100644
--- a/scripts/orbits/IMassive.cs
+++ b/scripts/orbits/IMassive.cs
@@ -3,5 +3,5 @@ using System;
public interface IMassive
{
- float Mass { get; }
+ double Mass { get; }
}
diff --git a/scripts/orbits/IMassiveExtensions.cs b/scripts/orbits/IMassiveExtensions.cs
new file mode 100644
index 0000000..6f4cf1a
--- /dev/null
+++ b/scripts/orbits/IMassiveExtensions.cs
@@ -0,0 +1,4 @@
+public static class IMassiveExtensions
+{
+ public static double U(this IMassive m) => Kepler.U(m.Mass);
+}
\ No newline at end of file
diff --git a/scripts/orbits/IOrbit.cs b/scripts/orbits/IOrbit.cs
new file mode 100644
index 0000000..20e8e0f
--- /dev/null
+++ b/scripts/orbits/IOrbit.cs
@@ -0,0 +1,6 @@
+using Vim.Math3d;
+
+public interface IOrbit : IEllipse
+{
+ DVector2 GetPosition(float t);
+}
\ No newline at end of file
diff --git a/scripts/orbits/IOrbitalElementsExtensions.cs b/scripts/orbits/IOrbitalElementsExtensions.cs
new file mode 100644
index 0000000..f8dbda1
--- /dev/null
+++ b/scripts/orbits/IOrbitalElementsExtensions.cs
@@ -0,0 +1,13 @@
+using Vim.Math3d;
+using static System.Math;
+
+public static class IOrbitalElementsExtensions
+{
+ ///
+ /// Get the Cartesian position on the orbit for a given mean anomaly.
+ ///
+ /// Mean anomaly
+ /// Position
+ public static DVector2 GetEccentricPosition2D(this OrbitalElements elements, double M) =>
+ Kepler.GetEccentricPosition2D(elements.e, elements.a, M);
+}
\ No newline at end of file
diff --git a/scripts/orbits/Orbit.cs b/scripts/orbits/Orbit.cs
index 3ebab6c..8be8617 100644
--- a/scripts/orbits/Orbit.cs
+++ b/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
- // orbit itself is presented as being centred on the primary, but
- // all our math is elliptical.
- private DVector2 PrimaryPosition => Ellipse.Focus0;
+ #region IEllipse
+ public double a => Ellipse.a;
+ public double b => Ellipse.b;
+ 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;
+
+ ///
+ /// An orbit is treated as one stationary body at one focus of the orbit's ellipse.
+ /// The satellite is treated as a massless point.
+ ///
+ /// Mass of the object being orbited
+ /// The orbiting body
+ 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
{
diff --git a/scripts/orbits/OrbitSystem.cs b/scripts/orbits/OrbitSystem.cs
index e60b2db..fb02ef1 100644
--- a/scripts/orbits/OrbitSystem.cs
+++ b/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] private NodePath _b;
- [Export] private NodePath _barycenter;
+ [Export]
+ private NodePath[] _pointMassPaths;
private double _semiMajorAxis;
[Export]
@@ -40,66 +40,37 @@ public class OrbitSystem : Node, IMassive, ILocation
}
}
- #region Point Masses
- private IPointMass Primary
+ private IPointMass[] _pointMasses;
+ public IPointMass[] Orbiters
{
get
{
- if (_pointMasses[0] == null)
+ if (_pointMasses == null)
{
- InitPointMasses();
+ _pointMasses = new IPointMass[_pointMassPaths.Length];
+ // give each body an orbit around ourselves
+ for (int i = 0; i < _pointMassPaths.Length; i++)
+ {
+ _pointMasses[i] = GetNode(_pointMassPaths[i]);
+ }
}
- return _pointMasses[0];
+ return _pointMasses;
}
}
- private IPointMass Secondary
+
+ 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];
- private void InitPointMasses()
- {
- var a = GetNode(_a);
- var b = GetNode(_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;
- }
+ public DVector3 Position { get; set; }
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);
}
}
diff --git a/scripts/orbits/OrbitalElements.cs b/scripts/orbits/OrbitalElements.cs
new file mode 100644
index 0000000..c501dc8
--- /dev/null
+++ b/scripts/orbits/OrbitalElements.cs
@@ -0,0 +1,46 @@
+public struct OrbitalElements
+{
+ private Ellipse _ellipse;
+ ///
+ /// Semi-major axis in AU
+ ///
+ public double a => _ellipse.a;
+ ///
+ /// Eccentricity
+ ///
+ public double e => _ellipse.e;
+ ///
+ /// Inclination in degrees
+ ///
+ public double i { get; private set; }
+ ///
+ /// Argument of ascending node in degrees
+ ///
+ public double W { get; private set; }
+ ///
+ /// Argument of periapsis in degrees
+ ///
+ public double p { get; private set; }
+ ///
+ /// Mean anomaly at epoch
+ ///
+ public double M0 { get; private set; }
+
+ ///
+ ///
+ ///
+ /// Semi-major axis in AU
+ /// eccentricity
+ /// Inclination in degrees
+ /// Argument of ascending node in degrees
+ /// Argument of periapsis in degrees
+ /// Mean anomaly at epoch
+ 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;
+ }
+}
\ No newline at end of file
diff --git a/scripts/orbits/Planet.cs b/scripts/orbits/Planet.cs
index b12c0bc..815a8f9 100644
--- a/scripts/orbits/Planet.cs
+++ b/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
diff --git a/scripts/orbits/math/Ellipse.cs b/scripts/orbits/math/Ellipse.cs
index e1fff2e..22292c1 100644
--- a/scripts/orbits/math/Ellipse.cs
+++ b/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)
+ };
}
///
@@ -29,8 +39,7 @@ public struct Ellipse
///
/// Semi-minor axis
///
- public double b => SemiMinorAxis;
- public double SemiMinorAxis => Math.Sqrt(Apisides.max * Apisides.min);
+ public double b { get; }
///
/// 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) =>
- new DVector2(Math.Cos(t) * a, Math.Sin(t) * b);
+ public DVector2[] Foci { get; }
- public DVector2 Focus0 => new DVector2(-GetFocusDistance(), 0);
- public DVector2 Focus1 => new DVector2(GetFocusDistance(), 0);
private double GetFocusDistance() => Math.Sqrt(a * a - b * b);
}
diff --git a/scripts/orbits/math/IEllipse.cs b/scripts/orbits/math/IEllipse.cs
new file mode 100644
index 0000000..c7b3d7b
--- /dev/null
+++ b/scripts/orbits/math/IEllipse.cs
@@ -0,0 +1,20 @@
+using Godot;
+using System;
+using Vim.Math3d;
+
+public interface IEllipse
+{
+ ///
+ /// Semi-major axis
+ ///
+ double a { get; }
+ double b { get; }
+ ///
+ /// Eccentricity
+ ///
+ double e { get; }
+ ///
+ /// Positions of the two foci of the ellipse
+ ///
+ DVector2[] Foci { get; }
+}
diff --git a/scripts/orbits/math/IEllipseExtensions.cs b/scripts/orbits/math/IEllipseExtensions.cs
new file mode 100644
index 0000000..40f1137
--- /dev/null
+++ b/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);
+}
\ No newline at end of file