Compare commits

...

3 Commits

Author SHA1 Message Date
Cat Flynn 6f79458d66 WIP - RESET THIS 2022-09-05 22:14:30 +01:00
Cat Flynn bc28d4d172 move according to kepler 2022-09-05 22:14:27 +01:00
ktyl f28670c083 implement ellipse
move planet around ellipse

adjust orbit speed
2022-09-05 22:14:24 +01:00
16 changed files with 499 additions and 112 deletions

View File

@ -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

View File

@ -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="."]

View File

@ -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

11
scenes/orbits/Planet.tscn Normal file
View File

@ -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

View File

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

View File

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

7
scripts/orbits/IOrbit.cs Normal file
View File

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

View File

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

View File

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

57
scripts/orbits/Orbit.cs Normal file
View File

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

View File

@ -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];
public Vector3 Position => Barycenter;
private double _semiMajorAxis;
[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;
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();
}
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
if (_pointMasses[0] == null)
{
x = Mathf.Sin(a),
z = Mathf.Cos(a)
});
InitPointMasses();
}
return _pointMasses[0];
}
_orbit.End();
}
private void InitGeometry()
private IPointMass Secondary
{
_orbit = new ImmediateGeometry();
var m = new SpatialMaterial();
m.VertexColorUseAsAlbedo = true;
m.FlagsUnshaded = true;
_orbit.MaterialOverride = m;
AddChild(_orbit);
get
{
if (_pointMasses[1] == null)
{
InitPointMasses();
}
return _pointMasses[1];
}
}
private readonly IPointMass[] _pointMasses = new IPointMass[2];
private void InitPointMasses()
{
var a = GetPointMass(_a);
var b = GetPointMass(_b);
var a = GetNode<IPointMass>(_a);
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);
//var massive = node.GetChild<IMassive>(0);
get => Barycenter;
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);
}
}

View File

@ -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;
}
}

View File

@ -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)
// {
//
// }
}
}

View File

@ -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;
}
}

View File

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

View File

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