mirror of
https://github.com/cosinekitty/astronomy.git
synced 2026-05-19 14:27:52 -04:00
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.
This commit is contained in:
Don Cross
@@ -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.
|
||||
*
|
||||
|
||||
@@ -106,12 +106,14 @@ these are used in function and type names.
|
||||
| Name | Summary |
|
||||
|---|---|
|
||||
| [angleFromSun](doc/angle-from-sun.md) | [jvm]<br>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)<br>Returns the angle between the given body and the Sun, as seen from the Earth. |
|
||||
| [baryState](doc/bary-state.md) | [jvm]<br>fun [baryState](doc/bary-state.md)(body: [Body](doc/-body/index.md), time: [Time](doc/-time/index.md)): [StateVector](doc/-state-vector/index.md)<br>Calculates barycentric position and velocity vectors for the given body. |
|
||||
| [constellation](doc/constellation.md) | [jvm]<br>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)<br>Determines the constellation that contains the given point in the sky. |
|
||||
| [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. |
|
||||
| [eclipticGeoMoon](doc/ecliptic-geo-moon.md) | [jvm]<br>fun [eclipticGeoMoon](doc/ecliptic-geo-moon.md)(time: [Time](doc/-time/index.md)): [Spherical](doc/-spherical/index.md)<br>Calculates spherical ecliptic geocentric position of the Moon. |
|
||||
| [eclipticLongitude](doc/ecliptic-longitude.md) | [jvm]<br>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)<br>Calculates heliocentric ecliptic longitude of a body based on the J2000 equinox. |
|
||||
| [equator](doc/equator.md) | [jvm]<br>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)<br>Calculates equatorial coordinates of a celestial body as seen by an observer on the Earth's surface. |
|
||||
| [equatorialToEcliptic](doc/equatorial-to-ecliptic.md) | [jvm]<br>fun [equatorialToEcliptic](doc/equatorial-to-ecliptic.md)(equ: [Vector](doc/-vector/index.md)): [Ecliptic](doc/-ecliptic/index.md)<br>Converts J2000 equatorial Cartesian coordinates to J2000 ecliptic coordinates. |
|
||||
| [geoEmbState](doc/geo-emb-state.md) | [jvm]<br>fun [geoEmbState](doc/geo-emb-state.md)(time: [Time](doc/-time/index.md)): [StateVector](doc/-state-vector/index.md)<br>Calculates the geocentric position and velocity of the Earth/Moon barycenter. |
|
||||
| [geoMoon](doc/geo-moon.md) | [jvm]<br>fun [geoMoon](doc/geo-moon.md)(time: [Time](doc/-time/index.md)): [Vector](doc/-vector/index.md)<br>Calculates equatorial geocentric position of the Moon at a given time. |
|
||||
| [geoMoonState](doc/geo-moon-state.md) | [jvm]<br>fun [geoMoonState](doc/geo-moon-state.md)(time: [Time](doc/-time/index.md)): [StateVector](doc/-state-vector/index.md)<br>Calculates equatorial geocentric position and velocity of the Moon at a given time. |
|
||||
| [geoVector](doc/geo-vector.md) | [jvm]<br>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)<br>Calculates geocentric Cartesian coordinates of a body in the J2000 equatorial system. |
|
||||
|
||||
23
source/kotlin/doc/bary-state.md
Normal file
23
source/kotlin/doc/bary-state.md
Normal file
@@ -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. |
|
||||
22
source/kotlin/doc/geo-emb-state.md
Normal file
22
source/kotlin/doc/geo-emb-state.md
Normal file
@@ -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. |
|
||||
@@ -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
|
||||
|
||||
|
||||
@@ -51,12 +51,14 @@
|
||||
| Name | Summary |
|
||||
|---|---|
|
||||
| [angleFromSun](angle-from-sun.md) | [jvm]<br>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)<br>Returns the angle between the given body and the Sun, as seen from the Earth. |
|
||||
| [baryState](bary-state.md) | [jvm]<br>fun [baryState](bary-state.md)(body: [Body](-body/index.md), time: [Time](-time/index.md)): [StateVector](-state-vector/index.md)<br>Calculates barycentric position and velocity vectors for the given body. |
|
||||
| [constellation](constellation.md) | [jvm]<br>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)<br>Determines the constellation that contains the given point in the sky. |
|
||||
| [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. |
|
||||
| [eclipticGeoMoon](ecliptic-geo-moon.md) | [jvm]<br>fun [eclipticGeoMoon](ecliptic-geo-moon.md)(time: [Time](-time/index.md)): [Spherical](-spherical/index.md)<br>Calculates spherical ecliptic geocentric position of the Moon. |
|
||||
| [eclipticLongitude](ecliptic-longitude.md) | [jvm]<br>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)<br>Calculates heliocentric ecliptic longitude of a body based on the J2000 equinox. |
|
||||
| [equator](equator.md) | [jvm]<br>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)<br>Calculates equatorial coordinates of a celestial body as seen by an observer on the Earth's surface. |
|
||||
| [equatorialToEcliptic](equatorial-to-ecliptic.md) | [jvm]<br>fun [equatorialToEcliptic](equatorial-to-ecliptic.md)(equ: [Vector](-vector/index.md)): [Ecliptic](-ecliptic/index.md)<br>Converts J2000 equatorial Cartesian coordinates to J2000 ecliptic coordinates. |
|
||||
| [geoEmbState](geo-emb-state.md) | [jvm]<br>fun [geoEmbState](geo-emb-state.md)(time: [Time](-time/index.md)): [StateVector](-state-vector/index.md)<br>Calculates the geocentric position and velocity of the Earth/Moon barycenter. |
|
||||
| [geoMoon](geo-moon.md) | [jvm]<br>fun [geoMoon](geo-moon.md)(time: [Time](-time/index.md)): [Vector](-vector/index.md)<br>Calculates equatorial geocentric position of the Moon at a given time. |
|
||||
| [geoMoonState](geo-moon-state.md) | [jvm]<br>fun [geoMoonState](geo-moon-state.md)(time: [Time](-time/index.md)): [StateVector](-state-vector/index.md)<br>Calculates equatorial geocentric position and velocity of the Moon at a given time. |
|
||||
| [geoVector](geo-vector.md) | [jvm]<br>fun [geoVector](geo-vector.md)(body: [Body](-body/index.md), time: [Time](-time/index.md), aberration: [Aberration](-aberration/index.md)): [Vector](-vector/index.md)<br>Calculates geocentric Cartesian coordinates of a body in the J2000 equatorial system. |
|
||||
|
||||
@@ -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.
|
||||
*
|
||||
|
||||
@@ -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<JplStateRecord> {
|
||||
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<JplStateRecord>()
|
||||
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) }
|
||||
}
|
||||
|
||||
//----------------------------------------------------------------------------------------
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user