From 1eeea5b46eb90e4bf426f84e7a30c57d11916f4e Mon Sep 17 00:00:00 2001 From: Don Cross Date: Tue, 29 Mar 2022 20:54:17 -0400 Subject: [PATCH] Kotlin: calculate Solar System Barycenter position. --- generate/template/astronomy.kt | 24 ++++++++++++++++++- source/kotlin/README.md | 2 +- source/kotlin/doc/index.md | 2 +- source/kotlin/doc/times.md | 2 ++ .../github/cosinekitty/astronomy/astronomy.kt | 24 ++++++++++++++++++- .../io/github/cosinekitty/astronomy/Tests.kt | 2 -- 6 files changed, 50 insertions(+), 6 deletions(-) diff --git a/generate/template/astronomy.kt b/generate/template/astronomy.kt index cbc82268..b9da73c2 100644 --- a/generate/template/astronomy.kt +++ b/generate/template/astronomy.kt @@ -633,6 +633,9 @@ data class AstroVector( } } +/** + * Multiply a scalar by a vector, yielding another vector. + */ operator fun Double.times(vec: AstroVector) = AstroVector(this*vec.x, this*vec.y, this*vec.z, vec.t) @@ -1998,6 +2001,9 @@ private fun vsopModel(body: Body) = else -> throw InvalidBodyException(body) } +private fun vsopHelioVector(body: Body, time: AstroTime) = + calcVsop(vsopModel(body), time) + private fun vsopHelioState(body: Body, tt: Double) = calcVsopPosVel(vsopModel(body), tt) @@ -3363,6 +3369,22 @@ object Astronomy { private fun helioEarth(time: AstroTime) = calcVsop(vsopTable[2], time) + private fun barycenterContrib(time: AstroTime, body: Body, planetGm: Double) = + (planetGm / (planetGm + SUN_GM)) * vsopHelioVector(body, time) + + private fun solarSystemBarycenter(time: AstroTime): AstroVector { + val j = barycenterContrib(time, Body.Jupiter, JUPITER_GM) + val s = barycenterContrib(time, Body.Saturn, SATURN_GM) + var u = barycenterContrib(time, Body.Uranus, URANUS_GM) + var n = barycenterContrib(time, Body.Neptune, NEPTUNE_GM) + return AstroVector( + j.x + s.x + u.x + n.x, + j.y + s.y + u.y + n.y, + j.z + s.z + u.z + n.z, + time + ) + } + /** * Calculates heliocentric Cartesian coordinates of a body in the J2000 equatorial system. * @@ -3400,7 +3422,7 @@ object Astronomy { Body.Pluto -> calcPluto(time, true).position() Body.Moon -> helioEarth(time) + geoMoon(time) Body.EMB -> helioEarth(time) + (geoMoon(time) / (1.0 + EARTH_MOON_MASS_RATIO)) - // FIXFIXFIX: add Solar System Barycenter + Body.SSB -> solarSystemBarycenter(time) else -> throw InvalidBodyException(body) } diff --git a/source/kotlin/README.md b/source/kotlin/README.md index ff4d63f8..60e6a64c 100644 --- a/source/kotlin/README.md +++ b/source/kotlin/README.md @@ -105,7 +105,7 @@ these are used in function and type names. |---|---| | [degreesToRadians](doc/degrees-to-radians.md) | [jvm]
fun [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html).[degreesToRadians](doc/degrees-to-radians.md)(): [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)
Convert an angle expressed in degrees to an angle expressed in radians. | | [radiansToDegrees](doc/radians-to-degrees.md) | [jvm]
fun [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html).[radiansToDegrees](doc/radians-to-degrees.md)(): [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)
Convert an angle expressed in radians to an angle expressed in degrees. | -| [times](doc/times.md) | [jvm]
operator fun [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html).[times](doc/times.md)(vec: [AstroVector](doc/-astro-vector/index.md)): [AstroVector](doc/-astro-vector/index.md) | +| [times](doc/times.md) | [jvm]
operator fun [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html).[times](doc/times.md)(vec: [AstroVector](doc/-astro-vector/index.md)): [AstroVector](doc/-astro-vector/index.md)
Multiply a scalar by a vector, yielding another vector. | ## Properties diff --git a/source/kotlin/doc/index.md b/source/kotlin/doc/index.md index ece5b0d4..4efebc47 100644 --- a/source/kotlin/doc/index.md +++ b/source/kotlin/doc/index.md @@ -50,7 +50,7 @@ |---|---| | [degreesToRadians](degrees-to-radians.md) | [jvm]
fun [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html).[degreesToRadians](degrees-to-radians.md)(): [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)
Convert an angle expressed in degrees to an angle expressed in radians. | | [radiansToDegrees](radians-to-degrees.md) | [jvm]
fun [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html).[radiansToDegrees](radians-to-degrees.md)(): [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)
Convert an angle expressed in radians to an angle expressed in degrees. | -| [times](times.md) | [jvm]
operator fun [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html).[times](times.md)(vec: [AstroVector](-astro-vector/index.md)): [AstroVector](-astro-vector/index.md) | +| [times](times.md) | [jvm]
operator fun [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html).[times](times.md)(vec: [AstroVector](-astro-vector/index.md)): [AstroVector](-astro-vector/index.md)
Multiply a scalar by a vector, yielding another vector. | ## Properties diff --git a/source/kotlin/doc/times.md b/source/kotlin/doc/times.md index 3db18a62..b3282bd9 100644 --- a/source/kotlin/doc/times.md +++ b/source/kotlin/doc/times.md @@ -4,3 +4,5 @@ [jvm]\ operator fun [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html).[times](times.md)(vec: [AstroVector](-astro-vector/index.md)): [AstroVector](-astro-vector/index.md) + +Multiply a scalar by a vector, yielding another vector. diff --git a/source/kotlin/src/main/kotlin/io/github/cosinekitty/astronomy/astronomy.kt b/source/kotlin/src/main/kotlin/io/github/cosinekitty/astronomy/astronomy.kt index 30dccc99..c95d965f 100644 --- a/source/kotlin/src/main/kotlin/io/github/cosinekitty/astronomy/astronomy.kt +++ b/source/kotlin/src/main/kotlin/io/github/cosinekitty/astronomy/astronomy.kt @@ -633,6 +633,9 @@ data class AstroVector( } } +/** + * Multiply a scalar by a vector, yielding another vector. + */ operator fun Double.times(vec: AstroVector) = AstroVector(this*vec.x, this*vec.y, this*vec.z, vec.t) @@ -1998,6 +2001,9 @@ private fun vsopModel(body: Body) = else -> throw InvalidBodyException(body) } +private fun vsopHelioVector(body: Body, time: AstroTime) = + calcVsop(vsopModel(body), time) + private fun vsopHelioState(body: Body, tt: Double) = calcVsopPosVel(vsopModel(body), tt) @@ -3363,6 +3369,22 @@ object Astronomy { private fun helioEarth(time: AstroTime) = calcVsop(vsopTable[2], time) + private fun barycenterContrib(time: AstroTime, body: Body, planetGm: Double) = + (planetGm / (planetGm + SUN_GM)) * vsopHelioVector(body, time) + + private fun solarSystemBarycenter(time: AstroTime): AstroVector { + val j = barycenterContrib(time, Body.Jupiter, JUPITER_GM) + val s = barycenterContrib(time, Body.Saturn, SATURN_GM) + var u = barycenterContrib(time, Body.Uranus, URANUS_GM) + var n = barycenterContrib(time, Body.Neptune, NEPTUNE_GM) + return AstroVector( + j.x + s.x + u.x + n.x, + j.y + s.y + u.y + n.y, + j.z + s.z + u.z + n.z, + time + ) + } + /** * Calculates heliocentric Cartesian coordinates of a body in the J2000 equatorial system. * @@ -3400,7 +3422,7 @@ object Astronomy { Body.Pluto -> calcPluto(time, true).position() Body.Moon -> helioEarth(time) + geoMoon(time) Body.EMB -> helioEarth(time) + (geoMoon(time) / (1.0 + EARTH_MOON_MASS_RATIO)) - // FIXFIXFIX: add Solar System Barycenter + Body.SSB -> solarSystemBarycenter(time) else -> throw InvalidBodyException(body) } diff --git a/source/kotlin/src/test/kotlin/io/github/cosinekitty/astronomy/Tests.kt b/source/kotlin/src/test/kotlin/io/github/cosinekitty/astronomy/Tests.kt index f17fdcee..b56fee24 100644 --- a/source/kotlin/src/test/kotlin/io/github/cosinekitty/astronomy/Tests.kt +++ b/source/kotlin/src/test/kotlin/io/github/cosinekitty/astronomy/Tests.kt @@ -373,14 +373,12 @@ class Tests { // 0.0020028944141634816, -0.015705734333781345, -0.006808363411687111 ) -/* verifyHelio( Body.SSB, time, 0.008844144038360908, -0.002316519421832873, -0.0012061496455500675 // 2.4572258206466927e-06, 8.124488621688556e-06, 3.38376990631319e-06 ) -*/ } private fun verifyHelio(