mirror of
https://github.com/cosinekitty/astronomy.git
synced 2026-05-19 06:17:03 -04:00
Implemented C# Lagrange point functions.
This commit is contained in:
Don Cross
@@ -2720,11 +2720,15 @@ namespace csharp_test
|
||||
return sqrt(diff_squared);
|
||||
}
|
||||
|
||||
interface IStateVectorFunc
|
||||
{
|
||||
StateVector Eval(AstroTime time);
|
||||
}
|
||||
|
||||
static int VerifyState(
|
||||
Func<Body, AstroTime, StateVector> func,
|
||||
IStateVectorFunc func,
|
||||
ref double max_rdiff,
|
||||
ref double max_vdiff,
|
||||
Body body,
|
||||
string filename,
|
||||
int lnum,
|
||||
AstroTime time,
|
||||
@@ -2733,7 +2737,7 @@ namespace csharp_test
|
||||
double r_thresh,
|
||||
double v_thresh)
|
||||
{
|
||||
StateVector state = func(body, time);
|
||||
StateVector state = func.Eval(time);
|
||||
|
||||
double rdiff = StateVectorDiff((r_thresh > 0.0), pos, state.x, state.y, state.z);
|
||||
if (rdiff > max_rdiff)
|
||||
@@ -2759,8 +2763,7 @@ namespace csharp_test
|
||||
}
|
||||
|
||||
static int VerifyStateBody(
|
||||
Func<Body, AstroTime, StateVector> func,
|
||||
Body body,
|
||||
IStateVectorFunc func,
|
||||
string filename,
|
||||
double r_thresh,
|
||||
double v_thresh)
|
||||
@@ -2777,7 +2780,7 @@ namespace csharp_test
|
||||
vel[0] = rec.state.vx;
|
||||
vel[1] = rec.state.vy;
|
||||
vel[2] = rec.state.vz;
|
||||
if (0 != VerifyState(func, ref max_rdiff, ref max_vdiff, body, filename, rec.lnum, rec.state.t, pos, vel, r_thresh, v_thresh))
|
||||
if (0 != VerifyState(func, ref max_rdiff, ref max_vdiff, filename, rec.lnum, rec.state.t, pos, vel, r_thresh, v_thresh))
|
||||
return 1;
|
||||
++count;
|
||||
}
|
||||
@@ -2788,89 +2791,124 @@ namespace csharp_test
|
||||
// Constants for use inside unit tests only; they doesn't make sense for public consumption.
|
||||
const Body Body_GeoMoon = (Body)(-100);
|
||||
const Body Body_Geo_EMB = (Body)(-101);
|
||||
static StateVector BaryState(Body body, AstroTime time)
|
||||
|
||||
class BaryStateFunc : IStateVectorFunc
|
||||
{
|
||||
if (body == Body_GeoMoon)
|
||||
return Astronomy.GeoMoonState(time);
|
||||
private Body body;
|
||||
|
||||
if (body == Body_Geo_EMB)
|
||||
return Astronomy.GeoEmbState(time);
|
||||
public BaryStateFunc(Body body)
|
||||
{
|
||||
this.body = body;
|
||||
}
|
||||
|
||||
return Astronomy.BaryState(body, time);
|
||||
public StateVector Eval(AstroTime time)
|
||||
{
|
||||
if (body == Body_GeoMoon)
|
||||
return Astronomy.GeoMoonState(time);
|
||||
|
||||
if (body == Body_Geo_EMB)
|
||||
return Astronomy.GeoEmbState(time);
|
||||
|
||||
return Astronomy.BaryState(body, time);
|
||||
}
|
||||
}
|
||||
|
||||
static int BaryStateTest()
|
||||
{
|
||||
if (0 != VerifyStateBody(BaryState, Body.Sun, "../../barystate/Sun.txt", -1.224e-05, -1.134e-07)) return 1;
|
||||
if (0 != VerifyStateBody(BaryState, Body.Mercury, "../../barystate/Mercury.txt", 1.672e-04, 2.698e-04)) return 1;
|
||||
if (0 != VerifyStateBody(BaryState, Body.Venus, "../../barystate/Venus.txt", 4.123e-05, 4.308e-05)) return 1;
|
||||
if (0 != VerifyStateBody(BaryState, Body.Earth, "../../barystate/Earth.txt", 2.296e-05, 6.359e-05)) return 1;
|
||||
if (0 != VerifyStateBody(BaryState, Body.Mars, "../../barystate/Mars.txt", 3.107e-05, 5.550e-05)) return 1;
|
||||
if (0 != VerifyStateBody(BaryState, Body.Jupiter, "../../barystate/Jupiter.txt", 7.389e-05, 2.471e-04)) return 1;
|
||||
if (0 != VerifyStateBody(BaryState, Body.Saturn, "../../barystate/Saturn.txt", 1.067e-04, 3.220e-04)) return 1;
|
||||
if (0 != VerifyStateBody(BaryState, Body.Uranus, "../../barystate/Uranus.txt", 9.035e-05, 2.519e-04)) return 1;
|
||||
if (0 != VerifyStateBody(BaryState, Body.Neptune, "../../barystate/Neptune.txt", 9.838e-05, 4.446e-04)) return 1;
|
||||
if (0 != VerifyStateBody(BaryState, Body.Pluto, "../../barystate/Pluto.txt", 4.259e-05, 7.827e-05)) return 1;
|
||||
if (0 != VerifyStateBody(BaryState, Body.Moon, "../../barystate/Moon.txt", 2.354e-05, 6.604e-05)) return 1;
|
||||
if (0 != VerifyStateBody(BaryState, Body.EMB, "../../barystate/EMB.txt", 2.353e-05, 6.511e-05)) return 1;
|
||||
if (0 != VerifyStateBody(BaryState, Body_GeoMoon, "../../barystate/GeoMoon.txt", 4.086e-05, 5.347e-05)) return 1;
|
||||
if (0 != VerifyStateBody(BaryState, Body_Geo_EMB, "../../barystate/GeoEMB.txt", 4.076e-05, 5.335e-05)) return 1;
|
||||
if (0 != VerifyStateBody(new BaryStateFunc(Body.Sun), "../../barystate/Sun.txt", -1.224e-05, -1.134e-07)) return 1;
|
||||
if (0 != VerifyStateBody(new BaryStateFunc(Body.Mercury), "../../barystate/Mercury.txt", 1.672e-04, 2.698e-04)) return 1;
|
||||
if (0 != VerifyStateBody(new BaryStateFunc(Body.Venus), "../../barystate/Venus.txt", 4.123e-05, 4.308e-05)) return 1;
|
||||
if (0 != VerifyStateBody(new BaryStateFunc(Body.Earth), "../../barystate/Earth.txt", 2.296e-05, 6.359e-05)) return 1;
|
||||
if (0 != VerifyStateBody(new BaryStateFunc(Body.Mars), "../../barystate/Mars.txt", 3.107e-05, 5.550e-05)) return 1;
|
||||
if (0 != VerifyStateBody(new BaryStateFunc(Body.Jupiter), "../../barystate/Jupiter.txt", 7.389e-05, 2.471e-04)) return 1;
|
||||
if (0 != VerifyStateBody(new BaryStateFunc(Body.Saturn), "../../barystate/Saturn.txt", 1.067e-04, 3.220e-04)) return 1;
|
||||
if (0 != VerifyStateBody(new BaryStateFunc(Body.Uranus), "../../barystate/Uranus.txt", 9.035e-05, 2.519e-04)) return 1;
|
||||
if (0 != VerifyStateBody(new BaryStateFunc(Body.Neptune), "../../barystate/Neptune.txt", 9.838e-05, 4.446e-04)) return 1;
|
||||
if (0 != VerifyStateBody(new BaryStateFunc(Body.Pluto), "../../barystate/Pluto.txt", 4.259e-05, 7.827e-05)) return 1;
|
||||
if (0 != VerifyStateBody(new BaryStateFunc(Body.Moon), "../../barystate/Moon.txt", 2.354e-05, 6.604e-05)) return 1;
|
||||
if (0 != VerifyStateBody(new BaryStateFunc(Body.EMB), "../../barystate/EMB.txt", 2.353e-05, 6.511e-05)) return 1;
|
||||
if (0 != VerifyStateBody(new BaryStateFunc(Body_GeoMoon), "../../barystate/GeoMoon.txt", 4.086e-05, 5.347e-05)) return 1;
|
||||
if (0 != VerifyStateBody(new BaryStateFunc(Body_Geo_EMB), "../../barystate/GeoEMB.txt", 4.076e-05, 5.335e-05)) return 1;
|
||||
Console.WriteLine("C# BaryStateTest: PASS");
|
||||
return 0;
|
||||
}
|
||||
|
||||
class HelioStateFunc : IStateVectorFunc
|
||||
{
|
||||
private Body body;
|
||||
|
||||
public HelioStateFunc(Body body)
|
||||
{
|
||||
this.body = body;
|
||||
}
|
||||
|
||||
public StateVector Eval(AstroTime time)
|
||||
{
|
||||
return Astronomy.HelioState(body, time);
|
||||
}
|
||||
}
|
||||
|
||||
static int HelioStateTest()
|
||||
{
|
||||
if (0 != VerifyStateBody(Astronomy.HelioState, Body.SSB, "../../heliostate/SSB.txt", -1.209e-05, -1.125e-07)) return 1;
|
||||
if (0 != VerifyStateBody(Astronomy.HelioState, Body.Mercury, "../../heliostate/Mercury.txt", 1.481e-04, 2.756e-04)) return 1;
|
||||
if (0 != VerifyStateBody(Astronomy.HelioState, Body.Venus, "../../heliostate/Venus.txt", 3.528e-05, 4.485e-05)) return 1;
|
||||
if (0 != VerifyStateBody(Astronomy.HelioState, Body.Earth, "../../heliostate/Earth.txt", 1.476e-05, 6.105e-05)) return 1;
|
||||
if (0 != VerifyStateBody(Astronomy.HelioState, Body.Mars, "../../heliostate/Mars.txt", 3.154e-05, 5.603e-05)) return 1;
|
||||
if (0 != VerifyStateBody(Astronomy.HelioState, Body.Jupiter, "../../heliostate/Jupiter.txt", 7.455e-05, 2.562e-04)) return 1;
|
||||
if (0 != VerifyStateBody(Astronomy.HelioState, Body.Saturn, "../../heliostate/Saturn.txt", 1.066e-04, 3.150e-04)) return 1;
|
||||
if (0 != VerifyStateBody(Astronomy.HelioState, Body.Uranus, "../../heliostate/Uranus.txt", 9.034e-05, 2.712e-04)) return 1;
|
||||
if (0 != VerifyStateBody(Astronomy.HelioState, Body.Neptune, "../../heliostate/Neptune.txt", 9.834e-05, 4.534e-04)) return 1;
|
||||
if (0 != VerifyStateBody(Astronomy.HelioState, Body.Pluto, "../../heliostate/Pluto.txt", 4.271e-05, 1.198e-04)) return 1;
|
||||
if (0 != VerifyStateBody(Astronomy.HelioState, Body.Moon, "../../heliostate/Moon.txt", 1.477e-05, 6.195e-05)) return 1;
|
||||
if (0 != VerifyStateBody(Astronomy.HelioState, Body.EMB, "../../heliostate/EMB.txt", 1.476e-05, 6.106e-05)) return 1;
|
||||
if (0 != VerifyStateBody(new HelioStateFunc(Body.SSB), "../../heliostate/SSB.txt", -1.209e-05, -1.125e-07)) return 1;
|
||||
if (0 != VerifyStateBody(new HelioStateFunc(Body.Mercury), "../../heliostate/Mercury.txt", 1.481e-04, 2.756e-04)) return 1;
|
||||
if (0 != VerifyStateBody(new HelioStateFunc(Body.Venus), "../../heliostate/Venus.txt", 3.528e-05, 4.485e-05)) return 1;
|
||||
if (0 != VerifyStateBody(new HelioStateFunc(Body.Earth), "../../heliostate/Earth.txt", 1.476e-05, 6.105e-05)) return 1;
|
||||
if (0 != VerifyStateBody(new HelioStateFunc(Body.Mars), "../../heliostate/Mars.txt", 3.154e-05, 5.603e-05)) return 1;
|
||||
if (0 != VerifyStateBody(new HelioStateFunc(Body.Jupiter), "../../heliostate/Jupiter.txt", 7.455e-05, 2.562e-04)) return 1;
|
||||
if (0 != VerifyStateBody(new HelioStateFunc(Body.Saturn), "../../heliostate/Saturn.txt", 1.066e-04, 3.150e-04)) return 1;
|
||||
if (0 != VerifyStateBody(new HelioStateFunc(Body.Uranus), "../../heliostate/Uranus.txt", 9.034e-05, 2.712e-04)) return 1;
|
||||
if (0 != VerifyStateBody(new HelioStateFunc(Body.Neptune), "../../heliostate/Neptune.txt", 9.834e-05, 4.534e-04)) return 1;
|
||||
if (0 != VerifyStateBody(new HelioStateFunc(Body.Pluto), "../../heliostate/Pluto.txt", 4.271e-05, 1.198e-04)) return 1;
|
||||
if (0 != VerifyStateBody(new HelioStateFunc(Body.Moon), "../../heliostate/Moon.txt", 1.477e-05, 6.195e-05)) return 1;
|
||||
if (0 != VerifyStateBody(new HelioStateFunc(Body.EMB), "../../heliostate/EMB.txt", 1.476e-05, 6.106e-05)) return 1;
|
||||
Console.WriteLine("C# HelioStateTest: PASS");
|
||||
return 0;
|
||||
}
|
||||
|
||||
static StateVector TopoStateFunc(Body body, AstroTime time)
|
||||
class TopoStateFunc : IStateVectorFunc
|
||||
{
|
||||
var observer = new Observer(30.0, -80.0, 1000.0);
|
||||
private Body body;
|
||||
|
||||
StateVector observer_state = Astronomy.ObserverState(time, observer, EquatorEpoch.J2000);
|
||||
StateVector state;
|
||||
if (body == Body_Geo_EMB)
|
||||
public TopoStateFunc(Body body)
|
||||
{
|
||||
state = Astronomy.GeoEmbState(time);
|
||||
}
|
||||
else if (body == Body.Earth)
|
||||
{
|
||||
state = new StateVector(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, time);
|
||||
}
|
||||
else
|
||||
{
|
||||
throw new ArgumentException($"C# TopoStateFunction: unsupported body {body}");
|
||||
this.body = body;
|
||||
}
|
||||
|
||||
state.x -= observer_state.x;
|
||||
state.y -= observer_state.y;
|
||||
state.z -= observer_state.z;
|
||||
state.vx -= observer_state.vx;
|
||||
state.vy -= observer_state.vy;
|
||||
state.vz -= observer_state.vz;
|
||||
public StateVector Eval(AstroTime time)
|
||||
{
|
||||
var observer = new Observer(30.0, -80.0, 1000.0);
|
||||
|
||||
return state;
|
||||
StateVector observer_state = Astronomy.ObserverState(time, observer, EquatorEpoch.J2000);
|
||||
StateVector state;
|
||||
if (body == Body_Geo_EMB)
|
||||
{
|
||||
state = Astronomy.GeoEmbState(time);
|
||||
}
|
||||
else if (body == Body.Earth)
|
||||
{
|
||||
state = new StateVector(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, time);
|
||||
}
|
||||
else
|
||||
{
|
||||
throw new ArgumentException($"C# TopoStateFunction: unsupported body {body}");
|
||||
}
|
||||
|
||||
state.x -= observer_state.x;
|
||||
state.y -= observer_state.y;
|
||||
state.z -= observer_state.z;
|
||||
state.vx -= observer_state.vx;
|
||||
state.vy -= observer_state.vy;
|
||||
state.vz -= observer_state.vz;
|
||||
|
||||
return state;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
static int TopoStateTest()
|
||||
{
|
||||
if (0 != VerifyStateBody(TopoStateFunc, Body.Earth, "../../topostate/Earth_N30_W80_1000m.txt", 2.108e-04, 2.430e-04)) return 1;
|
||||
if (0 != VerifyStateBody(TopoStateFunc, Body_Geo_EMB, "../../topostate/EMB_N30_W80_1000m.txt", 7.195e-04, 2.497e-04)) return 1;
|
||||
if (0 != VerifyStateBody(new TopoStateFunc(Body.Earth), "../../topostate/Earth_N30_W80_1000m.txt", 2.108e-04, 2.430e-04)) return 1;
|
||||
if (0 != VerifyStateBody(new TopoStateFunc(Body_Geo_EMB), "../../topostate/EMB_N30_W80_1000m.txt", 7.195e-04, 2.497e-04)) return 1;
|
||||
Console.WriteLine("C# TopoStateTest: PASS");
|
||||
return 0;
|
||||
}
|
||||
@@ -3298,8 +3336,54 @@ namespace csharp_test
|
||||
}
|
||||
}
|
||||
|
||||
//-----------------------------------------------------------------------------------------
|
||||
|
||||
class LagrangeFunc : IStateVectorFunc
|
||||
{
|
||||
private int point;
|
||||
private Body major_body;
|
||||
private Body minor_body;
|
||||
|
||||
public LagrangeFunc(int point, Body major_body, Body minor_body)
|
||||
{
|
||||
this.point = point;
|
||||
this.major_body = major_body;
|
||||
this.minor_body = minor_body;
|
||||
}
|
||||
|
||||
public StateVector Eval(AstroTime time)
|
||||
{
|
||||
return Astronomy.LagrangePoint(point, time, major_body, minor_body);
|
||||
}
|
||||
}
|
||||
|
||||
static int VerifyStateLagrange(
|
||||
Body major_body,
|
||||
Body minor_body,
|
||||
int point,
|
||||
string filename,
|
||||
double r_thresh,
|
||||
double v_thresh)
|
||||
{
|
||||
var func = new LagrangeFunc(point, major_body, minor_body);
|
||||
return VerifyStateBody(func, filename, r_thresh, v_thresh);
|
||||
}
|
||||
|
||||
static int LagrangeTest()
|
||||
{
|
||||
// Test Sun/EMB Lagrange points.
|
||||
if (0 != VerifyStateLagrange(Body.Sun, Body.EMB, 1, "../../lagrange/semb_L1.txt", 1.33e-5, 6.13e-5)) return 1;
|
||||
if (0 != VerifyStateLagrange(Body.Sun, Body.EMB, 2, "../../lagrange/semb_L2.txt", 1.33e-5, 6.13e-5)) return 1;
|
||||
if (0 != VerifyStateLagrange(Body.Sun, Body.EMB, 4, "../../lagrange/semb_L4.txt", 3.75e-5, 5.28e-5)) return 1;
|
||||
if (0 != VerifyStateLagrange(Body.Sun, Body.EMB, 5, "../../lagrange/semb_L5.txt", 3.75e-5, 5.28e-5)) return 1;
|
||||
|
||||
// Test Earth/Moon Lagrange points.
|
||||
if (0 != VerifyStateLagrange(Body.Earth, Body.Moon, 1, "../../lagrange/em_L1.txt", 3.79e-5, 5.06e-5)) return 1;
|
||||
if (0 != VerifyStateLagrange(Body.Earth, Body.Moon, 2, "../../lagrange/em_L2.txt", 3.79e-5, 5.06e-5)) return 1;
|
||||
if (0 != VerifyStateLagrange(Body.Earth, Body.Moon, 4, "../../lagrange/em_L4.txt", 3.79e-5, 1.59e-3)) return 1;
|
||||
if (0 != VerifyStateLagrange(Body.Earth, Body.Moon, 5, "../../lagrange/em_L5.txt", 3.79e-5, 1.59e-3)) return 1;
|
||||
|
||||
Console.WriteLine("C# LagrangeTest: PASS");
|
||||
return 0; // not yet implemented
|
||||
}
|
||||
|
||||
|
||||
@@ -3480,9 +3480,9 @@ astro_state_vector_t Astronomy_LagrangePoint(
|
||||
*
|
||||
* @param point A value 1..5 that selects which of the Lagrange points to calculate.
|
||||
* @param major_state The state vector of the major (more massive) of the pair of bodies.
|
||||
* @param major_mass The relative mass of the major body.
|
||||
* @param major_mass The mass product GM of the major body.
|
||||
* @param minor_state The state vector of the minor (less massive) of the pair of bodies.
|
||||
* @param minor_mass The relative mass of the minor body.
|
||||
* @param minor_mass The mass product GM of the minor body.
|
||||
* @return The position and velocity of the selected Lagrange point with respect to the major body's center.
|
||||
*/
|
||||
astro_state_vector_t Astronomy_LagrangePointFast(
|
||||
|
||||
@@ -2181,10 +2181,56 @@ $ASTRO_ADDSOL()
|
||||
the products GM[i] are known extremely accurately.
|
||||
*/
|
||||
private const double SUN_GM = 0.2959122082855911e-03;
|
||||
private const double MERCURY_GM = 0.4912547451450812e-10;
|
||||
private const double VENUS_GM = 0.7243452486162703e-09;
|
||||
private const double EARTH_GM = 0.8887692390113509e-09;
|
||||
private const double MARS_GM = 0.9549535105779258e-10;
|
||||
private const double JUPITER_GM = 0.2825345909524226e-06;
|
||||
private const double SATURN_GM = 0.8459715185680659e-07;
|
||||
private const double URANUS_GM = 0.1292024916781969e-07;
|
||||
private const double NEPTUNE_GM = 0.1524358900784276e-07;
|
||||
private const double PLUTO_GM = 0.2188699765425970e-11;
|
||||
|
||||
private const double MOON_GM = EARTH_GM / EARTH_MOON_MASS_RATIO;
|
||||
|
||||
/// <summary>
|
||||
/// Returns the product of mass and universal gravitational constant of a Solar System body.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// For problems involving the gravitational interactions of Solar System bodies,
|
||||
/// it is helpful to know the product GM, where G = the universal gravitational constant
|
||||
/// and M = the mass of the body. In practice, GM is known to a higher precision than
|
||||
/// either G or M alone, and thus using the product results in the most accurate results.
|
||||
/// This function returns the product GM in the units au^3/day^2.
|
||||
/// The values come from page 10 of a
|
||||
/// [JPL memorandum regarding the DE405/LE405 ephemeris](https://web.archive.org/web/20120220062549/http://iau-comm4.jpl.nasa.gov/de405iom/de405iom.pdf).
|
||||
/// </remarks>
|
||||
/// <param name="body">
|
||||
/// The body for which to find the GM product.
|
||||
/// Allowed to be the Sun, Moon, EMB (Earth/Moon Barycenter), or any planet.
|
||||
/// Any other value will cause an exception to be thrown.
|
||||
/// </param>
|
||||
/// <returns>The mass product of the given body in au^3/day^2.</returns>
|
||||
public static double MassProduct(Body body)
|
||||
{
|
||||
switch (body)
|
||||
{
|
||||
case Body.Sun: return SUN_GM;
|
||||
case Body.Mercury: return MERCURY_GM;
|
||||
case Body.Venus: return VENUS_GM;
|
||||
case Body.Earth: return EARTH_GM;
|
||||
case Body.Moon: return MOON_GM;
|
||||
case Body.EMB: return EARTH_GM + MOON_GM;
|
||||
case Body.Mars: return MARS_GM;
|
||||
case Body.Jupiter: return JUPITER_GM;
|
||||
case Body.Saturn: return SATURN_GM;
|
||||
case Body.Uranus: return URANUS_GM;
|
||||
case Body.Neptune: return NEPTUNE_GM;
|
||||
case Body.Pluto: return PLUTO_GM;
|
||||
default:
|
||||
throw new InvalidBodyException(body);
|
||||
}
|
||||
}
|
||||
|
||||
/// <summary>Counter used for performance testing.</summary>
|
||||
public static int CalcMoonCount;
|
||||
@@ -7090,6 +7136,275 @@ $ASTRO_IAU_DATA()
|
||||
throw new Exception("Peak magnitude search failed.");
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Calculates one of the 5 Lagrange points for a pair of co-orbiting bodies.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// Given a more massive "major" body and a much less massive "minor" body,
|
||||
/// calculates one of the five Lagrange points in relation to the minor body's
|
||||
/// orbit around the major body. The parameter `point` is an integer that
|
||||
/// selects the Lagrange point as follows:
|
||||
///
|
||||
/// 1 = the Lagrange point between the major body and minor body.
|
||||
/// 2 = the Lagrange point on the far side of the minor body.
|
||||
/// 3 = the Lagrange point on the far side of the major body.
|
||||
/// 4 = the Lagrange point 60 degrees ahead of the minor body's orbital position.
|
||||
/// 5 = the Lagrange point 60 degrees behind the minor body's orbital position.
|
||||
///
|
||||
/// The function returns the state vector for the selected Lagrange point
|
||||
/// in equatorial J2000 coordinates (EQJ), with respect to the center of the
|
||||
/// major body.
|
||||
///
|
||||
/// To calculate Sun/Earth Lagrange points, pass in `Body.Sun` for `major_body`
|
||||
/// and `Body.EMB` (Earth/Moon barycenter) for `minor_body`.
|
||||
/// For Lagrange points of the Sun and any other planet, pass in just that planet by itself,
|
||||
/// e.g. `Body.Jupiter` for `minor_body`.
|
||||
/// To calculate Earth/Moon Lagrange points, pass in `Body.Earth` and `Body.Moon`
|
||||
/// for the major and minor bodies respectively.
|
||||
///
|
||||
/// In some cases, it may be more efficient to call #Astronomy.LagrangePointFast,
|
||||
/// especially when the state vectors have already been calculated, or are needed
|
||||
/// for some other purpose.
|
||||
/// </remarks>
|
||||
/// <param name="point">A value 1..5 that selects which of the Lagrange points to calculate.</param>
|
||||
/// <param name="time">The time at which the Lagrange point is to be calculated.</param>
|
||||
/// <param name="major_body">The more massive of the co-orbiting bodies: `Body.Sun` or `Body.Earth`.</param>
|
||||
/// <param name="minor_body">The less massive of the co-orbiting bodies. See main remarks.</param>
|
||||
/// <returns></returns>
|
||||
public static StateVector LagrangePoint(
|
||||
int point,
|
||||
AstroTime time,
|
||||
Body major_body,
|
||||
Body minor_body)
|
||||
{
|
||||
double major_mass = MassProduct(major_body);
|
||||
double minor_mass = MassProduct(minor_body);
|
||||
|
||||
StateVector major_state;
|
||||
StateVector minor_state;
|
||||
|
||||
/* Calculate the state vectors for the major and minor bodies. */
|
||||
if (major_body == Body.Earth && minor_body == Body.Moon)
|
||||
{
|
||||
/* Use geocentric calculations for more precision. */
|
||||
|
||||
/* The Earth's geocentric state is trivial. */
|
||||
major_state.t = time;
|
||||
major_state.x = major_state.y = major_state.z = 0.0;
|
||||
major_state.vx = major_state.vy = major_state.vz = 0.0;
|
||||
|
||||
minor_state = GeoMoonState(time);
|
||||
}
|
||||
else
|
||||
{
|
||||
major_state = HelioState(major_body, time);
|
||||
minor_state = HelioState(minor_body, time);
|
||||
}
|
||||
|
||||
return LagrangePointFast(
|
||||
point,
|
||||
major_state,
|
||||
major_mass,
|
||||
minor_state,
|
||||
minor_mass
|
||||
);
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Calculates one of the 5 Lagrange points from body masses and state vectors.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// Given a more massive "major" body and a much less massive "minor" body,
|
||||
/// calculates one of the five Lagrange points in relation to the minor body's
|
||||
/// orbit around the major body. The parameter `point` is an integer that
|
||||
/// selects the Lagrange point as follows:
|
||||
///
|
||||
/// 1 = the Lagrange point between the major body and minor body.
|
||||
/// 2 = the Lagrange point on the far side of the minor body.
|
||||
/// 3 = the Lagrange point on the far side of the major body.
|
||||
/// 4 = the Lagrange point 60 degrees ahead of the minor body's orbital position.
|
||||
/// 5 = the Lagrange point 60 degrees behind the minor body's orbital position.
|
||||
///
|
||||
/// The caller passes in the state vector and mass for both bodies.
|
||||
/// The state vectors can be in any orientation and frame of reference.
|
||||
/// The body masses are expressed as GM products, where G = the universal
|
||||
/// gravitation constant and M = the body's mass. Thus the units for
|
||||
/// `major_mass` and `minor_mass` must be au^3/day^2.
|
||||
/// Use #Astronomy.MassProduct to obtain GM values for various solar system bodies.
|
||||
///
|
||||
/// The function returns the state vector for the selected Lagrange point
|
||||
/// using the same orientation as the state vector parameters `major_state` and `minor_state`,
|
||||
/// and the position and velocity components are with respect to the major body's center.
|
||||
///
|
||||
/// Consider calling #Astronomy.LagrangePoint, instead of this function, for simpler usage in most cases.
|
||||
/// </remarks>
|
||||
/// <param name="point">A value 1..5 that selects which of the Lagrange points to calculate.</param>
|
||||
/// <param name="major_state">The state vector of the major (more massive) of the pair of bodies.</param>
|
||||
/// <param name="major_mass">The mass product GM of the major body.</param>
|
||||
/// <param name="minor_state">The state vector of the minor (less massive) of the pair of bodies.</param>
|
||||
/// <param name="minor_mass">The mass product GM of the minor body.</param>
|
||||
/// <returns>The position and velocity of the selected Lagrange point with respect to the major body's center.</returns>
|
||||
public static StateVector LagrangePointFast(
|
||||
int point,
|
||||
StateVector major_state,
|
||||
double major_mass,
|
||||
StateVector minor_state,
|
||||
double minor_mass)
|
||||
{
|
||||
const double cos_60 = 0.5;
|
||||
const double sin_60 = 0.8660254037844386; /* sqrt(3) / 2 */
|
||||
|
||||
if (point < 1 || point > 5)
|
||||
throw new ArgumentException($"Invalid lagrange point {point}");
|
||||
|
||||
if (double.IsNaN(major_mass) || double.IsInfinity(major_mass) || major_mass <= 0.0)
|
||||
throw new ArgumentException("Major mass must be a positive number.");
|
||||
|
||||
if (double.IsNaN(minor_mass) || double.IsInfinity(minor_mass) || minor_mass <= 0.0)
|
||||
throw new ArgumentException("Minor mass must be a positive number.");
|
||||
|
||||
/* Find the relative position vector <dx, dy, dz>. */
|
||||
double dx = minor_state.x - major_state.x;
|
||||
double dy = minor_state.y - major_state.y;
|
||||
double dz = minor_state.z - major_state.z;
|
||||
double R2 = (dx*dx + dy*dy + dz*dz);
|
||||
|
||||
/* R = Total distance between the bodies. */
|
||||
double R = Math.Sqrt(R2);
|
||||
|
||||
/* Find the velocity vector <vx, vy, vz>. */
|
||||
double vx = minor_state.vx - major_state.vx;
|
||||
double vy = minor_state.vy - major_state.vy;
|
||||
double vz = minor_state.vz - major_state.vz;
|
||||
|
||||
StateVector p;
|
||||
if (point == 4 || point == 5)
|
||||
{
|
||||
/*
|
||||
For L4 and L5, we need to find points 60 degrees away from the
|
||||
line connecting the two bodies and in the instantaneous orbital plane.
|
||||
Define the instantaneous orbital plane as the unique plane that contains
|
||||
both the relative position vector and the relative velocity vector.
|
||||
*/
|
||||
|
||||
/* Take the cross product of position and velocity to find a normal vector <nx, ny, nz>. */
|
||||
double nx = dy*vz - dz*vy;
|
||||
double ny = dz*vx - dx*vz;
|
||||
double nz = dx*vy - dy*vx;
|
||||
|
||||
/* Take the cross product normal*position to get a tangential vector <ux, uy, uz>. */
|
||||
double ux = ny*dz - nz*dy;
|
||||
double uy = nz*dx - nx*dz;
|
||||
double uz = nx*dy - ny*dx;
|
||||
|
||||
/* Convert the tangential direction vector to a unit vector. */
|
||||
double U = Math.Sqrt(ux*ux + uy*uy + uz*uz);
|
||||
ux /= U;
|
||||
uy /= U;
|
||||
uz /= U;
|
||||
|
||||
/* Convert the relative position vector into a unit vector. */
|
||||
dx /= R;
|
||||
dy /= R;
|
||||
dz /= R;
|
||||
|
||||
/* Now we have two perpendicular unit vectors in the orbital plane: 'd' and 'u'. */
|
||||
|
||||
/* Create new unit vectors rotated (+/-)60 degrees from the radius/tangent directions. */
|
||||
double vert = (point == 4) ? +sin_60 : -sin_60;
|
||||
|
||||
/* Rotated radial vector */
|
||||
double Dx = cos_60*dx + vert*ux;
|
||||
double Dy = cos_60*dy + vert*uy;
|
||||
double Dz = cos_60*dz + vert*uz;
|
||||
|
||||
/* Rotated tangent vector */
|
||||
double Ux = cos_60*ux - vert*dx;
|
||||
double Uy = cos_60*uy - vert*dy;
|
||||
double Uz = cos_60*uz - vert*dz;
|
||||
|
||||
/* Calculate L4/L5 positions relative to the major body. */
|
||||
p.x = R * Dx;
|
||||
p.y = R * Dy;
|
||||
p.z = R * Dz;
|
||||
|
||||
/* Use dot products to find radial and tangential components of the relative velocity. */
|
||||
double vrad = vx*dx + vy*dy + vz*dz;
|
||||
double vtan = vx*ux + vy*uy + vz*uz;
|
||||
|
||||
/* Calculate L4/L5 velocities. */
|
||||
p.vx = vrad*Dx + vtan*Ux;
|
||||
p.vy = vrad*Dy + vtan*Uy;
|
||||
p.vz = vrad*Dz + vtan*Uz;
|
||||
}
|
||||
else
|
||||
{
|
||||
/*
|
||||
Calculate the distances of each body from their mutual barycenter.
|
||||
r1 = negative distance of major mass from barycenter (e.g. Sun to the left of barycenter)
|
||||
r2 = positive distance of minor mass from barycenter (e.g. Earth to the right of barycenter)
|
||||
*/
|
||||
double r1 = -R * (minor_mass / (major_mass + minor_mass));
|
||||
double r2 = +R * (major_mass / (major_mass + minor_mass));
|
||||
|
||||
/* Calculate the square of the angular orbital speed in [rad^2 / day^2]. */
|
||||
double omega2 = (major_mass + minor_mass) / (R2*R);
|
||||
|
||||
/*
|
||||
Use Newton's Method to numerically solve for the location where
|
||||
outward centrifugal acceleration in the rotating frame of reference
|
||||
is equal to net inward gravitational acceleration.
|
||||
First derive a good initial guess based on approximate analysis.
|
||||
*/
|
||||
double scale, numer1, numer2;
|
||||
if (point == 1 || point == 2)
|
||||
{
|
||||
scale = (major_mass / (major_mass + minor_mass)) * Math.Cbrt(minor_mass / (3.0 * major_mass));
|
||||
numer1 = -major_mass; /* The major mass is to the left of L1 and L2 */
|
||||
if (point == 1)
|
||||
{
|
||||
scale = 1.0 - scale;
|
||||
numer2 = +minor_mass; /* The minor mass is to the right of L1. */
|
||||
}
|
||||
else
|
||||
{
|
||||
scale = 1.0 + scale;
|
||||
numer2 = -minor_mass; /* The minor mass is to the left of L2. */
|
||||
}
|
||||
}
|
||||
else /* point == 3 */
|
||||
{
|
||||
scale = ((7.0/12.0)*minor_mass - major_mass) / (minor_mass + major_mass);
|
||||
numer1 = +major_mass; /* major mass is to the right of L3. */
|
||||
numer2 = +minor_mass; /* minor mass is to the right of L3. */
|
||||
}
|
||||
|
||||
/* Iterate Newton's Method until it converges. */
|
||||
double x = R*scale - r1;
|
||||
double deltax;
|
||||
do
|
||||
{
|
||||
double dr1 = x - r1;
|
||||
double dr2 = x - r2;
|
||||
double accel = omega2*x + numer1/(dr1*dr1) + numer2/(dr2*dr2);
|
||||
double deriv = omega2 - 2*numer1/(dr1*dr1*dr1) - 2*numer2/(dr2*dr2*dr2);
|
||||
deltax = accel/deriv;
|
||||
x -= deltax;
|
||||
}
|
||||
while (Math.Abs(deltax/R) > 1.0e-14);
|
||||
scale = (x - r1) / R;
|
||||
|
||||
p.x = scale * dx;
|
||||
p.y = scale * dy;
|
||||
p.z = scale * dz;
|
||||
p.vx = scale * vx;
|
||||
p.vy = scale * vy;
|
||||
p.vz = scale * vz;
|
||||
}
|
||||
p.t = major_state.t;
|
||||
return p;
|
||||
}
|
||||
|
||||
/// <summary>Calculates the inverse of a rotation matrix.</summary>
|
||||
/// <remarks>
|
||||
/// Given a rotation matrix that performs some coordinate transform,
|
||||
|
||||
@@ -1047,9 +1047,9 @@ Consider calling [`Astronomy_LagrangePoint`](#Astronomy_LagrangePoint), instead
|
||||
| --- | --- | --- |
|
||||
| `int` | `point` | A value 1..5 that selects which of the Lagrange points to calculate. |
|
||||
| [`astro_state_vector_t`](#astro_state_vector_t) | `major_state` | The state vector of the major (more massive) of the pair of bodies. |
|
||||
| `double` | `major_mass` | The relative mass of the major body. |
|
||||
| `double` | `major_mass` | The mass product GM of the major body. |
|
||||
| [`astro_state_vector_t`](#astro_state_vector_t) | `minor_state` | The state vector of the minor (less massive) of the pair of bodies. |
|
||||
| `double` | `minor_mass` | The relative mass of the minor body. |
|
||||
| `double` | `minor_mass` | The mass product GM of the minor body. |
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -4719,9 +4719,9 @@ astro_state_vector_t Astronomy_LagrangePoint(
|
||||
*
|
||||
* @param point A value 1..5 that selects which of the Lagrange points to calculate.
|
||||
* @param major_state The state vector of the major (more massive) of the pair of bodies.
|
||||
* @param major_mass The relative mass of the major body.
|
||||
* @param major_mass The mass product GM of the major body.
|
||||
* @param minor_state The state vector of the minor (less massive) of the pair of bodies.
|
||||
* @param minor_mass The relative mass of the minor body.
|
||||
* @param minor_mass The mass product GM of the minor body.
|
||||
* @return The position and velocity of the selected Lagrange point with respect to the major body's center.
|
||||
*/
|
||||
astro_state_vector_t Astronomy_LagrangePointFast(
|
||||
|
||||
@@ -755,6 +755,83 @@ to convert to geocentric positions.
|
||||
|
||||
**Returns:** Position and velocity vectors of Jupiter's largest 4 moons.
|
||||
|
||||
<a name="Astronomy.LagrangePoint"></a>
|
||||
### Astronomy.LagrangePoint(point, time, major_body, minor_body) ⇒ [`StateVector`](#StateVector)
|
||||
|
||||
**Calculates one of the 5 Lagrange points for a pair of co-orbiting bodies.**
|
||||
|
||||
Given a more massive "major" body and a much less massive "minor" body,
|
||||
calculates one of the five Lagrange points in relation to the minor body's
|
||||
orbit around the major body. The parameter `point` is an integer that
|
||||
selects the Lagrange point as follows:
|
||||
|
||||
1 = the Lagrange point between the major body and minor body.
|
||||
2 = the Lagrange point on the far side of the minor body.
|
||||
3 = the Lagrange point on the far side of the major body.
|
||||
4 = the Lagrange point 60 degrees ahead of the minor body's orbital position.
|
||||
5 = the Lagrange point 60 degrees behind the minor body's orbital position.
|
||||
|
||||
The function returns the state vector for the selected Lagrange point
|
||||
in equatorial J2000 coordinates (EQJ), with respect to the center of the
|
||||
major body.
|
||||
|
||||
To calculate Sun/Earth Lagrange points, pass in `Body.Sun` for `major_body`
|
||||
and `Body.EMB` (Earth/Moon barycenter) for `minor_body`.
|
||||
For Lagrange points of the Sun and any other planet, pass in just that planet by itself,
|
||||
e.g. `Body.Jupiter` for `minor_body`.
|
||||
To calculate Earth/Moon Lagrange points, pass in `Body.Earth` and `Body.Moon`
|
||||
for the major and minor bodies respectively.
|
||||
|
||||
In some cases, it may be more efficient to call [`Astronomy.LagrangePointFast`](#Astronomy.LagrangePointFast),
|
||||
especially when the state vectors have already been calculated, or are needed
|
||||
for some other purpose.
|
||||
|
||||
| Type | Parameter | Description |
|
||||
| --- | --- | --- |
|
||||
| `int` | `point` | A value 1..5 that selects which of the Lagrange points to calculate. |
|
||||
| [`AstroTime`](#AstroTime) | `time` | The time at which the Lagrange point is to be calculated. |
|
||||
| [`Body`](#Body) | `major_body` | The more massive of the co-orbiting bodies: `Body.Sun` or `Body.Earth`. |
|
||||
| [`Body`](#Body) | `minor_body` | The less massive of the co-orbiting bodies. See main remarks. |
|
||||
|
||||
<a name="Astronomy.LagrangePointFast"></a>
|
||||
### Astronomy.LagrangePointFast(point, major_state, major_mass, minor_state, minor_mass) ⇒ [`StateVector`](#StateVector)
|
||||
|
||||
**Calculates one of the 5 Lagrange points from body masses and state vectors.**
|
||||
|
||||
Given a more massive "major" body and a much less massive "minor" body,
|
||||
calculates one of the five Lagrange points in relation to the minor body's
|
||||
orbit around the major body. The parameter `point` is an integer that
|
||||
selects the Lagrange point as follows:
|
||||
|
||||
1 = the Lagrange point between the major body and minor body.
|
||||
2 = the Lagrange point on the far side of the minor body.
|
||||
3 = the Lagrange point on the far side of the major body.
|
||||
4 = the Lagrange point 60 degrees ahead of the minor body's orbital position.
|
||||
5 = the Lagrange point 60 degrees behind the minor body's orbital position.
|
||||
|
||||
The caller passes in the state vector and mass for both bodies.
|
||||
The state vectors can be in any orientation and frame of reference.
|
||||
The body masses are expressed as GM products, where G = the universal
|
||||
gravitation constant and M = the body's mass. Thus the units for
|
||||
`major_mass` and `minor_mass` must be au^3/day^2.
|
||||
Use [`Astronomy.MassProduct`](#Astronomy.MassProduct) to obtain GM values for various solar system bodies.
|
||||
|
||||
The function returns the state vector for the selected Lagrange point
|
||||
using the same orientation as the state vector parameters `major_state` and `minor_state`,
|
||||
and the position and velocity components are with respect to the major body's center.
|
||||
|
||||
Consider calling [`Astronomy.LagrangePoint`](#Astronomy.LagrangePoint), instead of this function, for simpler usage in most cases.
|
||||
|
||||
| Type | Parameter | Description |
|
||||
| --- | --- | --- |
|
||||
| `int` | `point` | A value 1..5 that selects which of the Lagrange points to calculate. |
|
||||
| [`StateVector`](#StateVector) | `major_state` | The state vector of the major (more massive) of the pair of bodies. |
|
||||
| `double` | `major_mass` | The mass product GM of the major body. |
|
||||
| [`StateVector`](#StateVector) | `minor_state` | The state vector of the minor (less massive) of the pair of bodies. |
|
||||
| `double` | `minor_mass` | The mass product GM of the minor body. |
|
||||
|
||||
**Returns:** The position and velocity of the selected Lagrange point with respect to the major body's center.
|
||||
|
||||
<a name="Astronomy.Libration"></a>
|
||||
### Astronomy.Libration(time) ⇒ [`LibrationInfo`](#LibrationInfo)
|
||||
|
||||
@@ -780,6 +857,25 @@ and the apparent angular diameter of the Moon `diam_deg`.
|
||||
|
||||
**Returns:** The Moon's ecliptic position and libration angles as seen from the Earth.
|
||||
|
||||
<a name="Astronomy.MassProduct"></a>
|
||||
### Astronomy.MassProduct(body) ⇒ `double`
|
||||
|
||||
**Returns the product of mass and universal gravitational constant of a Solar System body.**
|
||||
|
||||
For problems involving the gravitational interactions of Solar System bodies,
|
||||
it is helpful to know the product GM, where G = the universal gravitational constant
|
||||
and M = the mass of the body. In practice, GM is known to a higher precision than
|
||||
either G or M alone, and thus using the product results in the most accurate results.
|
||||
This function returns the product GM in the units au^3/day^2.
|
||||
The values come from page 10 of a
|
||||
[JPL memorandum regarding the DE405/LE405 ephemeris](https://web.archive.org/web/20120220062549/http://iau-comm4.jpl.nasa.gov/de405iom/de405iom.pdf).
|
||||
|
||||
| Type | Parameter | Description |
|
||||
| --- | --- | --- |
|
||||
| [`Body`](#Body) | `body` | The body for which to find the GM product. Allowed to be the Sun, Moon, EMB (Earth/Moon Barycenter), or any planet. Any other value will cause an exception to be thrown. |
|
||||
|
||||
**Returns:** The mass product of the given body in au^3/day^2.
|
||||
|
||||
<a name="Astronomy.MoonPhase"></a>
|
||||
### Astronomy.MoonPhase(time) ⇒ `double`
|
||||
|
||||
|
||||
@@ -2286,10 +2286,56 @@ namespace CosineKitty
|
||||
the products GM[i] are known extremely accurately.
|
||||
*/
|
||||
private const double SUN_GM = 0.2959122082855911e-03;
|
||||
private const double MERCURY_GM = 0.4912547451450812e-10;
|
||||
private const double VENUS_GM = 0.7243452486162703e-09;
|
||||
private const double EARTH_GM = 0.8887692390113509e-09;
|
||||
private const double MARS_GM = 0.9549535105779258e-10;
|
||||
private const double JUPITER_GM = 0.2825345909524226e-06;
|
||||
private const double SATURN_GM = 0.8459715185680659e-07;
|
||||
private const double URANUS_GM = 0.1292024916781969e-07;
|
||||
private const double NEPTUNE_GM = 0.1524358900784276e-07;
|
||||
private const double PLUTO_GM = 0.2188699765425970e-11;
|
||||
|
||||
private const double MOON_GM = EARTH_GM / EARTH_MOON_MASS_RATIO;
|
||||
|
||||
/// <summary>
|
||||
/// Returns the product of mass and universal gravitational constant of a Solar System body.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// For problems involving the gravitational interactions of Solar System bodies,
|
||||
/// it is helpful to know the product GM, where G = the universal gravitational constant
|
||||
/// and M = the mass of the body. In practice, GM is known to a higher precision than
|
||||
/// either G or M alone, and thus using the product results in the most accurate results.
|
||||
/// This function returns the product GM in the units au^3/day^2.
|
||||
/// The values come from page 10 of a
|
||||
/// [JPL memorandum regarding the DE405/LE405 ephemeris](https://web.archive.org/web/20120220062549/http://iau-comm4.jpl.nasa.gov/de405iom/de405iom.pdf).
|
||||
/// </remarks>
|
||||
/// <param name="body">
|
||||
/// The body for which to find the GM product.
|
||||
/// Allowed to be the Sun, Moon, EMB (Earth/Moon Barycenter), or any planet.
|
||||
/// Any other value will cause an exception to be thrown.
|
||||
/// </param>
|
||||
/// <returns>The mass product of the given body in au^3/day^2.</returns>
|
||||
public static double MassProduct(Body body)
|
||||
{
|
||||
switch (body)
|
||||
{
|
||||
case Body.Sun: return SUN_GM;
|
||||
case Body.Mercury: return MERCURY_GM;
|
||||
case Body.Venus: return VENUS_GM;
|
||||
case Body.Earth: return EARTH_GM;
|
||||
case Body.Moon: return MOON_GM;
|
||||
case Body.EMB: return EARTH_GM + MOON_GM;
|
||||
case Body.Mars: return MARS_GM;
|
||||
case Body.Jupiter: return JUPITER_GM;
|
||||
case Body.Saturn: return SATURN_GM;
|
||||
case Body.Uranus: return URANUS_GM;
|
||||
case Body.Neptune: return NEPTUNE_GM;
|
||||
case Body.Pluto: return PLUTO_GM;
|
||||
default:
|
||||
throw new InvalidBodyException(body);
|
||||
}
|
||||
}
|
||||
|
||||
/// <summary>Counter used for performance testing.</summary>
|
||||
public static int CalcMoonCount;
|
||||
@@ -8302,6 +8348,275 @@ namespace CosineKitty
|
||||
throw new Exception("Peak magnitude search failed.");
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Calculates one of the 5 Lagrange points for a pair of co-orbiting bodies.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// Given a more massive "major" body and a much less massive "minor" body,
|
||||
/// calculates one of the five Lagrange points in relation to the minor body's
|
||||
/// orbit around the major body. The parameter `point` is an integer that
|
||||
/// selects the Lagrange point as follows:
|
||||
///
|
||||
/// 1 = the Lagrange point between the major body and minor body.
|
||||
/// 2 = the Lagrange point on the far side of the minor body.
|
||||
/// 3 = the Lagrange point on the far side of the major body.
|
||||
/// 4 = the Lagrange point 60 degrees ahead of the minor body's orbital position.
|
||||
/// 5 = the Lagrange point 60 degrees behind the minor body's orbital position.
|
||||
///
|
||||
/// The function returns the state vector for the selected Lagrange point
|
||||
/// in equatorial J2000 coordinates (EQJ), with respect to the center of the
|
||||
/// major body.
|
||||
///
|
||||
/// To calculate Sun/Earth Lagrange points, pass in `Body.Sun` for `major_body`
|
||||
/// and `Body.EMB` (Earth/Moon barycenter) for `minor_body`.
|
||||
/// For Lagrange points of the Sun and any other planet, pass in just that planet by itself,
|
||||
/// e.g. `Body.Jupiter` for `minor_body`.
|
||||
/// To calculate Earth/Moon Lagrange points, pass in `Body.Earth` and `Body.Moon`
|
||||
/// for the major and minor bodies respectively.
|
||||
///
|
||||
/// In some cases, it may be more efficient to call #Astronomy.LagrangePointFast,
|
||||
/// especially when the state vectors have already been calculated, or are needed
|
||||
/// for some other purpose.
|
||||
/// </remarks>
|
||||
/// <param name="point">A value 1..5 that selects which of the Lagrange points to calculate.</param>
|
||||
/// <param name="time">The time at which the Lagrange point is to be calculated.</param>
|
||||
/// <param name="major_body">The more massive of the co-orbiting bodies: `Body.Sun` or `Body.Earth`.</param>
|
||||
/// <param name="minor_body">The less massive of the co-orbiting bodies. See main remarks.</param>
|
||||
/// <returns></returns>
|
||||
public static StateVector LagrangePoint(
|
||||
int point,
|
||||
AstroTime time,
|
||||
Body major_body,
|
||||
Body minor_body)
|
||||
{
|
||||
double major_mass = MassProduct(major_body);
|
||||
double minor_mass = MassProduct(minor_body);
|
||||
|
||||
StateVector major_state;
|
||||
StateVector minor_state;
|
||||
|
||||
/* Calculate the state vectors for the major and minor bodies. */
|
||||
if (major_body == Body.Earth && minor_body == Body.Moon)
|
||||
{
|
||||
/* Use geocentric calculations for more precision. */
|
||||
|
||||
/* The Earth's geocentric state is trivial. */
|
||||
major_state.t = time;
|
||||
major_state.x = major_state.y = major_state.z = 0.0;
|
||||
major_state.vx = major_state.vy = major_state.vz = 0.0;
|
||||
|
||||
minor_state = GeoMoonState(time);
|
||||
}
|
||||
else
|
||||
{
|
||||
major_state = HelioState(major_body, time);
|
||||
minor_state = HelioState(minor_body, time);
|
||||
}
|
||||
|
||||
return LagrangePointFast(
|
||||
point,
|
||||
major_state,
|
||||
major_mass,
|
||||
minor_state,
|
||||
minor_mass
|
||||
);
|
||||
}
|
||||
|
||||
/// <summary>
|
||||
/// Calculates one of the 5 Lagrange points from body masses and state vectors.
|
||||
/// </summary>
|
||||
/// <remarks>
|
||||
/// Given a more massive "major" body and a much less massive "minor" body,
|
||||
/// calculates one of the five Lagrange points in relation to the minor body's
|
||||
/// orbit around the major body. The parameter `point` is an integer that
|
||||
/// selects the Lagrange point as follows:
|
||||
///
|
||||
/// 1 = the Lagrange point between the major body and minor body.
|
||||
/// 2 = the Lagrange point on the far side of the minor body.
|
||||
/// 3 = the Lagrange point on the far side of the major body.
|
||||
/// 4 = the Lagrange point 60 degrees ahead of the minor body's orbital position.
|
||||
/// 5 = the Lagrange point 60 degrees behind the minor body's orbital position.
|
||||
///
|
||||
/// The caller passes in the state vector and mass for both bodies.
|
||||
/// The state vectors can be in any orientation and frame of reference.
|
||||
/// The body masses are expressed as GM products, where G = the universal
|
||||
/// gravitation constant and M = the body's mass. Thus the units for
|
||||
/// `major_mass` and `minor_mass` must be au^3/day^2.
|
||||
/// Use #Astronomy.MassProduct to obtain GM values for various solar system bodies.
|
||||
///
|
||||
/// The function returns the state vector for the selected Lagrange point
|
||||
/// using the same orientation as the state vector parameters `major_state` and `minor_state`,
|
||||
/// and the position and velocity components are with respect to the major body's center.
|
||||
///
|
||||
/// Consider calling #Astronomy.LagrangePoint, instead of this function, for simpler usage in most cases.
|
||||
/// </remarks>
|
||||
/// <param name="point">A value 1..5 that selects which of the Lagrange points to calculate.</param>
|
||||
/// <param name="major_state">The state vector of the major (more massive) of the pair of bodies.</param>
|
||||
/// <param name="major_mass">The mass product GM of the major body.</param>
|
||||
/// <param name="minor_state">The state vector of the minor (less massive) of the pair of bodies.</param>
|
||||
/// <param name="minor_mass">The mass product GM of the minor body.</param>
|
||||
/// <returns>The position and velocity of the selected Lagrange point with respect to the major body's center.</returns>
|
||||
public static StateVector LagrangePointFast(
|
||||
int point,
|
||||
StateVector major_state,
|
||||
double major_mass,
|
||||
StateVector minor_state,
|
||||
double minor_mass)
|
||||
{
|
||||
const double cos_60 = 0.5;
|
||||
const double sin_60 = 0.8660254037844386; /* sqrt(3) / 2 */
|
||||
|
||||
if (point < 1 || point > 5)
|
||||
throw new ArgumentException($"Invalid lagrange point {point}");
|
||||
|
||||
if (double.IsNaN(major_mass) || double.IsInfinity(major_mass) || major_mass <= 0.0)
|
||||
throw new ArgumentException("Major mass must be a positive number.");
|
||||
|
||||
if (double.IsNaN(minor_mass) || double.IsInfinity(minor_mass) || minor_mass <= 0.0)
|
||||
throw new ArgumentException("Minor mass must be a positive number.");
|
||||
|
||||
/* Find the relative position vector <dx, dy, dz>. */
|
||||
double dx = minor_state.x - major_state.x;
|
||||
double dy = minor_state.y - major_state.y;
|
||||
double dz = minor_state.z - major_state.z;
|
||||
double R2 = (dx*dx + dy*dy + dz*dz);
|
||||
|
||||
/* R = Total distance between the bodies. */
|
||||
double R = Math.Sqrt(R2);
|
||||
|
||||
/* Find the velocity vector <vx, vy, vz>. */
|
||||
double vx = minor_state.vx - major_state.vx;
|
||||
double vy = minor_state.vy - major_state.vy;
|
||||
double vz = minor_state.vz - major_state.vz;
|
||||
|
||||
StateVector p;
|
||||
if (point == 4 || point == 5)
|
||||
{
|
||||
/*
|
||||
For L4 and L5, we need to find points 60 degrees away from the
|
||||
line connecting the two bodies and in the instantaneous orbital plane.
|
||||
Define the instantaneous orbital plane as the unique plane that contains
|
||||
both the relative position vector and the relative velocity vector.
|
||||
*/
|
||||
|
||||
/* Take the cross product of position and velocity to find a normal vector <nx, ny, nz>. */
|
||||
double nx = dy*vz - dz*vy;
|
||||
double ny = dz*vx - dx*vz;
|
||||
double nz = dx*vy - dy*vx;
|
||||
|
||||
/* Take the cross product normal*position to get a tangential vector <ux, uy, uz>. */
|
||||
double ux = ny*dz - nz*dy;
|
||||
double uy = nz*dx - nx*dz;
|
||||
double uz = nx*dy - ny*dx;
|
||||
|
||||
/* Convert the tangential direction vector to a unit vector. */
|
||||
double U = Math.Sqrt(ux*ux + uy*uy + uz*uz);
|
||||
ux /= U;
|
||||
uy /= U;
|
||||
uz /= U;
|
||||
|
||||
/* Convert the relative position vector into a unit vector. */
|
||||
dx /= R;
|
||||
dy /= R;
|
||||
dz /= R;
|
||||
|
||||
/* Now we have two perpendicular unit vectors in the orbital plane: 'd' and 'u'. */
|
||||
|
||||
/* Create new unit vectors rotated (+/-)60 degrees from the radius/tangent directions. */
|
||||
double vert = (point == 4) ? +sin_60 : -sin_60;
|
||||
|
||||
/* Rotated radial vector */
|
||||
double Dx = cos_60*dx + vert*ux;
|
||||
double Dy = cos_60*dy + vert*uy;
|
||||
double Dz = cos_60*dz + vert*uz;
|
||||
|
||||
/* Rotated tangent vector */
|
||||
double Ux = cos_60*ux - vert*dx;
|
||||
double Uy = cos_60*uy - vert*dy;
|
||||
double Uz = cos_60*uz - vert*dz;
|
||||
|
||||
/* Calculate L4/L5 positions relative to the major body. */
|
||||
p.x = R * Dx;
|
||||
p.y = R * Dy;
|
||||
p.z = R * Dz;
|
||||
|
||||
/* Use dot products to find radial and tangential components of the relative velocity. */
|
||||
double vrad = vx*dx + vy*dy + vz*dz;
|
||||
double vtan = vx*ux + vy*uy + vz*uz;
|
||||
|
||||
/* Calculate L4/L5 velocities. */
|
||||
p.vx = vrad*Dx + vtan*Ux;
|
||||
p.vy = vrad*Dy + vtan*Uy;
|
||||
p.vz = vrad*Dz + vtan*Uz;
|
||||
}
|
||||
else
|
||||
{
|
||||
/*
|
||||
Calculate the distances of each body from their mutual barycenter.
|
||||
r1 = negative distance of major mass from barycenter (e.g. Sun to the left of barycenter)
|
||||
r2 = positive distance of minor mass from barycenter (e.g. Earth to the right of barycenter)
|
||||
*/
|
||||
double r1 = -R * (minor_mass / (major_mass + minor_mass));
|
||||
double r2 = +R * (major_mass / (major_mass + minor_mass));
|
||||
|
||||
/* Calculate the square of the angular orbital speed in [rad^2 / day^2]. */
|
||||
double omega2 = (major_mass + minor_mass) / (R2*R);
|
||||
|
||||
/*
|
||||
Use Newton's Method to numerically solve for the location where
|
||||
outward centrifugal acceleration in the rotating frame of reference
|
||||
is equal to net inward gravitational acceleration.
|
||||
First derive a good initial guess based on approximate analysis.
|
||||
*/
|
||||
double scale, numer1, numer2;
|
||||
if (point == 1 || point == 2)
|
||||
{
|
||||
scale = (major_mass / (major_mass + minor_mass)) * Math.Cbrt(minor_mass / (3.0 * major_mass));
|
||||
numer1 = -major_mass; /* The major mass is to the left of L1 and L2 */
|
||||
if (point == 1)
|
||||
{
|
||||
scale = 1.0 - scale;
|
||||
numer2 = +minor_mass; /* The minor mass is to the right of L1. */
|
||||
}
|
||||
else
|
||||
{
|
||||
scale = 1.0 + scale;
|
||||
numer2 = -minor_mass; /* The minor mass is to the left of L2. */
|
||||
}
|
||||
}
|
||||
else /* point == 3 */
|
||||
{
|
||||
scale = ((7.0/12.0)*minor_mass - major_mass) / (minor_mass + major_mass);
|
||||
numer1 = +major_mass; /* major mass is to the right of L3. */
|
||||
numer2 = +minor_mass; /* minor mass is to the right of L3. */
|
||||
}
|
||||
|
||||
/* Iterate Newton's Method until it converges. */
|
||||
double x = R*scale - r1;
|
||||
double deltax;
|
||||
do
|
||||
{
|
||||
double dr1 = x - r1;
|
||||
double dr2 = x - r2;
|
||||
double accel = omega2*x + numer1/(dr1*dr1) + numer2/(dr2*dr2);
|
||||
double deriv = omega2 - 2*numer1/(dr1*dr1*dr1) - 2*numer2/(dr2*dr2*dr2);
|
||||
deltax = accel/deriv;
|
||||
x -= deltax;
|
||||
}
|
||||
while (Math.Abs(deltax/R) > 1.0e-14);
|
||||
scale = (x - r1) / R;
|
||||
|
||||
p.x = scale * dx;
|
||||
p.y = scale * dy;
|
||||
p.z = scale * dz;
|
||||
p.vx = scale * vx;
|
||||
p.vy = scale * vy;
|
||||
p.vz = scale * vz;
|
||||
}
|
||||
p.t = major_state.t;
|
||||
return p;
|
||||
}
|
||||
|
||||
/// <summary>Calculates the inverse of a rotation matrix.</summary>
|
||||
/// <remarks>
|
||||
/// Given a rotation matrix that performs some coordinate transform,
|
||||
|
||||
Reference in New Issue
Block a user