Add sphere.py: calculate angle between points on the sphere

1 files changed, 65 insertions, 0 deletions

diff --git a/sphere.py b/sphere.py
new file mode 100644
index 0000000..f1b6460
--- /dev/null
+++ b/sphere.py
@@ -0,0 +1,65 @@
+# Copyright (c) 2017 Aaron LI
+# MIT license
+
+"""
+Spherical utilities.
+"""
+
+import numpy as np
+
+
+def central_angle(p0, points):
+    """
+    Calculate the central angle(s) between the points ``p0`` with respect to
+    the other point(s) ``points`` on the sphere.
+
+    Parameters
+    ----------
+    p0 : 2-element float tuple/list
+        (longitude/R.A., latitude/Dec.) coordinate of the reference point.
+        (Unit: deg)
+    points : 2-element float tuple/list, or 2-column float `~numpy` array
+        Coordinates of the other point(s)
+        (Unit: deg)
+
+    Returns
+    -------
+    angle : float, or 1D float `~numpy` array
+        Calculated central angle(s) (Unit: deg)
+
+    Algorithm
+    ---------
+    (radial, azimuthal, polar): (r, θ, φ)
+    central_angle: α
+    longitude/R.A.: λ = θ
+    latitude/Dec.: δ = π/2 - φ
+
+    Unit vector:
+        r1_vec = (cos(θ1)*sin(φ1), sin(θ1)*sin(φ1), cos(φ1))
+               = (cos(λ1)*cos(δ1), sin(λ1)*cos(δ1), sin(δ1))
+        r2_vec = (cos(θ2)*sin(φ2), sin(θ2)*sin(φ2), cos(φ2))
+               = (cos(λ2)*cos(δ2), sin(λ2)*cos(δ2), sin(δ2))
+
+    Therefore the angle (α) between r1_vec and r2_vec:
+        cos(α) = r1_vec * r2_vec
+               = cos(δ1)*cos(δ2)*cos(λ1-λ2) + sin(δ1)*sin(δ2)
+
+    References
+    ----------
+    [1] Spherical Coordinates - Wolfram MathWorld
+        http://mathworld.wolfram.com/SphericalCoordinates.html
+        Eq.(19)
+    [2] Great Circle - Wolfram MathWorld
+        http://mathworld.wolfram.com/GreatCircle.html
+        Eqs.(1,2,4)
+    """
+    lon0, lat0 = np.deg2rad(p0)
+    points = np.deg2rad(points)
+    try:
+        lon, lat = points[:, 0], points[:, 1]
+    except IndexError:
+        lon, lat = points  # single point
+    prod = (np.cos(lat0) * np.cos(lat) * np.cos(lon0-lon) +
+            np.sin(lat0) * np.sin(lat))
+    alpha = np.arccos(prod)
+    return np.rad2deg(alpha)