From c36f16e1bec053f6c7fae93c174941f346f509e6 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Thu, 2 Dec 2021 16:11:50 -0500 Subject: [PATCH] PY RotationAxis function. --- demo/browser/astronomy.browser.js | 2 + demo/nodejs/astronomy.js | 2 + demo/nodejs/calendar/astronomy.ts | 2 + demo/python/astronomy.py | 209 ++++++++++++++++++++++++++++++ generate/template/astronomy.cs | 3 + generate/template/astronomy.py | 209 ++++++++++++++++++++++++++++++ generate/template/astronomy.ts | 2 + generate/test.py | 51 ++++++++ source/c/README.md | 4 +- source/c/astronomy.h | 3 + source/csharp/README.md | 3 + source/csharp/astronomy.cs | 3 + source/js/README.md | 4 +- source/js/astronomy.browser.js | 2 + source/js/astronomy.d.ts | 2 + source/js/astronomy.js | 2 + source/js/astronomy.ts | 2 + source/js/esm/astronomy.js | 2 + source/python/README.md | 57 ++++++++ source/python/astronomy.py | 209 ++++++++++++++++++++++++++++++ 20 files changed, 771 insertions(+), 2 deletions(-) diff --git a/demo/browser/astronomy.browser.js b/demo/browser/astronomy.browser.js index bb91dfcd..197fa12f 100644 --- a/demo/browser/astronomy.browser.js +++ b/demo/browser/astronomy.browser.js @@ -8038,6 +8038,8 @@ exports.NextTransit = NextTransit; * The fields `ra`, `dec`, and `spin` correspond to the variables * α0, δ0, and W, respectively, from * [Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015](https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf). + * The field `north` is a unit vector pointing in the direction of the body's north pole. + * It is expressed in the equatorial J2000 system (EQJ). * * @property {number} ra * The J2000 right ascension of the body's north pole direction, in sidereal hours. diff --git a/demo/nodejs/astronomy.js b/demo/nodejs/astronomy.js index 77df22cc..23080d5c 100644 --- a/demo/nodejs/astronomy.js +++ b/demo/nodejs/astronomy.js @@ -8037,6 +8037,8 @@ exports.NextTransit = NextTransit; * The fields `ra`, `dec`, and `spin` correspond to the variables * α0, δ0, and W, respectively, from * [Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015](https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf). + * The field `north` is a unit vector pointing in the direction of the body's north pole. + * It is expressed in the equatorial J2000 system (EQJ). * * @property {number} ra * The J2000 right ascension of the body's north pole direction, in sidereal hours. diff --git a/demo/nodejs/calendar/astronomy.ts b/demo/nodejs/calendar/astronomy.ts index 6c830c53..7f5c8539 100644 --- a/demo/nodejs/calendar/astronomy.ts +++ b/demo/nodejs/calendar/astronomy.ts @@ -8338,6 +8338,8 @@ export function NextTransit(body: Body, prevTransitTime: AstroTime): TransitInfo * The fields `ra`, `dec`, and `spin` correspond to the variables * α0, δ0, and W, respectively, from * [Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015](https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf). + * The field `north` is a unit vector pointing in the direction of the body's north pole. + * It is expressed in the equatorial J2000 system (EQJ). * * @property {number} ra * The J2000 right ascension of the body's north pole direction, in sidereal hours. diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py index dfb9e76e..de9c8dd5 100644 --- a/demo/python/astronomy.py +++ b/demo/python/astronomy.py @@ -8925,3 +8925,212 @@ def Libration(time): bdash = math.degrees(bdash) bdash2 = sigma*math.cos(a) - rho*math.sin(a) return LibrationInfo(bdash + bdash2, ldash + ldash2, mlat, mlon, dist_km, diam_deg) + + +class AxisInfo: + """Information about a body's rotation axis at a given time. + + This structure is returned by #RotationAxis to report + the orientation of a body's rotation axis at a given moment in time. + The axis is specified by the direction in space that the body's north pole + points, using angular equatorial coordinates in the J2000 system (EQJ). + + Thus `ra` is the right ascension, and `dec` is the declination, of the + body's north pole vector at the given moment in time. The north pole + of a body is defined as the pole that lies on the north side of the + [Solar System's invariable plane](https://en.wikipedia.org/wiki/Invariable_plane), + regardless of the body's direction of rotation. + + The `spin` field indicates the angular position of a prime meridian + arbitrarily recommended for the body by the International Astronomical + Union (IAU). + + The fields `ra`, `dec`, and `spin` correspond to the variables + α0, δ0, and W, respectively, from + [Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015](https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf). + + The field `north` is a unit vector pointing in the direction of the body's north pole. + It is expressed in the equatorial J2000 system (EQJ). + + Attributes + ---------- + ra : float + The J2000 right ascension of the body's north pole direction, in sidereal hours. + dec : float + The J2000 declination of the body's north pole direction, in degrees. + spin : float + Rotation angle of the body's prime meridian, in degrees. + north : Vector + A J2000 dimensionless unit vector pointing in the direction of the body's north pole. + """ + def __init__(self, ra, dec, spin, north): + self.ra = ra + self.dec = dec + self.spin = spin + self.north = north + + +def _EarthRotationAxis(time): + # Unlike the other planets, we have a model of precession and nutation + # for the Earth's axis that provides a north pole vector. + # So calculate the vector first, then derive the (RA,DEC) angles from the vector. + + # Start with a north pole vector in equator-of-date coordinates: (0,0,1). + # Convert the vector into J2000 coordinates. + pos2 = _nutation([0, 0, 1], time, _PrecessDir.Into2000) + nvec = _precession(pos2, time, _PrecessDir.Into2000) + north = Vector(nvec[0], nvec[1], nvec[2], time) + + # Derive angular values: right ascension and declination. + equ = EquatorFromVector(north) + + # Use a modified version of the era() function that does not trim to 0..360 degrees. + # This expression is also corrected to give the correct angle at the J2000 epoch. + spin = 190.41375788700253 + (360.9856122880876 * time.ut) + + return AxisInfo(equ.ra, equ.dec, spin, north) + + +def RotationAxis(body, time): + """Calculates information about a body's rotation axis at a given time. + + Calculates the orientation of a body's rotation axis, along with + the rotation angle of its prime meridian, at a given moment in time. + + This function uses formulas standardized by the IAU Working Group + on Cartographics and Rotational Elements 2015 report, as described + in the following document: + + https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf + + See #AxisInfo for more detailed information. + + Parameters: + body : Body + One of the following values: + `Body.Sun`, `Body.Mercury`, `Body.Venus`, `Body.Earth`, `Body.Mars`, + `Body.Jupiter`, `Body.Saturn`, `Body.Uranus`, `Body.Neptune`, `Body.Pluto`. + + time : Time + The time at which to calculate the body's rotation axis. + + Returns + ------- + AxisInfo + The body's north pole direction and angle of its prime meridian. + """ + d = time.tt + T = d / 36525.0 + if body == Body.Sun: + ra = 286.13 + dec = 63.87 + w = 84.176 + (14.1844 * d) + elif body == Body.Mercury: + ra = 281.0103 - (0.0328 * T) + dec = 61.4155 - (0.0049 * T) + w = ( + 329.5988 + + (6.1385108 * d) + + (0.01067257 * math.sin(math.radians(174.7910857 + 4.092335*d))) + - (0.00112309 * math.sin(math.radians(349.5821714 + 8.184670*d))) + - (0.00011040 * math.sin(math.radians(164.3732571 + 12.277005*d))) + - (0.00002539 * math.sin(math.radians(339.1643429 + 16.369340*d))) + - (0.00000571 * math.sin(math.radians(153.9554286 + 20.461675*d))) + ) + elif body == Body.Venus: + ra = 272.76 + dec = 67.16 + w = 160.20 - (1.4813688 * d) + elif body == Body.Earth: + return _EarthRotationAxis(time) + elif body == Body.Mars: + ra = ( + 317.269202 - 0.10927547*T + + 0.000068 * math.sin(math.radians(198.991226 + 19139.4819985*T)) + + 0.000238 * math.sin(math.radians(226.292679 + 38280.8511281*T)) + + 0.000052 * math.sin(math.radians(249.663391 + 57420.7251593*T)) + + 0.000009 * math.sin(math.radians(266.183510 + 76560.6367950*T)) + + 0.419057 * math.sin(math.radians(79.398797 + 0.5042615*T)) + ) + + dec = ( + 54.432516 - 0.05827105*T + + 0.000051*math.cos(math.radians(122.433576 + 19139.9407476*T)) + + 0.000141*math.cos(math.radians(43.058401 + 38280.8753272*T)) + + 0.000031*math.cos(math.radians(57.663379 + 57420.7517205*T)) + + 0.000005*math.cos(math.radians(79.476401 + 76560.6495004*T)) + + 1.591274*math.cos(math.radians(166.325722 + 0.5042615*T)) + ) + + w = ( + 176.049863 + 350.891982443297*d + + 0.000145*math.sin(math.radians(129.071773 + 19140.0328244*T)) + + 0.000157*math.sin(math.radians(36.352167 + 38281.0473591*T)) + + 0.000040*math.sin(math.radians(56.668646 + 57420.9295360*T)) + + 0.000001*math.sin(math.radians(67.364003 + 76560.2552215*T)) + + 0.000001*math.sin(math.radians(104.792680 + 95700.4387578*T)) + + 0.584542*math.sin(math.radians(95.391654 + 0.5042615*T)) + ) + + elif body == Body.Jupiter: + Ja = math.radians(99.360714 + 4850.4046*T) + Jb = math.radians(175.895369 + 1191.9605*T) + Jc = math.radians(300.323162 + 262.5475*T) + Jd = math.radians(114.012305 + 6070.2476*T) + Je = math.radians(49.511251 + 64.3000*T) + + ra = ( + 268.056595 - 0.006499*T + + 0.000117*math.sin(Ja) + + 0.000938*math.sin(Jb) + + 0.001432*math.sin(Jc) + + 0.000030*math.sin(Jd) + + 0.002150*math.sin(Je) + ) + + dec = ( + 64.495303 + 0.002413*T + + 0.000050*math.cos(Ja) + + 0.000404*math.cos(Jb) + + 0.000617*math.cos(Jc) + - 0.000013*math.cos(Jd) + + 0.000926*math.cos(Je) + ) + + w = 284.95 + 870.536*d + + elif body == Body.Saturn: + ra = 40.589 - 0.036*T + dec = 83.537 - 0.004*T + w = 38.90 + 810.7939024*d + + elif body == Body.Uranus: + ra = 257.311 + dec = -15.175 + w = 203.81 - 501.1600928*d + + elif body == Body.Neptune: + N = math.radians(357.85 + 52.316*T) + ra = 299.36 + 0.70*math.sin(N) + dec = 43.46 - 0.51*math.cos(N) + w = 249.978 + 541.1397757*d - 0.48*math.sin(N) + + elif body == Body.Pluto: + ra = 132.993 + dec = -6.163 + w = 302.695 + 56.3625225*d + + else: + raise InvalidBodyError() + + # Calculate the north pole vector using the given angles. + radlat = math.radians(dec) + radlon = math.radians(ra) + rcoslat = math.cos(radlat) + north = Vector( + rcoslat * math.cos(radlon), + rcoslat * math.sin(radlon), + math.sin(radlat), + time + ) + return AxisInfo(ra/15.0, dec, w, north) diff --git a/generate/template/astronomy.cs b/generate/template/astronomy.cs index f7ed074a..e730b8ab 100644 --- a/generate/template/astronomy.cs +++ b/generate/template/astronomy.cs @@ -1436,6 +1436,9 @@ namespace CosineKitty /// The fields `ra`, `dec`, and `spin` correspond to the variables /// α0, δ0, and W, respectively, from /// [Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015](https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf). + /// + /// The field `north` is a unit vector pointing in the direction of the body's north pole. + /// It is expressed in the equatorial J2000 system (EQJ). /// public struct AxisInfo { diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index 36dd8d94..bc127331 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -6432,3 +6432,212 @@ def Libration(time): bdash = math.degrees(bdash) bdash2 = sigma*math.cos(a) - rho*math.sin(a) return LibrationInfo(bdash + bdash2, ldash + ldash2, mlat, mlon, dist_km, diam_deg) + + +class AxisInfo: + """Information about a body's rotation axis at a given time. + + This structure is returned by #RotationAxis to report + the orientation of a body's rotation axis at a given moment in time. + The axis is specified by the direction in space that the body's north pole + points, using angular equatorial coordinates in the J2000 system (EQJ). + + Thus `ra` is the right ascension, and `dec` is the declination, of the + body's north pole vector at the given moment in time. The north pole + of a body is defined as the pole that lies on the north side of the + [Solar System's invariable plane](https://en.wikipedia.org/wiki/Invariable_plane), + regardless of the body's direction of rotation. + + The `spin` field indicates the angular position of a prime meridian + arbitrarily recommended for the body by the International Astronomical + Union (IAU). + + The fields `ra`, `dec`, and `spin` correspond to the variables + α0, δ0, and W, respectively, from + [Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015](https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf). + + The field `north` is a unit vector pointing in the direction of the body's north pole. + It is expressed in the equatorial J2000 system (EQJ). + + Attributes + ---------- + ra : float + The J2000 right ascension of the body's north pole direction, in sidereal hours. + dec : float + The J2000 declination of the body's north pole direction, in degrees. + spin : float + Rotation angle of the body's prime meridian, in degrees. + north : Vector + A J2000 dimensionless unit vector pointing in the direction of the body's north pole. + """ + def __init__(self, ra, dec, spin, north): + self.ra = ra + self.dec = dec + self.spin = spin + self.north = north + + +def _EarthRotationAxis(time): + # Unlike the other planets, we have a model of precession and nutation + # for the Earth's axis that provides a north pole vector. + # So calculate the vector first, then derive the (RA,DEC) angles from the vector. + + # Start with a north pole vector in equator-of-date coordinates: (0,0,1). + # Convert the vector into J2000 coordinates. + pos2 = _nutation([0, 0, 1], time, _PrecessDir.Into2000) + nvec = _precession(pos2, time, _PrecessDir.Into2000) + north = Vector(nvec[0], nvec[1], nvec[2], time) + + # Derive angular values: right ascension and declination. + equ = EquatorFromVector(north) + + # Use a modified version of the era() function that does not trim to 0..360 degrees. + # This expression is also corrected to give the correct angle at the J2000 epoch. + spin = 190.41375788700253 + (360.9856122880876 * time.ut) + + return AxisInfo(equ.ra, equ.dec, spin, north) + + +def RotationAxis(body, time): + """Calculates information about a body's rotation axis at a given time. + + Calculates the orientation of a body's rotation axis, along with + the rotation angle of its prime meridian, at a given moment in time. + + This function uses formulas standardized by the IAU Working Group + on Cartographics and Rotational Elements 2015 report, as described + in the following document: + + https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf + + See #AxisInfo for more detailed information. + + Parameters: + body : Body + One of the following values: + `Body.Sun`, `Body.Mercury`, `Body.Venus`, `Body.Earth`, `Body.Mars`, + `Body.Jupiter`, `Body.Saturn`, `Body.Uranus`, `Body.Neptune`, `Body.Pluto`. + + time : Time + The time at which to calculate the body's rotation axis. + + Returns + ------- + AxisInfo + The body's north pole direction and angle of its prime meridian. + """ + d = time.tt + T = d / 36525.0 + if body == Body.Sun: + ra = 286.13 + dec = 63.87 + w = 84.176 + (14.1844 * d) + elif body == Body.Mercury: + ra = 281.0103 - (0.0328 * T) + dec = 61.4155 - (0.0049 * T) + w = ( + 329.5988 + + (6.1385108 * d) + + (0.01067257 * math.sin(math.radians(174.7910857 + 4.092335*d))) + - (0.00112309 * math.sin(math.radians(349.5821714 + 8.184670*d))) + - (0.00011040 * math.sin(math.radians(164.3732571 + 12.277005*d))) + - (0.00002539 * math.sin(math.radians(339.1643429 + 16.369340*d))) + - (0.00000571 * math.sin(math.radians(153.9554286 + 20.461675*d))) + ) + elif body == Body.Venus: + ra = 272.76 + dec = 67.16 + w = 160.20 - (1.4813688 * d) + elif body == Body.Earth: + return _EarthRotationAxis(time) + elif body == Body.Mars: + ra = ( + 317.269202 - 0.10927547*T + + 0.000068 * math.sin(math.radians(198.991226 + 19139.4819985*T)) + + 0.000238 * math.sin(math.radians(226.292679 + 38280.8511281*T)) + + 0.000052 * math.sin(math.radians(249.663391 + 57420.7251593*T)) + + 0.000009 * math.sin(math.radians(266.183510 + 76560.6367950*T)) + + 0.419057 * math.sin(math.radians(79.398797 + 0.5042615*T)) + ) + + dec = ( + 54.432516 - 0.05827105*T + + 0.000051*math.cos(math.radians(122.433576 + 19139.9407476*T)) + + 0.000141*math.cos(math.radians(43.058401 + 38280.8753272*T)) + + 0.000031*math.cos(math.radians(57.663379 + 57420.7517205*T)) + + 0.000005*math.cos(math.radians(79.476401 + 76560.6495004*T)) + + 1.591274*math.cos(math.radians(166.325722 + 0.5042615*T)) + ) + + w = ( + 176.049863 + 350.891982443297*d + + 0.000145*math.sin(math.radians(129.071773 + 19140.0328244*T)) + + 0.000157*math.sin(math.radians(36.352167 + 38281.0473591*T)) + + 0.000040*math.sin(math.radians(56.668646 + 57420.9295360*T)) + + 0.000001*math.sin(math.radians(67.364003 + 76560.2552215*T)) + + 0.000001*math.sin(math.radians(104.792680 + 95700.4387578*T)) + + 0.584542*math.sin(math.radians(95.391654 + 0.5042615*T)) + ) + + elif body == Body.Jupiter: + Ja = math.radians(99.360714 + 4850.4046*T) + Jb = math.radians(175.895369 + 1191.9605*T) + Jc = math.radians(300.323162 + 262.5475*T) + Jd = math.radians(114.012305 + 6070.2476*T) + Je = math.radians(49.511251 + 64.3000*T) + + ra = ( + 268.056595 - 0.006499*T + + 0.000117*math.sin(Ja) + + 0.000938*math.sin(Jb) + + 0.001432*math.sin(Jc) + + 0.000030*math.sin(Jd) + + 0.002150*math.sin(Je) + ) + + dec = ( + 64.495303 + 0.002413*T + + 0.000050*math.cos(Ja) + + 0.000404*math.cos(Jb) + + 0.000617*math.cos(Jc) + - 0.000013*math.cos(Jd) + + 0.000926*math.cos(Je) + ) + + w = 284.95 + 870.536*d + + elif body == Body.Saturn: + ra = 40.589 - 0.036*T + dec = 83.537 - 0.004*T + w = 38.90 + 810.7939024*d + + elif body == Body.Uranus: + ra = 257.311 + dec = -15.175 + w = 203.81 - 501.1600928*d + + elif body == Body.Neptune: + N = math.radians(357.85 + 52.316*T) + ra = 299.36 + 0.70*math.sin(N) + dec = 43.46 - 0.51*math.cos(N) + w = 249.978 + 541.1397757*d - 0.48*math.sin(N) + + elif body == Body.Pluto: + ra = 132.993 + dec = -6.163 + w = 302.695 + 56.3625225*d + + else: + raise InvalidBodyError() + + # Calculate the north pole vector using the given angles. + radlat = math.radians(dec) + radlon = math.radians(ra) + rcoslat = math.cos(radlat) + north = Vector( + rcoslat * math.cos(radlon), + rcoslat * math.sin(radlon), + math.sin(radlat), + time + ) + return AxisInfo(ra/15.0, dec, w, north) diff --git a/generate/template/astronomy.ts b/generate/template/astronomy.ts index e77edbc9..f64caab2 100644 --- a/generate/template/astronomy.ts +++ b/generate/template/astronomy.ts @@ -6881,6 +6881,8 @@ export function NextTransit(body: Body, prevTransitTime: AstroTime): TransitInfo * The fields `ra`, `dec`, and `spin` correspond to the variables * α0, δ0, and W, respectively, from * [Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015](https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf). + * The field `north` is a unit vector pointing in the direction of the body's north pole. + * It is expressed in the equatorial J2000 system (EQJ). * * @property {number} ra * The J2000 right ascension of the body's north pole direction, in sidereal hours. diff --git a/generate/test.py b/generate/test.py index a9899451..a0d390d8 100755 --- a/generate/test.py +++ b/generate/test.py @@ -2237,8 +2237,59 @@ def Libration(): #----------------------------------------------------------------------------------------------------------- +def Axis(): + if AxisTestBody(astronomy.Body.Sun, 'axis/Sun.txt', 0.0) : return 1 + if AxisTestBody(astronomy.Body.Mercury, 'axis/Mercury.txt', 0.074340) : return 1 + if AxisTestBody(astronomy.Body.Venus, 'axis/Venus.txt', 0.0) : return 1 + if AxisTestBody(astronomy.Body.Earth, 'axis/Earth.txt', 0.000591) : return 1 + if AxisTestBody(astronomy.Body.Mars, 'axis/Mars.txt', 0.075323) : return 1 + if AxisTestBody(astronomy.Body.Jupiter, 'axis/Jupiter.txt', 0.000324) : return 1 + if AxisTestBody(astronomy.Body.Saturn, 'axis/Saturn.txt', 0.000304) : return 1 + if AxisTestBody(astronomy.Body.Uranus, 'axis/Uranus.txt', 0.0) : return 1 + if AxisTestBody(astronomy.Body.Neptune, 'axis/Neptune.txt', 0.000464) : return 1 + if AxisTestBody(astronomy.Body.Pluto, 'axis/Pluto.txt', 0.0) : return 1 + print('PY AxisTest: PASS') + return 0 + +def AxisTestBody(body, filename, arcmin_tolerance): + max_arcmin = 0 + lnum = 0 + count = 0 + found_data = False + with open(filename, 'rt') as infile: + for line in infile: + line = line.strip() + lnum += 1 + if not found_data: + if line == '$$SOE': + found_data = True + else: + if line == '$$EOE': + break + token = line.split() + # [ '1970-Jan-01', '00:00', '2440587.500000000', '281.01954', '61.41577' ] + jd = float(token[2]) + ra = float(token[3]) + dec = float(token[4]) + time = astronomy.Time(jd - 2451545.0) + axis = astronomy.RotationAxis(body, time) + sphere = astronomy.Spherical(dec, ra, 1.0) + north = astronomy.VectorFromSphere(sphere, time) + arcmin = 60.0 * astronomy.AngleBetween(north, axis.north) + if arcmin > max_arcmin: + max_arcmin = arcmin + count += 1 + Debug('PY AxisTestBody({}): {} test cases, max arcmin error = {}.'.format(body, count, max_arcmin)) + if max_arcmin > arcmin_tolerance: + print('PY AxisTestBody({}): EXCESSIVE ERROR = {}'.format(body, max_arcmin)) + return 1 + return 0 + +#----------------------------------------------------------------------------------------------------------- + UnitTests = { 'aberration': Aberration, + 'axis': Axis, 'barystate': BaryState, 'constellation': Constellation, 'elongation': Elongation, diff --git a/source/c/README.md b/source/c/README.md index 9cf85877..15b38c83 100644 --- a/source/c/README.md +++ b/source/c/README.md @@ -3184,7 +3184,9 @@ Thus `ra` is the right ascension, and `dec` is the declination, of the body's no The `spin` field indicates the angular position of a prime meridian arbitrarily recommended for the body by the International Astronomical Union (IAU). -The fields `ra`, `dec`, and `spin` correspond to the variables α0, δ0, and W, respectively, from [Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015](https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf). +The fields `ra`, `dec`, and `spin` correspond to the variables α0, δ0, and W, respectively, from [Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015](https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf). + +The field `north` is a unit vector pointing in the direction of the body's north pole. It is expressed in the equatorial J2000 system (EQJ). | Type | Member | Description | | ---- | ------ | ----------- | diff --git a/source/c/astronomy.h b/source/c/astronomy.h index b5a262f6..587f054e 100644 --- a/source/c/astronomy.h +++ b/source/c/astronomy.h @@ -916,6 +916,9 @@ astro_libration_t; * The fields `ra`, `dec`, and `spin` correspond to the variables * α0, δ0, and W, respectively, from * [Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015](https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf). + * + * The field `north` is a unit vector pointing in the direction of the body's north pole. + * It is expressed in the equatorial J2000 system (EQJ). */ typedef struct { diff --git a/source/csharp/README.md b/source/csharp/README.md index fbaeceea..196467c6 100644 --- a/source/csharp/README.md +++ b/source/csharp/README.md @@ -2052,6 +2052,9 @@ The fields `ra`, `dec`, and `spin` correspond to the variables α0, δ0, and W, respectively, from [Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015](https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf). +The field `north` is a unit vector pointing in the direction of the body's north pole. +It is expressed in the equatorial J2000 system (EQJ). + | Type | Name | Description | | --- | --- | --- | | `double` | `ra` | The J2000 right ascension of the body's north pole direction, in sidereal hours. | diff --git a/source/csharp/astronomy.cs b/source/csharp/astronomy.cs index fef611c1..59cefb21 100644 --- a/source/csharp/astronomy.cs +++ b/source/csharp/astronomy.cs @@ -1436,6 +1436,9 @@ namespace CosineKitty /// The fields `ra`, `dec`, and `spin` correspond to the variables /// α0, δ0, and W, respectively, from /// [Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015](https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf). + /// + /// The field `north` is a unit vector pointing in the direction of the body's north pole. + /// It is expressed in the equatorial J2000 system (EQJ). /// public struct AxisInfo { diff --git a/source/js/README.md b/source/js/README.md index f4373b8b..7bbd78d5 100644 --- a/source/js/README.md +++ b/source/js/README.md @@ -774,7 +774,9 @@ Union (IAU). The fields `ra`, `dec`, and `spin` correspond to the variables α0, δ0, and W, respectively, from -[Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015](https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf). +[Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015](https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf). +The field `north` is a unit vector pointing in the direction of the body's north pole. +It is expressed in the equatorial J2000 system (EQJ). **Properties** | Name | Type | Description | diff --git a/source/js/astronomy.browser.js b/source/js/astronomy.browser.js index bb91dfcd..197fa12f 100644 --- a/source/js/astronomy.browser.js +++ b/source/js/astronomy.browser.js @@ -8038,6 +8038,8 @@ exports.NextTransit = NextTransit; * The fields `ra`, `dec`, and `spin` correspond to the variables * α0, δ0, and W, respectively, from * [Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015](https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf). + * The field `north` is a unit vector pointing in the direction of the body's north pole. + * It is expressed in the equatorial J2000 system (EQJ). * * @property {number} ra * The J2000 right ascension of the body's north pole direction, in sidereal hours. diff --git a/source/js/astronomy.d.ts b/source/js/astronomy.d.ts index ff4f3234..1a14e746 100644 --- a/source/js/astronomy.d.ts +++ b/source/js/astronomy.d.ts @@ -2614,6 +2614,8 @@ export declare function NextTransit(body: Body, prevTransitTime: AstroTime): Tra * The fields `ra`, `dec`, and `spin` correspond to the variables * α0, δ0, and W, respectively, from * [Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015](https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf). + * The field `north` is a unit vector pointing in the direction of the body's north pole. + * It is expressed in the equatorial J2000 system (EQJ). * * @property {number} ra * The J2000 right ascension of the body's north pole direction, in sidereal hours. diff --git a/source/js/astronomy.js b/source/js/astronomy.js index 77df22cc..23080d5c 100644 --- a/source/js/astronomy.js +++ b/source/js/astronomy.js @@ -8037,6 +8037,8 @@ exports.NextTransit = NextTransit; * The fields `ra`, `dec`, and `spin` correspond to the variables * α0, δ0, and W, respectively, from * [Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015](https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf). + * The field `north` is a unit vector pointing in the direction of the body's north pole. + * It is expressed in the equatorial J2000 system (EQJ). * * @property {number} ra * The J2000 right ascension of the body's north pole direction, in sidereal hours. diff --git a/source/js/astronomy.ts b/source/js/astronomy.ts index 6c830c53..7f5c8539 100644 --- a/source/js/astronomy.ts +++ b/source/js/astronomy.ts @@ -8338,6 +8338,8 @@ export function NextTransit(body: Body, prevTransitTime: AstroTime): TransitInfo * The fields `ra`, `dec`, and `spin` correspond to the variables * α0, δ0, and W, respectively, from * [Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015](https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf). + * The field `north` is a unit vector pointing in the direction of the body's north pole. + * It is expressed in the equatorial J2000 system (EQJ). * * @property {number} ra * The J2000 right ascension of the body's north pole direction, in sidereal hours. diff --git a/source/js/esm/astronomy.js b/source/js/esm/astronomy.js index c559be7d..8a7505dd 100644 --- a/source/js/esm/astronomy.js +++ b/source/js/esm/astronomy.js @@ -7928,6 +7928,8 @@ export function NextTransit(body, prevTransitTime) { * The fields `ra`, `dec`, and `spin` correspond to the variables * α0, δ0, and W, respectively, from * [Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015](https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf). + * The field `north` is a unit vector pointing in the direction of the body's north pole. + * It is expressed in the equatorial J2000 system (EQJ). * * @property {number} ra * The J2000 right ascension of the body's north pole direction, in sidereal hours. diff --git a/source/python/README.md b/source/python/README.md index cfc32772..0d694bd9 100644 --- a/source/python/README.md +++ b/source/python/README.md @@ -264,6 +264,38 @@ to iterate through consecutive alternating perigees and apogees. --- + +### class AxisInfo + +**Information about a body's rotation axis at a given time.** + +This structure is returned by [`RotationAxis`](#RotationAxis) to report +the orientation of a body's rotation axis at a given moment in time. +The axis is specified by the direction in space that the body's north pole +points, using angular equatorial coordinates in the J2000 system (EQJ). +Thus `ra` is the right ascension, and `dec` is the declination, of the +body's north pole vector at the given moment in time. The north pole +of a body is defined as the pole that lies on the north side of the +[Solar System's invariable plane](https://en.wikipedia.org/wiki/Invariable_plane), +regardless of the body's direction of rotation. +The `spin` field indicates the angular position of a prime meridian +arbitrarily recommended for the body by the International Astronomical +Union (IAU). +The fields `ra`, `dec`, and `spin` correspond to the variables +α0, δ0, and W, respectively, from +[Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015](https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf). +The field `north` is a unit vector pointing in the direction of the body's north pole. +It is expressed in the equatorial J2000 system (EQJ). + +| Type | Attribute | Description | +| --- | --- | --- | +| `float` | `ra` | The J2000 right ascension of the body's north pole direction, in sidereal hours. | +| `float` | `dec` | The J2000 declination of the body's north pole direction, in degrees. | +| `float` | `spin` | Rotation angle of the body's prime meridian, in degrees. | +| [`Vector`](#Vector) | `north` | A J2000 dimensionless unit vector pointing in the direction of the body's north pole. | + +--- + ### class ConstellationInfo @@ -1943,6 +1975,31 @@ A vector in the orientation specified by `rotation`. --- + +### RotationAxis(body, time) + +**Calculates information about a body's rotation axis at a given time.** + +Calculates the orientation of a body's rotation axis, along with +the rotation angle of its prime meridian, at a given moment in time. +This function uses formulas standardized by the IAU Working Group +on Cartographics and Rotational Elements 2015 report, as described +in the following document: +https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf +See [`AxisInfo`](#AxisInfo) for more detailed information. +Parameters: +body : Body + One of the following values: + `Body.Sun`, `Body.Mercury`, `Body.Venus`, `Body.Earth`, `Body.Mars`, + `Body.Jupiter`, `Body.Saturn`, `Body.Uranus`, `Body.Neptune`, `Body.Pluto`. +time : Time + The time at which to calculate the body's rotation axis. + +### Returns: [`AxisInfo`](#AxisInfo) +The body's north pole direction and angle of its prime meridian. + +--- + ### Rotation_ECL_EQD(time) diff --git a/source/python/astronomy.py b/source/python/astronomy.py index dfb9e76e..de9c8dd5 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -8925,3 +8925,212 @@ def Libration(time): bdash = math.degrees(bdash) bdash2 = sigma*math.cos(a) - rho*math.sin(a) return LibrationInfo(bdash + bdash2, ldash + ldash2, mlat, mlon, dist_km, diam_deg) + + +class AxisInfo: + """Information about a body's rotation axis at a given time. + + This structure is returned by #RotationAxis to report + the orientation of a body's rotation axis at a given moment in time. + The axis is specified by the direction in space that the body's north pole + points, using angular equatorial coordinates in the J2000 system (EQJ). + + Thus `ra` is the right ascension, and `dec` is the declination, of the + body's north pole vector at the given moment in time. The north pole + of a body is defined as the pole that lies on the north side of the + [Solar System's invariable plane](https://en.wikipedia.org/wiki/Invariable_plane), + regardless of the body's direction of rotation. + + The `spin` field indicates the angular position of a prime meridian + arbitrarily recommended for the body by the International Astronomical + Union (IAU). + + The fields `ra`, `dec`, and `spin` correspond to the variables + α0, δ0, and W, respectively, from + [Report of the IAU Working Group on Cartographic Coordinates and Rotational Elements: 2015](https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf). + + The field `north` is a unit vector pointing in the direction of the body's north pole. + It is expressed in the equatorial J2000 system (EQJ). + + Attributes + ---------- + ra : float + The J2000 right ascension of the body's north pole direction, in sidereal hours. + dec : float + The J2000 declination of the body's north pole direction, in degrees. + spin : float + Rotation angle of the body's prime meridian, in degrees. + north : Vector + A J2000 dimensionless unit vector pointing in the direction of the body's north pole. + """ + def __init__(self, ra, dec, spin, north): + self.ra = ra + self.dec = dec + self.spin = spin + self.north = north + + +def _EarthRotationAxis(time): + # Unlike the other planets, we have a model of precession and nutation + # for the Earth's axis that provides a north pole vector. + # So calculate the vector first, then derive the (RA,DEC) angles from the vector. + + # Start with a north pole vector in equator-of-date coordinates: (0,0,1). + # Convert the vector into J2000 coordinates. + pos2 = _nutation([0, 0, 1], time, _PrecessDir.Into2000) + nvec = _precession(pos2, time, _PrecessDir.Into2000) + north = Vector(nvec[0], nvec[1], nvec[2], time) + + # Derive angular values: right ascension and declination. + equ = EquatorFromVector(north) + + # Use a modified version of the era() function that does not trim to 0..360 degrees. + # This expression is also corrected to give the correct angle at the J2000 epoch. + spin = 190.41375788700253 + (360.9856122880876 * time.ut) + + return AxisInfo(equ.ra, equ.dec, spin, north) + + +def RotationAxis(body, time): + """Calculates information about a body's rotation axis at a given time. + + Calculates the orientation of a body's rotation axis, along with + the rotation angle of its prime meridian, at a given moment in time. + + This function uses formulas standardized by the IAU Working Group + on Cartographics and Rotational Elements 2015 report, as described + in the following document: + + https://astropedia.astrogeology.usgs.gov/download/Docs/WGCCRE/WGCCRE2015reprint.pdf + + See #AxisInfo for more detailed information. + + Parameters: + body : Body + One of the following values: + `Body.Sun`, `Body.Mercury`, `Body.Venus`, `Body.Earth`, `Body.Mars`, + `Body.Jupiter`, `Body.Saturn`, `Body.Uranus`, `Body.Neptune`, `Body.Pluto`. + + time : Time + The time at which to calculate the body's rotation axis. + + Returns + ------- + AxisInfo + The body's north pole direction and angle of its prime meridian. + """ + d = time.tt + T = d / 36525.0 + if body == Body.Sun: + ra = 286.13 + dec = 63.87 + w = 84.176 + (14.1844 * d) + elif body == Body.Mercury: + ra = 281.0103 - (0.0328 * T) + dec = 61.4155 - (0.0049 * T) + w = ( + 329.5988 + + (6.1385108 * d) + + (0.01067257 * math.sin(math.radians(174.7910857 + 4.092335*d))) + - (0.00112309 * math.sin(math.radians(349.5821714 + 8.184670*d))) + - (0.00011040 * math.sin(math.radians(164.3732571 + 12.277005*d))) + - (0.00002539 * math.sin(math.radians(339.1643429 + 16.369340*d))) + - (0.00000571 * math.sin(math.radians(153.9554286 + 20.461675*d))) + ) + elif body == Body.Venus: + ra = 272.76 + dec = 67.16 + w = 160.20 - (1.4813688 * d) + elif body == Body.Earth: + return _EarthRotationAxis(time) + elif body == Body.Mars: + ra = ( + 317.269202 - 0.10927547*T + + 0.000068 * math.sin(math.radians(198.991226 + 19139.4819985*T)) + + 0.000238 * math.sin(math.radians(226.292679 + 38280.8511281*T)) + + 0.000052 * math.sin(math.radians(249.663391 + 57420.7251593*T)) + + 0.000009 * math.sin(math.radians(266.183510 + 76560.6367950*T)) + + 0.419057 * math.sin(math.radians(79.398797 + 0.5042615*T)) + ) + + dec = ( + 54.432516 - 0.05827105*T + + 0.000051*math.cos(math.radians(122.433576 + 19139.9407476*T)) + + 0.000141*math.cos(math.radians(43.058401 + 38280.8753272*T)) + + 0.000031*math.cos(math.radians(57.663379 + 57420.7517205*T)) + + 0.000005*math.cos(math.radians(79.476401 + 76560.6495004*T)) + + 1.591274*math.cos(math.radians(166.325722 + 0.5042615*T)) + ) + + w = ( + 176.049863 + 350.891982443297*d + + 0.000145*math.sin(math.radians(129.071773 + 19140.0328244*T)) + + 0.000157*math.sin(math.radians(36.352167 + 38281.0473591*T)) + + 0.000040*math.sin(math.radians(56.668646 + 57420.9295360*T)) + + 0.000001*math.sin(math.radians(67.364003 + 76560.2552215*T)) + + 0.000001*math.sin(math.radians(104.792680 + 95700.4387578*T)) + + 0.584542*math.sin(math.radians(95.391654 + 0.5042615*T)) + ) + + elif body == Body.Jupiter: + Ja = math.radians(99.360714 + 4850.4046*T) + Jb = math.radians(175.895369 + 1191.9605*T) + Jc = math.radians(300.323162 + 262.5475*T) + Jd = math.radians(114.012305 + 6070.2476*T) + Je = math.radians(49.511251 + 64.3000*T) + + ra = ( + 268.056595 - 0.006499*T + + 0.000117*math.sin(Ja) + + 0.000938*math.sin(Jb) + + 0.001432*math.sin(Jc) + + 0.000030*math.sin(Jd) + + 0.002150*math.sin(Je) + ) + + dec = ( + 64.495303 + 0.002413*T + + 0.000050*math.cos(Ja) + + 0.000404*math.cos(Jb) + + 0.000617*math.cos(Jc) + - 0.000013*math.cos(Jd) + + 0.000926*math.cos(Je) + ) + + w = 284.95 + 870.536*d + + elif body == Body.Saturn: + ra = 40.589 - 0.036*T + dec = 83.537 - 0.004*T + w = 38.90 + 810.7939024*d + + elif body == Body.Uranus: + ra = 257.311 + dec = -15.175 + w = 203.81 - 501.1600928*d + + elif body == Body.Neptune: + N = math.radians(357.85 + 52.316*T) + ra = 299.36 + 0.70*math.sin(N) + dec = 43.46 - 0.51*math.cos(N) + w = 249.978 + 541.1397757*d - 0.48*math.sin(N) + + elif body == Body.Pluto: + ra = 132.993 + dec = -6.163 + w = 302.695 + 56.3625225*d + + else: + raise InvalidBodyError() + + # Calculate the north pole vector using the given angles. + radlat = math.radians(dec) + radlon = math.radians(ra) + rcoslat = math.cos(radlat) + north = Vector( + rcoslat * math.cos(radlon), + rcoslat * math.sin(radlon), + math.sin(radlat), + time + ) + return AxisInfo(ra/15.0, dec, w, north)