diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py index fa815f4f..92d25c2e 100644 --- a/demo/python/astronomy.py +++ b/demo/python/astronomy.py @@ -36,6 +36,11 @@ import datetime import enum import re +def _cbrt(x): + if x < 0.0: + return -((-x) ** (1.0 / 3.0)) + return x ** (1.0 / 3.0) + KM_PER_AU = 1.4959787069098932e+8 # The number of kilometers per astronomical unit. C_AUDAY = 173.1446326846693 # The speed of light expressed in astronomical units per day. @@ -105,10 +110,55 @@ _EARTH_MOON_MASS_RATIO = 81.30056 # the products GM[i] are known extremely accurately. # _SUN_GM = 0.2959122082855911e-03 +_MERCURY_GM = 0.4912547451450812e-10 +_VENUS_GM = 0.7243452486162703e-09 +_EARTH_GM = 0.8887692390113509e-09 +_MARS_GM = 0.9549535105779258e-10 _JUPITER_GM = 0.2825345909524226e-06 _SATURN_GM = 0.8459715185680659e-07 _URANUS_GM = 0.1292024916781969e-07 _NEPTUNE_GM = 0.1524358900784276e-07 +_PLUTO_GM = 0.2188699765425970e-11 + +_MOON_GM = _EARTH_GM / _EARTH_MOON_MASS_RATIO + + +def MassProduct(body): + """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). + + Parameters + ---------- + 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 + ------- + float + The mass product of the given body in au^3/day^2. + """ + if body == Body.Sun: return _SUN_GM + if body == Body.Mercury: return _MERCURY_GM + if body == Body.Venus: return _VENUS_GM + if body == Body.Earth: return _EARTH_GM + if body == Body.Moon: return _MOON_GM + if body == Body.EMB: return _EARTH_GM + _MOON_GM + if body == Body.Mars: return _MARS_GM + if body == Body.Jupiter: return _JUPITER_GM + if body == Body.Saturn: return _SATURN_GM + if body == Body.Uranus: return _URANUS_GM + if body == Body.Neptune: return _NEPTUNE_GM + if body == Body.Pluto: return _PLUTO_GM + raise InvalidBodyError() @enum.unique class _PrecessDir(enum.Enum): @@ -9193,3 +9243,244 @@ def RotationAxis(body, time): time ) return AxisInfo(ra/15.0, dec, w, north) + + +def LagrangePoint(point, time, major_body, minor_body): + """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 + (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 #LagrangePointFast, + especially when the state vectors have already been calculated, or are needed + for some other purpose. + + Parameters + ---------- + point : int + An integer 1..5 that selects which of the Lagrange points to calculate. + time : Time + The time for which the Lagrange point is to be calculated. + major_body : Body + The more massive of the co-orbiting bodies: `Body.Sun` or `Body.Earth`. + minor_body : Body + The less massive of the co-orbiting bodies. See main remarks. + + Returns + ------- + StateVector + The position and velocity of the selected Lagrange point with respect to the major body's center. + """ + major_mass = MassProduct(major_body) + minor_mass = MassProduct(minor_body) + + # Calculate the state vectors for the major and minor bodies. + if major_body == Body.Earth and minor_body == Body.Moon: + # Use geocentric calculations for more precision. + major_state = StateVector(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, time) + 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 + ) + + +def 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 #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 #LagrangePoint, instead of this function, for simpler usage in most cases. + + Parameters + ---------- + point : int + An integer 1..5 that selects which of the Lagrange points to calculate. + major_state : StateVector + The state vector of the major (more massive) of the pair of bodies. + major_mass : float + The mass product GM of the major body. + minor_state : StateVector + The state vector of the minor (less massive) of the pair of bodies. + minor_mass : float + The mass product GM of the minor body. + + Returns + ------- + StateVector + The position and velocity of the selected Lagrange point with respect to the major body's center. + """ + cos_60 = 0.5 + sin_60 = 0.8660254037844386 # sqrt(3) / 2 + + if point not in [1, 2, 3, 4, 5]: + raise Error('Invalid Lagrange point [{}]. Must be an integer 1..5.'.format(point)) + + if major_mass <= 0.0: + raise Error('Major mass must be a positive number.') + + if minor_mass <= 0.0: + raise Error('Minor mass must be a positive number.') + + # Find the relative position vector . + dx = minor_state.x - major_state.x + dy = minor_state.y - major_state.y + dz = minor_state.z - major_state.z + R2 = (dx*dx + dy*dy + dz*dz) + + # R = Total distance between the bodies. + R = math.sqrt(R2) + + # Find the velocity vector . + vx = minor_state.vx - major_state.vx + vy = minor_state.vy - major_state.vy + vz = minor_state.vz - major_state.vz + + if point in [4, 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 = dy*vz - dz*vy + ny = dz*vx - dx*vz + nz = dx*vy - dy*vx + + # Take the cross product normal*position to get a tangential vector . + ux = ny*dz - nz*dy + uy = nz*dx - nx*dz + uz = nx*dy - ny*dx + + # Convert the tangential direction vector to a unit vector. + 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. + if point == 4: + vert = +sin_60 + else: + vert = -sin_60 + + # Rotated radial vector + Dx = cos_60*dx + vert*ux + Dy = cos_60*dy + vert*uy + Dz = cos_60*dz + vert*uz + + # Rotated tangent vector + Ux = cos_60*ux - vert*dx + Uy = cos_60*uy - vert*dy + Uz = cos_60*uz - vert*dz + + # Calculate L4/L5 positions relative to the major body. + px = R * Dx + py = R * Dy + pz = R * Dz + + # Use dot products to find radial and tangential components of the relative velocity. + vrad = vx*dx + vy*dy + vz*dz + vtan = vx*ux + vy*uy + vz*uz + + # Calculate L4/L5 velocities. + pvx = vrad*Dx + vtan*Ux + pvy = vrad*Dy + vtan*Uy + pvz = vrad*Dz + vtan*Uz + + return StateVector(px, py, pz, pvx, pvy, pvz, major_state.t) + + # Handle Langrange points 1, 2, 3. + # 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) + r1 = -R * (minor_mass / (major_mass + minor_mass)) + r2 = +R * (major_mass / (major_mass + minor_mass)) + + # Calculate the square of the angular orbital speed in [rad^2 / day^2]. + 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. + if point in [1, 2]: + scale = (major_mass / (major_mass + minor_mass)) * _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. + x = R*scale - r1 + while True: + dr1 = x - r1 + dr2 = x - r2 + accel = omega2*x + numer1/(dr1*dr1) + numer2/(dr2*dr2) + deriv = omega2 - 2*numer1/(dr1*dr1*dr1) - 2*numer2/(dr2*dr2*dr2) + deltax = accel/deriv + x -= deltax + if abs(deltax/R) <= 1.0e-14: + break + scale = (x - r1) / R + return StateVector(scale*dx, scale*dy, scale*dz, scale*vx, scale*vy, scale*vz, major_state.t) diff --git a/generate/template/astronomy.cs b/generate/template/astronomy.cs index 6c4570b5..35f7d2c2 100644 --- a/generate/template/astronomy.cs +++ b/generate/template/astronomy.cs @@ -7166,8 +7166,8 @@ $ASTRO_IAU_DATA() /// 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. + /// An integer 1..5 that selects which of the Lagrange points to calculate. + /// The time for 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. /// The position and velocity of the selected Lagrange point with respect to the major body's center. @@ -7238,7 +7238,7 @@ $ASTRO_IAU_DATA() /// /// 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. + /// An integer 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. diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index 6a796a41..3f9aa4ab 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -36,6 +36,11 @@ import datetime import enum import re +def _cbrt(x): + if x < 0.0: + return -((-x) ** (1.0 / 3.0)) + return x ** (1.0 / 3.0) + KM_PER_AU = 1.4959787069098932e+8 # The number of kilometers per astronomical unit. C_AUDAY = 173.1446326846693 # The speed of light expressed in astronomical units per day. @@ -105,10 +110,55 @@ _EARTH_MOON_MASS_RATIO = 81.30056 # the products GM[i] are known extremely accurately. # _SUN_GM = 0.2959122082855911e-03 +_MERCURY_GM = 0.4912547451450812e-10 +_VENUS_GM = 0.7243452486162703e-09 +_EARTH_GM = 0.8887692390113509e-09 +_MARS_GM = 0.9549535105779258e-10 _JUPITER_GM = 0.2825345909524226e-06 _SATURN_GM = 0.8459715185680659e-07 _URANUS_GM = 0.1292024916781969e-07 _NEPTUNE_GM = 0.1524358900784276e-07 +_PLUTO_GM = 0.2188699765425970e-11 + +_MOON_GM = _EARTH_GM / _EARTH_MOON_MASS_RATIO + + +def MassProduct(body): + """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). + + Parameters + ---------- + 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 + ------- + float + The mass product of the given body in au^3/day^2. + """ + if body == Body.Sun: return _SUN_GM + if body == Body.Mercury: return _MERCURY_GM + if body == Body.Venus: return _VENUS_GM + if body == Body.Earth: return _EARTH_GM + if body == Body.Moon: return _MOON_GM + if body == Body.EMB: return _EARTH_GM + _MOON_GM + if body == Body.Mars: return _MARS_GM + if body == Body.Jupiter: return _JUPITER_GM + if body == Body.Saturn: return _SATURN_GM + if body == Body.Uranus: return _URANUS_GM + if body == Body.Neptune: return _NEPTUNE_GM + if body == Body.Pluto: return _PLUTO_GM + raise InvalidBodyError() @enum.unique class _PrecessDir(enum.Enum): @@ -6700,3 +6750,244 @@ def RotationAxis(body, time): time ) return AxisInfo(ra/15.0, dec, w, north) + + +def LagrangePoint(point, time, major_body, minor_body): + """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 + (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 #LagrangePointFast, + especially when the state vectors have already been calculated, or are needed + for some other purpose. + + Parameters + ---------- + point : int + An integer 1..5 that selects which of the Lagrange points to calculate. + time : Time + The time for which the Lagrange point is to be calculated. + major_body : Body + The more massive of the co-orbiting bodies: `Body.Sun` or `Body.Earth`. + minor_body : Body + The less massive of the co-orbiting bodies. See main remarks. + + Returns + ------- + StateVector + The position and velocity of the selected Lagrange point with respect to the major body's center. + """ + major_mass = MassProduct(major_body) + minor_mass = MassProduct(minor_body) + + # Calculate the state vectors for the major and minor bodies. + if major_body == Body.Earth and minor_body == Body.Moon: + # Use geocentric calculations for more precision. + major_state = StateVector(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, time) + 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 + ) + + +def 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 #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 #LagrangePoint, instead of this function, for simpler usage in most cases. + + Parameters + ---------- + point : int + An integer 1..5 that selects which of the Lagrange points to calculate. + major_state : StateVector + The state vector of the major (more massive) of the pair of bodies. + major_mass : float + The mass product GM of the major body. + minor_state : StateVector + The state vector of the minor (less massive) of the pair of bodies. + minor_mass : float + The mass product GM of the minor body. + + Returns + ------- + StateVector + The position and velocity of the selected Lagrange point with respect to the major body's center. + """ + cos_60 = 0.5 + sin_60 = 0.8660254037844386 # sqrt(3) / 2 + + if point not in [1, 2, 3, 4, 5]: + raise Error('Invalid Lagrange point [{}]. Must be an integer 1..5.'.format(point)) + + if major_mass <= 0.0: + raise Error('Major mass must be a positive number.') + + if minor_mass <= 0.0: + raise Error('Minor mass must be a positive number.') + + # Find the relative position vector . + dx = minor_state.x - major_state.x + dy = minor_state.y - major_state.y + dz = minor_state.z - major_state.z + R2 = (dx*dx + dy*dy + dz*dz) + + # R = Total distance between the bodies. + R = math.sqrt(R2) + + # Find the velocity vector . + vx = minor_state.vx - major_state.vx + vy = minor_state.vy - major_state.vy + vz = minor_state.vz - major_state.vz + + if point in [4, 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 = dy*vz - dz*vy + ny = dz*vx - dx*vz + nz = dx*vy - dy*vx + + # Take the cross product normal*position to get a tangential vector . + ux = ny*dz - nz*dy + uy = nz*dx - nx*dz + uz = nx*dy - ny*dx + + # Convert the tangential direction vector to a unit vector. + 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. + if point == 4: + vert = +sin_60 + else: + vert = -sin_60 + + # Rotated radial vector + Dx = cos_60*dx + vert*ux + Dy = cos_60*dy + vert*uy + Dz = cos_60*dz + vert*uz + + # Rotated tangent vector + Ux = cos_60*ux - vert*dx + Uy = cos_60*uy - vert*dy + Uz = cos_60*uz - vert*dz + + # Calculate L4/L5 positions relative to the major body. + px = R * Dx + py = R * Dy + pz = R * Dz + + # Use dot products to find radial and tangential components of the relative velocity. + vrad = vx*dx + vy*dy + vz*dz + vtan = vx*ux + vy*uy + vz*uz + + # Calculate L4/L5 velocities. + pvx = vrad*Dx + vtan*Ux + pvy = vrad*Dy + vtan*Uy + pvz = vrad*Dz + vtan*Uz + + return StateVector(px, py, pz, pvx, pvy, pvz, major_state.t) + + # Handle Langrange points 1, 2, 3. + # 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) + r1 = -R * (minor_mass / (major_mass + minor_mass)) + r2 = +R * (major_mass / (major_mass + minor_mass)) + + # Calculate the square of the angular orbital speed in [rad^2 / day^2]. + 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. + if point in [1, 2]: + scale = (major_mass / (major_mass + minor_mass)) * _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. + x = R*scale - r1 + while True: + dr1 = x - r1 + dr2 = x - r2 + accel = omega2*x + numer1/(dr1*dr1) + numer2/(dr2*dr2) + deriv = omega2 - 2*numer1/(dr1*dr1*dr1) - 2*numer2/(dr2*dr2*dr2) + deltax = accel/deriv + x -= deltax + if abs(deltax/R) <= 1.0e-14: + break + scale = (x - r1) / R + return StateVector(scale*dx, scale*dy, scale*dz, scale*vx, scale*vy, scale*vz, major_state.t) diff --git a/generate/test.py b/generate/test.py index 77f94405..1b66aaed 100755 --- a/generate/test.py +++ b/generate/test.py @@ -2302,6 +2302,39 @@ def AxisTestBody(body, filename, arcmin_tolerance): #----------------------------------------------------------------------------------------------------------- +class LagrangeFunc: + def __init__(self, point, major_body, minor_body): + self.point = point + self.major_body = major_body + self.minor_body = minor_body + + def Eval(self, time): + return astronomy.LagrangePoint(self.point, time, self.major_body, self.minor_body) + + +def VerifyStateLagrange(major_body, minor_body, point, filename, r_thresh, v_thresh): + func = LagrangeFunc(point, major_body, minor_body) + return VerifyStateBody(func, filename, r_thresh, v_thresh) + + +def Lagrange(): + # Test Sun/EMB Lagrange points. + if VerifyStateLagrange(astronomy.Body.Sun, astronomy.Body.EMB, 1, 'lagrange/semb_L1.txt', 1.33e-5, 6.13e-5): return 1 + if VerifyStateLagrange(astronomy.Body.Sun, astronomy.Body.EMB, 2, 'lagrange/semb_L2.txt', 1.33e-5, 6.13e-5): return 1 + if VerifyStateLagrange(astronomy.Body.Sun, astronomy.Body.EMB, 4, 'lagrange/semb_L4.txt', 3.75e-5, 5.28e-5): return 1 + if VerifyStateLagrange(astronomy.Body.Sun, astronomy.Body.EMB, 5, 'lagrange/semb_L5.txt', 3.75e-5, 5.28e-5): return 1 + + # Test Earth/Moon Lagrange points. + if VerifyStateLagrange(astronomy.Body.Earth, astronomy.Body.Moon, 1, 'lagrange/em_L1.txt', 3.79e-5, 5.06e-5): return 1 + if VerifyStateLagrange(astronomy.Body.Earth, astronomy.Body.Moon, 2, 'lagrange/em_L2.txt', 3.79e-5, 5.06e-5): return 1 + if VerifyStateLagrange(astronomy.Body.Earth, astronomy.Body.Moon, 4, 'lagrange/em_L4.txt', 3.79e-5, 1.59e-3): return 1 + if VerifyStateLagrange(astronomy.Body.Earth, astronomy.Body.Moon, 5, 'lagrange/em_L5.txt', 3.79e-5, 1.59e-3): return 1 + + print('PY Lagrange: PASS') + return 0 + +#----------------------------------------------------------------------------------------------------------- + UnitTests = { 'aberration': Aberration, 'axis': Axis, @@ -2313,6 +2346,7 @@ UnitTests = { 'heliostate': HelioState, 'issue_103': Issue103, 'jupiter_moons': JupiterMoons, + 'lagrange': Lagrange, 'libration': Libration, 'local_solar_eclipse': LocalSolarEclipse, 'lunar_apsis': LunarApsis, diff --git a/source/csharp/README.md b/source/csharp/README.md index 3fc73e87..410e7d5e 100644 --- a/source/csharp/README.md +++ b/source/csharp/README.md @@ -788,8 +788,8 @@ 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. | +| `int` | `point` | An integer 1..5 that selects which of the Lagrange points to calculate. | +| [`AstroTime`](#AstroTime) | `time` | The time for 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. | @@ -826,7 +826,7 @@ Consider calling [`Astronomy.LagrangePoint`](#Astronomy.LagrangePoint), instead | Type | Parameter | Description | | --- | --- | --- | -| `int` | `point` | A value 1..5 that selects which of the Lagrange points to calculate. | +| `int` | `point` | An integer 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. | diff --git a/source/csharp/astronomy.cs b/source/csharp/astronomy.cs index 98de9aed..ad9e703f 100644 --- a/source/csharp/astronomy.cs +++ b/source/csharp/astronomy.cs @@ -8378,8 +8378,8 @@ namespace CosineKitty /// 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. + /// An integer 1..5 that selects which of the Lagrange points to calculate. + /// The time for 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. /// The position and velocity of the selected Lagrange point with respect to the major body's center. @@ -8450,7 +8450,7 @@ namespace CosineKitty /// /// 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. + /// An integer 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. diff --git a/source/python/README.md b/source/python/README.md index e836026c..99bfeafb 100644 --- a/source/python/README.md +++ b/source/python/README.md @@ -1606,6 +1606,83 @@ The positions and velocities of Jupiter's 4 largest moons. --- + +### LagrangePoint(point, time, major_body, minor_body) + +**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 +(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 [`LagrangePointFast`](#LagrangePointFast), +especially when the state vectors have already been calculated, or are needed +for some other purpose. + +| Type | Parameter | Description | +| --- | --- | --- | +| `int` | `point` | An integer 1..5 that selects which of the Lagrange points to calculate. | +| [`Time`](#Time) | `time` | The time for 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. | + +### Returns: [`StateVector`](#StateVector) +The position and velocity of the selected Lagrange point with respect to the major body's center. + +--- + + +### 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 [`MassProduct`](#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 [`LagrangePoint`](#LagrangePoint), instead of this function, for simpler usage in most cases. + +| Type | Parameter | Description | +| --- | --- | --- | +| `int` | `point` | An integer 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. | +| `float` | `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. | +| `float` | `minor_mass` | The mass product GM of the minor body. | + +### Returns: [`StateVector`](#StateVector) +The position and velocity of the selected Lagrange point with respect to the major body's center. + +--- + ### Libration(time) @@ -1631,6 +1708,28 @@ and the apparent angular diameter of the Moon `diam_deg`. --- + +### MassProduct(body) + +**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: `float` +The mass product of the given body in au^3/day^2. + +--- + ### MoonPhase(time) diff --git a/source/python/astronomy.py b/source/python/astronomy.py index fa815f4f..92d25c2e 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -36,6 +36,11 @@ import datetime import enum import re +def _cbrt(x): + if x < 0.0: + return -((-x) ** (1.0 / 3.0)) + return x ** (1.0 / 3.0) + KM_PER_AU = 1.4959787069098932e+8 # The number of kilometers per astronomical unit. C_AUDAY = 173.1446326846693 # The speed of light expressed in astronomical units per day. @@ -105,10 +110,55 @@ _EARTH_MOON_MASS_RATIO = 81.30056 # the products GM[i] are known extremely accurately. # _SUN_GM = 0.2959122082855911e-03 +_MERCURY_GM = 0.4912547451450812e-10 +_VENUS_GM = 0.7243452486162703e-09 +_EARTH_GM = 0.8887692390113509e-09 +_MARS_GM = 0.9549535105779258e-10 _JUPITER_GM = 0.2825345909524226e-06 _SATURN_GM = 0.8459715185680659e-07 _URANUS_GM = 0.1292024916781969e-07 _NEPTUNE_GM = 0.1524358900784276e-07 +_PLUTO_GM = 0.2188699765425970e-11 + +_MOON_GM = _EARTH_GM / _EARTH_MOON_MASS_RATIO + + +def MassProduct(body): + """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). + + Parameters + ---------- + 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 + ------- + float + The mass product of the given body in au^3/day^2. + """ + if body == Body.Sun: return _SUN_GM + if body == Body.Mercury: return _MERCURY_GM + if body == Body.Venus: return _VENUS_GM + if body == Body.Earth: return _EARTH_GM + if body == Body.Moon: return _MOON_GM + if body == Body.EMB: return _EARTH_GM + _MOON_GM + if body == Body.Mars: return _MARS_GM + if body == Body.Jupiter: return _JUPITER_GM + if body == Body.Saturn: return _SATURN_GM + if body == Body.Uranus: return _URANUS_GM + if body == Body.Neptune: return _NEPTUNE_GM + if body == Body.Pluto: return _PLUTO_GM + raise InvalidBodyError() @enum.unique class _PrecessDir(enum.Enum): @@ -9193,3 +9243,244 @@ def RotationAxis(body, time): time ) return AxisInfo(ra/15.0, dec, w, north) + + +def LagrangePoint(point, time, major_body, minor_body): + """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 + (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 #LagrangePointFast, + especially when the state vectors have already been calculated, or are needed + for some other purpose. + + Parameters + ---------- + point : int + An integer 1..5 that selects which of the Lagrange points to calculate. + time : Time + The time for which the Lagrange point is to be calculated. + major_body : Body + The more massive of the co-orbiting bodies: `Body.Sun` or `Body.Earth`. + minor_body : Body + The less massive of the co-orbiting bodies. See main remarks. + + Returns + ------- + StateVector + The position and velocity of the selected Lagrange point with respect to the major body's center. + """ + major_mass = MassProduct(major_body) + minor_mass = MassProduct(minor_body) + + # Calculate the state vectors for the major and minor bodies. + if major_body == Body.Earth and minor_body == Body.Moon: + # Use geocentric calculations for more precision. + major_state = StateVector(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, time) + 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 + ) + + +def 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 #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 #LagrangePoint, instead of this function, for simpler usage in most cases. + + Parameters + ---------- + point : int + An integer 1..5 that selects which of the Lagrange points to calculate. + major_state : StateVector + The state vector of the major (more massive) of the pair of bodies. + major_mass : float + The mass product GM of the major body. + minor_state : StateVector + The state vector of the minor (less massive) of the pair of bodies. + minor_mass : float + The mass product GM of the minor body. + + Returns + ------- + StateVector + The position and velocity of the selected Lagrange point with respect to the major body's center. + """ + cos_60 = 0.5 + sin_60 = 0.8660254037844386 # sqrt(3) / 2 + + if point not in [1, 2, 3, 4, 5]: + raise Error('Invalid Lagrange point [{}]. Must be an integer 1..5.'.format(point)) + + if major_mass <= 0.0: + raise Error('Major mass must be a positive number.') + + if minor_mass <= 0.0: + raise Error('Minor mass must be a positive number.') + + # Find the relative position vector . + dx = minor_state.x - major_state.x + dy = minor_state.y - major_state.y + dz = minor_state.z - major_state.z + R2 = (dx*dx + dy*dy + dz*dz) + + # R = Total distance between the bodies. + R = math.sqrt(R2) + + # Find the velocity vector . + vx = minor_state.vx - major_state.vx + vy = minor_state.vy - major_state.vy + vz = minor_state.vz - major_state.vz + + if point in [4, 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 = dy*vz - dz*vy + ny = dz*vx - dx*vz + nz = dx*vy - dy*vx + + # Take the cross product normal*position to get a tangential vector . + ux = ny*dz - nz*dy + uy = nz*dx - nx*dz + uz = nx*dy - ny*dx + + # Convert the tangential direction vector to a unit vector. + 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. + if point == 4: + vert = +sin_60 + else: + vert = -sin_60 + + # Rotated radial vector + Dx = cos_60*dx + vert*ux + Dy = cos_60*dy + vert*uy + Dz = cos_60*dz + vert*uz + + # Rotated tangent vector + Ux = cos_60*ux - vert*dx + Uy = cos_60*uy - vert*dy + Uz = cos_60*uz - vert*dz + + # Calculate L4/L5 positions relative to the major body. + px = R * Dx + py = R * Dy + pz = R * Dz + + # Use dot products to find radial and tangential components of the relative velocity. + vrad = vx*dx + vy*dy + vz*dz + vtan = vx*ux + vy*uy + vz*uz + + # Calculate L4/L5 velocities. + pvx = vrad*Dx + vtan*Ux + pvy = vrad*Dy + vtan*Uy + pvz = vrad*Dz + vtan*Uz + + return StateVector(px, py, pz, pvx, pvy, pvz, major_state.t) + + # Handle Langrange points 1, 2, 3. + # 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) + r1 = -R * (minor_mass / (major_mass + minor_mass)) + r2 = +R * (major_mass / (major_mass + minor_mass)) + + # Calculate the square of the angular orbital speed in [rad^2 / day^2]. + 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. + if point in [1, 2]: + scale = (major_mass / (major_mass + minor_mass)) * _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. + x = R*scale - r1 + while True: + dr1 = x - r1 + dr2 = x - r2 + accel = omega2*x + numer1/(dr1*dr1) + numer2/(dr2*dr2) + deriv = omega2 - 2*numer1/(dr1*dr1*dr1) - 2*numer2/(dr2*dr2*dr2) + deltax = accel/deriv + x -= deltax + if abs(deltax/R) <= 1.0e-14: + break + scale = (x - r1) / R + return StateVector(scale*dx, scale*dy, scale*dz, scale*vx, scale*vy, scale*vz, major_state.t)