add move according to kepler

This commit is contained in:
Cat Flynn 2022-09-05 22:12:56 +01:00
parent 322952197f
commit 8c1bb4dadf
4 changed files with 144 additions and 6 deletions

View File

@ -14,4 +14,5 @@ transform = Transform( 0.515898, 0.606099, -0.605386, -0.393123, 0.795389, 0.461
transform = Transform( 0.856869, 0.3292, -0.396741, 0.0949996, 0.655565, 0.749139, 0.506706, -0.679604, 0.53046, -0.599122, 0, 0 ) 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 SemiMajorAxis = 6.166
Eccentricity = 0.239 Eccentricity = 0.239
Period = 1.0
_speed = 0.877 _speed = 0.877

View File

@ -5,8 +5,9 @@
[node name="Orbit System" type="Spatial"] [node name="Orbit System" type="Spatial"]
script = ExtResource( 1 ) script = ExtResource( 1 )
SemiMajorAxis = 5.994 SemiMajorAxis = 6.881
Eccentricity = 0.518 Eccentricity = 0.399
Period = 10.0
_a = NodePath("Planet A") _a = NodePath("Planet A")
_b = NodePath("Planet B") _b = NodePath("Planet B")
_barycenter = NodePath("Barycenter") _barycenter = NodePath("Barycenter")
@ -17,5 +18,5 @@ _barycenter = NodePath("Barycenter")
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, 0, 0, 0 )
[node name="Planet B" parent="." instance=ExtResource( 2 )] [node name="Planet B" parent="." instance=ExtResource( 2 )]
transform = Transform( 0.5, 0, 0, 0, 0.5, 0, 0, 0, 0.5, -6.57314, 0, 4.18169 ) transform = Transform( 0.5, 0, 0, 0, 0.5, 0, 0, 0, 0.5, -8.44634, 0, -3.53339 )
Mass = 0.125 Mass = 0.125

View File

@ -4,16 +4,25 @@ using Vim.Math3d;
public class Orbit public class Orbit
{ {
// the position of the primary body relative to the ellipse - the // The position of the primary body relative to the ellipse - the
// orbit itself is presented as being centred on the primary, but // orbit itself is presented as being centred on the primary, but
// all our math is ellipse stuff // all our math is elliptical.
private DVector2 PrimaryPosition => Ellipse.Focus0; private DVector2 PrimaryPosition => Ellipse.Focus0;
public Ellipse Ellipse { get; set; } public Ellipse Ellipse { get; set; }
public DVector3 GetPosition(float t) public DVector3 GetPosition(float t)
{ {
var p = PrimaryPosition + Ellipse.GetPosition(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); return new DVector3(p.X, 0, p.Y);
} }

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