From a97fc7da9c30028c790fc78cd809b7d413e5fe7e Mon Sep 17 00:00:00 2001 From: Don Cross Date: Sat, 27 Mar 2021 05:19:27 -0400 Subject: [PATCH] Ported IdentityMatrix, Pivot functions to Python. Added tests and camera demo. --- demo/csharp/camera/camera.cs | 14 ++++ demo/python/astronomy.py | 107 ++++++++++++++++++++++++++-- demo/python/camera.py | 97 +++++++++++++++++++++++++ demo/python/test/.gitignore | 1 + demo/python/test/camera_correct.txt | 6 ++ demo/python/test/test | 6 +- generate/template/astronomy.py | 107 ++++++++++++++++++++++++++-- generate/test.py | 56 +++++++++++++++ source/python/README.md | 41 +++++++++++ source/python/astronomy.py | 107 ++++++++++++++++++++++++++-- 10 files changed, 520 insertions(+), 22 deletions(-) create mode 100755 demo/python/camera.py create mode 100644 demo/python/test/camera_correct.txt diff --git a/demo/csharp/camera/camera.cs b/demo/csharp/camera/camera.cs index 66d31e1c..b31b8498 100644 --- a/demo/csharp/camera/camera.cs +++ b/demo/csharp/camera/camera.cs @@ -2,6 +2,20 @@ using demo_helper; using CosineKitty; +// camera.cs - by Don Cross - 2021-03-27 +// +// Example C# program for Astronomy Engine: +// https://github.com/cosinekitty/astronomy +// +// Suppose you want to photograph the Moon, +// and you want to know what it will look like in the photo. +// Given a location on the Earth, and a date/time, +// this program calculates the orientation of the sunlit +// side of the Moon with respect to the top of your +// photo image. It assumes the camera faces directly +// toward the Moon's azimuth and tilts upward to its +// altitude angle above the horizon. + namespace camera { class Program diff --git a/demo/python/astronomy.py b/demo/python/astronomy.py index fc7c6679..56978617 100644 --- a/demo/python/astronomy.py +++ b/demo/python/astronomy.py @@ -1323,13 +1323,20 @@ class Equatorial: Declination in degrees. dist : float Distance to the celestial body in AU. + vec : Vector + The equatorial coordinates in cartesian form, using AU distance units. + x = direction of the March equinox, + y = direction of the June solstice, + z = north. """ - def __init__(self, ra, dec, dist): + def __init__(self, ra, dec, dist, vec): self.ra = ra self.dec = dec self.dist = dist + self.vec = vec -def _vector2radec(pos): + +def _vector2radec(pos, time): xyproj = pos[0]*pos[0] + pos[1]*pos[1] dist = math.sqrt(xyproj + pos[2]*pos[2]) if xyproj == 0.0: @@ -1346,7 +1353,8 @@ def _vector2radec(pos): if ra < 0: ra += 24 dec = math.degrees(math.atan2(pos[2], math.sqrt(xyproj))) - return Equatorial(ra, dec, dist) + vec = Vector(pos[0], pos[1], pos[2], time) + return Equatorial(ra, dec, dist, vec) def _nutation_rot(time, direction): @@ -3765,10 +3773,10 @@ def Equator(body, time, observer, ofdate, aberration): gc.z - gc_observer[2] ] if not ofdate: - return _vector2radec(j2000) + return _vector2radec(j2000, time) temp = _precession(0, j2000, time.tt) datevect = _nutation(time, 0, temp) - return _vector2radec(datevect) + return _vector2radec(datevect, time) @enum.unique class Refraction(enum.Enum): @@ -5745,7 +5753,7 @@ def EquatorFromVector(vec): Angular coordinates expressed in the same equatorial system as `vec`. """ sphere = SphereFromVector(vec) - return Equatorial(sphere.lon / 15.0, sphere.lat, sphere.dist) + return Equatorial(sphere.lon / 15.0, sphere.lat, sphere.dist, vec) def SphereFromVector(vector): @@ -5918,6 +5926,91 @@ def CombineRotation(a, b): ]) +def IdentityMatrix(): + """Creates an identity rotation matrix. + + Returns a rotation matrix that has no effect on orientation. + This matrix can be the starting point for other operations, + such as using a series of calls to #Pivot to + create a custom rotation matrix. + + Returns + ------- + RotationMatrix + The identity rotation matrix. + """ + return RotationMatrix([ + [1, 0, 0], + [0, 1, 0], + [0, 0, 1] + ]) + + +def Pivot(rotation, axis, angle): + """Re-orients a rotation matrix by pivoting it by an angle around one of its axes. + + Given a rotation matrix, a selected coordinate axis, and an angle in degrees, + this function pivots the rotation matrix by that angle around that coordinate axis. + + For example, if you have rotation matrix that converts ecliptic coordinates (ECL) + to horizontal coordinates (HOR), but you really want to convert ECL to the orientation + of a telescope camera pointed at a given body, you can use `Pivot` twice: + (1) pivot around the zenith axis by the body's azimuth, then (2) pivot around the + western axis by the body's altitude angle. The resulting rotation matrix will then + reorient ECL coordinates to the orientation of your telescope camera. + + Parameters + ---------- + rotation : RotationMatrix + The input rotation matrix. + axis : int + An integer that selects which coordinate axis to rotate around: + 0 = x, 1 = y, 2 = z. Any other value will cause an exception. + angle : float + An angle in degrees indicating the amount of rotation around the specified axis. + Positive angles indicate rotation counterclockwise as seen from the positive + direction along that axis, looking towards the origin point of the orientation system. + Any finite number of degrees is allowed, but best precision will result from + keeping `angle` in the range [-360, +360]. + + Returns + ------- + RotationMatrix + A pivoted matrix object. + """ + # Check for an invalid coordinate axis. + if axis not in [0, 1, 2]: + raise Error('Invalid axis {}. Must be [0, 1, 2].'.format(axis)) + + radians = math.radians(angle) + c = math.cos(radians) + s = math.sin(radians) + + # We need to maintain the "right-hand" rule, no matter which + # axis was selected. That means we pick (i, j, k) axis order + # such that the following vector cross product is satisfied: + # i x j = k + i = (axis + 1) % 3 + j = (axis + 2) % 3 + k = axis + + rot = [[0, 0, 0], [0, 0, 0], [0, 0, 0]] + + rot[i][i] = c*rotation.rot[i][i] - s*rotation.rot[i][j] + rot[i][j] = s*rotation.rot[i][i] + c*rotation.rot[i][j] + rot[i][k] = rotation.rot[i][k] + + rot[j][i] = c*rotation.rot[j][i] - s*rotation.rot[j][j] + rot[j][j] = s*rotation.rot[j][i] + c*rotation.rot[j][j] + rot[j][k] = rotation.rot[j][k] + + rot[k][i] = c*rotation.rot[k][i] - s*rotation.rot[k][j] + rot[k][j] = s*rotation.rot[k][i] + c*rotation.rot[k][j] + rot[k][k] = rotation.rot[k][k] + + return RotationMatrix(rot) + + def RotateVector(rotation, vector): """Applies a rotation to a vector, yielding a rotated vector. @@ -6790,7 +6883,7 @@ def Constellation(ra, dec): _Epoch2000 = Time(0.0) # Convert coordinates from J2000 to B1875. - equ2000 = Equatorial(ra, dec, 1.0) + equ2000 = Equatorial(ra, dec, 1.0, None) # FIXFIXFIX - avoid this hack with unused vector vec2000 = VectorFromEquator(equ2000, _Epoch2000) vec1875 = RotateVector(_ConstelRot, vec2000) equ1875 = EquatorFromVector(vec1875) diff --git a/demo/python/camera.py b/demo/python/camera.py new file mode 100755 index 00000000..88fdf14d --- /dev/null +++ b/demo/python/camera.py @@ -0,0 +1,97 @@ +#!/usr/bin/env python3 +# +# camera.py - by Don Cross - 2021-03-27 +# +# Example Python program for Astronomy Engine: +# https://github.com/cosinekitty/astronomy +# +# Suppose you want to photograph the Moon, +# and you want to know what it will look like in the photo. +# Given a location on the Earth, and a date/time, +# this program calculates the orientation of the sunlit +# side of the Moon with respect to the top of your +# photo image. It assumes the camera faces directly +# toward the Moon's azimuth and tilts upward to its +# altitude angle above the horizon. +# +# To execute, run the command: +# +# python3 camera.py latitude longitude [yyyy-mm-ddThh:mm:ssZ] +# +import sys +import math +import astronomy +from astro_demo_common import ParseArgs + +def Camera(observer, time): + tolerance = 1.0e-15 + + # Calculate the topocentric equatorial coordinates of date for the Moon. + # Assume aberration does not matter because the Moon is so close and has such a small relative velocity. + moon_equ = astronomy.Equator(astronomy.Body.Moon, time, observer, True, False) + + # Also calculate the Sun's topocentric position in the same coordinate system. + sun_equ = astronomy.Equator(astronomy.Body.Sun, time, observer, True, False) + + # Get the Moon's horizontal coordinates, so we know how much to pivot azimuth and altitude. + moon_hor = astronomy.Horizon(time, observer, moon_equ.ra, moon_equ.dec, astronomy.Refraction.Airless) + print('Moon horizontal position: azimuth = {:0.3f}, altitude = {:0.3f}'.format(moon_hor.azimuth, moon_hor.altitude)) + + # Get the rotation matrix that converts equatorial to horizontal coordintes for this place and time. + rot = astronomy.Rotation_EQD_HOR(time, observer) + + # Modify the rotation matrix in two steps: + # First, rotate the orientation so we are facing the Moon's azimuth. + # We do this by pivoting around the zenith axis. + # Horizontal axes are: 0 = north, 1 = west, 2 = zenith. + # Tricky: because the pivot angle increases counterclockwise, and azimuth + # increases clockwise, we undo the azimuth by adding the positive value. + rot = astronomy.Pivot(rot, 2, moon_hor.azimuth) + + # Second, pivot around the leftward axis to bring the Moon to the camera's altitude level. + # From the point of view of the leftward axis, looking toward the camera, + # adding the angle is the correct sense for subtracting the altitude. + rot = astronomy.Pivot(rot, 1, moon_hor.altitude) + + # As a sanity check, apply this rotation to the Moon's equatorial (EQD) coordinates and verify x=0, y=0. + vec = astronomy.RotateVector(rot, moon_equ.vec) + + # Convert to unit vector. + radius = vec.Length() + vec.x /= radius + vec.y /= radius + vec.z /= radius + print('Moon check: {}'.format(vec)) + if abs(vec.x - 1.0) > tolerance: + print("Excessive error in moon check (x).") + return 1 + + if abs(vec.y) > tolerance: + print("Excessive error in moon check (y).") + return 1 + + if abs(vec.z) > tolerance: + print("Excessive error in moon check (z).") + return 1 + + # Apply the same rotation to the Sun's equatorial vector. + # The x- and y-coordinates now tell us which side appears sunlit in the camera! + + vec = astronomy.RotateVector(rot, sun_equ.vec) + + # Don't bother normalizing the Sun vector, because in AU it will be close to unit anyway. + print('Sun vector: {}'.format(vec)) + + # Calculate the tilt angle of the sunlit side, as seen by the camera. + # The x-axis is now pointing directly at the object, z is up in the camera image, y is to the left. + tilt = math.degrees(math.atan2(vec.z, vec.y)) + print('Tilt angle of sunlit side of the Moon = {:0.3f} degrees counterclockwise from up.'.format(tilt)) + illum = astronomy.Illumination(astronomy.Body.Moon, time) + print('Moon magnitude = {:0.2f}, phase angle = {:0.2f} degrees.'.format(illum.mag, illum.phase_angle)) + angle = astronomy.AngleFromSun(astronomy.Body.Moon, time) + print('Angle between Moon and Sun as seen from Earth = {:0.2f} degrees.'.format(angle)) + return 0 + +if __name__ == '__main__': + observer, time = ParseArgs(sys.argv) + sys.exit(Camera(observer, time)) diff --git a/demo/python/test/.gitignore b/demo/python/test/.gitignore index 73426264..4529900c 100644 --- a/demo/python/test/.gitignore +++ b/demo/python/test/.gitignore @@ -1,3 +1,4 @@ +camera.txt moonphase.txt positions.txt riseset.txt diff --git a/demo/python/test/camera_correct.txt b/demo/python/test/camera_correct.txt new file mode 100644 index 00000000..a76cf3db --- /dev/null +++ b/demo/python/test/camera_correct.txt @@ -0,0 +1,6 @@ +Moon horizontal position: azimuth = 274.486, altitude = 52.101 +Moon check: (1.0, -2.0725886704096896e-16, 1.2435532022458137e-16, 2021-03-22T02:45:00.000Z) +Sun vector: (-0.0885303206557285, -0.3193912784902719, -0.9396631499852888, 2021-03-22T02:45:00.000Z) +Tilt angle of sunlit side of the Moon = -108.773 degrees counterclockwise from up. +Moon magnitude = -10.27, phase angle = 84.22 degrees. +Angle between Moon and Sun as seen from Earth = 95.64 degrees. diff --git a/demo/python/test/test b/demo/python/test/test index 2617bf3b..c85ae2b0 100755 --- a/demo/python/test/test +++ b/demo/python/test/test @@ -5,7 +5,11 @@ Fail() exit 1 } -rm -f test/{moonphase,positions,riseset,seasons,culminate}.txt +rm -f test/{camera,moonphase,positions,riseset,seasons,culminate,horizon,lunar_eclipse}.txt + +echo "Testing example: camera.py" +./camera.py 29 -81 2021-03-22T02:45:00Z > test/camera.txt || Fail "Error testing camera.py." +diff test/camera.txt test/camera_correct.txt || Fail "Error comparing camera.py output." echo "Testing example: moonphase.py" ./moonphase.py 2019-06-15T09:15:32.987Z > test/moonphase.txt || Fail "Error running moonphase.py." diff --git a/generate/template/astronomy.py b/generate/template/astronomy.py index d7650d83..55762889 100644 --- a/generate/template/astronomy.py +++ b/generate/template/astronomy.py @@ -787,13 +787,20 @@ class Equatorial: Declination in degrees. dist : float Distance to the celestial body in AU. + vec : Vector + The equatorial coordinates in cartesian form, using AU distance units. + x = direction of the March equinox, + y = direction of the June solstice, + z = north. """ - def __init__(self, ra, dec, dist): + def __init__(self, ra, dec, dist, vec): self.ra = ra self.dec = dec self.dist = dist + self.vec = vec -def _vector2radec(pos): + +def _vector2radec(pos, time): xyproj = pos[0]*pos[0] + pos[1]*pos[1] dist = math.sqrt(xyproj + pos[2]*pos[2]) if xyproj == 0.0: @@ -810,7 +817,8 @@ def _vector2radec(pos): if ra < 0: ra += 24 dec = math.degrees(math.atan2(pos[2], math.sqrt(xyproj))) - return Equatorial(ra, dec, dist) + vec = Vector(pos[0], pos[1], pos[2], time) + return Equatorial(ra, dec, dist, vec) def _nutation_rot(time, direction): @@ -1868,10 +1876,10 @@ def Equator(body, time, observer, ofdate, aberration): gc.z - gc_observer[2] ] if not ofdate: - return _vector2radec(j2000) + return _vector2radec(j2000, time) temp = _precession(0, j2000, time.tt) datevect = _nutation(time, 0, temp) - return _vector2radec(datevect) + return _vector2radec(datevect, time) @enum.unique class Refraction(enum.Enum): @@ -3848,7 +3856,7 @@ def EquatorFromVector(vec): Angular coordinates expressed in the same equatorial system as `vec`. """ sphere = SphereFromVector(vec) - return Equatorial(sphere.lon / 15.0, sphere.lat, sphere.dist) + return Equatorial(sphere.lon / 15.0, sphere.lat, sphere.dist, vec) def SphereFromVector(vector): @@ -4021,6 +4029,91 @@ def CombineRotation(a, b): ]) +def IdentityMatrix(): + """Creates an identity rotation matrix. + + Returns a rotation matrix that has no effect on orientation. + This matrix can be the starting point for other operations, + such as using a series of calls to #Pivot to + create a custom rotation matrix. + + Returns + ------- + RotationMatrix + The identity rotation matrix. + """ + return RotationMatrix([ + [1, 0, 0], + [0, 1, 0], + [0, 0, 1] + ]) + + +def Pivot(rotation, axis, angle): + """Re-orients a rotation matrix by pivoting it by an angle around one of its axes. + + Given a rotation matrix, a selected coordinate axis, and an angle in degrees, + this function pivots the rotation matrix by that angle around that coordinate axis. + + For example, if you have rotation matrix that converts ecliptic coordinates (ECL) + to horizontal coordinates (HOR), but you really want to convert ECL to the orientation + of a telescope camera pointed at a given body, you can use `Pivot` twice: + (1) pivot around the zenith axis by the body's azimuth, then (2) pivot around the + western axis by the body's altitude angle. The resulting rotation matrix will then + reorient ECL coordinates to the orientation of your telescope camera. + + Parameters + ---------- + rotation : RotationMatrix + The input rotation matrix. + axis : int + An integer that selects which coordinate axis to rotate around: + 0 = x, 1 = y, 2 = z. Any other value will cause an exception. + angle : float + An angle in degrees indicating the amount of rotation around the specified axis. + Positive angles indicate rotation counterclockwise as seen from the positive + direction along that axis, looking towards the origin point of the orientation system. + Any finite number of degrees is allowed, but best precision will result from + keeping `angle` in the range [-360, +360]. + + Returns + ------- + RotationMatrix + A pivoted matrix object. + """ + # Check for an invalid coordinate axis. + if axis not in [0, 1, 2]: + raise Error('Invalid axis {}. Must be [0, 1, 2].'.format(axis)) + + radians = math.radians(angle) + c = math.cos(radians) + s = math.sin(radians) + + # We need to maintain the "right-hand" rule, no matter which + # axis was selected. That means we pick (i, j, k) axis order + # such that the following vector cross product is satisfied: + # i x j = k + i = (axis + 1) % 3 + j = (axis + 2) % 3 + k = axis + + rot = [[0, 0, 0], [0, 0, 0], [0, 0, 0]] + + rot[i][i] = c*rotation.rot[i][i] - s*rotation.rot[i][j] + rot[i][j] = s*rotation.rot[i][i] + c*rotation.rot[i][j] + rot[i][k] = rotation.rot[i][k] + + rot[j][i] = c*rotation.rot[j][i] - s*rotation.rot[j][j] + rot[j][j] = s*rotation.rot[j][i] + c*rotation.rot[j][j] + rot[j][k] = rotation.rot[j][k] + + rot[k][i] = c*rotation.rot[k][i] - s*rotation.rot[k][j] + rot[k][j] = s*rotation.rot[k][i] + c*rotation.rot[k][j] + rot[k][k] = rotation.rot[k][k] + + return RotationMatrix(rot) + + def RotateVector(rotation, vector): """Applies a rotation to a vector, yielding a rotated vector. @@ -4442,7 +4535,7 @@ def Constellation(ra, dec): _Epoch2000 = Time(0.0) # Convert coordinates from J2000 to B1875. - equ2000 = Equatorial(ra, dec, 1.0) + equ2000 = Equatorial(ra, dec, 1.0, None) # FIXFIXFIX - avoid this hack with unused vector vec2000 = VectorFromEquator(equ2000, _Epoch2000) vec1875 = RotateVector(_ConstelRot, vec2000) equ1875 = EquatorFromVector(vec1875) diff --git a/generate/test.py b/generate/test.py index 53aa6da9..0e3f0ecf 100755 --- a/generate/test.py +++ b/generate/test.py @@ -762,6 +762,23 @@ def CompareMatrices(caller, a, b, tolerance): sys.exit(1) +def CompareVectors(caller, a, b, tolerance): + diff = vabs(a.x - b.x) + if diff > tolerance: + print('PY CompareVectors ERROR({}): vector x = {}, expected {}, diff {}'.format(caller, a.x, b.x, diff)) + sys.exit(1) + + diff = vabs(a.y - b.y) + if diff > tolerance: + print('PY CompareVectors ERROR({}): vector y = {}, expected {}, diff {}'.format(caller, a.y, b.y, diff)) + sys.exit(1) + + diff = vabs(a.z - b.z) + if diff > tolerance: + print('PY CompareVectors ERROR({}): vector z = {}, expected {}, diff {}'.format(caller, a.z, b.z, diff)) + sys.exit(1) + + def Rotation_MatrixInverse(): a = astronomy.RotationMatrix([ [1, 4, 7], @@ -1010,9 +1027,48 @@ def Test_RotRoundTrip(): Debug('PY Test_RotRoundTrip: PASS') +def Rotation_Pivot(): + tolerance = 1.0e-15 + + # Start with an identity matrix. + ident = astronomy.IdentityMatrix() + + # Pivot 90 degrees counterclockwise around the z-axis. + r = astronomy.Pivot(ident, 2, +90.0) + + # Put the expected answer in 'a'. + a = astronomy.RotationMatrix([ + [ 0, +1, 0], + [-1, 0, 0], + [ 0, 0, +1], + ]) + + # Compare actual 'r' with expected 'a'. + CompareMatrices('Rotation_Pivot #1', r, a, tolerance) + + # Pivot again, -30 degrees around the x-axis. + r = astronomy.Pivot(r, 0, -30.0) + + # Pivot a third time, 180 degrees around the y-axis. + r = astronomy.Pivot(r, 1, +180.0) + + # Use the 'r' matrix to rotate a vector. + v1 = astronomy.Vector(1, 2, 3, astronomy.Time(0)) + v2 = astronomy.RotateVector(r, v1) + + # Initialize the expected vector 've'. + ve = astronomy.Vector(+2.0, +2.3660254037844390, -2.0980762113533156, v1.t) + + CompareVectors('Rotation_Pivot #2', v2, ve, tolerance) + + Debug('PY Rotation_Pivot: PASS') + + + def Rotation(): Rotation_MatrixInverse() Rotation_MatrixMultiply() + Rotation_Pivot() Test_EQJ_ECL() Test_EQJ_EQD(astronomy.Body.Mercury) Test_EQJ_EQD(astronomy.Body.Venus) diff --git a/source/python/README.md b/source/python/README.md index 242b2439..47634b85 100644 --- a/source/python/README.md +++ b/source/python/README.md @@ -265,6 +265,7 @@ equator projected onto the sky. | `float` | `ra` | Right ascension in sidereal hours. | | `float` | `dec` | Declination in degrees. | | `float` | `dist` | Distance to the celestial body in AU. | +| [`Vector`](#Vector) | `vec` | The equatorial coordinates in cartesian form, using AU distance units. x = direction of the March equinox, y = direction of the June solstice, z = north. | --- @@ -1197,6 +1198,21 @@ and is expressed in astronomical units (AU). --- + +### IdentityMatrix() + +**Creates an identity rotation matrix.** + +Returns a rotation matrix that has no effect on orientation. +This matrix can be the starting point for other operations, +such as using a series of calls to [`Pivot`](#Pivot) to +create a custom rotation matrix. + +### Returns: [`RotationMatrix`](#RotationMatrix) +The identity rotation matrix. + +--- + ### Illumination(body, time) @@ -1453,6 +1469,31 @@ Keep calling this function as many times as you want to keep finding more transi --- + +### Pivot(rotation, axis, angle) + +**Re-orients a rotation matrix by pivoting it by an angle around one of its axes.** + +Given a rotation matrix, a selected coordinate axis, and an angle in degrees, +this function pivots the rotation matrix by that angle around that coordinate axis. +For example, if you have rotation matrix that converts ecliptic coordinates (ECL) +to horizontal coordinates (HOR), but you really want to convert ECL to the orientation +of a telescope camera pointed at a given body, you can use `Pivot` twice: +(1) pivot around the zenith axis by the body's azimuth, then (2) pivot around the +western axis by the body's altitude angle. The resulting rotation matrix will then +reorient ECL coordinates to the orientation of your telescope camera. + +| Type | Parameter | Description | +| --- | --- | --- | +| [`RotationMatrix`](#RotationMatrix) | `rotation` | The input rotation matrix. | +| `int` | `axis` | An integer that selects which coordinate axis to rotate around: 0 = x, 1 = y, 2 = z. Any other value will cause an exception. | +| `float` | `angle` | An angle in degrees indicating the amount of rotation around the specified axis. Positive angles indicate rotation counterclockwise as seen from the positive direction along that axis, looking towards the origin point of the orientation system. Any finite number of degrees is allowed, but best precision will result from keeping `angle` in the range [-360, +360]. | + +### Returns: [`RotationMatrix`](#RotationMatrix) +A pivoted matrix object. + +--- + ### RefractionAngle(refraction, altitude) diff --git a/source/python/astronomy.py b/source/python/astronomy.py index fc7c6679..56978617 100644 --- a/source/python/astronomy.py +++ b/source/python/astronomy.py @@ -1323,13 +1323,20 @@ class Equatorial: Declination in degrees. dist : float Distance to the celestial body in AU. + vec : Vector + The equatorial coordinates in cartesian form, using AU distance units. + x = direction of the March equinox, + y = direction of the June solstice, + z = north. """ - def __init__(self, ra, dec, dist): + def __init__(self, ra, dec, dist, vec): self.ra = ra self.dec = dec self.dist = dist + self.vec = vec -def _vector2radec(pos): + +def _vector2radec(pos, time): xyproj = pos[0]*pos[0] + pos[1]*pos[1] dist = math.sqrt(xyproj + pos[2]*pos[2]) if xyproj == 0.0: @@ -1346,7 +1353,8 @@ def _vector2radec(pos): if ra < 0: ra += 24 dec = math.degrees(math.atan2(pos[2], math.sqrt(xyproj))) - return Equatorial(ra, dec, dist) + vec = Vector(pos[0], pos[1], pos[2], time) + return Equatorial(ra, dec, dist, vec) def _nutation_rot(time, direction): @@ -3765,10 +3773,10 @@ def Equator(body, time, observer, ofdate, aberration): gc.z - gc_observer[2] ] if not ofdate: - return _vector2radec(j2000) + return _vector2radec(j2000, time) temp = _precession(0, j2000, time.tt) datevect = _nutation(time, 0, temp) - return _vector2radec(datevect) + return _vector2radec(datevect, time) @enum.unique class Refraction(enum.Enum): @@ -5745,7 +5753,7 @@ def EquatorFromVector(vec): Angular coordinates expressed in the same equatorial system as `vec`. """ sphere = SphereFromVector(vec) - return Equatorial(sphere.lon / 15.0, sphere.lat, sphere.dist) + return Equatorial(sphere.lon / 15.0, sphere.lat, sphere.dist, vec) def SphereFromVector(vector): @@ -5918,6 +5926,91 @@ def CombineRotation(a, b): ]) +def IdentityMatrix(): + """Creates an identity rotation matrix. + + Returns a rotation matrix that has no effect on orientation. + This matrix can be the starting point for other operations, + such as using a series of calls to #Pivot to + create a custom rotation matrix. + + Returns + ------- + RotationMatrix + The identity rotation matrix. + """ + return RotationMatrix([ + [1, 0, 0], + [0, 1, 0], + [0, 0, 1] + ]) + + +def Pivot(rotation, axis, angle): + """Re-orients a rotation matrix by pivoting it by an angle around one of its axes. + + Given a rotation matrix, a selected coordinate axis, and an angle in degrees, + this function pivots the rotation matrix by that angle around that coordinate axis. + + For example, if you have rotation matrix that converts ecliptic coordinates (ECL) + to horizontal coordinates (HOR), but you really want to convert ECL to the orientation + of a telescope camera pointed at a given body, you can use `Pivot` twice: + (1) pivot around the zenith axis by the body's azimuth, then (2) pivot around the + western axis by the body's altitude angle. The resulting rotation matrix will then + reorient ECL coordinates to the orientation of your telescope camera. + + Parameters + ---------- + rotation : RotationMatrix + The input rotation matrix. + axis : int + An integer that selects which coordinate axis to rotate around: + 0 = x, 1 = y, 2 = z. Any other value will cause an exception. + angle : float + An angle in degrees indicating the amount of rotation around the specified axis. + Positive angles indicate rotation counterclockwise as seen from the positive + direction along that axis, looking towards the origin point of the orientation system. + Any finite number of degrees is allowed, but best precision will result from + keeping `angle` in the range [-360, +360]. + + Returns + ------- + RotationMatrix + A pivoted matrix object. + """ + # Check for an invalid coordinate axis. + if axis not in [0, 1, 2]: + raise Error('Invalid axis {}. Must be [0, 1, 2].'.format(axis)) + + radians = math.radians(angle) + c = math.cos(radians) + s = math.sin(radians) + + # We need to maintain the "right-hand" rule, no matter which + # axis was selected. That means we pick (i, j, k) axis order + # such that the following vector cross product is satisfied: + # i x j = k + i = (axis + 1) % 3 + j = (axis + 2) % 3 + k = axis + + rot = [[0, 0, 0], [0, 0, 0], [0, 0, 0]] + + rot[i][i] = c*rotation.rot[i][i] - s*rotation.rot[i][j] + rot[i][j] = s*rotation.rot[i][i] + c*rotation.rot[i][j] + rot[i][k] = rotation.rot[i][k] + + rot[j][i] = c*rotation.rot[j][i] - s*rotation.rot[j][j] + rot[j][j] = s*rotation.rot[j][i] + c*rotation.rot[j][j] + rot[j][k] = rotation.rot[j][k] + + rot[k][i] = c*rotation.rot[k][i] - s*rotation.rot[k][j] + rot[k][j] = s*rotation.rot[k][i] + c*rotation.rot[k][j] + rot[k][k] = rotation.rot[k][k] + + return RotationMatrix(rot) + + def RotateVector(rotation, vector): """Applies a rotation to a vector, yielding a rotated vector. @@ -6790,7 +6883,7 @@ def Constellation(ra, dec): _Epoch2000 = Time(0.0) # Convert coordinates from J2000 to B1875. - equ2000 = Equatorial(ra, dec, 1.0) + equ2000 = Equatorial(ra, dec, 1.0, None) # FIXFIXFIX - avoid this hack with unused vector vec2000 = VectorFromEquator(equ2000, _Epoch2000) vec1875 = RotateVector(_ConstelRot, vec2000) equ1875 = EquatorFromVector(vec1875)