Python's inverse formula - Vincenty does not converge (finding the distance between points on Earth)

Question

I am trying to implement the Vincenty inverse task as described in the wiki HERE

The problem is that lambda just doesn't fit. The value remains unchanged if I try to iterate over a sequence of formulas, and I really don't know why. Maybe I just stared at the obvious problem.

It should be noted that I am new to Python and still learning the language, so I'm not sure if it is abusing a language that might cause a problem, or if I have some errors in some of the calculations that I perform. I just can’t find the errors in the formulas.

Basically, I wrote the code in the format as close as possible to the wiki article, and the result is as follows:

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)

, , "" () . , .

?: -)

!

Answer 1

-, ( ):

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

// (int divisions) / (float divisions), lambda ( ):

2.5325205864224847
2.5325167509030906
2.532516759118641
2.532516759101044
2.5325167591010813
2.5325167591010813
2.5325167591010813

, , 10^(−12), , , .

(lambda ) , s.

: s .

Answer 2

, . ( .) , ; . , , . ( , , NGS, .)