# coord.py
"""
Coordinate transformations.
"""
import numpy as np
pi = np.pi
[docs]
def rotate_xzplane_to_source(
axis_zenith, source_zenith, source_azimuth, zenith, azimuth,
reverse=False, latlon=False):
"""
Rotate coordinates such that points in the x-z plane move to a source location.
Args:
axis_zenith (float): the destination source zenith
source_zenith (float or ndarray): the destination source zenith
source_azimuth (float or ndarray): the destination source azimuth
zenith (float or ndarray): the zenith of the map point(s)
azimuth (float or ndarray): the azimuth of the map point(s)
reverse (bool): if True, perform the inverse coordinate rotation
latlon (bool): if True, ``source_zenith`` and ``zenith`` are interpreted as
declinations rather than polar zeniths
Input coordinates are assumed to be standard spherical(in DC) or
equatorial.
"""
if latlon:
source_zenith = pi/2 - source_zenith
zenith = pi/2 - zenith
if reverse:
delta_zenith = source_zenith - axis_zenith
pre_delta_azimuth = -source_azimuth
post_delta_azimuth = 0
else:
delta_zenith = axis_zenith - source_zenith
pre_delta_azimuth = 0
post_delta_azimuth = source_azimuth
source_azimuth = source_azimuth + pre_delta_azimuth
azimuth = azimuth + pre_delta_azimuth
# cosines and sines
cz = np.cos(zenith)
sz = np.sin(zenith)
ca = np.cos(azimuth)
sa = np.sin(azimuth)
try:
len_source = len(source_zenith)
except:
len_source = 1
try:
len_point = len(zenith)
except:
len_point = 1
if len_source > 1 and len_point > 1 and len_source != len_point:
raise ValueError(
'input must be matching arrays, single source, single point, '
'or single source plus single point')
if len_source > 1 and len_point > 1:
coz = np.cos(delta_zenith)
soz = np.sin(delta_zenith)
filler_0 = np.zeros_like(coz)
filler_1 = np.ones_like(coz)
else:
filler_0 = np.zeros_like(cz)
filler_1 = np.ones_like(cz)
coz = filler_1 * np.cos(delta_zenith)
soz = filler_1 * np.sin(delta_zenith)
# coordinates(near x axis) and rotation matrix(about y axis)
xaxis_xyz = np.array([ca*sz, sa*sz, cz])
R_y = np.array([
[coz, filler_0, -soz ],
[filler_0, filler_1, filler_0 ],
[soz, filler_0, coz ]
])
# calculate R_y * xaxis_xyz manually(maybe could roll axes into place)
try:
# x = np.sum(R_y[0] * xaxis_xyz[:,:], axis=0)
# y = np.sum(R_y[1] * xaxis_xyz[:,:], axis=0)
# z = np.sum(R_y[2] * xaxis_xyz[:,:], axis=0)
Ryx, Ryy, Ryz = R_y
s = Ryx.shape
if len_source > 1 and len_point == 1:
Ryx, Ryy, Ryz = Ryx.T, Ryy.T, Ryz.T
x = np.sum(Ryx * xaxis_xyz, axis=0)
y = np.sum(Ryy * xaxis_xyz, axis=0)
z = np.sum(Ryz * xaxis_xyz, axis=0)
except:
raise
zenith_final = np.arccos(np.clip(z, -1, 1))
if latlon:
zenith_final = pi/2 - zenith_final
azimuth_final = np.arctan2 (y, x) + post_delta_azimuth
if zenith_final.shape:
azimuth_final[azimuth_final < 0] += 2*pi
azimuth_final[azimuth_final > 2*pi] -= 2*pi
else:
while azimuth_final > 2*pi:
azimuth_final -= 2*pi
while azimuth_final < 0:
azimuth_final += 2*pi
return zenith_final, azimuth_final
[docs]
def rotate_xaxis_to_source(source_zenith, source_azimuth, zenith, azimuth,
latlon=False):
"""
Rotate coordinates such that points on the xaxis move to a source location.
Args:
source_zenith (float or ndarray): the destination source zenith
source_azimuth (float or ndarray): the destination source azimuth
zenith (float or ndarray): the zenith of the map point(s)
azimuth (float or ndarray): the azimuth of the map point(s)
latlong (bool): if True, ``source_zenith`` and ``zenith`` are interpreted as
declinations rather than polar zeniths
Input coordinates are assumed to be standard spherical(either DC or CC).
"""
return rotate_xzplane_to_source(
pi/2, source_zenith, source_azimuth, zenith, azimuth,
reverse=False, latlon=latlon
)
[docs]
def rotate_source_to_xaxis(source_zenith, source_azimuth, zenith, azimuth,
latlon=False):
"""
Rotate coordinates such that points at the source location move to the x
axis.
Args:
source_zenith (float or ndarray): the destination source zenith
source_azimuth (float or ndarray): the destination source azimuth
zenith (float or ndarray): the zenith of the map point(s)
azimuth (float or ndarray): the azimuth of the map point(s)
latlon (bool): if True, ``source_zenith`` and ``zenith`` are interpreted as
declinations rather than polar zeniths
Input coordinates are assumed to be standard spherical (either DC or CC).
"""
return rotate_xzplane_to_source(
pi/2, source_zenith, source_azimuth, zenith, azimuth,
reverse=True, latlon=latlon
)
[docs]
def rotate_zaxis_to_source(source_zenith, source_azimuth, zenith, azimuth,
latlon=True):
"""
Rotate coordinates such that points on the zaxis move to a source location.
Args:
source_zenith (float or ndarray): the destination source zenith
source_azimuth (float or ndarray): the destination source azimuth
zenith (float or ndarray): the zenith of the map point(s)
azimuth (float or ndarray): the azimuth of the map point(s)
latlon (bool): if True, ``source_zenith`` and ``zenith`` are interpreted as
declinations rather than polar zeniths
Input coordinates are assumed to be standard spherical (either DC or CC).
"""
return rotate_xzplane_to_source(
0, source_zenith, source_azimuth, zenith, azimuth,
reverse=False, latlon=latlon
)
[docs]
def rotate_source_to_zaxis(source_zenith, source_azimuth, zenith, azimuth,
latlon=False):
"""
Rotate coordinates such that points at the source location move to the z
axis.
Args:
source_zenith (float or ndarray): the destination source zenith
source_azimuth (float or ndarray): the destination source azimuth
zenith (float or ndarray): the zenith of the map point(s)
azimuth (float or ndarray): the azimuth of the map point(s)
latlon (bool): if True, ``source_zenith`` and ``zenith`` are interpreted as
declinations rather than polar zeniths
Input coordinates are assumed to be standard spherical (either DC or CC).
"""
return rotate_xzplane_to_source(
0, source_zenith, source_azimuth, zenith, azimuth,
reverse=True, latlon=latlon
)
[docs]
def switch_ra_azimuth(phi_in, mjd):
"""Givin MJD, transform azimuth->RA **or** RA->azimuth."""
sidereal_day = 0.997269566 # sidereal day = length * solar day
sidereal_offset = 2.54199002505 # RA = offset + 2pi * (MJD/sidereal_length)%1 - azimuth
return (sidereal_offset + 2*pi*((mjd/sidereal_day)%1) - phi_in) % (2*pi)
[docs]
def delta_angle(zenith1, azimuth1, zenith2, azimuth2, cos=False, latlon=False):
"""Return the space angle between ``fit1`` and ``fit2``.
Args:
zenith1 (float or ndarray): zenith 1 or declination 1
azimuth1 (float or ndarray): azimuth 1
zenith2 (float or ndarray): zenith 2 or declination 2
azimuth2 (float or ndarray): azimuth 2
cos (bool): if True, return cos(delta_angle); otherwise return delta_angle
latlon (bool): if True, ``zenith1`` and ``zenith2`` are interpreted as
declinations rather than polar zeniths
Returns:
``numpy.ndarray``: of delta_angle or cos(delta_angle) values.
This function seems redundant but is useful in its treatment of nontrivial
numpy broadcasting use cases.
"""
return_cos = cos
sin, cos, arccos = np.sin, np.cos, np.arccos
z1, a1 = zenith1, azimuth1
z2, a2 = zenith2, azimuth2
if latlon:
sz1 = -cos(z1)
cz1 = -sin(z1)
sz2 = -cos(z2)
cz2 = -sin(z2)
else:
sz1 = sin(z1)
cz1 = cos(z1)
sz2 = sin(z2)
cz2 = cos(z2)
sa1 = sin(a1)
ca1 = cos(a1)
sa2 = sin(a2)
ca2 = cos(a2)
if np.shape(z1) != () and np.shape(z2) != () \
and np.shape(z1) != np.shape(z2):
vals = (
np.outer(sz1 * ca1, sz2 * ca2)
+ np.outer(sz1 * sa1, sz2 * sa2)
+ np.outer(cz1, cz2)
)
else:
vals = (
sz1 * ca1 * sz2 * ca2
+ sz1 * sa1 * sz2 * sa2
+ cz1 * cz2
)
if return_cos:
return vals
else:
# Account for numerical precision so acos doesn't throw warning
vals = np.clip(vals, -1, 1)
return arccos(vals)