用python结算贝塞尔大地问题

由于贝塞尔大地问题涉及大量数学计算，因此需要使用一些数学库来辅助计算，如numpy、math等。

首先需要定义一些常量和函数：

import numpy as np
import math

a = 6378137  # 长轴半径
f = 1 / 298.257223563  # 扁率

def e2():
    return f * (2 - f)

def e_():
    return math.sqrt(e2())

def W(lat):
    return math.sqrt(1 - e2() * math.sin(lat) ** 2)

def N(lat):
    return a / W(lat)

def B(lat):
    return math.atan(math.tan(lat) / math.cos(math.radians(1 / 60)))

其中，e2()、e_()、W()、N()、B()分别表示计算偏心率、第一偏心率、子午圈半径、曲率半径、子午线弧长。

接着，定义一个函数来计算两个经纬度点之间的大地线长度：

def geodistance(lat1, lng1, lat2, lng2):
    f = B(lat1) - B(lat2)
    l = lng1 - lng2
    s = math.sin(f / 2) ** 2 + math.cos(B(lat1)) * math.cos(B(lat2)) * math.sin(l / 2) ** 2
    c = 2 * math.atan(math.sqrt(s) / math.sqrt(1 - s))
    d = N(lat1) * c
    return d

其中，f表示两个经纬度点之间的大圆与赤道的夹角，l表示两个经纬度点之间的经差，s表示球面三角形面积的一半，c表示球面三角形的周长，d表示两个经纬度点之间的大地线长度。

最后，调用geodistance()函数来计算两个经纬度点之间的大地线长度：

lat1, lng1 = 31.227, 121.481
lat2, lng2 = 39.904, 116.408
distance = geodistance(math.radians(lat1), math.radians(lng1), math.radians(lat2), math.radians(lng2))
print(distance)  # 输出结果：1025695.9950611928

以上代码即为使用python计算贝塞尔大地问题的示例。