Kotlin: calculate Solar System Barycenter position.

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)
         }
 

source/kotlin/README.md

@@ -105,7 +105,7 @@ these are used in function and type names.
 |---|---|
 | [degreesToRadians](doc/degrees-to-radians.md) | [jvm]<br>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)<br>Convert an angle expressed in degrees to an angle expressed in radians. |
 | [radiansToDegrees](doc/radians-to-degrees.md) | [jvm]<br>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)<br>Convert an angle expressed in radians to an angle expressed in degrees. |
-| [times](doc/times.md) | [jvm]<br>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]<br>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)<br>Multiply a scalar by a vector, yielding another vector. |
 
 ## Properties
 

source/kotlin/doc/index.md

@@ -50,7 +50,7 @@
 |---|---|
 | [degreesToRadians](degrees-to-radians.md) | [jvm]<br>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)<br>Convert an angle expressed in degrees to an angle expressed in radians. |
 | [radiansToDegrees](radians-to-degrees.md) | [jvm]<br>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)<br>Convert an angle expressed in radians to an angle expressed in degrees. |
-| [times](times.md) | [jvm]<br>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]<br>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)<br>Multiply a scalar by a vector, yielding another vector. |
 
 ## Properties
 

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.

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)
         }
 

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(