From ab8d8ba295a0c620d6d8f2a9745812dcfff32936 Mon Sep 17 00:00:00 2001 From: Don Cross Date: Sun, 15 Dec 2019 11:34:16 -0500 Subject: [PATCH] Python: Starting to port new rotation code. Split _precession into _precession_rot/_precession. Split _nutation into _nutation_rot/_nutation. Removed trailing whitespace from readme. --- README.md | 12 +++---- generate/template/astronomy.py | 64 ++++++++++++++++++++++++---------- source/python/README.md | 12 +++++++ source/python/astronomy.py | 64 ++++++++++++++++++++++++---------- 4 files changed, 110 insertions(+), 42 deletions(-) diff --git a/README.md b/README.md index 942a73c6..ffc4eb8e 100644 --- a/README.md +++ b/README.md @@ -37,7 +37,7 @@ ### Overview -The Astronomy Engine is a suite of open source libraries for calculating positions of +The Astronomy Engine is a suite of open source libraries for calculating positions of the Sun, Moon, and planets, and for predicting interesting events like oppositions, conjunctions, rise and set times, lunar phases, and more. @@ -47,10 +47,10 @@ Function and type names are uniform across all the supported languages. The Astronomy Engine is designed to be small, fast, and accurate to within ±1 arcminute. It is based on the authoritative and well-tested models [VSOP87](https://en.wikipedia.org/wiki/VSOP_(planets)) -and +and [NOVAS C 3.1](https://aa.usno.navy.mil/software/novas/novas_c/novasc_info.php). -These libraries are rigorously unit-tested against NOVAS, +These libraries are rigorously unit-tested against NOVAS, [JPL Horizons](https://ssd.jpl.nasa.gov/horizons.cgi), and other reliable sources of ephemeris data. Calculations are also verified to be identical among all the supported programming languages. @@ -63,13 +63,13 @@ Calculations are also verified to be identical among all the supported programmi - Provides heliocentric and geocentric Cartesian vectors of all the above bodies. -- Determines apparent horizon-based positions for an observer anywhere on the Earth, - given that observer's latitude, longitude, and elevation in meters. +- Determines apparent horizon-based positions for an observer anywhere on the Earth, + given that observer's latitude, longitude, and elevation in meters. Optionally corrects for atmospheric refraction. - Calculates rise, set, and culmination times of Sun, Moon, and planets. -- Finds date and time of Moon phases: new, first quarter, full, third quarter +- Finds date and time of Moon phases: new, first quarter, full, third quarter (or anywhere in between as expressed in degrees of ecliptic longitude). - Predicts lunar apogee and perigee dates, times, and distances. diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index 49d20277..e4b5638e 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -435,6 +435,18 @@ class Observer: self.longitude = longitude self.height = height +class RotationMatrix: + """Contains a rotation matrix that can be used to transform one + coordinate system into another. + + Parameters + ---------- + rot : float[3][3] + A normalized 3x3 rotation matrix. + """ + def __init__(self, rot): + self.rot = rot + class _iau2000b: def __init__(self, time): t = time.tt / 36525.0 @@ -480,7 +492,7 @@ def _ecl2equ_vec(time, ecl): ecl[1]*sin_obl + ecl[2]*cos_obl ] -def _precession(tt1, pos1, tt2): +def _precession_rot(tt1, tt2): eps0 = 84381.406 if tt1 != 0 and tt2 != 0: raise Error('One of (tt1, tt2) must be zero.') @@ -531,17 +543,25 @@ def _precession(tt1, pos1, tt2): zz = -sc * cb * sa + cc * ca if tt2 == 0.0: # Perform rotation from other epoch to J2000.0. - return [ - xx * pos1[0] + xy * pos1[1] + xz * pos1[2], - yx * pos1[0] + yy * pos1[1] + yz * pos1[2], - zx * pos1[0] + zy * pos1[1] + zz * pos1[2] - ] + return RotationMatrix([ + [xx, yx, zx], + [xy, yy, zy], + [xz, yz, zz] + ]) # Perform rotation from J2000.0 to other epoch. + return RotationMatrix([ + [xx, xy, xz], + [yx, yy, yz], + [zx, zy, zz] + ]) + +def _precession(tt1, pos1, tt2): + r = _precession_rot(tt1, tt2) return [ - xx * pos1[0] + yx * pos1[1] + zx * pos1[2], - xy * pos1[0] + yy * pos1[1] + zy * pos1[2], - xz * pos1[0] + yz * pos1[1] + zz * pos1[2] + r.rot[0][0]*pos1[0] + r.rot[1][0]*pos1[1] + r.rot[2][0]*pos1[2], + r.rot[0][1]*pos1[0] + r.rot[1][1]*pos1[1] + r.rot[2][1]*pos1[2], + r.rot[0][2]*pos1[0] + r.rot[1][2]*pos1[1] + r.rot[2][2]*pos1[2] ] class Equatorial: @@ -586,7 +606,7 @@ def _vector2radec(pos): return Equatorial(ra, dec, dist) -def _nutation(time, direction, inpos): +def _nutation_rot(time, direction): tilt = time._etilt() oblm = math.radians(tilt.mobl) oblt = math.radians(tilt.tobl) @@ -610,17 +630,25 @@ def _nutation(time, direction, inpos): if direction == 0: # forward rotation - return [ - xx * inpos[0] + yx * inpos[1] + zx * inpos[2], - xy * inpos[0] + yy * inpos[1] + zy * inpos[2], - xz * inpos[0] + yz * inpos[1] + zz * inpos[2] - ] + return RotationMatrix([ + [xx, xy, xz], + [yx, yy, yz], + [zx, zy, zz] + ]) # inverse rotation + return RotationMatrix([ + [xx, yx, zx], + [xy, yy, zy], + [xz, yz, zz] + ]) + +def _nutation(time, direction, pos): + r = _nutation_rot(time, direction) return [ - xx * inpos[0] + xy * inpos[1] + xz * inpos[2], - yx * inpos[0] + yy * inpos[1] + yz * inpos[2], - zx * inpos[0] + zy * inpos[1] + zz * inpos[2] + r.rot[0][0]*pos[0] + r.rot[1][0]*pos[1] + r.rot[2][0]*pos[2], + r.rot[0][1]*pos[0] + r.rot[1][1]*pos[1] + r.rot[2][1]*pos[2], + r.rot[0][2]*pos[0] + r.rot[1][2]*pos[1] + r.rot[2][2]*pos[2] ] def _era(time): # Earth Rotation Angle diff --git a/source/python/README.md b/source/python/README.md index b687e8c1..4fb551c8 100644 --- a/source/python/README.md +++ b/source/python/README.md @@ -167,6 +167,18 @@ when the lunar quarter event happens. --- + +### class RotationMatrix + +Contains a rotation matrix that can be used to transform one +coordinate system into another. + +| Type | Parameter | Description | +| --- | --- | --- | +| `float[3][3]` | `rot` | A normalized 3x3 rotation matrix. | + +--- + ### class SeasonInfo diff --git a/source/python/astronomy.py b/source/python/astronomy.py index 7e38965e..293074e3 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -526,6 +526,18 @@ class Observer: self.longitude = longitude self.height = height +class RotationMatrix: + """Contains a rotation matrix that can be used to transform one + coordinate system into another. + + Parameters + ---------- + rot : float[3][3] + A normalized 3x3 rotation matrix. + """ + def __init__(self, rot): + self.rot = rot + class _iau2000b: def __init__(self, time): t = time.tt / 36525.0 @@ -1107,7 +1119,7 @@ def _ecl2equ_vec(time, ecl): ecl[1]*sin_obl + ecl[2]*cos_obl ] -def _precession(tt1, pos1, tt2): +def _precession_rot(tt1, tt2): eps0 = 84381.406 if tt1 != 0 and tt2 != 0: raise Error('One of (tt1, tt2) must be zero.') @@ -1158,17 +1170,25 @@ def _precession(tt1, pos1, tt2): zz = -sc * cb * sa + cc * ca if tt2 == 0.0: # Perform rotation from other epoch to J2000.0. - return [ - xx * pos1[0] + xy * pos1[1] + xz * pos1[2], - yx * pos1[0] + yy * pos1[1] + yz * pos1[2], - zx * pos1[0] + zy * pos1[1] + zz * pos1[2] - ] + return RotationMatrix([ + [xx, yx, zx], + [xy, yy, zy], + [xz, yz, zz] + ]) # Perform rotation from J2000.0 to other epoch. + return RotationMatrix([ + [xx, xy, xz], + [yx, yy, yz], + [zx, zy, zz] + ]) + +def _precession(tt1, pos1, tt2): + r = _precession_rot(tt1, tt2) return [ - xx * pos1[0] + yx * pos1[1] + zx * pos1[2], - xy * pos1[0] + yy * pos1[1] + zy * pos1[2], - xz * pos1[0] + yz * pos1[1] + zz * pos1[2] + r.rot[0][0]*pos1[0] + r.rot[1][0]*pos1[1] + r.rot[2][0]*pos1[2], + r.rot[0][1]*pos1[0] + r.rot[1][1]*pos1[1] + r.rot[2][1]*pos1[2], + r.rot[0][2]*pos1[0] + r.rot[1][2]*pos1[1] + r.rot[2][2]*pos1[2] ] class Equatorial: @@ -1213,7 +1233,7 @@ def _vector2radec(pos): return Equatorial(ra, dec, dist) -def _nutation(time, direction, inpos): +def _nutation_rot(time, direction): tilt = time._etilt() oblm = math.radians(tilt.mobl) oblt = math.radians(tilt.tobl) @@ -1237,17 +1257,25 @@ def _nutation(time, direction, inpos): if direction == 0: # forward rotation - return [ - xx * inpos[0] + yx * inpos[1] + zx * inpos[2], - xy * inpos[0] + yy * inpos[1] + zy * inpos[2], - xz * inpos[0] + yz * inpos[1] + zz * inpos[2] - ] + return RotationMatrix([ + [xx, xy, xz], + [yx, yy, yz], + [zx, zy, zz] + ]) # inverse rotation + return RotationMatrix([ + [xx, yx, zx], + [xy, yy, zy], + [xz, yz, zz] + ]) + +def _nutation(time, direction, pos): + r = _nutation_rot(time, direction) return [ - xx * inpos[0] + xy * inpos[1] + xz * inpos[2], - yx * inpos[0] + yy * inpos[1] + yz * inpos[2], - zx * inpos[0] + zy * inpos[1] + zz * inpos[2] + r.rot[0][0]*pos[0] + r.rot[1][0]*pos[1] + r.rot[2][0]*pos[2], + r.rot[0][1]*pos[0] + r.rot[1][1]*pos[1] + r.rot[2][1]*pos[2], + r.rot[0][2]*pos[0] + r.rot[1][2]*pos[1] + r.rot[2][2]*pos[2] ] def _era(time): # Earth Rotation Angle