Питон – обратная формула Винценти не сходится (нахождение расстояния между точками на Земле)

Я пытаюсь реализовать обратную задачу Винценти, как описано в wiki ЗДЕСЬ

Проблема в том, что лямбда просто не сходится. Значение остается неизменным, если я пытаюсь перебрать последовательность формул, и я действительно не знаю, почему. Возможно, я просто уставился на очевидную проблему.

Следует отметить, что я новичок в 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) 

Как уже упоминалось, проблема в том, что значение «ягненок» (лямбда) не станет меньше. Я даже пытался сравнить свой код с другими реализациями, но они выглядели примерно одинаково.

Что я здесь делаю неправильно? 🙂

Спасибо вам всем!

Во-первых, вы также должны преобразовать ваши широты в радианах (вы уже делаете это для своих долгот):

 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 .

Примечание. Вы можете проверить свое окончательное значение здесь .

Даже если он будет правильно реализован, алгоритм Винценти не сможет сходиться для некоторых точек. (Эта проблема была отмечена Винценти.) Я даю алгоритм, который гарантированно сходится в алгоритмах для геодезических ; здесь доступна реализация python. Наконец, вы можете найти дополнительную информацию о проблеме на странице Википедии, Геодезические на эллипсоиде . (На странице обсуждения есть примеры пар точек, для которых Винценти, как реализовано NGS, не сходится.)