diff --git a/demo/browser/astronomy.browser.js b/demo/browser/astronomy.browser.js index d6ddf743..8dfb1edd 100644 --- a/demo/browser/astronomy.browser.js +++ b/demo/browser/astronomy.browser.js @@ -5902,7 +5902,7 @@ exports.Refraction = Refraction; * @brief Calculates the inverse of an atmospheric refraction angle. * * Given an observed altitude angle that includes atmospheric refraction, - * calculate the negative angular correction to obtain the unrefracted + * calculates the negative angular correction to obtain the unrefracted * altitude. This is useful for cases where observed horizontal * coordinates are to be converted to another orientation system, * but refraction first must be removed from the observed position. diff --git a/demo/nodejs/astronomy.js b/demo/nodejs/astronomy.js index 55a90de7..6a50104a 100644 --- a/demo/nodejs/astronomy.js +++ b/demo/nodejs/astronomy.js @@ -5901,7 +5901,7 @@ exports.Refraction = Refraction; * @brief Calculates the inverse of an atmospheric refraction angle. * * Given an observed altitude angle that includes atmospheric refraction, - * calculate the negative angular correction to obtain the unrefracted + * calculates the negative angular correction to obtain the unrefracted * altitude. This is useful for cases where observed horizontal * coordinates are to be converted to another orientation system, * but refraction first must be removed from the observed position. diff --git a/demo/nodejs/calendar/astronomy.ts b/demo/nodejs/calendar/astronomy.ts index 0bf5ca01..bd0741a6 100644 --- a/demo/nodejs/calendar/astronomy.ts +++ b/demo/nodejs/calendar/astronomy.ts @@ -6510,7 +6510,7 @@ export function Refraction(refraction: string, altitude: number): number { * @brief Calculates the inverse of an atmospheric refraction angle. * * Given an observed altitude angle that includes atmospheric refraction, - * calculate the negative angular correction to obtain the unrefracted + * calculates the negative angular correction to obtain the unrefracted * altitude. This is useful for cases where observed horizontal * coordinates are to be converted to another orientation system, * but refraction first must be removed from the observed position. diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py index d7713bc3..7f7c87ea 100644 --- a/demo/python/astronomy.py +++ b/demo/python/astronomy.py @@ -5061,7 +5061,7 @@ def InverseRefractionAngle(refraction, bent_altitude): """Calculates the inverse of an atmospheric refraction angle. Given an observed altitude angle that includes atmospheric refraction, - calculate the negative angular correction to obtain the unrefracted + calculates the negative angular correction to obtain the unrefracted altitude. This is useful for cases where observed horizontal coordinates are to be converted to another orientation system, but refraction first must be removed from the observed position. diff --git a/generate/template/astronomy.c b/generate/template/astronomy.c index a014011a..454f8bac 100644 --- a/generate/template/astronomy.c +++ b/generate/template/astronomy.c @@ -7248,7 +7248,7 @@ double Astronomy_Refraction(astro_refraction_t refraction, double altitude) * Calculates the inverse of an atmospheric refraction angle. * * Given an observed altitude angle that includes atmospheric refraction, - * calculate the negative angular correction to obtain the unrefracted + * calculates the negative angular correction to obtain the unrefracted * altitude. This is useful for cases where observed horizontal * coordinates are to be converted to another orientation system, * but refraction first must be removed from the observed position. diff --git a/generate/template/astronomy.cs b/generate/template/astronomy.cs index 1a02b8b9..2880dc3e 100644 --- a/generate/template/astronomy.cs +++ b/generate/template/astronomy.cs @@ -8037,7 +8037,7 @@ $ASTRO_IAU_DATA() /// /// /// Given an observed altitude angle that includes atmospheric refraction, - /// calculate the negative angular correction to obtain the unrefracted + /// calculates the negative angular correction to obtain the unrefracted /// altitude. This is useful for cases where observed horizontal /// coordinates are to be converted to another orientation system, /// but refraction first must be removed from the observed position. diff --git a/generate/template/astronomy.kt b/generate/template/astronomy.kt index 6a656022..9a10029e 100644 --- a/generate/template/astronomy.kt +++ b/generate/template/astronomy.kt @@ -30,10 +30,12 @@ package io.github.cosinekitty.astronomy import java.text.SimpleDateFormat import java.util.* import kotlin.math.absoluteValue -import kotlin.math.roundToLong +import kotlin.math.atan2 import kotlin.math.cos +import kotlin.math.roundToLong import kotlin.math.sin import kotlin.math.sqrt +import kotlin.math.tan /** @@ -41,15 +43,34 @@ import kotlin.math.sqrt */ const val DEG2RAD = 0.017453292519943296 +/** + * Convert an angle expressed in degrees to an angle expressed in radians. + */ +private fun Double.degreesToRadians() = this * DEG2RAD + /** * The factor to convert radians to degrees = 180/pi. */ const val RAD2DEG = 57.295779513082321 +/** + * Convert an angle expressed in radians to an angle expressed in degrees. + */ +private fun Double.radiansToDegrees() = this * RAD2DEG private val TimeZoneUtc = TimeZone.getTimeZone("UTC") private const val DAYS_PER_TROPICAL_YEAR = 365.24217 +private fun Double.withMinDegreeValue(min: Double): Double { + var deg = this + while (deg < min) + deg += 360.0 + while (deg >= min + 360.0) + deg -= 360.0 + return deg +} + +private fun toggleAzimuthDirection(az: Double) = (360.0 - az).withMinDegreeValue(0.0) /** * The enumeration of celestial bodies supported by Astronomy Engine. @@ -312,6 +333,7 @@ internal data class TerseVector(val x: Double, val y: Double, val z: Double) { } } + data class AstroVector( /** * A Cartesian x-coordinate expressed in AU. @@ -361,6 +383,73 @@ data class AstroVector( operator fun div(denom: Double) = AstroVector(x/denom, y/denom, z/denom, t) + + /** + * Converts Cartesian coordinates to spherical coordinates. + * + * Given a Cartesian vector, returns latitude, longitude, and distance. + */ + fun toSpherical(): Spherical { + val xyproj = x*x + y*y + val dist = sqrt(xyproj + z*z) + var lat: Double + var lon: Double + if (xyproj == 0.0) { + if (z == 0.0) { + // Indeterminate coordinates; pos vector has zero length. + throw IllegalArgumentException("Cannot convert zero-length vector to spherical coordinates.") + } + lon = 0.0 + lat = if (z < 0.0) -90.0 else +90.0 + } else { + lon = atan2(y, x).radiansToDegrees().withMinDegreeValue(0.0) + lat = atan2(z, sqrt(xyproj)).radiansToDegrees() + } + return Spherical(lat, lon, dist) + } + + /** + * Given an equatorial vector, calculates equatorial angular coordinates. + */ + fun toEquatorial(): Equatorial { + val sphere = toSpherical() + return Equatorial(sphere.lon / 15.0, sphere.lat, sphere.dist, this) + } + + /** + * Converts Cartesian coordinates to horizontal coordinates. + * + * Given a horizontal Cartesian vector, returns horizontal azimuth and altitude. + * *IMPORTANT:* This function differs from #AstroVector.toSpherical in two ways: + * - `toSpherical` returns a `lon` value that represents azimuth defined counterclockwise + * from north (e.g., west = +90), but this function represents a clockwise rotation + * (e.g., east = +90). The difference is because `toSpherical` is intended + * to preserve the vector "right-hand rule", while this function defines azimuth in a more + * traditional way as used in navigation and cartography. + * - This function optionally corrects for atmospheric refraction, while `toSpherical` + * does not. + * + * The returned object contains the azimuth in `lon`. + * It is measured in degrees clockwise from north: east = +90 degrees, west = +270 degrees. + * + * The altitude is stored in `lat`. + * + * The distance to the observed object is stored in `dist`, + * and is expressed in astronomical units (AU). + * + * @param refraction + * `Refraction.None`: no atmospheric refraction correction is performed. + * `Refraction.Normal`: correct altitude for atmospheric refraction. + * `Refraction.JplHor`: for JPL Horizons compatibility testing only; not recommended for normal use. + */ + fun toHorizontal(refraction: Refraction): Spherical { + val sphere = toSpherical() + return Spherical( + sphere.lat + Astronomy.refractionAngle(refraction, sphere.lat), + toggleAzimuthDirection(sphere.lon), + sphere.dist + ) + } } operator fun Double.times(vec: AstroVector) = @@ -456,7 +545,7 @@ class RotationMatrix( val rot: Array ) { init { - if (rot.size != 3 || rot[0].size != 3 || rot[1].size != 3 || rot[2].size != 3) + if (rot.size != 3 || rot.any { it.size != 3 }) throw IllegalArgumentException("Rotation matrix must be a 3x3 array.") } @@ -566,7 +655,7 @@ class RotationMatrix( if (!angle.isFinite()) throw IllegalArgumentException("Angle must be a finite number.") - val radians = angle * DEG2RAD + val radians = angle.degreesToRadians() val c = cos(radians) val s = sin(radians) @@ -628,7 +717,56 @@ data class Spherical( * Distance in AU. */ val dist: Double -) +) { + /** + * Converts spherical coordinates to Cartesian coordinates. + * + * Given spherical coordinates and a time at which they are valid, + * returns a vector of Cartesian coordinates. The returned value + * includes the time, as required by the type #AstroVector. + * + * @param time + * The time that should be included in the return value. + */ + fun toVector(time: AstroTime): AstroVector { + val radlat = lat.degreesToRadians() + val radlon = lon.degreesToRadians() + val rcoslat = dist * cos(radlat) + return AstroVector( + rcoslat * cos(radlon), + rcoslat * sin(radlon), + dist * sin(radlat), + time + ) + } + + /** + * Given apparent angular horizontal coordinates, calculate the unrefracted horizontal vector. + * + * Assumes `this` contains apparent horizontal coordinates: + * `lat` holds the refracted azimuth angle, + * `lon` holds the azimuth in degrees clockwise from north, + * and `dist` holds the distance from the observer to the object in AU. + * + * @param time + * The date and time of the observation. This is needed because the returned + * #AstroVector requires a valid time value when passed to certain other functions. + * + * @param refraction + * The refraction option used to model atmospheric lensing. See #Astronomy.refractionAngle. + * This specifies how refraction is to be removed from the altitude stored in `this.lat`. + * + * @returns + * A vector in the horizontal system: `x` = north, `y` = west, and `z` = zenith (up). + */ + fun toVectorFromHorizon(time: AstroTime, refraction: Refraction): AstroVector = + Spherical( + lat + Astronomy.inverseRefractionAngle(refraction, lat), + toggleAzimuthDirection(lon), + dist + ) + .toVector(time) +} /** * The location of an observer on (or near) the surface of the Earth. @@ -900,7 +1038,7 @@ object Astronomy { } internal fun universalTime(tt: Double): Double { - // This is the inverse function of TerrestrialTime. + // This is the inverse function of terrestrialTime. // This is an iterative numerical solver, but because // the relationship between UT and TT is almost perfectly linear, // it converges extremely fast (never more than 3 iterations). @@ -915,4 +1053,91 @@ object Astronomy { dt += err } } + + /** + * Calculates the amount of "lift" to an altitude angle caused by atmospheric refraction. + * + * Given an altitude angle and a refraction option, calculates + * the amount of "lift" caused by atmospheric refraction. + * This is the number of degrees higher in the sky an object appears + * due to the lensing of the Earth's atmosphere. + * + * @param refraction + * The option selecting which refraction correction to use. + * If `Refraction.Normal`, uses a well-behaved refraction model that works well for + * all valid values (-90 to +90) of `altitude`. + * If `Refraction.JplHor`, this function returns a compatible value with the JPL Horizons tool. + * If any other value, including `Refraction.None`, this function returns 0. + * + * @param altitude + * An altitude angle in a horizontal coordinate system. Must be a value between -90 and +90. + */ + fun refractionAngle(refraction: Refraction, altitude: Double): Double { + if (altitude < -90.0 || altitude > +90.0) + return 0.0 // no attempt to correct an invalid altitude + + var angle: Double + if (refraction == Refraction.Normal || refraction == Refraction.JplHor) { + // http://extras.springer.com/1999/978-1-4471-0555-8/chap4/horizons/horizons.pdf + // JPL Horizons says it uses refraction algorithm from + // Meeus "Astronomical Algorithms", 1991, p. 101-102. + // I found the following Go implementation: + // https://github.com/soniakeys/meeus/blob/master/v3/refraction/refract.go + // This is a translation from the function "Saemundsson" there. + // I found experimentally that JPL Horizons clamps the angle to 1 degree below the horizon. + // This is important because the 'refr' formula below goes crazy near hd = -5.11. + var hd = altitude + if (hd < -1.0) + hd = -1.0 + + angle = (1.02 / tan((hd+10.3/(hd+5.11))*DEG2RAD)) / 60.0 + + if (refraction == Refraction.Normal && altitude < -1.0) { + // In "normal" mode we gradually reduce refraction toward the nadir + // so that we never get an altitude angle less than -90 degrees. + // When horizon angle is -1 degrees, the factor is exactly 1. + // As altitude approaches -90 (the nadir), the fraction approaches 0 linearly. + angle *= (altitude + 90.0) / 89.0 + } + } else { + // No refraction, or the refraction option is invalid. + angle = 0.0 + } + return angle + } + + /** + * Calculates the inverse of an atmospheric refraction angle. + * + * Given an observed altitude angle that includes atmospheric refraction, + * calculates the negative angular correction to obtain the unrefracted + * altitude. This is useful for cases where observed horizontal + * coordinates are to be converted to another orientation system, + * but refraction first must be removed from the observed position. + * + * @param refraction + * The option selecting which refraction correction to use. + * + * @param bentAltitude + * The apparent altitude that includes atmospheric refraction. + * + * @param + * The angular adjustment in degrees to be added to the + * altitude angle to remove atmospheric lensing. + * This will be less than or equal to zero. + */ + fun inverseRefractionAngle(refraction: Refraction, bentAltitude: Double): Double { + if (bentAltitude < -90.0 || bentAltitude > +90.0) + return 0.0 // no attempt to correct an invalid altitude + + // Find the pre-adjusted altitude whose refraction correction leads to 'altitude'. + var altitude = bentAltitude - refractionAngle(refraction, bentAltitude) + while (true) { + // See how close we got. + var diff = (altitude + refractionAngle(refraction, altitude)) - bentAltitude + if (diff.absoluteValue < 1.0e-14) + return altitude - bentAltitude + altitude -= diff + } + } } diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index 126d24a8..d3f74f2d 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -3019,7 +3019,7 @@ def InverseRefractionAngle(refraction, bent_altitude): """Calculates the inverse of an atmospheric refraction angle. Given an observed altitude angle that includes atmospheric refraction, - calculate the negative angular correction to obtain the unrefracted + calculates the negative angular correction to obtain the unrefracted altitude. This is useful for cases where observed horizontal coordinates are to be converted to another orientation system, but refraction first must be removed from the observed position. diff --git a/generate/template/astronomy.ts b/generate/template/astronomy.ts index 08b24c33..37ebc58a 100644 --- a/generate/template/astronomy.ts +++ b/generate/template/astronomy.ts @@ -5504,7 +5504,7 @@ export function Refraction(refraction: string, altitude: number): number { * @brief Calculates the inverse of an atmospheric refraction angle. * * Given an observed altitude angle that includes atmospheric refraction, - * calculate the negative angular correction to obtain the unrefracted + * calculates the negative angular correction to obtain the unrefracted * altitude. This is useful for cases where observed horizontal * coordinates are to be converted to another orientation system, * but refraction first must be removed from the observed position. diff --git a/source/c/README.md b/source/c/README.md index f883f81e..d5e43594 100644 --- a/source/c/README.md +++ b/source/c/README.md @@ -947,7 +947,7 @@ When the body is Saturn, the returned structure contains a field `ring_tilt` tha -Given an observed altitude angle that includes atmospheric refraction, calculate the negative angular correction to obtain the unrefracted altitude. This is useful for cases where observed horizontal coordinates are to be converted to another orientation system, but refraction first must be removed from the observed position. +Given an observed altitude angle that includes atmospheric refraction, calculates the negative angular correction to obtain the unrefracted altitude. This is useful for cases where observed horizontal coordinates are to be converted to another orientation system, but refraction first must be removed from the observed position. diff --git a/source/c/astronomy.c b/source/c/astronomy.c index 1c076a1b..c064173c 100644 --- a/source/c/astronomy.c +++ b/source/c/astronomy.c @@ -8487,7 +8487,7 @@ double Astronomy_Refraction(astro_refraction_t refraction, double altitude) * Calculates the inverse of an atmospheric refraction angle. * * Given an observed altitude angle that includes atmospheric refraction, - * calculate the negative angular correction to obtain the unrefracted + * calculates the negative angular correction to obtain the unrefracted * altitude. This is useful for cases where observed horizontal * coordinates are to be converted to another orientation system, * but refraction first must be removed from the observed position. diff --git a/source/csharp/README.md b/source/csharp/README.md index dd96e671..90a18fb9 100644 --- a/source/csharp/README.md +++ b/source/csharp/README.md @@ -747,7 +747,7 @@ the rings appear edge-on, and are thus nearly invisible from the Earth. The `rin **Calculates the inverse of an atmospheric refraction angle.** Given an observed altitude angle that includes atmospheric refraction, -calculate the negative angular correction to obtain the unrefracted +calculates the negative angular correction to obtain the unrefracted altitude. This is useful for cases where observed horizontal coordinates are to be converted to another orientation system, but refraction first must be removed from the observed position. diff --git a/source/csharp/astronomy.cs b/source/csharp/astronomy.cs index 17ba1f68..df3dd3f7 100644 --- a/source/csharp/astronomy.cs +++ b/source/csharp/astronomy.cs @@ -9249,7 +9249,7 @@ namespace CosineKitty /// /// /// Given an observed altitude angle that includes atmospheric refraction, - /// calculate the negative angular correction to obtain the unrefracted + /// calculates the negative angular correction to obtain the unrefracted /// altitude. This is useful for cases where observed horizontal /// coordinates are to be converted to another orientation system, /// but refraction first must be removed from the observed position. diff --git a/source/js/README.md b/source/js/README.md index 997c0a66..90d5e65f 100644 --- a/source/js/README.md +++ b/source/js/README.md @@ -2216,7 +2216,7 @@ due to the lensing of the Earth's atmosphere. **Brief**: Calculates the inverse of an atmospheric refraction angle. Given an observed altitude angle that includes atmospheric refraction, -calculate the negative angular correction to obtain the unrefracted +calculates the negative angular correction to obtain the unrefracted altitude. This is useful for cases where observed horizontal coordinates are to be converted to another orientation system, but refraction first must be removed from the observed position. diff --git a/source/js/astronomy.browser.js b/source/js/astronomy.browser.js index d6ddf743..8dfb1edd 100644 --- a/source/js/astronomy.browser.js +++ b/source/js/astronomy.browser.js @@ -5902,7 +5902,7 @@ exports.Refraction = Refraction; * @brief Calculates the inverse of an atmospheric refraction angle. * * Given an observed altitude angle that includes atmospheric refraction, - * calculate the negative angular correction to obtain the unrefracted + * calculates the negative angular correction to obtain the unrefracted * altitude. This is useful for cases where observed horizontal * coordinates are to be converted to another orientation system, * but refraction first must be removed from the observed position. diff --git a/source/js/astronomy.d.ts b/source/js/astronomy.d.ts index 852fb182..9ebb42cb 100644 --- a/source/js/astronomy.d.ts +++ b/source/js/astronomy.d.ts @@ -2002,7 +2002,7 @@ export declare function Refraction(refraction: string, altitude: number): number * @brief Calculates the inverse of an atmospheric refraction angle. * * Given an observed altitude angle that includes atmospheric refraction, - * calculate the negative angular correction to obtain the unrefracted + * calculates the negative angular correction to obtain the unrefracted * altitude. This is useful for cases where observed horizontal * coordinates are to be converted to another orientation system, * but refraction first must be removed from the observed position. diff --git a/source/js/astronomy.js b/source/js/astronomy.js index 55a90de7..6a50104a 100644 --- a/source/js/astronomy.js +++ b/source/js/astronomy.js @@ -5901,7 +5901,7 @@ exports.Refraction = Refraction; * @brief Calculates the inverse of an atmospheric refraction angle. * * Given an observed altitude angle that includes atmospheric refraction, - * calculate the negative angular correction to obtain the unrefracted + * calculates the negative angular correction to obtain the unrefracted * altitude. This is useful for cases where observed horizontal * coordinates are to be converted to another orientation system, * but refraction first must be removed from the observed position. diff --git a/source/js/astronomy.ts b/source/js/astronomy.ts index 0bf5ca01..bd0741a6 100644 --- a/source/js/astronomy.ts +++ b/source/js/astronomy.ts @@ -6510,7 +6510,7 @@ export function Refraction(refraction: string, altitude: number): number { * @brief Calculates the inverse of an atmospheric refraction angle. * * Given an observed altitude angle that includes atmospheric refraction, - * calculate the negative angular correction to obtain the unrefracted + * calculates the negative angular correction to obtain the unrefracted * altitude. This is useful for cases where observed horizontal * coordinates are to be converted to another orientation system, * but refraction first must be removed from the observed position. diff --git a/source/js/esm/astronomy.js b/source/js/esm/astronomy.js index 663a7e57..2f7c0591 100644 --- a/source/js/esm/astronomy.js +++ b/source/js/esm/astronomy.js @@ -5821,7 +5821,7 @@ export function Refraction(refraction, altitude) { * @brief Calculates the inverse of an atmospheric refraction angle. * * Given an observed altitude angle that includes atmospheric refraction, - * calculate the negative angular correction to obtain the unrefracted + * calculates the negative angular correction to obtain the unrefracted * altitude. This is useful for cases where observed horizontal * coordinates are to be converted to another orientation system, * but refraction first must be removed from the observed position. 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 6a656022..9a10029e 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 @@ -30,10 +30,12 @@ package io.github.cosinekitty.astronomy import java.text.SimpleDateFormat import java.util.* import kotlin.math.absoluteValue -import kotlin.math.roundToLong +import kotlin.math.atan2 import kotlin.math.cos +import kotlin.math.roundToLong import kotlin.math.sin import kotlin.math.sqrt +import kotlin.math.tan /** @@ -41,15 +43,34 @@ import kotlin.math.sqrt */ const val DEG2RAD = 0.017453292519943296 +/** + * Convert an angle expressed in degrees to an angle expressed in radians. + */ +private fun Double.degreesToRadians() = this * DEG2RAD + /** * The factor to convert radians to degrees = 180/pi. */ const val RAD2DEG = 57.295779513082321 +/** + * Convert an angle expressed in radians to an angle expressed in degrees. + */ +private fun Double.radiansToDegrees() = this * RAD2DEG private val TimeZoneUtc = TimeZone.getTimeZone("UTC") private const val DAYS_PER_TROPICAL_YEAR = 365.24217 +private fun Double.withMinDegreeValue(min: Double): Double { + var deg = this + while (deg < min) + deg += 360.0 + while (deg >= min + 360.0) + deg -= 360.0 + return deg +} + +private fun toggleAzimuthDirection(az: Double) = (360.0 - az).withMinDegreeValue(0.0) /** * The enumeration of celestial bodies supported by Astronomy Engine. @@ -312,6 +333,7 @@ internal data class TerseVector(val x: Double, val y: Double, val z: Double) { } } + data class AstroVector( /** * A Cartesian x-coordinate expressed in AU. @@ -361,6 +383,73 @@ data class AstroVector( operator fun div(denom: Double) = AstroVector(x/denom, y/denom, z/denom, t) + + /** + * Converts Cartesian coordinates to spherical coordinates. + * + * Given a Cartesian vector, returns latitude, longitude, and distance. + */ + fun toSpherical(): Spherical { + val xyproj = x*x + y*y + val dist = sqrt(xyproj + z*z) + var lat: Double + var lon: Double + if (xyproj == 0.0) { + if (z == 0.0) { + // Indeterminate coordinates; pos vector has zero length. + throw IllegalArgumentException("Cannot convert zero-length vector to spherical coordinates.") + } + lon = 0.0 + lat = if (z < 0.0) -90.0 else +90.0 + } else { + lon = atan2(y, x).radiansToDegrees().withMinDegreeValue(0.0) + lat = atan2(z, sqrt(xyproj)).radiansToDegrees() + } + return Spherical(lat, lon, dist) + } + + /** + * Given an equatorial vector, calculates equatorial angular coordinates. + */ + fun toEquatorial(): Equatorial { + val sphere = toSpherical() + return Equatorial(sphere.lon / 15.0, sphere.lat, sphere.dist, this) + } + + /** + * Converts Cartesian coordinates to horizontal coordinates. + * + * Given a horizontal Cartesian vector, returns horizontal azimuth and altitude. + * *IMPORTANT:* This function differs from #AstroVector.toSpherical in two ways: + * - `toSpherical` returns a `lon` value that represents azimuth defined counterclockwise + * from north (e.g., west = +90), but this function represents a clockwise rotation + * (e.g., east = +90). The difference is because `toSpherical` is intended + * to preserve the vector "right-hand rule", while this function defines azimuth in a more + * traditional way as used in navigation and cartography. + * - This function optionally corrects for atmospheric refraction, while `toSpherical` + * does not. + * + * The returned object contains the azimuth in `lon`. + * It is measured in degrees clockwise from north: east = +90 degrees, west = +270 degrees. + * + * The altitude is stored in `lat`. + * + * The distance to the observed object is stored in `dist`, + * and is expressed in astronomical units (AU). + * + * @param refraction + * `Refraction.None`: no atmospheric refraction correction is performed. + * `Refraction.Normal`: correct altitude for atmospheric refraction. + * `Refraction.JplHor`: for JPL Horizons compatibility testing only; not recommended for normal use. + */ + fun toHorizontal(refraction: Refraction): Spherical { + val sphere = toSpherical() + return Spherical( + sphere.lat + Astronomy.refractionAngle(refraction, sphere.lat), + toggleAzimuthDirection(sphere.lon), + sphere.dist + ) + } } operator fun Double.times(vec: AstroVector) = @@ -456,7 +545,7 @@ class RotationMatrix( val rot: Array ) { init { - if (rot.size != 3 || rot[0].size != 3 || rot[1].size != 3 || rot[2].size != 3) + if (rot.size != 3 || rot.any { it.size != 3 }) throw IllegalArgumentException("Rotation matrix must be a 3x3 array.") } @@ -566,7 +655,7 @@ class RotationMatrix( if (!angle.isFinite()) throw IllegalArgumentException("Angle must be a finite number.") - val radians = angle * DEG2RAD + val radians = angle.degreesToRadians() val c = cos(radians) val s = sin(radians) @@ -628,7 +717,56 @@ data class Spherical( * Distance in AU. */ val dist: Double -) +) { + /** + * Converts spherical coordinates to Cartesian coordinates. + * + * Given spherical coordinates and a time at which they are valid, + * returns a vector of Cartesian coordinates. The returned value + * includes the time, as required by the type #AstroVector. + * + * @param time + * The time that should be included in the return value. + */ + fun toVector(time: AstroTime): AstroVector { + val radlat = lat.degreesToRadians() + val radlon = lon.degreesToRadians() + val rcoslat = dist * cos(radlat) + return AstroVector( + rcoslat * cos(radlon), + rcoslat * sin(radlon), + dist * sin(radlat), + time + ) + } + + /** + * Given apparent angular horizontal coordinates, calculate the unrefracted horizontal vector. + * + * Assumes `this` contains apparent horizontal coordinates: + * `lat` holds the refracted azimuth angle, + * `lon` holds the azimuth in degrees clockwise from north, + * and `dist` holds the distance from the observer to the object in AU. + * + * @param time + * The date and time of the observation. This is needed because the returned + * #AstroVector requires a valid time value when passed to certain other functions. + * + * @param refraction + * The refraction option used to model atmospheric lensing. See #Astronomy.refractionAngle. + * This specifies how refraction is to be removed from the altitude stored in `this.lat`. + * + * @returns + * A vector in the horizontal system: `x` = north, `y` = west, and `z` = zenith (up). + */ + fun toVectorFromHorizon(time: AstroTime, refraction: Refraction): AstroVector = + Spherical( + lat + Astronomy.inverseRefractionAngle(refraction, lat), + toggleAzimuthDirection(lon), + dist + ) + .toVector(time) +} /** * The location of an observer on (or near) the surface of the Earth. @@ -900,7 +1038,7 @@ object Astronomy { } internal fun universalTime(tt: Double): Double { - // This is the inverse function of TerrestrialTime. + // This is the inverse function of terrestrialTime. // This is an iterative numerical solver, but because // the relationship between UT and TT is almost perfectly linear, // it converges extremely fast (never more than 3 iterations). @@ -915,4 +1053,91 @@ object Astronomy { dt += err } } + + /** + * Calculates the amount of "lift" to an altitude angle caused by atmospheric refraction. + * + * Given an altitude angle and a refraction option, calculates + * the amount of "lift" caused by atmospheric refraction. + * This is the number of degrees higher in the sky an object appears + * due to the lensing of the Earth's atmosphere. + * + * @param refraction + * The option selecting which refraction correction to use. + * If `Refraction.Normal`, uses a well-behaved refraction model that works well for + * all valid values (-90 to +90) of `altitude`. + * If `Refraction.JplHor`, this function returns a compatible value with the JPL Horizons tool. + * If any other value, including `Refraction.None`, this function returns 0. + * + * @param altitude + * An altitude angle in a horizontal coordinate system. Must be a value between -90 and +90. + */ + fun refractionAngle(refraction: Refraction, altitude: Double): Double { + if (altitude < -90.0 || altitude > +90.0) + return 0.0 // no attempt to correct an invalid altitude + + var angle: Double + if (refraction == Refraction.Normal || refraction == Refraction.JplHor) { + // http://extras.springer.com/1999/978-1-4471-0555-8/chap4/horizons/horizons.pdf + // JPL Horizons says it uses refraction algorithm from + // Meeus "Astronomical Algorithms", 1991, p. 101-102. + // I found the following Go implementation: + // https://github.com/soniakeys/meeus/blob/master/v3/refraction/refract.go + // This is a translation from the function "Saemundsson" there. + // I found experimentally that JPL Horizons clamps the angle to 1 degree below the horizon. + // This is important because the 'refr' formula below goes crazy near hd = -5.11. + var hd = altitude + if (hd < -1.0) + hd = -1.0 + + angle = (1.02 / tan((hd+10.3/(hd+5.11))*DEG2RAD)) / 60.0 + + if (refraction == Refraction.Normal && altitude < -1.0) { + // In "normal" mode we gradually reduce refraction toward the nadir + // so that we never get an altitude angle less than -90 degrees. + // When horizon angle is -1 degrees, the factor is exactly 1. + // As altitude approaches -90 (the nadir), the fraction approaches 0 linearly. + angle *= (altitude + 90.0) / 89.0 + } + } else { + // No refraction, or the refraction option is invalid. + angle = 0.0 + } + return angle + } + + /** + * Calculates the inverse of an atmospheric refraction angle. + * + * Given an observed altitude angle that includes atmospheric refraction, + * calculates the negative angular correction to obtain the unrefracted + * altitude. This is useful for cases where observed horizontal + * coordinates are to be converted to another orientation system, + * but refraction first must be removed from the observed position. + * + * @param refraction + * The option selecting which refraction correction to use. + * + * @param bentAltitude + * The apparent altitude that includes atmospheric refraction. + * + * @param + * The angular adjustment in degrees to be added to the + * altitude angle to remove atmospheric lensing. + * This will be less than or equal to zero. + */ + fun inverseRefractionAngle(refraction: Refraction, bentAltitude: Double): Double { + if (bentAltitude < -90.0 || bentAltitude > +90.0) + return 0.0 // no attempt to correct an invalid altitude + + // Find the pre-adjusted altitude whose refraction correction leads to 'altitude'. + var altitude = bentAltitude - refractionAngle(refraction, bentAltitude) + while (true) { + // See how close we got. + var diff = (altitude + refractionAngle(refraction, altitude)) - bentAltitude + if (diff.absoluteValue < 1.0e-14) + return altitude - bentAltitude + altitude -= diff + } + } } diff --git a/source/python/README.md b/source/python/README.md index d6955834..45b38039 100644 --- a/source/python/README.md +++ b/source/python/README.md @@ -1618,7 +1618,7 @@ the rings appear edge-on, and are thus nearly invisible from the Earth. The `rin **Calculates the inverse of an atmospheric refraction angle.** Given an observed altitude angle that includes atmospheric refraction, -calculate the negative angular correction to obtain the unrefracted +calculates the negative angular correction to obtain the unrefracted altitude. This is useful for cases where observed horizontal coordinates are to be converted to another orientation system, but refraction first must be removed from the observed position. diff --git a/source/python/astronomy/astronomy.py b/source/python/astronomy/astronomy.py index d7713bc3..7f7c87ea 100644 --- a/source/python/astronomy/astronomy.py +++ b/source/python/astronomy/astronomy.py @@ -5061,7 +5061,7 @@ def InverseRefractionAngle(refraction, bent_altitude): """Calculates the inverse of an atmospheric refraction angle. Given an observed altitude angle that includes atmospheric refraction, - calculate the negative angular correction to obtain the unrefracted + calculates the negative angular correction to obtain the unrefracted altitude. This is useful for cases where observed horizontal coordinates are to be converted to another orientation system, but refraction first must be removed from the observed position.