17298437009

Committed 28 Aug 2025 02:10PM UTC coverage: 55.908% (+0.04%) from 55.866%

Build # 17298437009

Build Type

push

github

Committed by

mhasself

Commit Message

quat: {rotation,decompose}_lonlat support azel option

So azel=True will flip the sign of longitude; useful for horizon
coordinates.

Run Details

6 of 6 new or added lines in 1 file covered. (100.0%)

2 existing lines in 1 file now uncovered.

1339 of 2395 relevant lines covered (55.91%)

0.56 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

95.65

/python/proj/quat.py

import numpy as np

from spt3g.core import quat, G3VectorQuat


"""We are using the spt3g quaternion containers,
i.e. cu3g.G3VectorQuat and cu3g.quat.  One way these are nice is that
they have accelerated quaternion arithmetic with operator overloading.
The component ordering is (1,i,j,k).  We are restricted to 0d or 1d,
and that's fine."""

DEG = np.pi/180


def euler(axis, angle):
    """
    The quaternion representing of an Euler rotation.
    
    For example, if axis=2 the computed quaternion(s) will have
    components:

        q = (cos(angle/2), 0, 0, sin(angle/2))
    
    Parameters
    ----------
    axis : {0, 1, 2}
        The index of the cartesian axis of the rotation (x, y, z).
    angle : float or 1-d float array
        Angle of rotation, in radians.
    
    Returns
    -------
    quat or G3VectorQuat, depending on ndim(angle).
    """
    # Either angle or axis or both can be vectors.
    angle = np.asarray(angle)
    shape = np.broadcast(axis, angle).shape + (4,)
    c, s = np.cos(angle/2), np.sin(angle/2)
    q = np.zeros(shape)
    q[..., 0] = c
    q[..., axis+1] = s
    if len(shape) == 1:
        return quat(*q)
    return G3VectorQuat(q)


def rotation_iso(theta, phi, psi=None):
    """Returns the quaternion that composes the Euler rotations:
   
        Qz(phi) Qy(theta) Qz(psi)

    Note arguments are in radians.
    """
    output = euler(2, phi) * euler(1, theta)
    if psi is None:
        return output
    return output * euler(2, psi)


def rotation_lonlat(lon, lat, psi=0., azel=False):
    """Returns the quaternion that composes the Euler rotations:
        
        Qz(lon) Qy(pi/2 - lat) Qz(psi)
    
    Note the three angle arguments are in radians.

    If azel is True, then the sign of lon is flipped (as though lon
    and lat were azimuth and elevation).

    """
    if azel:
        return rotation_iso(np.pi/2 - lat, -lon, psi)
    return rotation_iso(np.pi/2 - lat, lon, psi)

def rotation_xieta(xi, eta, gamma=0):
    """Returns the quaternion that rotates the center of focal plane to
    (xi, eta, gamma).  This is equivalent to composed Euler rotations:

        Qz(phi) Qy(theta) Qz (psi)

    where

        xi = - sin(theta) * sin(phi)
        eta = - sin(theta) * cos(phi)
        gamma = psi + phi

    Note arguments are in radians.

    """
    phi = np.arctan2(-xi, -eta)
    theta = np.arcsin(np.sqrt(xi**2 + eta**2))
    psi = gamma - phi
    return rotation_iso(theta, phi, psi)

def decompose_iso(q):
    """Decomposes the rotation encoded by q into the product of Euler
    rotations:
    
        q = Qz(phi) Qy(theta) Qz(psi)
    
    Parameters
    ----------
    q : quat or G3VectorQuat
        The quaternion(s) to be decomposed.
    
    Returns
    -------
    (theta, phi, psi) : tuple of floats or of 1-d arrays
        The rotation angles, in radians.
    """

    if isinstance(q, quat):
        a,b,c,d = q.a, q.b, q.c, q.d
    else:
        a,b,c,d = np.transpose(q)

    psi = np.arctan2(a*b+c*d, a*c-b*d)
    phi = np.arctan2(c*d-a*b, a*c+b*d)

    # There are many ways to get theta; this is probably the most
    # expensive but the most robust.
    theta = 2 * np.arctan2((b**2 + c**2)**.5, (a**2 + d**2)**.5)

    return (theta, phi, psi)

def decompose_lonlat(q, azel=False):
    """Like decompose_iso, but returns (lon, lat, psi) assuming that the
    input quaternion(s) were constructed as rotation_lonlat(lon, lat,
    psi).

    If azel is True, returns (-lon, lat, psi), so that first two
    coords are easily interpreted as an azimuth and elevation.

    """
    theta, phi, psi = decompose_iso(q)
    if azel:
        return (-phi, np.pi/2 - theta, psi)
    return (phi, np.pi/2 - theta, psi)

def decompose_xieta(q):
    """Like decompose_iso, but returns (xi, eta, gamma) assuming that the
    input quaternion(s) were constructed as rotation_xieta(xi, eta,
    gamma).

    """
    if isinstance(q, quat):
        a,b,c,d = q.a, q.b, q.c, q.d
    else:
        a,b,c,d = np.transpose(q)
    return (
        2 * (a * b - c * d),
        -2 * (a * c + b * d),
        np.arctan2(2 * a * d, a * a - d * d)
    )

