diff --git a/generate/dotnet/csharp_test/csharp_test.cs b/generate/dotnet/csharp_test/csharp_test.cs index e148a8fc..27acabd4 100644 --- a/generate/dotnet/csharp_test/csharp_test.cs +++ b/generate/dotnet/csharp_test/csharp_test.cs @@ -2720,11 +2720,15 @@ namespace csharp_test return sqrt(diff_squared); } + interface IStateVectorFunc + { + StateVector Eval(AstroTime time); + } + static int VerifyState( - Func 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 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 } diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index ffa52646..88b4a9a7 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -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( diff --git a/generate/template/astronomy.cs b/generate/template/astronomy.cs index 34d15e03..db16d54c 100644 --- a/generate/template/astronomy.cs +++ b/generate/template/astronomy.cs @@ -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; + + /// + /// 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). + /// + /// + /// 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. + /// + /// The mass product of the given body in au^3/day^2. + 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); + } + } /// Counter used for performance testing. public static int CalcMoonCount; @@ -7090,6 +7136,275 @@ $ASTRO_IAU_DATA() throw new Exception("Peak magnitude search failed."); } + /// + /// 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, + /// especially when the state vectors have already been calculated, or are needed + /// for some other purpose. + /// + /// A value 1..5 that selects which of the Lagrange points to calculate. + /// The time at which the Lagrange point is to be calculated. + /// The more massive of the co-orbiting bodies: `Body.Sun` or `Body.Earth`. + /// The less massive of the co-orbiting bodies. See main remarks. + /// + 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 + ); + } + + /// + /// 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 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. + /// + /// A value 1..5 that selects which of the Lagrange points to calculate. + /// The state vector of the major (more massive) of the pair of bodies. + /// The mass product GM of the major body. + /// The state vector of the minor (less massive) of the pair of bodies. + /// The mass product GM of the minor body. + /// The position and velocity of the selected Lagrange point with respect to the major body's center. + 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 . */ + 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 . */ + 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 . */ + 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 . */ + 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; + } + /// Calculates the inverse of a rotation matrix. /// /// Given a rotation matrix that performs some coordinate transform, diff --git a/source/c/README.md b/source/c/README.md index f97767f0..19eabd48 100644 --- a/source/c/README.md +++ b/source/c/README.md @@ -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. | diff --git a/source/c/astronomy.c b/source/c/astronomy.c index feab020b..663ab606 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -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( diff --git a/source/csharp/README.md b/source/csharp/README.md index 17ee53c7..99b93b3d 100644 --- a/source/csharp/README.md +++ b/source/csharp/README.md @@ -755,6 +755,83 @@ to convert to geocentric positions. **Returns:** Position and velocity vectors of Jupiter's largest 4 moons. + +### 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. | + + +### 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. + ### 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. + +### 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. + ### Astronomy.MoonPhase(time) ⇒ `double` diff --git a/source/csharp/astronomy.cs b/source/csharp/astronomy.cs index 360c5115..94b0e61b 100644 --- a/source/csharp/astronomy.cs +++ b/source/csharp/astronomy.cs @@ -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; + + /// + /// 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). + /// + /// + /// 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. + /// + /// The mass product of the given body in au^3/day^2. + 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); + } + } /// Counter used for performance testing. public static int CalcMoonCount; @@ -8302,6 +8348,275 @@ namespace CosineKitty throw new Exception("Peak magnitude search failed."); } + /// + /// 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, + /// especially when the state vectors have already been calculated, or are needed + /// for some other purpose. + /// + /// A value 1..5 that selects which of the Lagrange points to calculate. + /// The time at which the Lagrange point is to be calculated. + /// The more massive of the co-orbiting bodies: `Body.Sun` or `Body.Earth`. + /// The less massive of the co-orbiting bodies. See main remarks. + /// + 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 + ); + } + + /// + /// 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 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. + /// + /// A value 1..5 that selects which of the Lagrange points to calculate. + /// The state vector of the major (more massive) of the pair of bodies. + /// The mass product GM of the major body. + /// The state vector of the minor (less massive) of the pair of bodies. + /// The mass product GM of the minor body. + /// The position and velocity of the selected Lagrange point with respect to the major body's center. + 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 . */ + 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 . */ + 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 . */ + 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 . */ + 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; + } + /// Calculates the inverse of a rotation matrix. /// /// Given a rotation matrix that performs some coordinate transform,