1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
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)
|
