python两点之间的距离公式_Python-Vincenty的逆公式不收敛（查找地球上两点之间的距离）...

I'm attempting to implement Vincenty's inverse problem as described on wiki HERE

The problem is that lambda is simply not converging. The value stays the same if I try to iterate over the sequence of formulas, and I'm really not sure why. Perhaps I've just stared myself blind on an obvious problem.

It should be noted that I'm new to Python and still learning the language, so I'm not sure if it's misuse of the language that might cause the problem, or if I do have some mistakes in some of the calculations that I perform. I just can't seem to find any mistakes in the formulas.

Basically, I've written in the code in as close of a format as I could to the wiki article, and the result is this:

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)

As mentioned the problem is that the value "lamb" (lambda) will not become smaller. I've even tried to compare my code to other implementations, but they looked just about the same.

What am I doing wrong here? :-)

Thank you all!

解决方案

First, you should convert you latitudes in radians too (you already do this for your longitudes):

u1 = math.atan( (1 - f) * math.tan(math.radians(la1)) )

u2 = math.atan( (1 - f) * math.tan(math.radians(la2)) )

L = math.radians((lo2 - lo1)) # better than * 0.0174532925

Once you do this and get rid of // (int divisions) and replace them by / (float divisions), lambda stops repeating the same value through your iterations and starts following this path (based on your example coordinates):

2.5325205864224847

2.5325167509030906

2.532516759118641

2.532516759101044

2.5325167591010813

2.5325167591010813

2.5325167591010813

As you seem to expect a convergence precision of 10^(−12), it seems to make the point.

You can now exit the loop (lambda having converged) and keep going until you compute the desired geodesic distance s.