mirror of
https://github.com/cosinekitty/astronomy.git
synced 2026-05-19 14:27:52 -04:00
Kotlin: calculate Lagrange points.
Calculate the position and velocity of any
of the five Lagrange points for a pair of
Solar System bodies. Added the functions:
lagrangePoint
lagrangePointFast
This commit is contained in:
Don Cross
@@ -41,6 +41,7 @@ import kotlin.math.hypot
|
||||
import kotlin.math.log10
|
||||
import kotlin.math.min
|
||||
import kotlin.math.PI
|
||||
import kotlin.math.pow
|
||||
import kotlin.math.round
|
||||
import kotlin.math.roundToLong
|
||||
import kotlin.math.sin
|
||||
@@ -202,6 +203,11 @@ private const val NEPTUNE_GM = 0.1524358900784276e-07
|
||||
private const val PLUTO_GM = 0.2188699765425970e-11
|
||||
private const val MOON_GM = EARTH_GM / EARTH_MOON_MASS_RATIO
|
||||
|
||||
private fun cbrt(x: Double): Double = // cube root of any real number
|
||||
if (x < 0.0)
|
||||
-cbrt(-x)
|
||||
else
|
||||
x.pow(1.0 / 3.0)
|
||||
|
||||
private fun Double.withMinDegreeValue(min: Double): Double {
|
||||
var deg = this
|
||||
@@ -6556,6 +6562,273 @@ fun observerGravity(latitude: Double, height: Double): Double {
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* 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 `majorBody`
|
||||
* and `Body.EMB` (Earth/Moon barycenter) for `minorBody`.
|
||||
* For Lagrange points of the Sun and any other planet, pass in just that planet
|
||||
* (e.g. `Body.Jupiter`) for `minorBody`.
|
||||
* 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.
|
||||
*
|
||||
* @param point
|
||||
* An integer 1..5 that selects which of the Lagrange points to calculate.
|
||||
*
|
||||
* @param time
|
||||
* The time for which the Lagrange point is to be calculated.
|
||||
*
|
||||
* @param majorBody
|
||||
* The more massive of the co-orbiting bodies: `Body.Sun` or `Body.Earth`.
|
||||
*
|
||||
* @param minorBody
|
||||
* The less massive of the co-orbiting bodies. See main remarks.
|
||||
*
|
||||
* @return
|
||||
* The position and velocity of the selected Lagrange point with respect to the major body's center.
|
||||
*/
|
||||
fun lagrangePoint(point: Int, time: Time, majorBody: Body, minorBody: Body): StateVector {
|
||||
val majorMass = massProduct(majorBody)
|
||||
val minorMass = massProduct(minorBody)
|
||||
var majorState: StateVector
|
||||
var minorState: StateVector
|
||||
|
||||
// Calculate the state vectors for the major and minor bodies.
|
||||
if (majorBody == Body.Earth && minorBody == Body.Moon) {
|
||||
// Use geocentric calculations for more precision.
|
||||
majorState = StateVector(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, time)
|
||||
minorState = geoMoonState(time)
|
||||
} else {
|
||||
majorState = helioState(majorBody, time)
|
||||
minorState = helioState(minorBody, time)
|
||||
}
|
||||
|
||||
return lagrangePointFast(
|
||||
point,
|
||||
majorState,
|
||||
majorMass,
|
||||
minorState,
|
||||
minorMass
|
||||
)
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* 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 `majorState` and `minorState`,
|
||||
* 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.
|
||||
*
|
||||
* @param point
|
||||
* An integer 1..5 that selects which of the Lagrange points to calculate.
|
||||
*
|
||||
* @param majorState
|
||||
* The state vector of the major (more massive) of the pair of bodies.
|
||||
*
|
||||
* @param majorMass
|
||||
* The mass product GM of the major body.
|
||||
*
|
||||
* @param minorState
|
||||
* The state vector of the minor (less massive) of the pair of bodies.
|
||||
*
|
||||
* @param minorMass
|
||||
* 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.
|
||||
*/
|
||||
fun lagrangePointFast(
|
||||
point: Int,
|
||||
majorState: StateVector,
|
||||
majorMass: Double,
|
||||
minorState: StateVector,
|
||||
minorMass: Double
|
||||
): StateVector {
|
||||
val cos_60 = 0.5
|
||||
val sin_60 = 0.8660254037844386 // sqrt(3) / 2
|
||||
|
||||
if (point < 1 || point > 5)
|
||||
throw IllegalArgumentException("Invalid lagrange point $point")
|
||||
|
||||
if (!majorMass.isFinite() || majorMass <= 0.0)
|
||||
throw IllegalArgumentException("Major mass must be a positive number.")
|
||||
|
||||
if (!minorMass.isFinite() || minorMass <= 0.0)
|
||||
throw IllegalArgumentException("Minor mass must be a positive number.")
|
||||
|
||||
verifyIdenticalTimes(majorState.t, minorState.t)
|
||||
|
||||
// Find the relative position vector <dx, dy, dz>.
|
||||
var dx = minorState.x - majorState.x
|
||||
var dy = minorState.y - majorState.y
|
||||
var dz = minorState.z - majorState.z
|
||||
val R2 = dx*dx + dy*dy + dz*dz
|
||||
|
||||
// R = Total distance between the bodies.
|
||||
val R = sqrt(R2)
|
||||
|
||||
// Find the velocity vector <vx, vy, vz>.
|
||||
val vx = minorState.vx - majorState.vx
|
||||
val vy = minorState.vy - majorState.vy
|
||||
val vz = minorState.vz - majorState.vz
|
||||
|
||||
if (point == 4 || point == 5) {
|
||||
// For L4 and L5, we need to find points 60 degrees away from the
|
||||
// line connecting the two bodies and in the instantaneous orbital plane.
|
||||
// Define the instantaneous orbital plane as the unique plane that contains
|
||||
// both the relative position vector and the relative velocity vector.
|
||||
|
||||
// Take the cross product of position and velocity to find a normal vector <nx, ny, nz>.
|
||||
val nx = dy*vz - dz*vy
|
||||
val ny = dz*vx - dx*vz
|
||||
val nz = dx*vy - dy*vx
|
||||
|
||||
// Take the cross product normal*position to get a tangential vector <ux, uy, uz>.
|
||||
var ux = ny*dz - nz*dy
|
||||
var uy = nz*dx - nx*dz
|
||||
var uz = nx*dy - ny*dx
|
||||
|
||||
// Convert the tangential direction vector to a unit vector.
|
||||
val U = 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.
|
||||
val vert = if (point == 4) +sin_60 else -sin_60
|
||||
|
||||
// Rotated radial vector
|
||||
val Dx = cos_60*dx + vert*ux
|
||||
val Dy = cos_60*dy + vert*uy
|
||||
val Dz = cos_60*dz + vert*uz
|
||||
|
||||
// Rotated tangent vector
|
||||
val Ux = cos_60*ux - vert*dx
|
||||
val Uy = cos_60*uy - vert*dy
|
||||
val Uz = cos_60*uz - vert*dz
|
||||
|
||||
|
||||
// Use dot products to find radial and tangential components of the relative velocity.
|
||||
val vrad = vx*dx + vy*dy + vz*dz
|
||||
val vtan = vx*ux + vy*uy + vz*uz
|
||||
|
||||
return StateVector(
|
||||
R * Dx,
|
||||
R * Dy,
|
||||
R * Dz,
|
||||
vrad*Dx + vtan*Ux,
|
||||
vrad*Dy + vtan*Uy,
|
||||
vrad*Dz + vtan*Uz,
|
||||
majorState.t
|
||||
)
|
||||
}
|
||||
|
||||
// point = 1..3; all these lie on a straight line between the centers of the two bodies.
|
||||
// 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)
|
||||
val r1 = -R * (minorMass / (majorMass + minorMass))
|
||||
val r2 = +R * (majorMass / (majorMass + minorMass))
|
||||
|
||||
// Calculate the square of the angular orbital speed in [rad^2 / day^2].
|
||||
val omega2 = (majorMass + minorMass) / (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.
|
||||
var scale: Double
|
||||
var numer1: Double
|
||||
var numer2: Double
|
||||
|
||||
if (point == 1 || point == 2) {
|
||||
scale = (majorMass / (majorMass + minorMass)) * cbrt(minorMass / (3.0 * majorMass))
|
||||
numer1 = -majorMass // The major mass is to the left of L1 and L2
|
||||
if (point == 1) {
|
||||
scale = 1.0 - scale
|
||||
numer2 = +minorMass // The minor mass is to the right of L1.
|
||||
} else {
|
||||
scale = 1.0 + scale
|
||||
numer2 = -minorMass // The minor mass is to the left of L2.
|
||||
}
|
||||
} else { // point == 3
|
||||
scale = ((7.0/12.0)*minorMass - majorMass) / (minorMass + majorMass)
|
||||
numer1 = +majorMass // major mass is to the right of L3.
|
||||
numer2 = +minorMass // minor mass is to the right of L3.
|
||||
}
|
||||
|
||||
// Iterate Newton's Method until it converges.
|
||||
var x = R*scale - r1
|
||||
var deltax: Double
|
||||
do {
|
||||
val dr1 = x - r1
|
||||
val dr2 = x - r2
|
||||
val accel = omega2*x + numer1/(dr1*dr1) + numer2/(dr2*dr2)
|
||||
val deriv = omega2 - 2*numer1/(dr1*dr1*dr1) - 2*numer2/(dr2*dr2*dr2)
|
||||
deltax = accel/deriv
|
||||
x -= deltax
|
||||
} while (abs(deltax/R) > 1.0e-14)
|
||||
scale = (x - r1) / R
|
||||
|
||||
return StateVector(
|
||||
scale * dx,
|
||||
scale * dy,
|
||||
scale * dz,
|
||||
scale * vx,
|
||||
scale * vy,
|
||||
scale * vz,
|
||||
majorState.t
|
||||
)
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Calculates a rotation matrix from equatorial J2000 (EQJ) to ecliptic J2000 (ECL).
|
||||
*
|
||||
|
||||
@@ -125,6 +125,8 @@ these are used in function and type names.
|
||||
| [illumination](doc/illumination.md) | [jvm]<br>fun [illumination](doc/illumination.md)(body: [Body](doc/-body/index.md), time: [Time](doc/-time/index.md)): [IlluminationInfo](doc/-illumination-info/index.md)<br>Finds visual magnitude, phase angle, and other illumination information about a celestial body. |
|
||||
| [inverseRefractionAngle](doc/inverse-refraction-angle.md) | [jvm]<br>fun [inverseRefractionAngle](doc/inverse-refraction-angle.md)(refraction: [Refraction](doc/-refraction/index.md), bentAltitude: [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)): [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)<br>Calculates the inverse of an atmospheric refraction angle. |
|
||||
| [jupiterMoons](doc/jupiter-moons.md) | [jvm]<br>fun [jupiterMoons](doc/jupiter-moons.md)(time: [Time](doc/-time/index.md)): [JupiterMoonsInfo](doc/-jupiter-moons-info/index.md)<br>Calculates jovicentric positions and velocities of Jupiter's largest 4 moons. |
|
||||
| [lagrangePoint](doc/lagrange-point.md) | [jvm]<br>fun [lagrangePoint](doc/lagrange-point.md)(point: [Int](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-int/index.html), time: [Time](doc/-time/index.md), majorBody: [Body](doc/-body/index.md), minorBody: [Body](doc/-body/index.md)): [StateVector](doc/-state-vector/index.md)<br>Calculates one of the 5 Lagrange points for a pair of co-orbiting bodies. |
|
||||
| [lagrangePointFast](doc/lagrange-point-fast.md) | [jvm]<br>fun [lagrangePointFast](doc/lagrange-point-fast.md)(point: [Int](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-int/index.html), majorState: [StateVector](doc/-state-vector/index.md), majorMass: [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html), minorState: [StateVector](doc/-state-vector/index.md), minorMass: [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)): [StateVector](doc/-state-vector/index.md)<br>Calculates one of the 5 Lagrange points from body masses and state vectors. |
|
||||
| [massProduct](doc/mass-product.md) | [jvm]<br>fun [massProduct](doc/mass-product.md)(body: [Body](doc/-body/index.md)): [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)<br>Returns the product of mass and universal gravitational constant of a Solar System body. |
|
||||
| [moonPhase](doc/moon-phase.md) | [jvm]<br>fun [moonPhase](doc/moon-phase.md)(time: [Time](doc/-time/index.md)): [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)<br>Returns the Moon's phase as an angle from 0 to 360 degrees. |
|
||||
| [nextGlobalSolarEclipse](doc/next-global-solar-eclipse.md) | [jvm]<br>fun [nextGlobalSolarEclipse](doc/next-global-solar-eclipse.md)(prevEclipseTime: [Time](doc/-time/index.md)): [GlobalSolarEclipseInfo](doc/-global-solar-eclipse-info/index.md)<br>Searches for the next global solar eclipse in a series. |
|
||||
|
||||
@@ -70,6 +70,8 @@
|
||||
| [illumination](illumination.md) | [jvm]<br>fun [illumination](illumination.md)(body: [Body](-body/index.md), time: [Time](-time/index.md)): [IlluminationInfo](-illumination-info/index.md)<br>Finds visual magnitude, phase angle, and other illumination information about a celestial body. |
|
||||
| [inverseRefractionAngle](inverse-refraction-angle.md) | [jvm]<br>fun [inverseRefractionAngle](inverse-refraction-angle.md)(refraction: [Refraction](-refraction/index.md), bentAltitude: [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)): [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)<br>Calculates the inverse of an atmospheric refraction angle. |
|
||||
| [jupiterMoons](jupiter-moons.md) | [jvm]<br>fun [jupiterMoons](jupiter-moons.md)(time: [Time](-time/index.md)): [JupiterMoonsInfo](-jupiter-moons-info/index.md)<br>Calculates jovicentric positions and velocities of Jupiter's largest 4 moons. |
|
||||
| [lagrangePoint](lagrange-point.md) | [jvm]<br>fun [lagrangePoint](lagrange-point.md)(point: [Int](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-int/index.html), time: [Time](-time/index.md), majorBody: [Body](-body/index.md), minorBody: [Body](-body/index.md)): [StateVector](-state-vector/index.md)<br>Calculates one of the 5 Lagrange points for a pair of co-orbiting bodies. |
|
||||
| [lagrangePointFast](lagrange-point-fast.md) | [jvm]<br>fun [lagrangePointFast](lagrange-point-fast.md)(point: [Int](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-int/index.html), majorState: [StateVector](-state-vector/index.md), majorMass: [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html), minorState: [StateVector](-state-vector/index.md), minorMass: [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)): [StateVector](-state-vector/index.md)<br>Calculates one of the 5 Lagrange points from body masses and state vectors. |
|
||||
| [massProduct](mass-product.md) | [jvm]<br>fun [massProduct](mass-product.md)(body: [Body](-body/index.md)): [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)<br>Returns the product of mass and universal gravitational constant of a Solar System body. |
|
||||
| [moonPhase](moon-phase.md) | [jvm]<br>fun [moonPhase](moon-phase.md)(time: [Time](-time/index.md)): [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)<br>Returns the Moon's phase as an angle from 0 to 360 degrees. |
|
||||
| [nextGlobalSolarEclipse](next-global-solar-eclipse.md) | [jvm]<br>fun [nextGlobalSolarEclipse](next-global-solar-eclipse.md)(prevEclipseTime: [Time](-time/index.md)): [GlobalSolarEclipseInfo](-global-solar-eclipse-info/index.md)<br>Searches for the next global solar eclipse in a series. |
|
||||
|
||||
34
source/kotlin/doc/lagrange-point-fast.md
Normal file
34
source/kotlin/doc/lagrange-point-fast.md
Normal file
@@ -0,0 +1,34 @@
|
||||
//[astronomy](../../index.md)/[io.github.cosinekitty.astronomy](index.md)/[lagrangePointFast](lagrange-point-fast.md)
|
||||
|
||||
# lagrangePointFast
|
||||
|
||||
[jvm]\
|
||||
fun [lagrangePointFast](lagrange-point-fast.md)(point: [Int](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-int/index.html), majorState: [StateVector](-state-vector/index.md), majorMass: [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html), minorState: [StateVector](-state-vector/index.md), minorMass: [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)): [StateVector](-state-vector/index.md)
|
||||
|
||||
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](mass-product.md) 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 majorState and minorState, and the position and velocity components are with respect to the major body's center.
|
||||
|
||||
Consider calling [lagrangePoint](lagrange-point.md), instead of this function, for simpler usage in most cases.
|
||||
|
||||
#### Return
|
||||
|
||||
The position and velocity of the selected Lagrange point with respect to the major body's center.
|
||||
|
||||
## Parameters
|
||||
|
||||
jvm
|
||||
|
||||
| | |
|
||||
|---|---|
|
||||
| point | An integer 1..5 that selects which of the Lagrange points to calculate. |
|
||||
| majorState | The state vector of the major (more massive) of the pair of bodies. |
|
||||
| majorMass | The mass product GM of the major body. |
|
||||
| minorState | The state vector of the minor (less massive) of the pair of bodies. |
|
||||
| minorMass | The mass product GM of the minor body. |
|
||||
33
source/kotlin/doc/lagrange-point.md
Normal file
33
source/kotlin/doc/lagrange-point.md
Normal file
@@ -0,0 +1,33 @@
|
||||
//[astronomy](../../index.md)/[io.github.cosinekitty.astronomy](index.md)/[lagrangePoint](lagrange-point.md)
|
||||
|
||||
# lagrangePoint
|
||||
|
||||
[jvm]\
|
||||
fun [lagrangePoint](lagrange-point.md)(point: [Int](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-int/index.html), time: [Time](-time/index.md), majorBody: [Body](-body/index.md), minorBody: [Body](-body/index.md)): [StateVector](-state-vector/index.md)
|
||||
|
||||
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 majorBody and Body.EMB (Earth/Moon barycenter) for minorBody. For Lagrange points of the Sun and any other planet, pass in just that planet (e.g. Body.Jupiter) for minorBody. 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](lagrange-point-fast.md), especially when the state vectors have already been calculated, or are needed for some other purpose.
|
||||
|
||||
#### Return
|
||||
|
||||
The position and velocity of the selected Lagrange point with respect to the major body's center.
|
||||
|
||||
## Parameters
|
||||
|
||||
jvm
|
||||
|
||||
| | |
|
||||
|---|---|
|
||||
| point | An integer 1..5 that selects which of the Lagrange points to calculate. |
|
||||
| time | The time for which the Lagrange point is to be calculated. |
|
||||
| majorBody | The more massive of the co-orbiting bodies: `Body.Sun` or `Body.Earth`. |
|
||||
| minorBody | The less massive of the co-orbiting bodies. See main remarks. |
|
||||
@@ -41,6 +41,7 @@ import kotlin.math.hypot
|
||||
import kotlin.math.log10
|
||||
import kotlin.math.min
|
||||
import kotlin.math.PI
|
||||
import kotlin.math.pow
|
||||
import kotlin.math.round
|
||||
import kotlin.math.roundToLong
|
||||
import kotlin.math.sin
|
||||
@@ -202,6 +203,11 @@ private const val NEPTUNE_GM = 0.1524358900784276e-07
|
||||
private const val PLUTO_GM = 0.2188699765425970e-11
|
||||
private const val MOON_GM = EARTH_GM / EARTH_MOON_MASS_RATIO
|
||||
|
||||
private fun cbrt(x: Double): Double = // cube root of any real number
|
||||
if (x < 0.0)
|
||||
-cbrt(-x)
|
||||
else
|
||||
x.pow(1.0 / 3.0)
|
||||
|
||||
private fun Double.withMinDegreeValue(min: Double): Double {
|
||||
var deg = this
|
||||
@@ -6556,6 +6562,273 @@ fun observerGravity(latitude: Double, height: Double): Double {
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* 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 `majorBody`
|
||||
* and `Body.EMB` (Earth/Moon barycenter) for `minorBody`.
|
||||
* For Lagrange points of the Sun and any other planet, pass in just that planet
|
||||
* (e.g. `Body.Jupiter`) for `minorBody`.
|
||||
* 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.
|
||||
*
|
||||
* @param point
|
||||
* An integer 1..5 that selects which of the Lagrange points to calculate.
|
||||
*
|
||||
* @param time
|
||||
* The time for which the Lagrange point is to be calculated.
|
||||
*
|
||||
* @param majorBody
|
||||
* The more massive of the co-orbiting bodies: `Body.Sun` or `Body.Earth`.
|
||||
*
|
||||
* @param minorBody
|
||||
* The less massive of the co-orbiting bodies. See main remarks.
|
||||
*
|
||||
* @return
|
||||
* The position and velocity of the selected Lagrange point with respect to the major body's center.
|
||||
*/
|
||||
fun lagrangePoint(point: Int, time: Time, majorBody: Body, minorBody: Body): StateVector {
|
||||
val majorMass = massProduct(majorBody)
|
||||
val minorMass = massProduct(minorBody)
|
||||
var majorState: StateVector
|
||||
var minorState: StateVector
|
||||
|
||||
// Calculate the state vectors for the major and minor bodies.
|
||||
if (majorBody == Body.Earth && minorBody == Body.Moon) {
|
||||
// Use geocentric calculations for more precision.
|
||||
majorState = StateVector(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, time)
|
||||
minorState = geoMoonState(time)
|
||||
} else {
|
||||
majorState = helioState(majorBody, time)
|
||||
minorState = helioState(minorBody, time)
|
||||
}
|
||||
|
||||
return lagrangePointFast(
|
||||
point,
|
||||
majorState,
|
||||
majorMass,
|
||||
minorState,
|
||||
minorMass
|
||||
)
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* 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 `majorState` and `minorState`,
|
||||
* 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.
|
||||
*
|
||||
* @param point
|
||||
* An integer 1..5 that selects which of the Lagrange points to calculate.
|
||||
*
|
||||
* @param majorState
|
||||
* The state vector of the major (more massive) of the pair of bodies.
|
||||
*
|
||||
* @param majorMass
|
||||
* The mass product GM of the major body.
|
||||
*
|
||||
* @param minorState
|
||||
* The state vector of the minor (less massive) of the pair of bodies.
|
||||
*
|
||||
* @param minorMass
|
||||
* 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.
|
||||
*/
|
||||
fun lagrangePointFast(
|
||||
point: Int,
|
||||
majorState: StateVector,
|
||||
majorMass: Double,
|
||||
minorState: StateVector,
|
||||
minorMass: Double
|
||||
): StateVector {
|
||||
val cos_60 = 0.5
|
||||
val sin_60 = 0.8660254037844386 // sqrt(3) / 2
|
||||
|
||||
if (point < 1 || point > 5)
|
||||
throw IllegalArgumentException("Invalid lagrange point $point")
|
||||
|
||||
if (!majorMass.isFinite() || majorMass <= 0.0)
|
||||
throw IllegalArgumentException("Major mass must be a positive number.")
|
||||
|
||||
if (!minorMass.isFinite() || minorMass <= 0.0)
|
||||
throw IllegalArgumentException("Minor mass must be a positive number.")
|
||||
|
||||
verifyIdenticalTimes(majorState.t, minorState.t)
|
||||
|
||||
// Find the relative position vector <dx, dy, dz>.
|
||||
var dx = minorState.x - majorState.x
|
||||
var dy = minorState.y - majorState.y
|
||||
var dz = minorState.z - majorState.z
|
||||
val R2 = dx*dx + dy*dy + dz*dz
|
||||
|
||||
// R = Total distance between the bodies.
|
||||
val R = sqrt(R2)
|
||||
|
||||
// Find the velocity vector <vx, vy, vz>.
|
||||
val vx = minorState.vx - majorState.vx
|
||||
val vy = minorState.vy - majorState.vy
|
||||
val vz = minorState.vz - majorState.vz
|
||||
|
||||
if (point == 4 || point == 5) {
|
||||
// For L4 and L5, we need to find points 60 degrees away from the
|
||||
// line connecting the two bodies and in the instantaneous orbital plane.
|
||||
// Define the instantaneous orbital plane as the unique plane that contains
|
||||
// both the relative position vector and the relative velocity vector.
|
||||
|
||||
// Take the cross product of position and velocity to find a normal vector <nx, ny, nz>.
|
||||
val nx = dy*vz - dz*vy
|
||||
val ny = dz*vx - dx*vz
|
||||
val nz = dx*vy - dy*vx
|
||||
|
||||
// Take the cross product normal*position to get a tangential vector <ux, uy, uz>.
|
||||
var ux = ny*dz - nz*dy
|
||||
var uy = nz*dx - nx*dz
|
||||
var uz = nx*dy - ny*dx
|
||||
|
||||
// Convert the tangential direction vector to a unit vector.
|
||||
val U = 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.
|
||||
val vert = if (point == 4) +sin_60 else -sin_60
|
||||
|
||||
// Rotated radial vector
|
||||
val Dx = cos_60*dx + vert*ux
|
||||
val Dy = cos_60*dy + vert*uy
|
||||
val Dz = cos_60*dz + vert*uz
|
||||
|
||||
// Rotated tangent vector
|
||||
val Ux = cos_60*ux - vert*dx
|
||||
val Uy = cos_60*uy - vert*dy
|
||||
val Uz = cos_60*uz - vert*dz
|
||||
|
||||
|
||||
// Use dot products to find radial and tangential components of the relative velocity.
|
||||
val vrad = vx*dx + vy*dy + vz*dz
|
||||
val vtan = vx*ux + vy*uy + vz*uz
|
||||
|
||||
return StateVector(
|
||||
R * Dx,
|
||||
R * Dy,
|
||||
R * Dz,
|
||||
vrad*Dx + vtan*Ux,
|
||||
vrad*Dy + vtan*Uy,
|
||||
vrad*Dz + vtan*Uz,
|
||||
majorState.t
|
||||
)
|
||||
}
|
||||
|
||||
// point = 1..3; all these lie on a straight line between the centers of the two bodies.
|
||||
// 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)
|
||||
val r1 = -R * (minorMass / (majorMass + minorMass))
|
||||
val r2 = +R * (majorMass / (majorMass + minorMass))
|
||||
|
||||
// Calculate the square of the angular orbital speed in [rad^2 / day^2].
|
||||
val omega2 = (majorMass + minorMass) / (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.
|
||||
var scale: Double
|
||||
var numer1: Double
|
||||
var numer2: Double
|
||||
|
||||
if (point == 1 || point == 2) {
|
||||
scale = (majorMass / (majorMass + minorMass)) * cbrt(minorMass / (3.0 * majorMass))
|
||||
numer1 = -majorMass // The major mass is to the left of L1 and L2
|
||||
if (point == 1) {
|
||||
scale = 1.0 - scale
|
||||
numer2 = +minorMass // The minor mass is to the right of L1.
|
||||
} else {
|
||||
scale = 1.0 + scale
|
||||
numer2 = -minorMass // The minor mass is to the left of L2.
|
||||
}
|
||||
} else { // point == 3
|
||||
scale = ((7.0/12.0)*minorMass - majorMass) / (minorMass + majorMass)
|
||||
numer1 = +majorMass // major mass is to the right of L3.
|
||||
numer2 = +minorMass // minor mass is to the right of L3.
|
||||
}
|
||||
|
||||
// Iterate Newton's Method until it converges.
|
||||
var x = R*scale - r1
|
||||
var deltax: Double
|
||||
do {
|
||||
val dr1 = x - r1
|
||||
val dr2 = x - r2
|
||||
val accel = omega2*x + numer1/(dr1*dr1) + numer2/(dr2*dr2)
|
||||
val deriv = omega2 - 2*numer1/(dr1*dr1*dr1) - 2*numer2/(dr2*dr2*dr2)
|
||||
deltax = accel/deriv
|
||||
x -= deltax
|
||||
} while (abs(deltax/R) > 1.0e-14)
|
||||
scale = (x - r1) / R
|
||||
|
||||
return StateVector(
|
||||
scale * dx,
|
||||
scale * dy,
|
||||
scale * dz,
|
||||
scale * vx,
|
||||
scale * vy,
|
||||
scale * vz,
|
||||
majorState.t
|
||||
)
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* Calculates a rotation matrix from equatorial J2000 (EQJ) to ecliptic J2000 (ECL).
|
||||
*
|
||||
|
||||
@@ -2037,4 +2037,34 @@ class Tests {
|
||||
}
|
||||
|
||||
//----------------------------------------------------------------------------------------
|
||||
|
||||
@Test
|
||||
fun `Lagrange Points`() {
|
||||
// Test Sun/EMB Lagrange points.
|
||||
verifyStateLagrange(Body.Sun, Body.EMB, 1, "lagrange/semb_L1.txt", 1.33e-5, 6.13e-5)
|
||||
verifyStateLagrange(Body.Sun, Body.EMB, 2, "lagrange/semb_L2.txt", 1.33e-5, 6.13e-5)
|
||||
verifyStateLagrange(Body.Sun, Body.EMB, 4, "lagrange/semb_L4.txt", 3.75e-5, 5.28e-5)
|
||||
verifyStateLagrange(Body.Sun, Body.EMB, 5, "lagrange/semb_L5.txt", 3.75e-5, 5.28e-5)
|
||||
|
||||
// Test Earth/Moon Lagrange points.
|
||||
verifyStateLagrange(Body.Earth, Body.Moon, 1, "lagrange/em_L1.txt", 3.79e-5, 5.06e-5)
|
||||
verifyStateLagrange(Body.Earth, Body.Moon, 2, "lagrange/em_L2.txt", 3.79e-5, 5.06e-5)
|
||||
verifyStateLagrange(Body.Earth, Body.Moon, 4, "lagrange/em_L4.txt", 3.79e-5, 1.59e-3)
|
||||
verifyStateLagrange(Body.Earth, Body.Moon, 5, "lagrange/em_L5.txt", 3.79e-5, 1.59e-3)
|
||||
}
|
||||
|
||||
private fun verifyStateLagrange(
|
||||
majorBody: Body,
|
||||
minorBody: Body,
|
||||
point: Int,
|
||||
relativeFileName: String,
|
||||
rThresh: Double,
|
||||
vThresh: Double
|
||||
) {
|
||||
verifyStateBody(relativeFileName, rThresh, vThresh) { time ->
|
||||
lagrangePoint(point, time, majorBody, minorBody)
|
||||
}
|
||||
}
|
||||
|
||||
//----------------------------------------------------------------------------------------
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user