From 0ed52f3cea037d817dae2ea6cbcb808fe43cbbe7 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Sat, 16 Apr 2022 15:31:25 -0400 Subject: [PATCH] Kotlin: barycentric states, geocentric moon/EMB Added geoEmbState: calculates the geocentric state vector of the Earth/Moon Barycenter (EMB). Added baryState: calculates the barycentric state vector of a body. Unit tests now support parsing JPL Horizons state vector text files. --- generate/template/astronomy.kt | 82 ++++++++++++ source/kotlin/README.md | 2 + source/kotlin/doc/bary-state.md | 23 ++++ source/kotlin/doc/geo-emb-state.md | 22 ++++ source/kotlin/doc/helio-state.md | 2 +- source/kotlin/doc/index.md | 2 + .../github/cosinekitty/astronomy/astronomy.kt | 82 ++++++++++++ .../io/github/cosinekitty/astronomy/Tests.kt | 118 ++++++++++++++++++ 8 files changed, 332 insertions(+), 1 deletion(-) create mode 100644 source/kotlin/doc/bary-state.md create mode 100644 source/kotlin/doc/geo-emb-state.md diff --git a/generate/template/astronomy.kt b/generate/template/astronomy.kt index 6e4246e7..d6e47667 100644 --- a/generate/template/astronomy.kt +++ b/generate/template/astronomy.kt @@ -2797,6 +2797,15 @@ internal class BodyState( } } + +private fun exportState(bodyState: BodyState, time: Time) = + StateVector( + bodyState.r.x, bodyState.r.y, bodyState.r.z, + bodyState.v.x, bodyState.v.y, bodyState.v.z, + time + ) + + private fun calcVsopPosVel(model: VsopModel, tt: Double): BodyState { val t = tt / DAYS_PER_MILLENNIUM @@ -4358,6 +4367,22 @@ fun geoMoonState(time: Time): StateVector { ) } +/** + * Calculates the geocentric position and velocity of the Earth/Moon barycenter. + * + * Given a time of observation, calculates the geocentric position and velocity vectors + * of the Earth/Moon barycenter (EMB). + * The position (x, y, z) components are expressed in AU (astronomical units). + * The velocity (vx, vy, vz) components are expressed in AU/day. + * + * @param time + * The date and time for which to calculate the EMB vectors. + * + * @return The EMB's position and velocity vectors in geocentric J2000 equatorial coordinates. + */ +fun geoEmbState(time: Time): StateVector = + geoMoonState(time) / (1.0 + EARTH_MOON_MASS_RATIO) + private fun helioEarthPos(time: Time) = calcVsop(vsopModel(Body.Earth), time) @@ -4508,6 +4533,63 @@ fun helioState(body: Body, time: Time): StateVector = else -> throw InvalidBodyException(body) } + +/** + * Calculates barycentric position and velocity vectors for the given body. + * + * Given a body and a time, calculates the barycentric position and velocity + * vectors for the center of that body at that time. + * The vectors are expressed in equatorial J2000 coordinates (EQJ). + * + * @param body + * The celestial body whose barycentric state vector is to be calculated. + * Supported values are `Body.Sun`, `Body.Moon`, `Body.EMB`, `Body.SSB`, and all planets: + * `Body.Mercury`, `Body.Venus`, `Body.Earth`, `Body.Mars`, `Body.Jupiter`, + * `Body.Saturn`, `Body.Uranus`, `Body.Neptune`, `Body.Pluto`. + * + * @param time + * The date and time for which to calculate position and velocity. + * + * @return The barycentric position and velocity vectors of the body. + */ +fun baryState(body: Body, time: Time): StateVector { + // Trival case: the solar system barycenter itself: + if (body == Body.SSB) + return StateVector(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, time) + + if (body == Body.Pluto) + return calcPluto(time, false) + + // Find the barycentric positions and velocities for the 5 major bodies. + val bary = majorBodyBary(time.tt) + + return when (body) { + // If the caller is asking for one of the major bodies, we already have the answer. + Body.Sun -> exportState(bary.sun, time) + Body.Jupiter -> exportState(bary.jupiter, time) + Body.Saturn -> exportState(bary.saturn, time) + Body.Uranus -> exportState(bary.uranus, time) + Body.Neptune -> exportState(bary.neptune, time) + + // The Moon and EMB require calculating both the heliocentric Earth and geocentric Moon. + Body.Moon -> exportState(bary.sun, time) + helioEarthState(time) + geoMoonState(time) + Body.EMB -> exportState(bary.sun, time) + helioEarthState(time) + geoEmbState(time) + + else -> { + val planet: BodyState = calcVsopPosVel(vsopModel(body), time.tt) + StateVector( + bary.sun.r.x + planet.r.x, + bary.sun.r.y + planet.r.y, + bary.sun.r.z + planet.r.z, + bary.sun.v.x + planet.v.x, + bary.sun.v.y + planet.v.y, + bary.sun.v.z + planet.v.z, + time + ) + } + } +} + /** * Calculates geocentric Cartesian coordinates of a body in the J2000 equatorial system. * diff --git a/source/kotlin/README.md b/source/kotlin/README.md index 0a6b48dd..6fda4b84 100644 --- a/source/kotlin/README.md +++ b/source/kotlin/README.md @@ -106,12 +106,14 @@ these are used in function and type names. | Name | Summary | |---|---| | [angleFromSun](doc/angle-from-sun.md) | [jvm]
fun [angleFromSun](doc/angle-from-sun.md)(body: [Body](doc/-body/index.md), time: [Time](doc/-time/index.md)): [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)
Returns the angle between the given body and the Sun, as seen from the Earth. | +| [baryState](doc/bary-state.md) | [jvm]
fun [baryState](doc/bary-state.md)(body: [Body](doc/-body/index.md), time: [Time](doc/-time/index.md)): [StateVector](doc/-state-vector/index.md)
Calculates barycentric position and velocity vectors for the given body. | | [constellation](doc/constellation.md) | [jvm]
fun [constellation](doc/constellation.md)(ra: [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html), dec: [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)): [ConstellationInfo](doc/-constellation-info/index.md)
Determines the constellation that contains the given point in the sky. | | [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. | | [eclipticGeoMoon](doc/ecliptic-geo-moon.md) | [jvm]
fun [eclipticGeoMoon](doc/ecliptic-geo-moon.md)(time: [Time](doc/-time/index.md)): [Spherical](doc/-spherical/index.md)
Calculates spherical ecliptic geocentric position of the Moon. | | [eclipticLongitude](doc/ecliptic-longitude.md) | [jvm]
fun [eclipticLongitude](doc/ecliptic-longitude.md)(body: [Body](doc/-body/index.md), time: [Time](doc/-time/index.md)): [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)
Calculates heliocentric ecliptic longitude of a body based on the J2000 equinox. | | [equator](doc/equator.md) | [jvm]
fun [equator](doc/equator.md)(body: [Body](doc/-body/index.md), time: [Time](doc/-time/index.md), observer: [Observer](doc/-observer/index.md), equdate: [EquatorEpoch](doc/-equator-epoch/index.md), aberration: [Aberration](doc/-aberration/index.md)): [Equatorial](doc/-equatorial/index.md)
Calculates equatorial coordinates of a celestial body as seen by an observer on the Earth's surface. | | [equatorialToEcliptic](doc/equatorial-to-ecliptic.md) | [jvm]
fun [equatorialToEcliptic](doc/equatorial-to-ecliptic.md)(equ: [Vector](doc/-vector/index.md)): [Ecliptic](doc/-ecliptic/index.md)
Converts J2000 equatorial Cartesian coordinates to J2000 ecliptic coordinates. | +| [geoEmbState](doc/geo-emb-state.md) | [jvm]
fun [geoEmbState](doc/geo-emb-state.md)(time: [Time](doc/-time/index.md)): [StateVector](doc/-state-vector/index.md)
Calculates the geocentric position and velocity of the Earth/Moon barycenter. | | [geoMoon](doc/geo-moon.md) | [jvm]
fun [geoMoon](doc/geo-moon.md)(time: [Time](doc/-time/index.md)): [Vector](doc/-vector/index.md)
Calculates equatorial geocentric position of the Moon at a given time. | | [geoMoonState](doc/geo-moon-state.md) | [jvm]
fun [geoMoonState](doc/geo-moon-state.md)(time: [Time](doc/-time/index.md)): [StateVector](doc/-state-vector/index.md)
Calculates equatorial geocentric position and velocity of the Moon at a given time. | | [geoVector](doc/geo-vector.md) | [jvm]
fun [geoVector](doc/geo-vector.md)(body: [Body](doc/-body/index.md), time: [Time](doc/-time/index.md), aberration: [Aberration](doc/-aberration/index.md)): [Vector](doc/-vector/index.md)
Calculates geocentric Cartesian coordinates of a body in the J2000 equatorial system. | diff --git a/source/kotlin/doc/bary-state.md b/source/kotlin/doc/bary-state.md new file mode 100644 index 00000000..ff1264b0 --- /dev/null +++ b/source/kotlin/doc/bary-state.md @@ -0,0 +1,23 @@ +//[astronomy](../../index.md)/[io.github.cosinekitty.astronomy](index.md)/[baryState](bary-state.md) + +# baryState + +[jvm]\ +fun [baryState](bary-state.md)(body: [Body](-body/index.md), time: [Time](-time/index.md)): [StateVector](-state-vector/index.md) + +Calculates barycentric position and velocity vectors for the given body. + +Given a body and a time, calculates the barycentric position and velocity vectors for the center of that body at that time. The vectors are expressed in equatorial J2000 coordinates (EQJ). + +#### Return + +The barycentric position and velocity vectors of the body. + +## Parameters + +jvm + +| | | +|---|---| +| body | The celestial body whose barycentric state vector is to be calculated. Supported values are `Body.Sun`, `Body.Moon`, `Body.EMB`, `Body.SSB`, and all planets: `Body.Mercury`, `Body.Venus`, `Body.Earth`, `Body.Mars`, `Body.Jupiter`, `Body.Saturn`, `Body.Uranus`, `Body.Neptune`, `Body.Pluto`. | +| time | The date and time for which to calculate position and velocity. | diff --git a/source/kotlin/doc/geo-emb-state.md b/source/kotlin/doc/geo-emb-state.md new file mode 100644 index 00000000..f09a3051 --- /dev/null +++ b/source/kotlin/doc/geo-emb-state.md @@ -0,0 +1,22 @@ +//[astronomy](../../index.md)/[io.github.cosinekitty.astronomy](index.md)/[geoEmbState](geo-emb-state.md) + +# geoEmbState + +[jvm]\ +fun [geoEmbState](geo-emb-state.md)(time: [Time](-time/index.md)): [StateVector](-state-vector/index.md) + +Calculates the geocentric position and velocity of the Earth/Moon barycenter. + +Given a time of observation, calculates the geocentric position and velocity vectors of the Earth/Moon barycenter (EMB). The position (x, y, z) components are expressed in AU (astronomical units). The velocity (vx, vy, vz) components are expressed in AU/day. + +#### Return + +The EMB's position and velocity vectors in geocentric J2000 equatorial coordinates. + +## Parameters + +jvm + +| | | +|---|---| +| time | The date and time for which to calculate the EMB vectors. | diff --git a/source/kotlin/doc/helio-state.md b/source/kotlin/doc/helio-state.md index b3167cdf..c3a499a7 100644 --- a/source/kotlin/doc/helio-state.md +++ b/source/kotlin/doc/helio-state.md @@ -7,7 +7,7 @@ fun [helioState](helio-state.md)(body: [Body](-body/index.md), time: [Time](-tim Calculates heliocentric position and velocity vectors for the given body. -Given a body and a time, calculates the position and velocity vectors for the center of that body at that time, relative to the center of the Sun. The vectors are expressed in equatorial J2000 coordinates (EQJ). If you need the position vector only, it is more efficient to call [helioVector](helio-vector.md). The Sun's center is a non-inertial frame of reference. In other words, the Sun experiences acceleration due to gravitational forces, mostly from the larger planets (Jupiter, Saturn, Uranus, and Neptune). If you want to calculate momentum, kinetic energy, or other quantities that require a non-accelerating frame of reference, consider using baryState instead. +Given a body and a time, calculates the position and velocity vectors for the center of that body at that time, relative to the center of the Sun. The vectors are expressed in equatorial J2000 coordinates (EQJ). If you need the position vector only, it is more efficient to call [helioVector](helio-vector.md). The Sun's center is a non-inertial frame of reference. In other words, the Sun experiences acceleration due to gravitational forces, mostly from the larger planets (Jupiter, Saturn, Uranus, and Neptune). If you want to calculate momentum, kinetic energy, or other quantities that require a non-accelerating frame of reference, consider using [baryState](bary-state.md) instead. #### Return diff --git a/source/kotlin/doc/index.md b/source/kotlin/doc/index.md index 5874760e..6bc3c28e 100644 --- a/source/kotlin/doc/index.md +++ b/source/kotlin/doc/index.md @@ -51,12 +51,14 @@ | Name | Summary | |---|---| | [angleFromSun](angle-from-sun.md) | [jvm]
fun [angleFromSun](angle-from-sun.md)(body: [Body](-body/index.md), time: [Time](-time/index.md)): [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)
Returns the angle between the given body and the Sun, as seen from the Earth. | +| [baryState](bary-state.md) | [jvm]
fun [baryState](bary-state.md)(body: [Body](-body/index.md), time: [Time](-time/index.md)): [StateVector](-state-vector/index.md)
Calculates barycentric position and velocity vectors for the given body. | | [constellation](constellation.md) | [jvm]
fun [constellation](constellation.md)(ra: [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html), dec: [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)): [ConstellationInfo](-constellation-info/index.md)
Determines the constellation that contains the given point in the sky. | | [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. | | [eclipticGeoMoon](ecliptic-geo-moon.md) | [jvm]
fun [eclipticGeoMoon](ecliptic-geo-moon.md)(time: [Time](-time/index.md)): [Spherical](-spherical/index.md)
Calculates spherical ecliptic geocentric position of the Moon. | | [eclipticLongitude](ecliptic-longitude.md) | [jvm]
fun [eclipticLongitude](ecliptic-longitude.md)(body: [Body](-body/index.md), time: [Time](-time/index.md)): [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)
Calculates heliocentric ecliptic longitude of a body based on the J2000 equinox. | | [equator](equator.md) | [jvm]
fun [equator](equator.md)(body: [Body](-body/index.md), time: [Time](-time/index.md), observer: [Observer](-observer/index.md), equdate: [EquatorEpoch](-equator-epoch/index.md), aberration: [Aberration](-aberration/index.md)): [Equatorial](-equatorial/index.md)
Calculates equatorial coordinates of a celestial body as seen by an observer on the Earth's surface. | | [equatorialToEcliptic](equatorial-to-ecliptic.md) | [jvm]
fun [equatorialToEcliptic](equatorial-to-ecliptic.md)(equ: [Vector](-vector/index.md)): [Ecliptic](-ecliptic/index.md)
Converts J2000 equatorial Cartesian coordinates to J2000 ecliptic coordinates. | +| [geoEmbState](geo-emb-state.md) | [jvm]
fun [geoEmbState](geo-emb-state.md)(time: [Time](-time/index.md)): [StateVector](-state-vector/index.md)
Calculates the geocentric position and velocity of the Earth/Moon barycenter. | | [geoMoon](geo-moon.md) | [jvm]
fun [geoMoon](geo-moon.md)(time: [Time](-time/index.md)): [Vector](-vector/index.md)
Calculates equatorial geocentric position of the Moon at a given time. | | [geoMoonState](geo-moon-state.md) | [jvm]
fun [geoMoonState](geo-moon-state.md)(time: [Time](-time/index.md)): [StateVector](-state-vector/index.md)
Calculates equatorial geocentric position and velocity of the Moon at a given time. | | [geoVector](geo-vector.md) | [jvm]
fun [geoVector](geo-vector.md)(body: [Body](-body/index.md), time: [Time](-time/index.md), aberration: [Aberration](-aberration/index.md)): [Vector](-vector/index.md)
Calculates geocentric Cartesian coordinates of a body in the J2000 equatorial system. | 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 2827d918..905f5226 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 @@ -2797,6 +2797,15 @@ internal class BodyState( } } + +private fun exportState(bodyState: BodyState, time: Time) = + StateVector( + bodyState.r.x, bodyState.r.y, bodyState.r.z, + bodyState.v.x, bodyState.v.y, bodyState.v.z, + time + ) + + private fun calcVsopPosVel(model: VsopModel, tt: Double): BodyState { val t = tt / DAYS_PER_MILLENNIUM @@ -4358,6 +4367,22 @@ fun geoMoonState(time: Time): StateVector { ) } +/** + * Calculates the geocentric position and velocity of the Earth/Moon barycenter. + * + * Given a time of observation, calculates the geocentric position and velocity vectors + * of the Earth/Moon barycenter (EMB). + * The position (x, y, z) components are expressed in AU (astronomical units). + * The velocity (vx, vy, vz) components are expressed in AU/day. + * + * @param time + * The date and time for which to calculate the EMB vectors. + * + * @return The EMB's position and velocity vectors in geocentric J2000 equatorial coordinates. + */ +fun geoEmbState(time: Time): StateVector = + geoMoonState(time) / (1.0 + EARTH_MOON_MASS_RATIO) + private fun helioEarthPos(time: Time) = calcVsop(vsopModel(Body.Earth), time) @@ -4508,6 +4533,63 @@ fun helioState(body: Body, time: Time): StateVector = else -> throw InvalidBodyException(body) } + +/** + * Calculates barycentric position and velocity vectors for the given body. + * + * Given a body and a time, calculates the barycentric position and velocity + * vectors for the center of that body at that time. + * The vectors are expressed in equatorial J2000 coordinates (EQJ). + * + * @param body + * The celestial body whose barycentric state vector is to be calculated. + * Supported values are `Body.Sun`, `Body.Moon`, `Body.EMB`, `Body.SSB`, and all planets: + * `Body.Mercury`, `Body.Venus`, `Body.Earth`, `Body.Mars`, `Body.Jupiter`, + * `Body.Saturn`, `Body.Uranus`, `Body.Neptune`, `Body.Pluto`. + * + * @param time + * The date and time for which to calculate position and velocity. + * + * @return The barycentric position and velocity vectors of the body. + */ +fun baryState(body: Body, time: Time): StateVector { + // Trival case: the solar system barycenter itself: + if (body == Body.SSB) + return StateVector(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, time) + + if (body == Body.Pluto) + return calcPluto(time, false) + + // Find the barycentric positions and velocities for the 5 major bodies. + val bary = majorBodyBary(time.tt) + + return when (body) { + // If the caller is asking for one of the major bodies, we already have the answer. + Body.Sun -> exportState(bary.sun, time) + Body.Jupiter -> exportState(bary.jupiter, time) + Body.Saturn -> exportState(bary.saturn, time) + Body.Uranus -> exportState(bary.uranus, time) + Body.Neptune -> exportState(bary.neptune, time) + + // The Moon and EMB require calculating both the heliocentric Earth and geocentric Moon. + Body.Moon -> exportState(bary.sun, time) + helioEarthState(time) + geoMoonState(time) + Body.EMB -> exportState(bary.sun, time) + helioEarthState(time) + geoEmbState(time) + + else -> { + val planet: BodyState = calcVsopPosVel(vsopModel(body), time.tt) + StateVector( + bary.sun.r.x + planet.r.x, + bary.sun.r.y + planet.r.y, + bary.sun.r.z + planet.r.z, + bary.sun.v.x + planet.v.x, + bary.sun.v.y + planet.v.y, + bary.sun.v.z + planet.v.z, + time + ) + } + } +} + /** * Calculates geocentric Cartesian coordinates of a body in the J2000 equatorial system. * 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 6084c48e..c7375393 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 @@ -1429,5 +1429,123 @@ class Tests { assertTrue(diff < 1.0e-16, "excessive ecliptic longitude error for $body: $diff degrees") } + //---------------------------------------------------------------------------------------- + // Generic utility function for loading a series of state vectors + // from a JPL Horizons output file. + + private class JplStateRecord( + val lnum: Int, // the line number where the state vector ends in the JPL Horizons text file + val state: StateVector // the state vector itself: position, velocity, and time + ) + + private fun jplHorizonsStateVectors(filename: String): List { + val infile = File(filename) + var lnum = 0 + var foundBegin = false + var part = 0 + var time: Time? = null + val pos = arrayOf(0.0, 0.0, 0.0) + val vel = arrayOf(0.0, 0.0, 0.0) + var list = ArrayList() + var match: MatchResult? + val regexPos = Regex("""\s*X =\s*(\S+) Y =\s*(\S+) Z =\s*(\S+)""") + val regexVel = Regex("""\s*VX=\s*(\S+) VY=\s*(\S+) VZ=\s*(\S+)""") + for (line in infile.readLines()) { + ++lnum + if (!foundBegin) { + if (line == "\$\$SOE") + foundBegin = true + } else { + // Input comes in triplets of lines: + // + // 2444249.500000000 = A.D. 1980-Jan-11 00:00:00.0000 TDB + // X =-3.314860345089456E-01 Y = 8.463418210972562E-01 Z = 3.667227830514760E-01 + // VX=-1.642704711077836E-02 VY=-5.494770742558920E-03 VZ=-2.383170237527642E-03 + // + // Track which of these 3 cases we are in using the 'part' variable... + when (part) { + 0 -> { + if (line == "\$\$EOE") // end-of-data marker: the list is complete. + return list + + // 2444249.500000000 = A.D. 1980-Jan-11 00:00:00.0000 TDB + // Convert JD to J2000 TT. + val tt = tokenize(line)[0].toDouble() - 2451545.0 + time = Time.fromTerrestrialTime(tt) + } + 1 -> { + match = regexPos.matchEntire(line) + assertTrue(match != null, "$filename line $lnum: cannot parse position vector") + pos[0] = groupDouble(match, 1, filename, lnum) + pos[1] = groupDouble(match, 2, filename, lnum) + pos[2] = groupDouble(match, 3, filename, lnum) + } + 2 -> { + match = regexVel.matchEntire(line) + assertTrue(match != null, "$filename line $lnum: cannot parse velocity vector") + vel[0] = groupDouble(match, 1, filename, lnum) + vel[1] = groupDouble(match, 2, filename, lnum) + vel[2] = groupDouble(match, 3, filename, lnum) + assertTrue(time != null, "$filename line $lnum: time value was never parsed") + val state = StateVector(pos[0], pos[1], pos[2], vel[0], vel[1], vel[2], time) + list.add(JplStateRecord(lnum, state)) + time = null + } + else -> fail("$filename line $lnum: invalid part=$part") + } + part = (part + 1) % 3 + } + } + fail("$filename: never found end-of-data marker") + } + + private fun stateVectorDiff(relative: Boolean, expected: Vector, calculated: Vector): Double { + val diff = (expected - calculated).length() + return if (relative) + diff / expected.length() + else + diff + } + + private fun verifyStateBody( + relativeFileName: String, + rThresh: Double, + vThresh: Double, + func: (Time) -> StateVector + ) { + var filename = dataRootDir + relativeFileName + for (rec in jplHorizonsStateVectors(filename)) { + val state = func(rec.state.t) + val rdiff = stateVectorDiff(rThresh > 0.0, rec.state.position(), state.position()) + val vdiff = stateVectorDiff(vThresh > 0.0, rec.state.velocity(), state.velocity()) + assertTrue(rdiff < abs(rThresh), "$filename line ${rec.lnum}: excessive position error = $rdiff") + assertTrue(vdiff < abs(vThresh), "$filename line ${rec.lnum}: excessive velocity error = $vdiff") + } + } + + //---------------------------------------------------------------------------------------- + + @Test + fun `Barycentric state vectors`() { + verifyStateBody("barystate/Sun.txt", -1.224e-05, -1.134e-07) { time -> baryState(Body.Sun, time) } + verifyStateBody("barystate/Mercury.txt", 1.672e-04, 2.698e-04) { time -> baryState(Body.Mercury, time) } + verifyStateBody("barystate/Venus.txt", 4.123e-05, 4.308e-05) { time -> baryState(Body.Venus, time) } + verifyStateBody("barystate/Earth.txt", 2.296e-05, 6.359e-05) { time -> baryState(Body.Earth, time) } + verifyStateBody("barystate/Mars.txt", 3.107e-05, 5.550e-05) { time -> baryState(Body.Mars, time) } + verifyStateBody("barystate/Jupiter.txt", 7.389e-05, 2.471e-04) { time -> baryState(Body.Jupiter, time) } + verifyStateBody("barystate/Saturn.txt", 1.067e-04, 3.220e-04) { time -> baryState(Body.Saturn, time) } + verifyStateBody("barystate/Uranus.txt", 9.035e-05, 2.519e-04) { time -> baryState(Body.Uranus, time) } + verifyStateBody("barystate/Neptune.txt", 9.838e-05, 4.446e-04) { time -> baryState(Body.Neptune, time) } + verifyStateBody("barystate/Pluto.txt", 4.259e-05, 7.827e-05) { time -> baryState(Body.Pluto, time) } + verifyStateBody("barystate/Moon.txt", 2.354e-05, 6.604e-05) { time -> baryState(Body.Moon, time) } + verifyStateBody("barystate/EMB.txt", 2.353e-05, 6.511e-05) { time -> baryState(Body.EMB, time) } + } + + @Test + fun `Geocentric state vectors`() { + verifyStateBody("barystate/GeoMoon.txt", 4.086e-05, 5.347e-05) { time -> geoMoonState(time) } + verifyStateBody("barystate/GeoEMB.txt", 4.076e-05, 5.335e-05) { time -> geoEmbState(time) } + } + //---------------------------------------------------------------------------------------- }