diff options
| -rw-r--r-- | sphere.py | 65 |
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) |