十行代码复现 1052.7

正文第 ④ 节给出的计算链条——从两地经纬度构造单位方向向量，通过点积求解圆心角，再乘以地球平均半径求得大圆弧长——在计算机中没有任何黑箱。用 Python 与 NumPy，十行代码即可完整复现：

# 球面大圆弧距离：经纬度 → 单位向量 → 点积 → 夹角 → R·θ
import numpy as np

def unit_vec(B, L):
    B, L = np.radians(B), np.radians(L)
    return np.array([np.cos(B)*np.cos(L), np.cos(B)*np.sin(L), np.sin(B)])

R  = 6371.0
v1 = unit_vec(30.5928, 114.3055)   # 武汉
v2 = unit_vec(39.9042, 116.4074)   # 北京

theta = np.arccos(np.clip(v1 @ v2, -1, 1))
print(f"{theta:.5f} rad = {np.degrees(theta):.2f}°")
print(f"{R*theta:.1f} km")           # → 1052.7 km

代码逻辑与数学推导完全一一对应。其中唯一需要注意的工程实现细节是 np.clip(..., -1, 1)：在计算机浮点运算中，当两点重合或极为接近时，微小的舍入误差可能让点积结果变成类似 1.0000000000000002 的数值；若直接传入反余弦函数 np.arccos，会导致定义域溢出报错或返回 nan。在实际工程中用 clip 将点积强行收束在 $[-1,1]$ 内，是数值计算实务中必写的一笔防御。

完整的 Python 练习脚本（含边界对跖点测试与更多反直觉点对）收录在课程练习集 exercises/01-ellipsoid/01-01-great-circle 中。