我在 Python 3 中从事天体物理学项目。我需要计算天空中一组点的最小外接圆（如赤经和赤纬所述）。我在这里找到了一个代码（这里是实时演示）。
中实现了一种算法，而不是简单的算法，这对我来说非常有趣，因为我有大量数据集。但它是写在欧几里得平面上的。我已经尽可能地对其进行了转换，以匹配天空上的球面三角，但有一部分我不知道如何转换。
我当前的代码是球面几何和平面近似的混合。不幸的是，虽然我的物体在天空中非常接近（它们是紧凑群中的星系），但这种近似还不够好，我需要全球面计算。
此后我修改过的代码，尽可能多地添加注释（Original code、Spherical code等）来描述我做了什么以及我不知道该做什么的部分。它使用导入numpy as np和import math.
# Circumscribed circle - PARTIALLY ADAPTED TO SPHERICAL GEOMETRY!
# From https://www.nayuki.io/res/smallest-enclosing-circle/smallestenclosingcircle.py
# Demo on https://www.nayuki.io/page/smallest-enclosing-circle
def calc_sep1(a,b):
(phi1,theta1)=(a[0],a[1])
(phi2,theta2)=(b[0],b[1])
return math.acos(min(1,math.sin(theta1)*math.sin(theta2) +
math.cos(theta1)*math.cos(theta2)*math.cos(phi2 - phi1)) )
#Defines barycenter of a triangle
def calc_bary1(a,b,c):
phi = [a[0],b[0],c[0]]
theta = [a[1],b[1],c[1]]
x = np.cos(theta) * np.cos(phi)
y = np.cos(theta) * np.sin(phi)
z = np.sin(theta)
x_O, y_O, z_O = np.mean(x), np.mean(y), np.mean(z)
l = np.sqrt(x_O**2 + y_O**2 + z_O**2)
(x_M, y_M, z_M) = (x_O, y_O, z_O)/l
Dec_M = np.arcsin(z_M)
if (y_M <0):
RA_M = 2.0*np.pi - np.arccos(x_M / np.cos(Dec_M))
else:
RA_M = np.arccos(x_M / np.cos(Dec_M))
return (RA_M,Dec_M)
def midpoint(a,b):
(phi1,theta1)=(a[0],a[1])
(phi2,theta2)=(b[0],b[1])
Bx = math.cos(theta2)*math.cos(phi2-phi1)
By = math.cos(theta2)*math.sin(phi2-phi1)
theta_m = math.atan2(math.sin(theta1) + math.sin(theta2), math.sqrt((math.cos(theta1)+Bx)**2 + By**2) )
phi_m = phi1 + math.atan2(By, math.cos(theta1)+Bx)
return (phi_m,theta_m)
# Data conventions: A point is a pair of floats (x, y).
# A circle is a triple of floats (center x, center y, radius).
# Returns the smallest circle that encloses all the given points. Runs in expected O(n) time, randomized.
# Input: A sequence of pairs of floats or ints, e.g. [(0,5), (3.1,-2.7)].
# Output: A triple of floats representing a circle.
# Note: If 0 points are given, None is returned. If 1 point is given, a circle of radius 0 is returned.
#
# Initially: No boundary points known
def make_circle(points):
# Randomize order
# Carfull: inpuot and outpu in degrees, inside: radiants
shuffled = [(d2r(row['RA']), d2r(row['Dec'])) for index, row in points.iterrows()]
random.shuffle(shuffled)
# Progressively add points to circle or recompute circle
c = None
for (i, p) in enumerate(shuffled):
if c is None or not is_in_circle(c, p):
c = _make_circle_one_point(shuffled[ : i + 1], p)
return (r2d(c[0]),r2d(c[1]),c[2])
# One boundary point known
def _make_circle_one_point(points, p):
c = (p[0], p[1], 0.0)
for (i, q) in enumerate(points):
if not is_in_circle(c, q):
if c[2] == 0.0:
c = make_diameter(p, q)
else:
c = _make_circle_two_points(points[ : i + 1], p, q)
return c
# Two boundary points known
def _make_circle_two_points(points, p, q):
circ = make_diameter(p, q)
left = None
right = None
px, py = p
qx, qy = q
# For each point not in the two-point circle
for r in points:
if is_in_circle(circ, r):
continue
# Form a circumcircle and classify it on left or right side
cross = _cross_product(px, py, qx, qy, r[0], r[1])
c = make_circumcircle(p, q, r)
if c is None:
continue
elif cross > 0.0 and (left is None or _cross_product(px, py, qx, qy, c[0], c[1]) >
_cross_product(px, py, qx, qy, left[0], left[1])):
left = c
elif cross < 0.0 and (right is None or _cross_product(px, py, qx, qy, c[0], c[1]) <
_cross_product(px, py, qx, qy, right[0], right[1])):
right = c
# Select which circle to return
if left is None and right is None:
return circ
elif left is None:
return right
elif right is None:
return left
else:
return left if (left[2] <= right[2]) else right
def make_diameter(a, b):
# Original code:
# cx = (a[0] + b[0]) / 2.0
# cy = (a[1] + b[1]) / 2.0
# r0 = math.hypot(cx - a[0], cy - a[1])
# r1 = math.hypot(cx - b[0], cy - b[1])
# return (cx, cy, max(r0, r1))
# Sperical code:
(phi_m,theta_m) = midpoint(a,b)
return (phi_m,theta_m,calc_sep1(a,b)*.5)
def make_circumcircle(a, b, c):
# Mathematical algorithm from Wikipedia: Circumscribed circle
#---------------------------------------------------------------------------
# I don't know how to sphericalize this part
ox = (min(a[0], b[0], c[0]) + max(a[0], b[0], c[0])) / 2.0
oy = (min(a[1], b[1], c[1]) + max(a[1], b[1], c[1])) / 2.0
ax = a[0] - ox; ay = a[1] - oy
bx = b[0] - ox; by = b[1] - oy
cx = c[0] - ox; cy = c[1] - oy
d = (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by)) * 2.0
if d == 0.0:
return None
x = ox + ((ax * ax + ay * ay) * (by - cy) +
(bx * bx + by * by) * (cy - ay) +
(cx * cx + cy * cy) * (ay - by)) / d
y = oy + ((ax * ax + ay * ay) * (cx - bx) +
(bx * bx + by * by) * (ax - cx) +
(cx * cx + cy * cy) * (bx - ax)) / d
#---------------------------------------------------------------------------
#---------------
# Original code:
# ra = math.hypot(x - a[0], y - a[1])
# rb = math.hypot(x - b[0], y - b[1])
# rc = math.hypot(x - c[0], y - c[1])
# Spherical code:
ra = calc_sep1((x,y),a)
rb = calc_sep1((x,y),b)
rc = calc_sep1((x,y),c)
#---------------
return (x, y, max(ra, rb, rc))
_MULTIPLICATIVE_EPSILON = 1 + 1e-14
def is_in_circle(c, p):
# return c is not None and math.hypot(p[0] - c[0], p[1] - c[1]) <= c[2] * _MULTIPLICATIVE_EPSILON
return c is not None and calc_sep1(p,(c[0],c[1])) <= c[2] * _MULTIPLICATIVE_EPSILON
# Returns twice the signed area of the triangle defined by (x0, y0), (x1, y1), (x2, y2).
def _cross_product(x0, y0, x1, y1, x2, y2):
# Original code:
#return (x1 - x0) * (y2 - y0) - (y1 - y0) * (x2 - x0)
# Spherical code:
# (Uses http://mathforum.org/library/drmath/view/65316.html
# then adds sign accordingly to ordinary cross product)
side1 = calc_sep1((x0, y0), (x1, y1))
side2 = calc_sep1((x1, y1), (x2, y2))
side3 = calc_sep1((x2, y2), (x0, y0))
s = (side1+side2+side3)/2
Area = 4* math.atan(math.sqrt(math.tan(s/2)*
math.tan((s-side1)/2)*
math.tan((s-side2)/2)*
math.tan((s-side3)/2)))
Eucl_cross_prod = (x1 - x0) * (y2 - y0) - (y1 - y0) * (x2 - x0)
cross_prod_sign = 1 if (Eucl_cross_prod >=0) else -1
return 2*cross_prod_sign*Area