mirror/astronomy
1
0
Fork 0
mirror of https://github.com/cosinekitty/astronomy.git synced 2026-04-04 14:55:38 -04:00
Code Issues Packages Projects Releases Wiki Activity

Ported IdentityMatrix, Pivot functions to Python. Added tests and camera demo.

Browse Source
This commit is contained in:
Don Cross
2021-03-27 05:19:27 -04:00
parent a4d61c872a
commit a97fc7da9c
10 changed files with 520 additions and 22 deletions
Download Patch File Download Diff File Expand all files Collapse all files

41
source/python/README.md
View File

@@ -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).
---
<a name="IdentityMatrix"></a>
### 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.
---
<a name="Illumination"></a>
### Illumination(body, time)
@@ -1453,6 +1469,31 @@ Keep calling this function as many times as you want to keep finding more transi
---
<a name="Pivot"></a>
### 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.
---
<a name="RefractionAngle"></a>
### RefractionAngle(refraction, altitude)

107
source/python/astronomy.py
View File

@@ -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)