Python - Vincenty's inverse formula not converging (Finding distance between points on Earth)

Question

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!

Answer 1

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.

Note: you can test your final value s here.

Answer 2

Even if it is correctly implemented, Vincenty's algorithm will fail to converge for some points. (This problem was noted by Vincenty.) I give an algorithm which is guaranteed to converge in Algorithms for geodesics; there's a python implementation available here. Finally, you can find more information on the problem at the Wikipedia page, Geodesics on an ellipsoid. (The talk page has examples of pairs of points for which Vincenty, as implemented by the NGS, fails to converge.)