feat(interferometry): Import from ch_util and driftscan

This imports two general-purpose angle-related functions from other
modules using caput, with the intention of re-exporting them back
to those modules:

* `rotate_ypr`: rotates a 3-basis by a yaw, pitch, and roll.
  Originally from drifscan:
  https://github.com/radiocosmology/driftscan/blob/8def16f7d1f9d35f729305dd8998c40a7cb7aceb/drift/telescope/cylbeam.py#L44
* `sphdist`: calculates the angular distance betwee two points
  on the sphere.  Originally from ch_util:
  https://github.com/chime-experiment/ch_util/blob/c460b37b27f5bcaa1dfbe7897f3ade88e6e6775b/ch_util/ephemeris.py#L358

I have made some minor changes to the docstrings.

caput/interferometry.py

@@ -2,9 +2,11 @@
 
 Coordinates
 -----------
+- :py:meth:`sphdist`
 - :py:meth:`sph_to_ground`
 - :py:meth:`ground_to_sph`
 - :py:meth:`project_distance`
+- :py:meth:`rotate_ypr`
 
 Interferometry
 --------------
@@ -14,6 +16,40 @@
 import numpy as np
 
 
+def sphdist(long1, lat1, long2, lat2):
+    """Return the angular distance between two coordinates on the sphere.
+
+    Parameters
+    ----------
+    long1, lat1 : Skyfield Angle objects
+        longitude and latitude of the first coordinate. Each should be the
+        same length; can be one or longer.
+
+    long2, lat2 : Skyfield Angle objects
+        longitude and latitude of the second coordinate. Each should be the
+        same length. If long1, lat1 have length longer than 1, long2 and
+        lat2 should either have the same length as coordinate 1 or length 1.
+
+    Returns
+    -------
+    dist : Skyfield Angle object
+        angle between the two coordinates
+    """
+    from skyfield.units import Angle
+
+    dsinb = np.sin((lat1.radians - lat2.radians) / 2.0) ** 2
+
+    dsinl = (
+        np.cos(lat1.radians)
+        * np.cos(lat2.radians)
+        * (np.sin((long1.radians - long2.radians) / 2.0)) ** 2
+    )
+
+    dist = np.arcsin(np.sqrt(dsinl + dsinb))
+
+    return Angle(radians=2 * dist)
+
+
 def sph_to_ground(ha, lat, dec):
     """Get the ground based XYZ coordinates.
 
@@ -103,6 +139,45 @@ def projected_distance(ha, lat, dec, x, y, z=0.0):
     return dist
 
 
+def rotate_ypr(rot, xhat, yhat, zhat):
+    """Rotate a basis by a yaw, pitch and roll.
+
+    Parameters
+    ----------
+    rot : [yaw, pitch, roll]
+        Angles of rotation, in radians.
+    xhat: np.ndarray
+        X-component of the basis.  X is the axis of rotation for pitch.
+    yhat: np.ndarray
+        Y-component of the basis.  Y is the axis of rotation for roll.
+    zhat: np.ndarray
+        Z-component of the basis.  Z is the axis of rotation for yaw.
+
+    Returns
+    -------
+    xhat, yhat, zhat : np.ndarray[3]
+        New basis vectors.
+    """
+    yaw, pitch, roll = rot
+
+    # Yaw rotation
+    xhat1 = np.cos(yaw) * xhat - np.sin(yaw) * yhat
+    yhat1 = np.sin(yaw) * xhat + np.cos(yaw) * yhat
+    zhat1 = zhat
+
+    # Pitch rotation
+    xhat2 = xhat1
+    yhat2 = np.cos(pitch) * yhat1 + np.sin(pitch) * zhat1
+    zhat2 = -np.sin(pitch) * yhat1 + np.cos(pitch) * zhat1
+
+    # Roll rotation
+    xhat3 = np.cos(roll) * xhat2 - np.sin(roll) * zhat2
+    yhat3 = yhat2
+    zhat3 = np.sin(roll) * xhat2 + np.cos(roll) * zhat2
+
+    return xhat3, yhat3, zhat3
+
+
 def fringestop_phase(ha, lat, dec, u, v, w=0.0):
     """Return the phase required to fringestop. All angle inputs are radians.
 

caput/tests/test_interferometry.py

@@ -0,0 +1,61 @@
+"""Test interferometry routines."""
+
+import pytest
+from math import pi, sqrt
+
+from caput import interferometry
+
+
+def test_sphdist():
+    from skyfield.units import Angle
+
+    # 90 degrees
+    assert interferometry.sphdist(
+        Angle(radians=0), Angle(radians=0), Angle(radians=pi / 2), Angle(radians=0)
+    ).radians == pytest.approx(pi / 2)
+    assert interferometry.sphdist(
+        Angle(radians=0), Angle(radians=0), Angle(radians=0), Angle(radians=pi / 2)
+    ).radians == pytest.approx(pi / 2)
+
+    # 60 degrees
+    assert interferometry.sphdist(
+        Angle(radians=0), Angle(radians=0), Angle(radians=pi / 4), Angle(radians=pi / 4)
+    ).radians == pytest.approx(pi / 3)
+    assert interferometry.sphdist(
+        Angle(radians=pi / 4), Angle(radians=pi / 4), Angle(radians=0), Angle(radians=0)
+    ).radians == pytest.approx(pi / 3)
+
+
+def test_rotate_ypr():
+
+    # No rotation
+    basis = interferometry.rotate_ypr([0, 0, 0], 1, 2, 3)
+    assert basis == pytest.approx([1, 2, 3])
+
+    # Rotate into +Y two ways
+    basis = interferometry.rotate_ypr([pi / 2, 0, 0], 1, 0, 0)
+    assert basis == pytest.approx([0, 1, 0])
+
+    basis = interferometry.rotate_ypr([0, pi / 2, 0], 0, 0, 1)
+    assert basis == pytest.approx([0, 1, 0])
+
+    # General rotation
+    x0 = sqrt(2) / 2
+    y0 = 0.5
+    z0 = 0.5
+    basis = interferometry.rotate_ypr([pi / 2, pi / 3, pi / 4], x0, y0, z0)
+
+    # The calculation goes like this ("v" denotes a square root):
+
+    # After yawing by 90 degrees:
+    # x1 = -y0               y1 = x0                     z1 = z0
+
+    # Then pitching by 60 degrees:
+    # x2 = x1                y2 = 0.5 * y1 + v3/2 z1     z2 = -v3/2 y1 + 0.5 z1
+
+    # Then rolling by 45 degrees:
+    # x3 = v2/2 (x2 - z2)    y3 = y2                     z3 = v2/2 (x2 + z2)
+
+    assert basis[0] == pytest.approx(sqrt(2) / 2 * (-y0 + sqrt(3) * x0 / 2 - z0 / 2))
+    assert basis[1] == pytest.approx(x0 / 2 + sqrt(3) * z0 / 2)
+    assert basis[2] == pytest.approx(sqrt(2) / 2 * (-y0 - sqrt(3) * x0 / 2 + z0 / 2))