using static System.Math;
using Vim.Math3d;
public static class Kepler
{
///
/// 2 * PI
///
private static double TAU => PI * 2.0;
///
/// Universal gravitational constant
///
private static double G => 6.677e-11;
///
/// Get the mean anomaly of an orbit given at a given time, given
/// the period of the orbit.
///
/// Time of mean anomaly
/// Orbital period
///
public static double GetMeanAnomaly(double time, double P)
{
var M = ((time % P / P) * TAU) % TAU;
return M % TAU;
}
///
///
///
/// Orbital period
///
public static double GetMeanAngularVelocity(double P) => TAU / P;
///
///
///
/// Semimajor axis in AU
/// Standard gravitational parameter
///
public static double GetPeriod(double a, double U)
{
var a3 = a * a * a;
var T = TAU * Sqrt(a3 / U);
return T;
}
///
/// Get the standard gravitational parameter of a mass
///
/// Mass in kilograms
/// Standard gravitational parameter
public static double U(double M) => M * G;
///
/// Get the standard gravitational parameter of a mass
///
/// Mass in kilograms
/// Standard gravitational parameter
public static double GetStandardGravitationalParameter(double M) => M * G;
///
/// Get the eccentric anomaly for a given mean anomaly.
///
/// Eccentricity
/// Mean anomaly
/// Eccentric anomaly
public static double GetEccentricAnomaly(double e, double M) =>
GetEccentricAnomalyBS(e, M);
///
/// Get the eccentric anomaly for a given mean anomaly
/// using a binary search.
///
/// Eccentricity
/// Mean anomaly
/// Eccentric anomaly
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;
}
///
/// Get the Cartesian position on an orbit given its eccentricity,
/// its semi-major axis and a given mean anomaly.
///
/// Eccentricity
/// Semi-major axis
/// Mean anomaly
/// Position on the orbit
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);
}
}