我试图实现Vincenty的反问题,如wiki HERE所述
问题是lambda根本没有收敛。如果我尝试在一系列公式上迭代,值会保持不变,我真的不知道为什么。也许我只是对一个明显的问题视而不见。在
需要注意的是,我是Python新手,还在学习这门语言,所以我不确定是语言的误用导致了这个问题,还是我在执行的一些计算中确实有一些错误。我只是在公式中找不到任何错误。在
基本上,我用尽可能接近wiki文章的格式编写了代码,结果是:
import math
# Length of radius at equator of the ellipsoid
a = 6378137.0
# Flattening of the ellipsoid
f = 1/298.257223563
# Length of radius at the poles of the ellipsoid
b = (1 - f) * a
# Latitude points
la1, la2 = 10, 60
# Longitude points
lo1, lo2 = 5, 150
# For the inverse problem, we calculate U1, U2 and L.
# We set the initial value of lamb = L
u1 = math.atan( (1 - f) * math.tan(la1) )
u2 = math.atan( (1 - f) * math.tan(la2) )
L = (lo2 - lo1) * 0.0174532925
lamb = L
while True:
sinArc = math.sqrt( math.pow(math.cos(u2) * math.sin(lamb),2) + math.pow(math.cos(u1) * math.sin(u2) - math.sin(u1) * math.cos(u2) * math.cos(lamb),2) )
cosArc = math.sin(u1) * math.sin(u2) + math.cos(u1) * math.cos(u2) * math.cos(lamb)
arc = math.atan2(sinArc, cosArc)
sinAzimuth = ( math.cos(u1) * math.cos(u2) * math.sin(lamb) ) // ( sinArc )
cosAzimuthSqr = 1 - math.pow(sinAzimuth, 2)
cosProduct = cosArc - ((2 * math.sin(u1) * math.sin(u2) ) // (cosAzimuthSqr))
C = (f//16) * cosAzimuthSqr * (4 + f * (4 - 3 * cosAzimuthSqr))
lamb = L + (1 - C) * f * sinAzimuth * ( arc + C * sinArc * ( cosProduct + C * cosArc * (-1 + 2 * math.pow(cosProduct, 2))))
print(lamb)
如前所述,问题是“lamb”(lambda)值不会变小。我甚至试图将我的代码与其他实现进行比较,但它们看起来几乎一样。在
我做错什么了?:-)
谢谢大家!在
即使它被正确地实现了,Vincenty的算法也将失败 收敛于某些点。(Vincenty注意到了这个问题。) 我给出了一个保证 在Algorithms for geodesics中聚合;有一个python 实现可用here。最后,你可以找到更多 维基百科页面上关于这个问题的信息, Geodesics on an ellipsoid。(谈话页面有一些例子 由国家统计局实施的Vincenty, 无法聚合。)
首先,你也应该用弧度来转换纬度(你已经为你的经度这样做了):
一旦您这样做并去掉
^{pr2}$//
(int
)并将其替换为/
(float
个divisions),lambda
在迭代过程中停止重复相同的值,并开始遵循以下路径(基于示例坐标):正如您期望的收敛精度是
10^(−12)
,它似乎说明了这一点。在现在可以退出循环(
lambda
已经收敛)并继续下去,直到计算出所需的测地距离s
。在注意:you can test your final value ^{} here 。在
相关问题 更多 >
编程相关推荐