From 044ecfb77499f4e1d4f166b822f7df38d867ca6c Mon Sep 17 00:00:00 2001 From: Don Cross Date: Mon, 28 Mar 2022 12:33:01 -0400 Subject: [PATCH] Kotlin: VSOP87 position, velocity vectors. Implemented the VSOP87 calculation functions for heliocentric position vectors, heliocentric velocity vectors, and heliocentric distances. Implemented Astronomy.helioVector for everything except Pluto and SSB. Corrected small errors in C# documentation. --- generate/template/astronomy.cs | 11 +- generate/template/astronomy.kt | 226 ++++++++++++++++-- source/csharp/README.md | 6 +- source/csharp/astronomy.cs | 11 +- .../kotlin/doc/-astronomy/helio-distance.md | 19 ++ source/kotlin/doc/-astronomy/helio-vector.md | 23 ++ source/kotlin/doc/-astronomy/index.md | 2 + .../github/cosinekitty/astronomy/astronomy.kt | 226 ++++++++++++++++-- .../io/github/cosinekitty/astronomy/Tests.kt | 161 ++++++++++++- 9 files changed, 639 insertions(+), 46 deletions(-) create mode 100644 source/kotlin/doc/-astronomy/helio-distance.md create mode 100644 source/kotlin/doc/-astronomy/helio-vector.md diff --git a/generate/template/astronomy.cs b/generate/template/astronomy.cs index 973fec77..e4e813fc 100644 --- a/generate/template/astronomy.cs +++ b/generate/template/astronomy.cs @@ -3812,9 +3812,12 @@ $ASTRO_IAU_DATA() /// The position is not corrected for light travel time or aberration. /// This is different from the behavior of #Astronomy.GeoVector. /// - /// If given an invalid value for `body`, this function will throw an `ArgumentException`. + /// If given an invalid value for `body`, this function will throw an #InvalidBodyException. /// - /// A body for which to calculate a heliocentric position: the Sun, Moon, EMB, SSB, or any of the planets. + /// + /// A body for which to calculate a heliocentric position: + /// the Sun, Moon, EMB, SSB, or any of the planets. + /// /// The date and time for which to calculate the position. /// A heliocentric position vector of the center of the given body. public static AstroVector HelioVector(Body body, AstroTime time) @@ -3863,14 +3866,14 @@ $ASTRO_IAU_DATA() /// /// /// Given a date and time, this function calculates the distance between - /// the center of `body` and the center of the Sun. + /// the center of `body` and the center of the Sun, expressed in AU. /// For the planets Mercury through Neptune, this function is significantly /// more efficient than calling #Astronomy.HelioVector followed by taking the length /// of the resulting vector. /// /// /// A body for which to calculate a heliocentric distance: - /// the Sun, Moon, or any of the planets. + /// the Sun, Moon, EMB, SSB, or any of the planets. /// /// /// The date and time for which to calculate the heliocentric distance. diff --git a/generate/template/astronomy.kt b/generate/template/astronomy.kt index 3036b0dc..6d4d94a6 100644 --- a/generate/template/astronomy.kt +++ b/generate/template/astronomy.kt @@ -131,6 +131,7 @@ const val KM_PER_AU = 1.4959787069098932e+8 private val TimeZoneUtc = TimeZone.getTimeZone("UTC") private const val DAYS_PER_TROPICAL_YEAR = 365.24217 +private const val DAYS_PER_MILLENNIUM = 365250.0 private const val ASEC360 = 1296000.0 private const val ASEC2RAD = 4.848136811095359935899141e-6 @@ -416,6 +417,7 @@ class AstroTime private constructor( fun addDays(days: Double) = AstroTime(ut + days) internal fun julianCenturies() = tt / 36525.0 + internal fun julianMillennia() = tt / DAYS_PER_MILLENNIUM companion object { private val origin = GregorianCalendar(TimeZoneUtc).also { @@ -464,7 +466,7 @@ internal data class TerseVector(val x: Double, val y: Double, val z: Double) { operator fun div(other: Double) = TerseVector(x / other, y / other, z / other) - val quadrature get() = x * x + y * y + z * z + val quadrature get() = (x * x) + (y * y) + (z * z) val magnitude get() = sqrt(quadrature) companion object { @@ -497,8 +499,7 @@ data class AstroVector( /** * The total distance in AU represented by this vector. */ - fun length() = - sqrt(x*x + y*y + z*z) + fun length() = sqrt((x * x) + (y * y) + (z * z)) private fun verifyIdenticalTimes(otherTime: AstroTime): AstroTime { if (t.tt != otherTime.tt) @@ -1818,20 +1819,139 @@ private class VsopModel( val rad: VsopFormula ) -private fun vsopTableIndex(body: Body): Int = - when (body) { - Body.Mercury -> 0 - Body.Venus -> 1 - Body.Earth -> 2 - Body.Mars -> 3 - Body.Jupiter -> 4 - Body.Saturn -> 5 - Body.Uranus -> 6 - Body.Neptune -> 7 - else -> throw InvalidBodyException(body) +private fun vsopFormulaCalc(formula: VsopFormula, t: Double, clampAngle: Boolean): Double { + var coord = 0.0 + var tpower = 1.0 + for (series in formula.series) { + var sum = 0.0 + for (term in series.term) + sum += term.amplitude * cos(term.phase + (t * term.frequency)) + coord += + if (clampAngle) + (tpower * sum) % PI2 // improve precision: longitude angles can be hundreds of radians + else + tpower * sum + tpower *= t } + return coord +} + +private fun vsopDistance(model: VsopModel, time: AstroTime) = + vsopFormulaCalc(model.rad, time.julianMillennia(), false) + +private fun vsopRotate(eclip: TerseVector) = + TerseVector( + eclip.x + 0.000000440360*eclip.y - 0.000000190919*eclip.z, + -0.000000479966*eclip.x + 0.917482137087*eclip.y - 0.397776982902*eclip.z, + 0.397776982902*eclip.y + 0.917482137087*eclip.z + ) + +private fun vsopSphereToRect(lon: Double, lat: Double, radius: Double): TerseVector { + val rCosLat = radius * cos(lat) + return TerseVector ( + rCosLat * cos(lon), + rCosLat * sin(lon), + radius * sin(lat) + ) +} + +private fun calcVsop(model: VsopModel, time: AstroTime): AstroVector { + val t = time.julianMillennia() + + // Calculate the VSOP "B" trigonometric series to obtain ecliptic spherical coordinates. + val lon = vsopFormulaCalc(model.lon, t, true) + val lat = vsopFormulaCalc(model.lat, t, false) + val rad = vsopFormulaCalc(model.rad, t, false) + + // Convert ecliptic spherical coordinates to ecliptic Cartesian coordinates. + val eclip = vsopSphereToRect(lon, lat, rad) + + // Convert ecliptic Cartesian coordinates to equatorial Cartesian coordinates. + // Also convert from TerseVector (coordinates only) to AstroVector (coordinates + time). + return vsopRotate(eclip).toAstroVector(time) +} + +private fun vsopDerivCalc(formula: VsopFormula, t: Double): Double { + var tpower = 1.0 // t^s + var dpower = 0.0 // t^(s-1) + var deriv = 0.0 + var s: Int = 0 + for (series in formula.series) { + var sinSum = 0.0 + var cosSum = 0.0 + for (term in series.term) { + val angle = term.phase + (term.frequency * t) + sinSum += term.amplitude * (term.frequency * sin(angle)) + if (s > 0) + cosSum += term.amplitude * cos(angle) + } + deriv += (s * dpower * cosSum) - (tpower * sinSum) + dpower = tpower + tpower *= t + ++s + } + return deriv +} + +private class BodyState( + val tt: Double, + val r: TerseVector, + val v: TerseVector +) + +private fun calcVsopPosVel(model: VsopModel, tt: Double): BodyState { + val t = tt / DAYS_PER_MILLENNIUM + + // Calculate the VSOP "B" trigonometric series to obtain ecliptic spherical coordinates. + val lon = vsopFormulaCalc(model.lon, t, true) + val lat = vsopFormulaCalc(model.lat, t, false) + val rad = vsopFormulaCalc(model.rad, t, false) + + val eclipPos = vsopSphereToRect(lon, lat, rad) + + // Calculate derivatives of spherical coordinates with respect to time. + val dlon = vsopDerivCalc(model.lon, t); // dlon = d(lon) / dt + val dlat = vsopDerivCalc(model.lat, t); // dlat = d(lat) / dt + val drad = vsopDerivCalc(model.rad, t); // drad = d(rad) / dt + + // Use spherical coords and spherical derivatives to calculate + // the velocity vector in rectangular coordinates. + + val coslon = cos(lon); + val sinlon = sin(lon); + val coslat = cos(lat); + val sinlat = sin(lat); + + val vx = ( + + (drad * coslat * coslon) + - (rad * sinlat * coslon * dlat) + - (rad * coslat * sinlon * dlon) + ) + + val vy = ( + + (drad * coslat * sinlon) + - (rad * sinlat * sinlon * dlat) + + (rad * coslat * coslon * dlon) + ) + + val vz = ( + + (drad * sinlat) + + (rad * coslat * dlat) + ) + + // Convert speed units from [AU/millennium] to [AU/day]. + val eclipVel = TerseVector( + vx / DAYS_PER_MILLENNIUM, + vy / DAYS_PER_MILLENNIUM, + vz / DAYS_PER_MILLENNIUM + ) + + // Rotate the vectors from ecliptic to equatorial coordinates. + val equPos = vsopRotate(eclipPos) + val equVel = vsopRotate(eclipVel) + return BodyState(tt, equPos, equVel) +} -private fun vsopModel(body: Body) = vsopTable[vsopTableIndex(body)] //--------------------------------------------------------------------------------------- // Geocentric Moon @@ -2935,6 +3055,82 @@ object Astronomy { val equVec = eclipticToEquatorial(eclVec) return precession(equVec, PrecessDirection.Into2000) } + + private fun helioEarth(time: AstroTime) = calcVsop(vsopTable[2], time) + + /** + * Calculates heliocentric Cartesian coordinates of a body in the J2000 equatorial system. + * + * This function calculates the position of the given celestial body as a vector, + * using the center of the Sun as the origin. The result is expressed as a Cartesian + * vector in the J2000 equatorial system: the coordinates are based on the mean equator + * of the Earth at noon UTC on 1 January 2000. + * + * The position is not corrected for light travel time or aberration. + * This is different from the behavior of #Astronomy.geoVector. + * + * If given an invalid value for `body`, this function will throw an #InvalidBodyException. + * + * @param body + * A body for which to calculate a heliocentric position: + * the Sun, Moon, EMB, SSB, or any of the planets. + * + * @param time + * The date and time for which to calculate the position. + * + * @returns + * The heliocentric position vector of the center of the given body. + */ + fun helioVector(body: Body, time: AstroTime): AstroVector = + when (body) { + Body.Sun -> AstroVector(0.0, 0.0, 0.0, time) + Body.Mercury -> calcVsop(vsopTable[0], time) + Body.Venus -> calcVsop(vsopTable[1], time) + Body.Earth -> calcVsop(vsopTable[2], time) + Body.Mars -> calcVsop(vsopTable[3], time) + Body.Jupiter -> calcVsop(vsopTable[4], time) + Body.Saturn -> calcVsop(vsopTable[5], time) + Body.Uranus -> calcVsop(vsopTable[6], time) + Body.Neptune -> calcVsop(vsopTable[7], time) + Body.Moon -> helioEarth(time) + geoMoon(time) + Body.EMB -> helioEarth(time) + (geoMoon(time) / (1.0 + EARTH_MOON_MASS_RATIO)) + // FIXFIXFIX: add Pluto + // FIXFIXFIX: add Solar System Barycenter + else -> throw InvalidBodyException(body) + } + + /** + * Calculates the distance between a body and the Sun at a given time. + * + * Given a date and time, this function calculates the distance between + * the center of `body` and the center of the Sun, expressed in AU. + * For the planets Mercury through Neptune, this function is significantly + * more efficient than calling #Astronomy.helioVector followed by taking the length + * of the resulting vector. + * + * @param body + * A body for which to calculate a heliocentric distance: + * the Sun, Moon, EMB, SSB, or any of the planets. + * + * @param time + * The date and time for which to calculate the distance. + * + * @returns + * The heliocentric distance in AU. + */ + fun helioDistance(body: Body, time: AstroTime): Double = + when (body) { + Body.Sun -> 0.0 + Body.Mercury -> vsopDistance(vsopTable[0], time) + Body.Venus -> vsopDistance(vsopTable[1], time) + Body.Earth -> vsopDistance(vsopTable[2], time) + Body.Mars -> vsopDistance(vsopTable[3], time) + Body.Jupiter -> vsopDistance(vsopTable[4], time) + Body.Saturn -> vsopDistance(vsopTable[5], time) + Body.Uranus -> vsopDistance(vsopTable[6], time) + Body.Neptune -> vsopDistance(vsopTable[7], time) + else -> helioVector(body, time).length() + } } //======================================================================================= diff --git a/source/csharp/README.md b/source/csharp/README.md index 0d588cb2..70c53873 100644 --- a/source/csharp/README.md +++ b/source/csharp/README.md @@ -574,14 +574,14 @@ movement of the Earth with respect to the rays of light coming from that body. **Calculates the distance between a body and the Sun at a given time.** Given a date and time, this function calculates the distance between -the center of `body` and the center of the Sun. +the center of `body` and the center of the Sun, expressed in AU. For the planets Mercury through Neptune, this function is significantly more efficient than calling [`Astronomy.HelioVector`](#Astronomy.HelioVector) followed by taking the length of the resulting vector. | Type | Parameter | Description | | --- | --- | --- | -| [`Body`](#Body) | `body` | A body for which to calculate a heliocentric distance: the Sun, Moon, or any of the planets. | +| [`Body`](#Body) | `body` | A body for which to calculate a heliocentric distance: the Sun, Moon, EMB, SSB, or any of the planets. | | [`AstroTime`](#AstroTime) | `time` | The date and time for which to calculate the heliocentric distance. | **Returns:** The heliocentric distance in AU. @@ -621,7 +621,7 @@ of the Earth at noon UTC on 1 January 2000. The position is not corrected for light travel time or aberration. This is different from the behavior of [`Astronomy.GeoVector`](#Astronomy.GeoVector). -If given an invalid value for `body`, this function will throw an `ArgumentException`. +If given an invalid value for `body`, this function will throw an [`InvalidBodyException`](#InvalidBodyException). | Type | Parameter | Description | | --- | --- | --- | diff --git a/source/csharp/astronomy.cs b/source/csharp/astronomy.cs index 32651332..75f70e23 100644 --- a/source/csharp/astronomy.cs +++ b/source/csharp/astronomy.cs @@ -5024,9 +5024,12 @@ namespace CosineKitty /// The position is not corrected for light travel time or aberration. /// This is different from the behavior of #Astronomy.GeoVector. /// - /// If given an invalid value for `body`, this function will throw an `ArgumentException`. + /// If given an invalid value for `body`, this function will throw an #InvalidBodyException. /// - /// A body for which to calculate a heliocentric position: the Sun, Moon, EMB, SSB, or any of the planets. + /// + /// A body for which to calculate a heliocentric position: + /// the Sun, Moon, EMB, SSB, or any of the planets. + /// /// The date and time for which to calculate the position. /// A heliocentric position vector of the center of the given body. public static AstroVector HelioVector(Body body, AstroTime time) @@ -5075,14 +5078,14 @@ namespace CosineKitty /// /// /// Given a date and time, this function calculates the distance between - /// the center of `body` and the center of the Sun. + /// the center of `body` and the center of the Sun, expressed in AU. /// For the planets Mercury through Neptune, this function is significantly /// more efficient than calling #Astronomy.HelioVector followed by taking the length /// of the resulting vector. /// /// /// A body for which to calculate a heliocentric distance: - /// the Sun, Moon, or any of the planets. + /// the Sun, Moon, EMB, SSB, or any of the planets. /// /// /// The date and time for which to calculate the heliocentric distance. diff --git a/source/kotlin/doc/-astronomy/helio-distance.md b/source/kotlin/doc/-astronomy/helio-distance.md new file mode 100644 index 00000000..b7c191e6 --- /dev/null +++ b/source/kotlin/doc/-astronomy/helio-distance.md @@ -0,0 +1,19 @@ +//[astronomy](../../../index.md)/[io.github.cosinekitty.astronomy](../index.md)/[Astronomy](index.md)/[helioDistance](helio-distance.md) + +# helioDistance + +[jvm]\ +fun [helioDistance](helio-distance.md)(body: [Body](../-body/index.md), time: [AstroTime](../-astro-time/index.md)): [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html) + +Calculates the distance between a body and the Sun at a given time. + +Given a date and time, this function calculates the distance between the center of body and the center of the Sun, expressed in AU. For the planets Mercury through Neptune, this function is significantly more efficient than calling #Astronomy.helioVector followed by taking the length of the resulting vector. + +## Parameters + +jvm + +| | | +|---|---| +| body | A body for which to calculate a heliocentric distance: the Sun, Moon, EMB, SSB, or any of the planets. | +| time | The date and time for which to calculate the distance. | diff --git a/source/kotlin/doc/-astronomy/helio-vector.md b/source/kotlin/doc/-astronomy/helio-vector.md new file mode 100644 index 00000000..19abeeb5 --- /dev/null +++ b/source/kotlin/doc/-astronomy/helio-vector.md @@ -0,0 +1,23 @@ +//[astronomy](../../../index.md)/[io.github.cosinekitty.astronomy](../index.md)/[Astronomy](index.md)/[helioVector](helio-vector.md) + +# helioVector + +[jvm]\ +fun [helioVector](helio-vector.md)(body: [Body](../-body/index.md), time: [AstroTime](../-astro-time/index.md)): [AstroVector](../-astro-vector/index.md) + +Calculates heliocentric Cartesian coordinates of a body in the J2000 equatorial system. + +This function calculates the position of the given celestial body as a vector, using the center of the Sun as the origin. The result is expressed as a Cartesian vector in the J2000 equatorial system: the coordinates are based on the mean equator of the Earth at noon UTC on 1 January 2000. + +The position is not corrected for light travel time or aberration. This is different from the behavior of #Astronomy.geoVector. + +If given an invalid value for body, this function will throw an #InvalidBodyException. + +## Parameters + +jvm + +| | | +|---|---| +| body | A body for which to calculate a heliocentric position: the Sun, Moon, EMB, SSB, or any of the planets. | +| time | The date and time for which to calculate the position. | diff --git a/source/kotlin/doc/-astronomy/index.md b/source/kotlin/doc/-astronomy/index.md index fa1c289e..51cda8b1 100644 --- a/source/kotlin/doc/-astronomy/index.md +++ b/source/kotlin/doc/-astronomy/index.md @@ -14,6 +14,8 @@ The main container of astronomy calculation functions. | [eclipticGeoMoon](ecliptic-geo-moon.md) | [jvm]
fun [eclipticGeoMoon](ecliptic-geo-moon.md)(time: [AstroTime](../-astro-time/index.md)): [Spherical](../-spherical/index.md)
Calculates spherical ecliptic geocentric position of the Moon. | | [equatorFromVector](equator-from-vector.md) | [jvm]
fun [equatorFromVector](equator-from-vector.md)(vector: [AstroVector](../-astro-vector/index.md)): [Equatorial](../-equatorial/index.md)
Given an equatorial vector, calculates equatorial angular coordinates. | | [geoMoon](geo-moon.md) | [jvm]
fun [geoMoon](geo-moon.md)(time: [AstroTime](../-astro-time/index.md)): [AstroVector](../-astro-vector/index.md)
Calculates equatorial geocentric position of the Moon at a given time. | +| [helioDistance](helio-distance.md) | [jvm]
fun [helioDistance](helio-distance.md)(body: [Body](../-body/index.md), time: [AstroTime](../-astro-time/index.md)): [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)
Calculates the distance between a body and the Sun at a given time. | +| [helioVector](helio-vector.md) | [jvm]
fun [helioVector](helio-vector.md)(body: [Body](../-body/index.md), time: [AstroTime](../-astro-time/index.md)): [AstroVector](../-astro-vector/index.md)
Calculates heliocentric Cartesian coordinates of a body in the J2000 equatorial system. | | [inverseRefractionAngle](inverse-refraction-angle.md) | [jvm]
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)
Calculates the inverse of an atmospheric refraction angle. | | [massProduct](mass-product.md) | [jvm]
fun [massProduct](mass-product.md)(body: [Body](../-body/index.md)): [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)
Returns the product of mass and universal gravitational constant of a Solar System body. | | [refractionAngle](refraction-angle.md) | [jvm]
fun [refractionAngle](refraction-angle.md)(refraction: [Refraction](../-refraction/index.md), altitude: [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)): [Double](https://kotlinlang.org/api/latest/jvm/stdlib/kotlin/-double/index.html)
Calculates the amount of "lift" to an altitude angle caused by atmospheric refraction. | 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 a918c1da..a053a5be 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 @@ -131,6 +131,7 @@ const val KM_PER_AU = 1.4959787069098932e+8 private val TimeZoneUtc = TimeZone.getTimeZone("UTC") private const val DAYS_PER_TROPICAL_YEAR = 365.24217 +private const val DAYS_PER_MILLENNIUM = 365250.0 private const val ASEC360 = 1296000.0 private const val ASEC2RAD = 4.848136811095359935899141e-6 @@ -416,6 +417,7 @@ class AstroTime private constructor( fun addDays(days: Double) = AstroTime(ut + days) internal fun julianCenturies() = tt / 36525.0 + internal fun julianMillennia() = tt / DAYS_PER_MILLENNIUM companion object { private val origin = GregorianCalendar(TimeZoneUtc).also { @@ -464,7 +466,7 @@ internal data class TerseVector(val x: Double, val y: Double, val z: Double) { operator fun div(other: Double) = TerseVector(x / other, y / other, z / other) - val quadrature get() = x * x + y * y + z * z + val quadrature get() = (x * x) + (y * y) + (z * z) val magnitude get() = sqrt(quadrature) companion object { @@ -497,8 +499,7 @@ data class AstroVector( /** * The total distance in AU represented by this vector. */ - fun length() = - sqrt(x*x + y*y + z*z) + fun length() = sqrt((x * x) + (y * y) + (z * z)) private fun verifyIdenticalTimes(otherTime: AstroTime): AstroTime { if (t.tt != otherTime.tt) @@ -1818,20 +1819,139 @@ private class VsopModel( val rad: VsopFormula ) -private fun vsopTableIndex(body: Body): Int = - when (body) { - Body.Mercury -> 0 - Body.Venus -> 1 - Body.Earth -> 2 - Body.Mars -> 3 - Body.Jupiter -> 4 - Body.Saturn -> 5 - Body.Uranus -> 6 - Body.Neptune -> 7 - else -> throw InvalidBodyException(body) +private fun vsopFormulaCalc(formula: VsopFormula, t: Double, clampAngle: Boolean): Double { + var coord = 0.0 + var tpower = 1.0 + for (series in formula.series) { + var sum = 0.0 + for (term in series.term) + sum += term.amplitude * cos(term.phase + (t * term.frequency)) + coord += + if (clampAngle) + (tpower * sum) % PI2 // improve precision: longitude angles can be hundreds of radians + else + tpower * sum + tpower *= t } + return coord +} + +private fun vsopDistance(model: VsopModel, time: AstroTime) = + vsopFormulaCalc(model.rad, time.julianMillennia(), false) + +private fun vsopRotate(eclip: TerseVector) = + TerseVector( + eclip.x + 0.000000440360*eclip.y - 0.000000190919*eclip.z, + -0.000000479966*eclip.x + 0.917482137087*eclip.y - 0.397776982902*eclip.z, + 0.397776982902*eclip.y + 0.917482137087*eclip.z + ) + +private fun vsopSphereToRect(lon: Double, lat: Double, radius: Double): TerseVector { + val rCosLat = radius * cos(lat) + return TerseVector ( + rCosLat * cos(lon), + rCosLat * sin(lon), + radius * sin(lat) + ) +} + +private fun calcVsop(model: VsopModel, time: AstroTime): AstroVector { + val t = time.julianMillennia() + + // Calculate the VSOP "B" trigonometric series to obtain ecliptic spherical coordinates. + val lon = vsopFormulaCalc(model.lon, t, true) + val lat = vsopFormulaCalc(model.lat, t, false) + val rad = vsopFormulaCalc(model.rad, t, false) + + // Convert ecliptic spherical coordinates to ecliptic Cartesian coordinates. + val eclip = vsopSphereToRect(lon, lat, rad) + + // Convert ecliptic Cartesian coordinates to equatorial Cartesian coordinates. + // Also convert from TerseVector (coordinates only) to AstroVector (coordinates + time). + return vsopRotate(eclip).toAstroVector(time) +} + +private fun vsopDerivCalc(formula: VsopFormula, t: Double): Double { + var tpower = 1.0 // t^s + var dpower = 0.0 // t^(s-1) + var deriv = 0.0 + var s: Int = 0 + for (series in formula.series) { + var sinSum = 0.0 + var cosSum = 0.0 + for (term in series.term) { + val angle = term.phase + (term.frequency * t) + sinSum += term.amplitude * (term.frequency * sin(angle)) + if (s > 0) + cosSum += term.amplitude * cos(angle) + } + deriv += (s * dpower * cosSum) - (tpower * sinSum) + dpower = tpower + tpower *= t + ++s + } + return deriv +} + +private class BodyState( + val tt: Double, + val r: TerseVector, + val v: TerseVector +) + +private fun calcVsopPosVel(model: VsopModel, tt: Double): BodyState { + val t = tt / DAYS_PER_MILLENNIUM + + // Calculate the VSOP "B" trigonometric series to obtain ecliptic spherical coordinates. + val lon = vsopFormulaCalc(model.lon, t, true) + val lat = vsopFormulaCalc(model.lat, t, false) + val rad = vsopFormulaCalc(model.rad, t, false) + + val eclipPos = vsopSphereToRect(lon, lat, rad) + + // Calculate derivatives of spherical coordinates with respect to time. + val dlon = vsopDerivCalc(model.lon, t); // dlon = d(lon) / dt + val dlat = vsopDerivCalc(model.lat, t); // dlat = d(lat) / dt + val drad = vsopDerivCalc(model.rad, t); // drad = d(rad) / dt + + // Use spherical coords and spherical derivatives to calculate + // the velocity vector in rectangular coordinates. + + val coslon = cos(lon); + val sinlon = sin(lon); + val coslat = cos(lat); + val sinlat = sin(lat); + + val vx = ( + + (drad * coslat * coslon) + - (rad * sinlat * coslon * dlat) + - (rad * coslat * sinlon * dlon) + ) + + val vy = ( + + (drad * coslat * sinlon) + - (rad * sinlat * sinlon * dlat) + + (rad * coslat * coslon * dlon) + ) + + val vz = ( + + (drad * sinlat) + + (rad * coslat * dlat) + ) + + // Convert speed units from [AU/millennium] to [AU/day]. + val eclipVel = TerseVector( + vx / DAYS_PER_MILLENNIUM, + vy / DAYS_PER_MILLENNIUM, + vz / DAYS_PER_MILLENNIUM + ) + + // Rotate the vectors from ecliptic to equatorial coordinates. + val equPos = vsopRotate(eclipPos) + val equVel = vsopRotate(eclipVel) + return BodyState(tt, equPos, equVel) +} -private fun vsopModel(body: Body) = vsopTable[vsopTableIndex(body)] //--------------------------------------------------------------------------------------- // Geocentric Moon @@ -2935,6 +3055,82 @@ object Astronomy { val equVec = eclipticToEquatorial(eclVec) return precession(equVec, PrecessDirection.Into2000) } + + private fun helioEarth(time: AstroTime) = calcVsop(vsopTable[2], time) + + /** + * Calculates heliocentric Cartesian coordinates of a body in the J2000 equatorial system. + * + * This function calculates the position of the given celestial body as a vector, + * using the center of the Sun as the origin. The result is expressed as a Cartesian + * vector in the J2000 equatorial system: the coordinates are based on the mean equator + * of the Earth at noon UTC on 1 January 2000. + * + * The position is not corrected for light travel time or aberration. + * This is different from the behavior of #Astronomy.geoVector. + * + * If given an invalid value for `body`, this function will throw an #InvalidBodyException. + * + * @param body + * A body for which to calculate a heliocentric position: + * the Sun, Moon, EMB, SSB, or any of the planets. + * + * @param time + * The date and time for which to calculate the position. + * + * @returns + * The heliocentric position vector of the center of the given body. + */ + fun helioVector(body: Body, time: AstroTime): AstroVector = + when (body) { + Body.Sun -> AstroVector(0.0, 0.0, 0.0, time) + Body.Mercury -> calcVsop(vsopTable[0], time) + Body.Venus -> calcVsop(vsopTable[1], time) + Body.Earth -> calcVsop(vsopTable[2], time) + Body.Mars -> calcVsop(vsopTable[3], time) + Body.Jupiter -> calcVsop(vsopTable[4], time) + Body.Saturn -> calcVsop(vsopTable[5], time) + Body.Uranus -> calcVsop(vsopTable[6], time) + Body.Neptune -> calcVsop(vsopTable[7], time) + Body.Moon -> helioEarth(time) + geoMoon(time) + Body.EMB -> helioEarth(time) + (geoMoon(time) / (1.0 + EARTH_MOON_MASS_RATIO)) + // FIXFIXFIX: add Pluto + // FIXFIXFIX: add Solar System Barycenter + else -> throw InvalidBodyException(body) + } + + /** + * Calculates the distance between a body and the Sun at a given time. + * + * Given a date and time, this function calculates the distance between + * the center of `body` and the center of the Sun, expressed in AU. + * For the planets Mercury through Neptune, this function is significantly + * more efficient than calling #Astronomy.helioVector followed by taking the length + * of the resulting vector. + * + * @param body + * A body for which to calculate a heliocentric distance: + * the Sun, Moon, EMB, SSB, or any of the planets. + * + * @param time + * The date and time for which to calculate the distance. + * + * @returns + * The heliocentric distance in AU. + */ + fun helioDistance(body: Body, time: AstroTime): Double = + when (body) { + Body.Sun -> 0.0 + Body.Mercury -> vsopDistance(vsopTable[0], time) + Body.Venus -> vsopDistance(vsopTable[1], time) + Body.Earth -> vsopDistance(vsopTable[2], time) + Body.Mars -> vsopDistance(vsopTable[3], time) + Body.Jupiter -> vsopDistance(vsopTable[4], time) + Body.Saturn -> vsopDistance(vsopTable[5], time) + Body.Uranus -> vsopDistance(vsopTable[6], time) + Body.Neptune -> vsopDistance(vsopTable[7], time) + else -> helioVector(body, time).length() + } } //======================================================================================= 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 a7ef09b3..f01f6a32 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 @@ -17,6 +17,35 @@ private fun tokenize(s: String): List { } class Tests { + private fun checkVector( + vec: AstroVector, + x: Double, y: Double, z: Double, + tolerance: Double, + message: String + ) { + assertTrue(vec.x.isFinite(), "$message: x is not a finite value") + assertTrue(vec.y.isFinite(), "$message: y is not a finite value") + assertTrue(vec.z.isFinite(), "$message: z is not a finite value") + val dx = vec.x - x + val dy = vec.y - y + val dz = vec.z - z + val diff = sqrt(dx*dx + dy*dy + dz*dz) + assertTrue( + diff < tolerance, + "$message: EXCESSIVE ERROR: calc=(${vec.x}, ${vec.y}, ${vec.z}), expected=($x, $y, $z), diff=$diff" + ) + } + + private fun checkScalar( + calculated: Double, + expected: Double, + tolerance: Double, + message: String + ) { + assertTrue(calculated.isFinite(), "$message: scalar is not finite") + val diff = abs(calculated - expected) + assertTrue(diff < tolerance, "$message: EXCESSIVE ERROR: calc=$calculated, expected=$expected, diff=$diff") + } @Test fun `test deltaT calculation`() { @@ -242,11 +271,133 @@ class Tests { assertTrue(abs(drad) < 1.0e-17, "eclipticGeoMoon: excessive distance error $drad") val equVec = Astronomy.geoMoon(time) - val dx = equVec.x - (+0.002674037026701135) - val dy = equVec.y - (-0.0001531610316600666) - val dz = equVec.z - (-0.0003150159927069429) - val diff = sqrt(dx*dx + dy*dy + dz*dz) - assertTrue(diff < 5.43e-20, "geoMoon excessive error = $diff") + checkVector( + equVec, + +0.002674037026701135, -0.0001531610316600666, -0.0003150159927069429, + 5.43e-20, + "geoMoon" + ) + } + + //---------------------------------------------------------------------------------------- + + @Test + fun `Heliocentric vectors and distances`() { + val time = AstroTime(2022, 3, 28, 15, 21, 41.0) + + verifyHelio( + Body.Sun, + time, + 0.0, 0.0, 0.0 + // 0.0, 0.0, 0.0 + ) + + // I used the Python version to calculate reference state vectors. + // This is just a sanity check. + // More exhaustive tests will be added later. + + verifyHelio( + Body.Mercury, + time, + 0.3601225052845418, -0.05702761412960074, -0.06778598705660548 + // 0.0005993053269706367, 0.025455638276077892, 0.013536525517872652 + ) + + verifyHelio( + Body.Venus, + time, + -0.41031775375890733, -0.553797703820094, -0.2232212930914451 + // 0.01652666689705999, -0.010155161581682825, -0.005615571915795434 + ) + + verifyHelio( + Body.Earth, + time, + -0.989323089873633, -0.12148064775730126, -0.05265719320840268 + // 0.001998233870734938, -0.01571097012164241, -0.00681060160402074 + ) + + verifyHelio( + Body.Mars, + time, + 0.33465883914714956, -1.2596467270049816, -0.5868032433717035 + // 0.014131155779203354, 0.0042119182043996745, 0.0015503545422968985 + ) + + verifyHelio( + Body.Jupiter, + time, + 4.843202980566162, -1.0038257126069707, -0.5481395157798546 + // 0.0016395952415457383, 0.007099637074264691, 0.0030031754891337277 + ) + + verifyHelio( + Body.Saturn, + time, + 7.2700610693589, -6.096752520060888, -2.8318384086244084 + // 0.0034786432886192487, 0.003838295762587206, 0.00143581155860176 + ) + + verifyHelio( + Body.Uranus, + time, + 14.159758260480487, 12.631575584016622, 5.332054883582301 + // -0.002763132653154497, 0.0024129846977416014, 0.0010960328116071437 + ) + + verifyHelio( + Body.Neptune, + time, + 29.66860888675136, -3.261664651496572, -2.0721405654341702 + // 0.00038196533635370553, 0.0029105280185167605, 0.0011819337466936254 + ) + +/* + verifyHelio( + Body.Pluto, + time, + 15.377665594383952, -27.85223298004125, -13.32288901996256 + // 0.0028887491551056084, 0.0010313534744949828, -0.0005469417328328084 + ) +*/ + + verifyHelio( + Body.Moon, + time, + -0.9873532819712121, -0.12279866865704744, -0.05347300317660862 + // 0.002381799204808397, -0.015280061848635642, -0.006626397121575363 + ) + + verifyHelio( + Body.EMB, + time, + -0.9892991555540802, -0.1214966624830789, -0.05266710577726255 + // 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( + body: Body, + time: AstroTime, + x: Double, y: Double, z: Double + // , vx: Double, vy: Double, vz: Double + ) { + val tolerance = 1.0e-11 + val calcPos = Astronomy.helioVector(body, time) + checkVector(calcPos, x, y, z, tolerance, "helioVector($body)") + + val calcDist = Astronomy.helioDistance(body, time) + val expectedDist = sqrt(x*x + y*y + z*z) + checkScalar(calcDist, expectedDist, tolerance, "helioDistance($body)") } //----------------------------------------------------------------------------------------