r""" Module for the sky mapping.
A sky map is sometimes needed to display the instrument FOV in an instrument-center frame.
This is identical to the azimuth-altitude system.
https://en.wikipedia.org/wiki/File:Azimut_altitude.svg
To avoid the confusion in the term "altitude", frequently used to describe the distance from the surface,
we instead use the "elevation" here.
- Elevation angle |alpha| or El: The angle from the "horizon".
- Polar angle |theta| or T: 90 degrees - elevation angle.
- Azimuthal angle |beta| or Az: North is 0, East is 90, South is +-180, and West is -90.
Usually, the Az shall be defined in the range ``(-180, 180]``.
- Primary axis $\hat{p}$ or P: The primary axis defines the zenith, i.e., elevation = 90. (Frequently $z$ axis)
- Secondary axis $\hat{s}$ or S: The secondary axis defines the north, frequenty $x$ axis.
More precisely speaking, the meridan defined by the primary and secondary axes contains the north.
- Third axis $\hat{t}$ or T: The third axis. This defines the west (frequently $y$ axis).
This axis is calculated from the primary and secondary axes.
- Distance, d or D, is also important for space physics, as a distance to the object of interest.
The interface always accept the angles in the unit of "degrees".
.. note::
Traditionally, several variants of the azimuth angle have been defined, so it sometimes confusing.
Here irfpy follow the defintition described in
https://www.sciencedirect.com/science/article/abs/pii/S0038092X0300450X
See also https://en.wikipedia.org/wiki/Solar_azimuth_angle.
"""
import numpy as _np
from irfpy.util import intersection as _ix
from irfpy.util import sphere as _sphere
[docs]class SkymapFrame:
def __init__(self, primary_axis, secondary_axis):
""" Define the skymap frame.
A skymap frame is defined, following the defintion above.
:param primary_axis: The primary axis, (3,) shape numpy array.
:param secondary_axis: The secondary axis, (3,) shape numpy array.
:raises ValueError: Raised when the given axes do not provide appropriate cartesian.
>>> skymap_frame = SkymapFrame([0, 0, 1], [1, 0, 0])
>>> print(skymap_frame) # doctest: +NORMALIZE_WHITESPACE
Skymap object: P= 0.0000 S= 1.0000 T= 0.0000
0.0000 0.0000 1.0000
1.0000 0.0000 0.0000
This skymap is relatively simple example, namely, z is the zenith and x is "north".
In this case, the vector (0, 0, 100) will be converted to d=100, el=90 and az=(N.A.).
>>> d, el, az = skymap_frame.vec2dea(0, 0, 100)
>>> print(d)
100.0
>>> print(el)
90.0
The vector (1, 0, 0), looking north will be d=1, el=0, az=0.
>>> print(skymap_frame.vec2dea(1, 0, 0))
(1.0, 0.0, -0.0)
The vector (0, 1, 0), looking west will be d=1, el=0, az=-90. Here note the sign of az.
>>> print(skymap_frame.vec2dea(0, 1, 0))
(1.0, 0.0, -90.0)
:meth:`vec2dea` can eat a list as arguments.
>>> print(skymap_frame.vec2dea([1, 0], [0, 1], [0, 0])) # doctest: +NORMALIZE_WHITESPACE
(array([1., 1.]), array([0., 0.]), array([ -0., -90.]))
"""
self._Pax = _np.array(primary_axis)
self._Sax = _np.array(secondary_axis)
if self._Pax.shape != (3,):
raise ValueError('Given primary axis shape must be (3,). Given={}'.format(str(self._Pax.shape)))
if self._Sax.shape != (3,):
raise ValueError('Given secondary axis shape must be!= (3,). Given={}'.format(str(self._Sax.shape)))
self._Tax = _np.cross(self._Pax, self._Sax)
lT = _np.linalg.norm(self._Tax)
if lT == 0:
raise ValueError('Given primary and secondary axes looks parallel '
'or one of them have zero-length')
self._Sax = _np.cross(self._Tax, self._Pax)
self._Pax = self._Pax / _np.linalg.norm(self._Pax)
self._Sax = self._Sax / _np.linalg.norm(self._Sax)
self._Tax = self._Tax / _np.linalg.norm(self._Tax)
def __str__(self):
p = self._Pax
s = self._Sax
t = self._Tax
val = (f"Skymap object: P={p[0]:7.4f} S={s[0]:7.4f} T={t[0]:7.4f}\n" +
f" {p[1]:7.4f} {s[1]:7.4f} {t[1]:7.4f}\n" +
f" {p[2]:7.4f} {s[2]:7.4f} {t[2]:7.4f}")
return val
[docs] def vec2dea(self, x, y, z):
r""" Vector list ([x_i], [y_i], [z_i]) is converted to distance-elevation-azimuth.
A vector $\vec{v}$ is converted to (d, e, a).
"""
vec = _np.array([x, y, z]) # (3, ...) shape
xx = _np.tensordot(self._Sax, vec, axes=(-1, 0)) # (...) shape, or float
yy = _np.tensordot(self._Tax, vec, axes=(-1, 0)) # (...) shape, or float
zz = _np.tensordot(self._Pax, vec, axes=(-1, 0)) # (...) shape, or float
d = _np.sqrt(xx ** 2 + yy ** 2 + zz ** 2)
theta = _np.degrees(_np.arccos(_np.clip(zz / d, -1, 1))) # clip is needed as zz/d seldomly gives >1 or <-1 value due to numerical error
elevation = _np.clip(90 - theta, -90, 90)
phi = _np.degrees(_np.arctan2(yy, xx)) # Phi, the normal azimuthal angle but is "opposite-singed" definition
azimuth = -phi
return (d, elevation, azimuth)
@property
def primary_axis(self):
""" The primary axis, the 90 degrees of elevation.
The primary axis, or equivalently the z axis. (Zenithward)
"""
return self._Pax.copy()
@property
def secondary_axis(self):
""" The secondary axis, the 0 elevation and 0 azimuth.
The secondary axis, or equivalently the x axis. (Northward)
"""
return self._Sax.copy()
@property
def third_axis(self):
""" The third axis, the 0 elevation and 90 azimuth
The third axis, or eqiuvalently the y axis. (Eastward)
"""
return self._Tax.copy()
[docs] def limb_of_sphere(self, position: _np.ndarray, radius: float, resolution: int=16) -> tuple[_np.ndarray, _np.ndarray]:
""" Return a limb of sphere coordinates mapped on the sky.
:param position: Position of the sphere. (3,) shaped numpy array.
:param radius: Radius of the sphere.
:returns: A tuple, (elevation_array, azimuth_array) in degrees.
>>> sk = SkymapFrame([0, 0, 1], [1, 0, 0])
>>> e, a = sk.limb_of_sphere([0, 0, 3], 1)
>>> print(f'{e.max():.2f} {e.min():.2f}') # Must be 70.52877937 for both.
70.53 70.53
>>> print(f'{a[1]:.2f}')
-24.00
>>> print(a) # doctest: +SKIP
[0 -24 -48 -72 -96 -120 -144 -168 168 144 120 96 72 48 24 0]
>>> e, a = sk.limb_of_sphere([0, 3, 0], 1, resolution=360)
>>> print(f'{e.max():.2f} {e.min():.2f}')
19.47 -19.47
>>> print(a) # doctest: +SKIP
[ -90. -90.35451559 -90.70889547 -91.06300404 -91.41670596
-91.76986629 -92.1223506 -92.47402508 -92.82475668
...
-86.82558674 -87.17524332 -87.52597492 -87.8776494 -88.23013371
-88.58329404 -88.93699596 -89.29110453 -89.64548441 -90. ]
>>> print(f'{a.min():.2f} {a.max():.2f}')
-109.47 -70.53
"""
# First calculate the theta, phi.
_t, _p = _ix.sphere_limb_plane_projection(position, radius, resolution=resolution)
_x = _np.sin(_t) * _np.cos(_p)
_y = _np.sin(_t) * _np.sin(_p)
_z = _np.cos(_t)
d_limb, e_limb, a_limb = self.vec2dea(_x, _y, _z)
return (e_limb, a_limb)
[docs] def parallel_on_sphere(self, sphere_position, sphere_radius, sphere_primary_axis, angle,
resolution=16, mask_invisible=True):
""" Obtain elevation and azimuth angle list of parallel line on the sphere
Suppose we have a sphere (e.g., Venus), and want to draw equator.
This method do that.
:param sphere_position: The central position of the sphere
:param sphere_radius: The radius of sphere (in the dimension of length, e.g., km)
:param sphere_primary_axis: The primary axis direction measured from the central position of the sphere.
:param angle: The angle measured from the sphere_primary_axis.
:keyword resolution: Resolution.
:keyword mask_invisible: Mask the elevation, azimuth angle if that location cannot
be visible from the observer. Default, True (i.e., masked)
"""
sphere = _sphere.Sphere(sphere_radius) # Sphere is first defined, centered at the origin.
parallel_distance = sphere_radius * _np.cos(_np.radians(angle))
cuts = sphere.cut(sphere_primary_axis, parallel_distance, resolution=resolution)
cuts = _np.array(cuts) # (N, 3) shape. Vectors mesaured from the sphere center.
cuts_from_observer = cuts + sphere_position[_np.newaxis, :] # (N, 3) shape
d, e, a = self.vec2dea(cuts_from_observer[:, 0], cuts_from_observer[:, 1], cuts_from_observer[:, 2])
# Visibility check?
if mask_invisible:
mask = [not _ix.isvisible([0, 0, 0], vec, sphere_position, sphere_radius * 0.9999) for vec in cuts_from_observer] # True for invisible
e = _np.ma.masked_array(e, mask=mask)
a = _np.ma.masked_array(a, mask=mask)
# print(mask)
return (e, a)