1	import numpy as np	1✔
2
3	from spt3g.core import quat, G3VectorQuat	1✔
4
5
6	"""We are using the spt3g quaternion containers,	1✔
7	i.e. cu3g.G3VectorQuat and cu3g.quat. One way these are nice is that
8	they have accelerated quaternion arithmetic with operator overloading.
9	The component ordering is (1,i,j,k). We are restricted to 0d or 1d,
10	and that's fine."""
11
12	DEG = np.pi/180	1✔
13
14
15	def euler(axis, angle):	1✔
16	"""
17	The quaternion representing of an Euler rotation.
18
19	For example, if axis=2 the computed quaternion(s) will have
20	components:
21
22	q = (cos(angle/2), 0, 0, sin(angle/2))
23
24	Parameters
25	----------
26	axis : {0, 1, 2}
27	The index of the cartesian axis of the rotation (x, y, z).
28	angle : float or 1-d float array
29	Angle of rotation, in radians.
30
31	Returns
32	-------
33	quat or G3VectorQuat, depending on ndim(angle).
34	"""
35	# Either angle or axis or both can be vectors.
36	angle = np.asarray(angle)	1✔
37	shape = np.broadcast(axis, angle).shape + (4,)	1✔
38	c, s = np.cos(angle/2), np.sin(angle/2)	1✔
39	q = np.zeros(shape)	1✔
40	q[..., 0] = c	1✔
41	q[..., axis+1] = s	1✔
42	if len(shape) == 1:	1✔
43	return quat(*q)	1✔
44	return G3VectorQuat(q)	1✔
45
46
47	def rotation_iso(theta, phi, psi=None):	1✔
48	"""Returns the quaternion that composes the Euler rotations:
49
50	Qz(phi) Qy(theta) Qz(psi)
51
52	Note arguments are in radians.
53	"""
54	output = euler(2, phi) * euler(1, theta)	1✔
55	if psi is None:	1✔
UNCOV 56	return output	×
57	return output * euler(2, psi)	1✔
58
59
60	def rotation_lonlat(lon, lat, psi=0., azel=False):	1✔
61	"""Returns the quaternion that composes the Euler rotations:
62
63	Qz(lon) Qy(pi/2 - lat) Qz(psi)
64
65	Note the three angle arguments are in radians.
66
67	If azel is True, then the sign of lon is flipped (as though lon
68	and lat were azimuth and elevation).
69
70	"""
71	if azel:	1✔
72	return rotation_iso(np.pi/2 - lat, -lon, psi)	1✔
73	return rotation_iso(np.pi/2 - lat, lon, psi)	1✔
74
75	def rotation_xieta(xi, eta, gamma=0):	1✔
76	"""Returns the quaternion that rotates the center of focal plane to
77	(xi, eta, gamma). This is equivalent to composed Euler rotations:
78
79	Qz(phi) Qy(theta) Qz (psi)
80
81	where
82
83	xi = - sin(theta) * sin(phi)
84	eta = - sin(theta) * cos(phi)
85	gamma = psi + phi
86
87	Note arguments are in radians.
88
89	"""
90	phi = np.arctan2(-xi, -eta)	1✔
91	theta = np.arcsin(np.sqrt(xi2 + eta2))	1✔
92	psi = gamma - phi	1✔
93	return rotation_iso(theta, phi, psi)	1✔
94
95	def decompose_iso(q):	1✔
96	"""Decomposes the rotation encoded by q into the product of Euler
97	rotations:
98
99	q = Qz(phi) Qy(theta) Qz(psi)
100
101	Parameters
102	----------
103	q : quat or G3VectorQuat
104	The quaternion(s) to be decomposed.
105
106	Returns
107	-------
108	(theta, phi, psi) : tuple of floats or of 1-d arrays
109	The rotation angles, in radians.
110	"""
111
112	if isinstance(q, quat):	1✔
113	a,b,c,d = q.a, q.b, q.c, q.d	1✔
114	else:
UNCOV 115	a,b,c,d = np.transpose(q)	×
116
117	psi = np.arctan2(ab+cd, ac-bd)	1✔
118	phi = np.arctan2(cd-ab, ac+bd)	1✔
119
120	# There are many ways to get theta; this is probably the most
121	# expensive but the most robust.
122	theta = 2 * np.arctan2((b2 + c2).5, (a2 + d2).5)	1✔
123
124	return (theta, phi, psi)	1✔
125
126	def decompose_lonlat(q, azel=False):	1✔
127	"""Like decompose_iso, but returns (lon, lat, psi) assuming that the
128	input quaternion(s) were constructed as rotation_lonlat(lon, lat,
129	psi).
130
131	If azel is True, returns (-lon, lat, psi), so that first two
132	coords are easily interpreted as an azimuth and elevation.
133
134	"""
135	theta, phi, psi = decompose_iso(q)	1✔
136	if azel:	1✔
137	return (-phi, np.pi/2 - theta, psi)	1✔
138	return (phi, np.pi/2 - theta, psi)	1✔
139
140	def decompose_xieta(q):	1✔
141	"""Like decompose_iso, but returns (xi, eta, gamma) assuming that the
142	input quaternion(s) were constructed as rotation_xieta(xi, eta,
143	gamma).
144
145	"""
146	if isinstance(q, quat):	1✔
147	a,b,c,d = q.a, q.b, q.c, q.d	1✔
148	else:
149	a,b,c,d = np.transpose(q)	1✔
150	return (	1✔
151	2 * (a * b - c * d),
152	-2 * (a * c + b * d),
153	np.arctan2(2 * a * d, a * a - d * d)
154	)

simonsobs / so3g / 17298437009

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous