Изобретаю свой сложный способ поиска координат точки пересечения двух линий / forpes.ru

Главная
Изобретаю свой сложный способ поиска координат точки пересечения двух линий

Изобретаю свой сложный способ поиска координат точки пересечения двух линий +3

10.09.2024 12:30

velon 25 1800 Источник

Начну с громкого заявления: я придумал способ найти точку пересечения двух отрезков, заданных координатами концов. Придумал давно, лет 7 назад, в 2017 году, примерно, да, путь к этой публикации был долог, в основном из‑за лени.

И да, я его придумал потому что не смог нагуглить, может он где‑то и без меня описан был, может за эти 7 лет кто‑то написал что‑то похожее, а может я придумал какую‑то фигню, которую умные люди изобретать не станут...

Вы, возможно, скажете, что нет тут ничего сложного, открываешь Википедию: пересечение прямых и смотришь. Хм, даже на Хабр об этом уже писали: Нахождение точки пересечения двух прямых (и отрезков). Но не всё так однозначно... недаром на КДПВ изображена Земля. Интересует меня именно пересечение прямых не на плоскости, а на поверхности Земли. Если уж пошла такая пьянка, то прямые вовсе не прямые, а ортодромия. Об этом тоже писали на Хабре и не раз, например вот: Занимательная геодезия.

В общем, Земля у нас не плоская! Была бы плоская, уже бы закончили статью и не мучились, и не начинали бы.

Поскольку Земля у нас в лучшем случае круглая, а может даже сферическая в вакууме, то применим Сферическую Тригонометрию. Здесь я ничего нового не придумывал — здесь уже всё придумано до нас, просто опишу, то что удалось нагуглить. Такое длительное вступление поможет лучше проникнуться проблемой.

Вообще уже на этом шаге будет уже много мороки, количество формул и тригонометрии в публикации зашкаливает, даже я заколебался писать их в TeX, спасибо@NewTechAuditи публикации Как красиво писать формулы c LaTeX, — без этой статьи я бы не справился.

Но на сферической тригонометрии мы не остановимся, она конечно позволит найти аналитическое решение для ситуации, когда полярный и экваториальный радиусы Земли равны (возможно для Ваших задач такое допущение приемлемо, — я не настаиваю), но мы то знаем, что Земля не только не плоская, и имеет форму шара, так она ещё и сплюснута на полюсах — к‑к-к‑комбо.

Усложнили себе жизнь на ровном месте, простите за плоскоземельный каламбур.

В общем, аналитического решения с формулами, для Сплюснутой Земли — нет. С этим вопросом я вышел в интернет... и ничего не нашёл, кроме решения главных геодезических задач — это близко, но не оно. Я вышел из интернета. Однако придумал способ применить главные геодезические задачи для решения своей проблемы, — вот его и опишу, помимо прочего.

На начнём всё же со Сферической Земли, без сплющивания. Как я уже писал, возможно, для Вашей задачи этого будет достаточно.

В процессе чтения этой статьи могут возникнуть МатАн'овские или ВышМат'овские флешбеки, — мужайтесь, но примеры кода приложу, на C#. Вероятно код будет проще воспринимать, чем формулы, кстати, вот он, на GitHub: Geodesic.

А теперь, приступим.

Расчёты
Итоги

Расчёты

И так, нам нужно найти координаты точки пересечения. Решаем задачу на сфероиде, у которого полярный и экваториальный радиусы равны.

Начнём с того, что на Википедии есть формула широты промежуточной точки как функция долготы:

$\varphi=\arctan\left(\frac{\tan\varphi_1*\sin\left(\lambda_2-\lambda\right)}{\sin\left(\lambda_2-\lambda_1\right)}+\frac{\tan\varphi_2*\sin\left(\lambda-\lambda_1\right)}{\sin\left(\lambda_2-\lambda_1\right)}\right)$

Здесь:

$\varphi_1$ и $\lambda_1$ — широта и долгота начальной точки отрезка

$\varphi_2$ и $\lambda_2$ — широта и долгота конечной точки отрезка

$\varphi$ и $\lambda$ — широта и долгота промежуточной точки

Т.е. формула позволяет найти широту $\varphi$ промежуточной точки, в интересующей нас долготе $\lambda$ . Начальные и конечные точки отрезка считаем известными.

Скрытый текст

private double GetLatitudeSpheroid(double longitude, Point coord1, Point coord2)
{
    return Math.Atan(Math.Tan(coord1.LatR) / Math.Sin(coord2.LonR - coord1.LonR) *
                         Math.Sin(coord2.LonR - longitude) +
                         Math.Tan(coord2.LatR) / Math.Sin(coord2.LonR - coord1.LonR) *
                         Math.Sin(longitude - coord1.LonR)) * 180 / Math.PI;
}

Метод на GitHub: IntermediatePointService.cs

Наоборот — найти долготу по широте, не нашёл формулу, но если не лень, можем её вывести сами.

Вывод формулы долготы промежуточной точки как функции широты

Для начала введём обозначения для удобства:

$\begin{align} A_1=\frac{\tan\varphi_1}{\sin\left(\lambda_2-\lambda_1\right)} \\ A_2=\frac{\tan\varphi_2}{\sin\left(\lambda_2-\lambda_1\right)} \end{align}$

Представление на C#:

double a1 = Math.Tan(coord1.LatR) / Math.Sin(coord2.LonR - coord1.LonR);
double a2 = Math.Tan(coord2.LatR) / Math.Sin(coord2.LonR - coord1.LonR);

Здесь все известные нам величины. Формулу широты промежуточной точки теперь можно переписать вот так (будем искать тангенс широты, чтобы избавиться от арктангенса):

$\tan\varphi=A_1\sin\left(\lambda_2-\lambda\right)+A_2\sin\left(\lambda-\lambda_1\right)$

Теперь мы можем раскрыть синус разности, по формуле:

$\sin\left(\alpha\pm\beta\right)=\sin\alpha*\cos\beta\pm\cos\alpha*\sin\beta$

Получим:

$\tan\varphi=A_1\left(\sin\lambda_2\cos\lambda-\cos\lambda_2\sin\lambda\right)+A_2\left(\sin\lambda\cos\lambda_1-\cos\lambda\sin\lambda_1\right)$

Раскроем скобки:

$\tan\varphi=A_1\sin\lambda_2\cos\lambda-A_1\cos\lambda_2\sin\lambda+A_2\sin\lambda\cos\lambda_1-A_2\cos\lambda\sin\lambda_1$

Приведём подобные и снова занесём всё в скобки:

$\tan\varphi=\cos\lambda\left(A_1\sin\lambda_2-A_2\sin\lambda_1\right)+\sin\lambda\left(A_2\cos\lambda_1-A_1\cos\lambda_2\right)$

Для вычисления той части выражения, что оказалась внутри скобок, у нас есть всё необходимое — все величины известны, поэтому можем ввести для них буквенные обозначения:

$\begin{align} B=A_1\sin\lambda_2-A_2\sin\lambda_1 \\ C=A_2\cos\lambda_1-A_1\cos\lambda_2 \\ \end{align}$

Оно же в виде программного кода:

double b = a1 * Math.Sin(coord2.LonR) - a2 * Math.Sin(coord1.LonR);
double c = a2 * Math.Cos(coord1.LonR) - a1 * Math.Cos(coord2.LonR);

Подставим эти обозначения в формулу, получим:

$\tan\varphi=B\cos\lambda+C\sin\lambda$

Существует такое представление:

$X\sin\alpha+Y\cos\alpha=\sqrt{X^2+Y^2}*\sin\left(\alpha+\phi\right)$

где угол $\phi$ находится из соотношения:

$\left\{ \begin{array}{1} \sin\phi=\frac{Y}{\sqrt{X^2+Y^2}} \\ \cos\phi=\frac{X}{\sqrt{X^2+Y^2}} \\ \end{array} \right.$

Её тоже можно посмотреть на Википедии, в разделе суммы, только А и B, заменил на X и Y, что бы не было путаницы с ведёнными нами обозначениями.

С принятыми нами обозначениями получим:

$\begin{align} \tan\varphi=\sin\left(\lambda+\phi\right)\sqrt{C^2+B^2} \\ \sin\left(\lambda+\phi\right)=\frac{\tan\varphi}{\sqrt{C^2+B^2}} \\ \end{align}$

или

$\lambda+\phi=\arcsin\left(\frac{\tan\varphi}{\sqrt{C^2+B^2}}\right)$

Соответствующая строка в программе:

double lonfet = Math.Asin(Math.Tan(latitude) / Math.Sqrt(b * b + c * c));

где $\lambda$ — долгота промежуточной точки, а $\phi$ получается из соотношения:

$\left\{ \begin{array}{1} \sin\phi=\frac{B}{\sqrt{C^2+B^2}} \\ \cos\phi=\frac{C}{\sqrt{C^2+B^2}} \\ \end{array} \right.$

Теперь программная реализация.

Скрытый текст

private double GetLongitudeSpheroid(double latitude, Point coord1, Point coord2)
{
    double a1 = Math.Tan(coord1.LatR) / Math.Sin(coord2.LonR - coord1.LonR);
    double a2 = Math.Tan(coord2.LatR) / Math.Sin(coord2.LonR - coord1.LonR);
    double b = a1 * Math.Sin(coord2.LonR) - a2 * Math.Sin(coord1.LonR);
    double c = a2 * Math.Cos(coord1.LonR) - a1 * Math.Cos(coord2.LonR);

    double lonfet = Math.Asin(Math.Tan(latitude) / Math.Sqrt(b * b + c * c));

    // В зависимости от значений c, b и угла - разные преобразования ответа

    var angle = Math.Atan((coord2.LatR - coord1.LatR) / (coord2.LonR - coord1.LonR)) * 180 / Math.PI;

    if (c < 0 && b < 0 && angle > 0)
        return (lonfet + Math.Atan(c / b)) * 180 / Math.PI + 90;
    if (c > 0 && b > 0 && angle > 0)
        return (lonfet + Math.Atan(c / b)) * 180 / Math.PI - 90;
    if (c < 0 && b < 0 && angle < 0)
        return -90 - (lonfet - Math.Atan(c / b)) * 180 / Math.PI;
    if (c > 0 && b > 0 && angle < 0)
        return 90 - (lonfet - Math.Atan(c / b)) * 180 / Math.PI;

    if (Math.Abs(b) < 0.00000001 && angle > 0)
        return -90 + (lonfet + Math.Atan(c / b)) * 180 / Math.PI;
    if (Math.Abs(b) < 0.00000001 && angle < 0)
        return 90 - (lonfet - Math.Atan(c / b)) * 180 / Math.PI;

    if (c * angle > 0)
    {
        if (c > 0)
            return (lonfet + Math.Atan(c / b)) * 180 / Math.PI + 90; 
        return 90 - (lonfet - Math.Atan(c / b)) * 180 / Math.PI;
    }
    else
    {
        if (c > 0)
            return -90 - (lonfet - Math.Atan(c / b)) * 180 / Math.PI;
        return -90 + (lonfet + Math.Atan(c / b)) * 180 / Math.PI;
    }
}

Метод на GutHub: IntermediatePointService.cs

Но вернёмся к поиску пересечения.

Поиск координат точки пересечения линий на сфероиде

Нам же нужно найти точку пересечения двух линий, логично, что в точке пересечения значения широт двух линий должны совпасть, как и значения долготы, иначе это не пересечение вовсе. Вооружившись сим размышлением и формулой расчёта широты, можем вывести равенство:

$\begin{align} \frac{\tan\varphi_1*\sin\left(\lambda_2-\lambda\right)+\tan\varphi_2*\sin\left(\lambda-\lambda_1\right)}{\sin\left(\lambda_2-\lambda_1\right)}=\\ \frac{\tan\varphi_3*\sin\left(\lambda_4-\lambda\right)+\tan\varphi_4*\sin\left(\lambda-\lambda_3\right)}{\sin\left(\lambda_4-\lambda_3\right)} \end{align}$

Здесь:

$\{\varphi_1;\lambda_1\}$ и $\{\varphi_2;\lambda_2\}$ — широты и долготы граничных точек первой линии.

$\{\varphi_3;\lambda_3\}$ и $\{\varphi_4;\lambda_4\}$ — широты и долготы граничных точек второй линии.

$\lambda$ — долгота точки пересечения

В этом уравнении у нас одна неизвестная — долгота точки пересечения $\lambda$ , остальные должны быть известны изначально. Если мы её сможем вывести, то подставив в формулу широты, мы получим долготу и широту точки пересечения.

Выведем из получившегося равенства долготу.

Для удобства введём следующие обозначения:

$\begin{align} A=\frac{\tan\varphi_1}{\sin\left(\lambda_2-\lambda_1\right)} \\ B=\frac{\tan\varphi_2}{\sin\left(\lambda_2-\lambda_1\right)} \\ C=\frac{\tan\varphi_3}{\sin\left(\lambda_4-\lambda_3\right)} \\ D=\frac{\tan\varphi_4}{\sin\left(\lambda_4-\lambda_3\right)} \\ \end{align}$

Т.е. введём буквенные обозначения для известной части уравнения, которые мы можем вычислить просто зная координаты границ отрезков (погодите, пока рано это делать), теперь равенство можно записать следующим образом:

$A*\sin\left(\lambda_2-\lambda\right)+B*\sin\left(\lambda-\lambda_1\right)=C*\sin\left(\lambda_4-\lambda\right)+D*\sin\left(\lambda-\lambda_3\right)$

Или перенесём всё в левую часть:

$A*\sin\left(\lambda_2-\lambda\right)+B*\sin\left(\lambda-\lambda_1\right)-C*\sin\left(\lambda_4-\lambda\right)-D*\sin\left(\lambda-\lambda_3\right)=0$

Преобразуем тригонометрические функцию разности двух аргументов, воспользовавшись формулами сложения (эта формула уже использовалась выше, но приведу её ещё раз):

$\sin\left(\alpha\pm\beta\right)=\sin\alpha*\cos\beta\pm\cos\alpha*\sin\beta$

Получим:

$\begin{align} &A*\left(\sin\lambda_2\cos\lambda-\cos\lambda_2\sin\lambda\right) \\ &+B*\left(\sin\lambda\cos\lambda_1-\cos\lambda\sin\lambda_1\right) \\ &-C*\left(\sin\lambda_4\cos\lambda-\cos\lambda_4\sin\lambda\right) \\ &-D*\left(\sin\lambda\cos\lambda_3-\cos\lambda\sin\lambda_3\right)=0 \\ \end{align}$

Мы можем вновь раскрыть скобки, например (для части выражения):

$\begin{align} &A*\left(\sin\lambda_2\cos\lambda-\cos\lambda_2\sin\lambda\right) \\ &=A\sin\lambda_2*A\cos\lambda-A\cos\lambda_2*A\sin\lambda \end{align}$

И ввести буквенные обозначения для тех частей уравнения, в которых у нас изначально есть всё для вычислений т.е. нет неизвестной, как например для $A\sin\lambda_2$ .

$\begin{align} &A_1=A*\sin\lambda_2=\frac{\tan\varphi_1}{\sin\left(\lambda_2-\lambda_1\right)}*\sin\lambda_2 \\ &A_2=A*\cos\lambda_2=\frac{\tan\varphi_1}{\sin\left(\lambda_2-\lambda_1\right)}*\cos\lambda_2 \\ &B_1=B*\cos\lambda_1=\frac{\tan\varphi_2}{\sin\left(\lambda_2-\lambda_1\right)}*\cos\lambda_1 \\ &B_2=B*\sin\lambda_1=\frac{\tan\varphi_2}{\sin\left(\lambda_2-\lambda_1\right)}*\sin\lambda_1 \\ &C_1=C*\sin\lambda_4=\frac{\tan\varphi_3}{\sin\left(\lambda_4-\lambda_3\right)}*\sin\lambda_4 \\ &C_2=C*\cos\lambda_4=\frac{\tan\varphi_3}{\sin\left(\lambda_4-\lambda_3\right)}*\cos\lambda_4 \\ &D_1=D*\cos\lambda_3=\frac{\tan\varphi_4}{\sin\left(\lambda_4-\lambda_3\right)}*\cos\lambda_3 \\ &D_2=D*\sin\lambda_3=\frac{\tan\varphi_4}{\sin\left(\lambda_4-\lambda_3\right)}*\sin\lambda_3 \\ \end{align}$

Всё это в виде программного кода:

double dl1 = Math.Sin(coord12.LonR - coord11.LonR);
double dl2 = Math.Sin(coord22.LonR - coord21.LonR);

// Первая ортодрома
double a1S = Math.Tan(coord11.LatR) / dl1 * Math.Sin(coord12.LonR);
double a2S = Math.Tan(coord12.LatR) / dl1 * Math.Cos(coord11.LonR);
double a1Ss = Math.Tan(coord11.LatR) / dl1 * Math.Cos(coord12.LonR);
double a2Ss = Math.Tan(coord12.LatR) / dl1 * Math.Sin(coord11.LonR);

// Вторая ортодрома
double a3S = Math.Tan(coord21.LatR) / dl2 * Math.Sin(coord22.LonR);
double a4S = Math.Tan(coord22.LatR) / dl2 * Math.Cos(coord21.LonR);
double a3Ss = Math.Tan(coord21.LatR) / dl2 * Math.Cos(coord22.LonR);
double a4Ss = Math.Tan(coord22.LatR) / dl2 * Math.Sin(coord21.LonR);

Выражение приобретёт следующий вид:

$\begin{align} &A_1\cos\lambda-A_2\sin\lambda+B_1\sin\lambda-B_2\cos\lambda \\ &-C_1\cos\lambda+C_2\sin\lambda-D_1\sin\lambda+D_2\cos\lambda=0 \\ \end{align}$ $\left(A_1-B_2-C_1+D_2\right)\cos\lambda + \left(B_1-A_2+C_2-D_1\right)\sin\lambda = 0$

поделим на $\cos\lambda$ :

$\left(B_1-A_2+C_2-D_1\right)\tan\lambda=B_2+C_1-D_2-A_1$ $\tan\lambda=\frac{B_2+C_1-D_2-A_1}{B_1-A_2+C_2-D_1}$

или

$\lambda=\arctan\left(\frac{B_2+C_1-D_2-A_1}{B_1-A_2+C_2-D_1}\right)$

Соответствует коду:

double b1 = a1S - a2Ss - a3S + a4Ss;
double b2 = a2S - a1Ss + a3Ss - a4S;

var answer = Math.Atan(-b1 / b2) * 180 / Math.PI;

Код метода целиком:

Скрытый текст

private double SpheroidLongitude(Point coord11, Point coord12, Point coord21, Point coord22)
{
    double dl1 = Math.Sin(coord12.LonR - coord11.LonR);
    double dl2 = Math.Sin(coord22.LonR - coord21.LonR);

    // Первая ортодрома
    double a1S = Math.Tan(coord11.LatR) / dl1 * Math.Sin(coord12.LonR);
    double a2S = Math.Tan(coord12.LatR) / dl1 * Math.Cos(coord11.LonR);
    double a1Ss = Math.Tan(coord11.LatR) / dl1 * Math.Cos(coord12.LonR);
    double a2Ss = Math.Tan(coord12.LatR) / dl1 * Math.Sin(coord11.LonR);

    // Вторая ортодрома
    double a3S = Math.Tan(coord21.LatR) / dl2 * Math.Sin(coord22.LonR);
    double a4S = Math.Tan(coord22.LatR) / dl2 * Math.Cos(coord21.LonR);
    double a3Ss = Math.Tan(coord21.LatR) / dl2 * Math.Cos(coord22.LonR);
    double a4Ss = Math.Tan(coord22.LatR) / dl2 * Math.Sin(coord21.LonR);

    double b1 = a1S - a2Ss - a3S + a4Ss;
    double b2 = a2S - a1Ss + a3Ss - a4S;

    var answer = Math.Atan(-b1 / b2) * 180 / Math.PI;
    if (answer < 0 && coord11.Longitude > 0 && coord12.Longitude > 0)
        return 180 + answer;
    if (answer > 0 && coord11.Longitude < 0 && coord12.Longitude < 0)
        return -180 + answer;

    return answer;
}

Метод на GitHub: IntersectService.cs

Таким образом мы можем найти долготу точки пересечения, найти широту можно подставив значение в функцию широты. У нас есть две линии, в функции вычисления широты по долготе, используется одна, поэтому в качестве $\varphi_1$ , $\lambda_1$ и $\varphi_2$ , $\lambda_2$ можно использовать координаты любой из двух пересекающихся линий - значение совпадёт. А можно использовать вторую пару точек для контрольной проверки.

var inLat1 = _interPoint.GetLatitude(inLon, coord11, coord12);
var inLat2 = _interPoint.GetLatitude(inLon, coord21, coord22);

Вот так мы вычислили долготу и широту точки пересечения двух линий на сфероиде - сферическая точка пересечения найдена.

Со сфероидом всё хорошо, есть аналитическое решение, осталось только определиться с эллипсоидом вращения. И вот тут аналитического решения не будет, - придётся изобретать алгоритм.

Решая задачу для сфероида, мы начали с того, что у нас была формула широты промежуточной точки по известному значению долготы, да эллипсоида так же начнём с этой задачи.

Определение широты или долготы промежуточной точки на эллипсоиде

Для этого предлагаю использовать прямую и обратную геодезические задачи:

Прямая геодезическая задача: по известной координате начала, направлению и расстоянию, позволяет найти координаты конца отрезка.
Обратная геодезическая задача: по известным координатам двух точек вычисляется расстояние и угол.

В нашем случае как раз известны координаты начала $\left(\varphi_1,\lambda_1\right)$ и координаты конца отрезка $\left(\varphi_2,\lambda_2\right)$ . И нам нужно либо найти на этой линии долготу $\lambda$ по известной широте $\varphi$ , либо найти широту $\varphi$ , по известной долготе $\lambda$ . Давайте, для примера, по известной долготе, найдём значение широты, алгоритм таков:

Используя обратную геодезическую задачу, определим расстояние между двумя известными точками и угол.
Полученное расстояние разделим пополам, и, применив прямую геодезическую задачу, по найденному ранее углу, и полученному расстоянию, найдём координаты средней точки.
ЕСЛИ долгота этой точки равна целевой долготе, пусть даже с некоторой погрешностью, ТО и её широта — наш ответ, алгоритм ЗАВЕРШАЕТСЯ.
ИНАЧЕ мы проверяем: целевая долгота ближе в началу отрезка, или ближе к концу:
1. ЕСЛИ ближе к началу, ТО алгоритм рекурсивно продолжается с ШАГА 1, где в качестве начальной точки используется начальная точка, а в качестве конечной — средняя.
2. ЕСЛИ ближе к концу, ТО алгоритм рекурсивно продолжается с ШАГА 1, где в качестве начальной точки используется средняя точка, а в качестве конечной — конечная.

Этот алгоритм будет рекурсивно продолжаться пока мы не достигнем на ШАГЕ 3 заданной точности.

Собственно для этого подхода не важно, ищем ли мы широту по долготе, либо долготу по широте. Подход итеративный, т.е. на каждом шаге мы асимптотически приближаемся к ответу, однако точный ответ может быть недостижим, поэтому алгоритм стоит прервать при достижении результата с некоторой погрешностью.

И так, мы научились на эллипсоиде находить промежуточную точку ортодромии (в границах двух координат), с помощью итеративного алгоритма и геодезических задач, однако мы хотели найти точку пересечения двух ортодромий, для решения нам опять же понадобится итеративный алгоритм.

Поиск координат точки пересечения линий на эллипсоиде

И так нам известны координаты первой точки $\left(\varphi_1,\lambda_1\right), \left(\varphi_2,\lambda_2\right)$ и координаты второй точки $\left(\varphi_3,\lambda_3\right), \left(\varphi_4,\lambda_4\right)$ . Требуется найти точку пересечения.

Отсортируем полученные 4 точки, например по долготе, и возьмём две средних, т.е. те, которые в отсортированном списке оказались на втором и третьем местах. Значения их долготы обозначим $\lambda_f$ и $\lambda_s$

По $\lambda_f$ и $\lambda_s$ найдём среднюю арифметическую и обозначим $\lambda_m$ .

У нас есть две линии, и выше мы уже описали алгоритм нахождения значения широты, по известному значению долготы, поэтому вычислим для обеих линий значения широты соответствующих значениям долготы $\lambda_f$ , $\lambda_m$ , $\lambda_s$ , получим такой набор точек для первой и для второй линии соответственно:

$\begin{align} \left\{\varphi_{f1},\lambda_f\right\},\left\{\varphi_{m1},\lambda_m\right\},\left\{\varphi_{s1},\lambda_s\right\} \\ \left\{\varphi_{f2},\lambda_f\right\},\left\{\varphi_{m2},\lambda_m\right\},\left\{\varphi_{s2},\lambda_s\right\} \\ \end{align}$

ЕСЛИ какая-либо пара координат первой ортодромии совпала (опять же с некоторой погрешностью) с какой-либо парой координат второй, ТО это и есть точка пересечения. И алгоритм ЗАВЕРШАЕТСЯ. ИНАЧЕ переходим на следующий шаг.
Возьмём попарно широты из двух разных линий отложенные на одной долготе: $\varphi_{f1}$ и $\varphi_{f2}$ , $\varphi_{m1}$ и $\varphi_{m2}$ , $\varphi_{s1}$ и $\varphi_{s2}$ , вычитаем их попарно друг из друга:

$\begin{align} &\Delta\varphi_f=\varphi_{f1}-\varphi_{f2} \\ &\Delta\varphi_m=\varphi_{m1}-\varphi_{m2} \\ &\Delta\varphi_s=\varphi_{s1}-\varphi_{s2} \\ \end{align}$

Теперь попарно сравним $\Delta\varphi_f$ с $\Delta\varphi_m$ и $\Delta\varphi_m$ с $\Delta\varphi_s$
1. ЕСЛИ у $\Delta\varphi_f$ и $\Delta\varphi_m$ разные знаки, ТО $\lambda_s=\lambda_m$
2. ЕСЛИ у $\Delta\varphi_m$ и $\Delta\varphi_s$ разные знаки, ТО $\lambda_f=\lambda_m$
  
  Различие знаков говорит о том, что пересечение произошло на этом участке, между соответствующими значениями долготы. Алгоритм рекурсивно продолжается, переходим на ШАГ 2.

$\left\{ \begin{array}{1} \lambda_s=\lambda_m\text{, если }&\Delta\varphi_f*\Delta\varphi_m<0 \\ \lambda_f=\lambda_m\text{, если }&\Delta\varphi_m*\Delta\varphi_s<0 \\ \end{array} \right.$

Алгоритм так же на каждой итерации приближается к ответу, рано или поздно, с некоторой точностью он найдёт решение, и это будут координаты точки пересечения двух линий.

Ура, мы нашли решение, которое позволяет найти точку пересечения линий для эллипсоида. Это собственно и было изначальной задачей.

Для поиска пересечений мы активно применяем прямую и обратную геодезические задачи, однако за скобками осталось само решение этих геодезических задач. Для полноты изложения необходимо привести и их.

Решение геодезических задач на сфероиде

Решение геодезической задачи на сфероиде вероятно можно загуглить, так я и сделал когда‑то, и нашёл решение на сайте geogr.msu.ru в виде pdf файла, однако со временем ссылка стала битой, а веб‑архив ничего не сохранил, аналогичного решения сходу найти не получилось, поэтому опишу вычисления так, как я их помню и могу воспроизвести.

С этим вопросом я выходил в интернет, тогда интернет помог, а теперь уже - нет

Прямая геодезическая задача

Как выше уже говорилось, нам известны координаты начальной точки $\left(\varphi_1,\lambda_1\right)$ , азимут $\alpha_1$ , т.е. направление к конечной точке, и расстояние до конечной точки. Необходимо найти эту конечную точку $\left(\varphi_2,\lambda_2\right)$ , и обратный азимут $\alpha_2$ .

Примем за радиус сферы.

Широту определим используя теоремы сферической тригонометрии (теорему косинусов).

$\sin\varphi_2=\sin\varphi_1\cos\sigma+\cos\varphi_1\sin\sigma\cos\alpha_1$

где $\sigma$ — длинна , выраженная в долях радиуса сферы $\sigma=s/R$

var sigma = s / _ellipsoid.PolarRadius;

Теперь определим долготу.

Из теоремы котангенсов:

$\cot\sigma\cos\varphi_1=\sin\varphi_1\cos\alpha_1+\sin\alpha_1\cot\left(\lambda_2-\lambda_1\right)$

Выполним некоторые преобразования:

$\sin\alpha_1\cot\left(\lambda_2-\lambda_1\right)=\cot\sigma\cos\varphi_1-\sin\varphi_1\cos\alpha_1$ $\cot\left(\lambda_2-\lambda_1\right)=\frac{\cot\sigma\cos\varphi_1-\sin\varphi_1\cos\alpha_1}{\sin\alpha_1}$ $\tan\left(\lambda_2-\lambda_1\right)=\frac{\sin\alpha_1}{\cot\sigma\cos\varphi_1-\sin\varphi_1\cos\alpha_1}$

Домножим всё на $\sin\sigma$ чтобы избавиться от всяких котангенсов в знаменателе:

$\tan\left(\lambda_2-\lambda_1\right)=\frac{\sin\alpha_1\sin\sigma}{\cos\sigma\cos\varphi_1-\sin\sigma\sin\varphi_1\cos\alpha_1}$

В коде:

var lambda = Math.Atan(Math.Sin(sigma) * Math.Sin(a1) /
                                   (Math.Cos(sigma) * Math.Cos(lat1) - Math.Sin(sigma) * Math.Sin(lat1) * Math.Cos(a1)));

Здесь $\lambda_2-\lambda_1$ это разность долгот, второй и первой точек, зная долготу первой, мы можем теперь вычислить долготу второй.

var lon2 = -lambda + coord.LonR;

Осталось вычислить обратный азимут $\alpha_2$ , по теореме котангенсов:

$\tan\varphi_1\sin\sigma=\cos\sigma\cos\alpha_1-\sin\alpha_1\cot\alpha_2$

Выполним всё те же некоторые преобразования, аналогичные тем что были выше, я уже не буду их расписывать, получим:

$\tan\alpha_2=\frac{\cos\varphi_1\sin\alpha_1}{\cos\varphi_1\cos\sigma\cos\alpha_1-\sin\varphi_1\sin\sigma}$

Соответствующий код:

 var a2 = -Math.Atan(Math.Cos(lat1) * Math.Sin(a1) /
                                (Math.Cos(lat1) * Math.Cos(sigma) * Math.Cos(a1) - Math.Sin(lat1) * Math.Sin(sigma))) *
                     180 / Math.PI;

Ну и метод полностью:

Скрытый текст

private DirectProblemAnswer DirectProblemSpheroid(Point coord, double a1, double s)
{
    var aDegree = a1;
    a1 = a1 * Math.PI / 180;
    var lat1 = coord.LatR;

    var sigma = s / _ellipsoid.PolarRadius;
    var lat2 = Math.Asin(Math.Sin(lat1) * Math.Cos(sigma) + Math.Cos(lat1) * Math.Sin(sigma) * Math.Cos(a1));

    var lambda = Math.Atan(Math.Sin(sigma) * Math.Sin(a1) /
                           (Math.Cos(sigma) * Math.Cos(lat1) - Math.Sin(sigma) * Math.Sin(lat1) * Math.Cos(a1)));
    var lon2 = -lambda + coord.LonR;

    var a2 = -Math.Atan(Math.Cos(lat1) * Math.Sin(a1) /
                        (Math.Cos(lat1) * Math.Cos(sigma) * Math.Cos(a1) - Math.Sin(lat1) * Math.Sin(sigma))) *
             180 / Math.PI;
    

    lon2 = lon2 * 180 / Math.PI;
    lat2 = lat2 * 180 / Math.PI;

    var point1 = IsPolisIntersect(coord, aDegree, s)
        ? new Point(lon2 + (lon2 < 0 ? 180 : -180), lat2)
        : new Point(lon2, lat2);

    a2 = Azimuth.AzimuthRecovery(point1, new Point(coord.Longitude, lat1 * 180 / Math.PI), a2);

    return  new DirectProblemAnswer(point1, a2);
}

Метод на GitHub: DirectProblemService.cs

Мы привели все формулы, необходимые для решения прямой геодезической задачи на сфере, применяя сферическую тригонометрию.

Обратная геодезическая задача

Обратная геодезическая задача всё ещё заключается в том, чтобы по координатам двух точек, определить длину и азимутов, как прямого $\alpha_1$ , так и обратного $\alpha_2$ .

Для сферы радиуса , решение следующее, опять же по теореме косинусов:

$\cos\sigma=\sin\varphi_1\sin\varphi_2+\cos\varphi_1\cos\varphi_2\cos\omega$

где $\omega=\lambda_2-\lambda_1$ , $\sigma$ — расстояние, выраженное в долях радиуса сферы. Само расстояние определяется как:

$\left\{ \begin{array}{1} &s=R\sigma&\text{, при }\cos\sigma\geq0 \\ &s=R\left(\pi-|\sigma|\right)&\text{, при }\cos\sigma<0 \\ \end{array} \right.$

Немного кода:

var ω = coord2.LonR - coord1.LonR;

var cosSigma = Math.Sin(coord1.LatR) * Math.Sin(coord2.LatR) + Math.Cos(coord1.LatR) * Math.Cos(coord2.LatR) * Math.Cos(ω);

double s;
if (cosSigma < 0)
    s = _ellipsoid.PolarRadius * (Math.PI - Math.Abs(Math.Acos(cosSigma)));
else if (Math.Abs(cosSigma - 1) < TOLERANCE)
    s = 0;
else
    s = _ellipsoid.PolarRadius * Math.Acos(cosSigma);

Прямой $\alpha_1$ и обратный $\alpha_2$ азимуты, можно выразить применяя формулы котангенсов:

$\tan\alpha_1=\frac{\cos\varphi_2\sin\omega}{\cos\varphi_1\sin\varphi_2-\sin\varphi_1\cos\varphi_2\cos\omega}$ $\tan\alpha_2=\frac{\cos\varphi_1\sin\omega}{\cos\varphi_1\sin\varphi_2\cos\omega-\sin\varphi_1\cos\varphi_2}$

Ещё немного кода:

var a1 =
    Math.Atan(Math.Cos(coord2.LatR) * Math.Sin(ω) /
              (Math.Cos(coord1.LatR) * Math.Sin(coord2.LatR) - Math.Sin(coord1.LatR) * Math.Cos(coord2.LatR) * Math.Cos(ω))) * 180 / Math.PI;
var a2 =
    Math.Atan(Math.Cos(coord1.LatR) * Math.Sin(ω) /
              (Math.Cos(coord1.LatR) * Math.Sin(coord2.LatR) * Math.Cos(ω) - Math.Sin(coord1.LatR) * Math.Cos(coord2.LatR))) * 180 / Math.PI;

Ну и метод полностью:

Скрытый текст

private InverseProblemAnswer OrthodromicSpheroidDistance(Point coord1, Point coord2)
{
    var ω = coord2.LonR - coord1.LonR;

    var cosSigma = Math.Sin(coord1.LatR) * Math.Sin(coord2.LatR) + Math.Cos(coord1.LatR) * Math.Cos(coord2.LatR) * Math.Cos(ω);

    double s;
    if (cosSigma < 0)
        s = _ellipsoid.PolarRadius * (Math.PI - Math.Abs(Math.Acos(cosSigma)));
    else if (Math.Abs(cosSigma - 1) < TOLERANCE)
        s = 0;
    else
        s = _ellipsoid.PolarRadius * Math.Acos(cosSigma);

    var a1 =
        Math.Atan(Math.Cos(coord2.LatR) * Math.Sin(ω) /
                  (Math.Cos(coord1.LatR) * Math.Sin(coord2.LatR) - Math.Sin(coord1.LatR) * Math.Cos(coord2.LatR) * Math.Cos(ω))) * 180 /
        Math.PI;
    var a2 =
        Math.Atan(Math.Cos(coord1.LatR) * Math.Sin(ω) /
                  (Math.Cos(coord1.LatR) * Math.Sin(coord2.LatR) * Math.Cos(ω) - Math.Sin(coord1.LatR) * Math.Cos(coord2.LatR))) * 180 /
        Math.PI;

    a1 = Azimuth.AzimuthRecovery(coord1, coord2, a1);
    a2 = Azimuth.AzimuthRecovery(coord2, coord1, a2);

    return new InverseProblemAnswer(a1, a2, s);
}

Метод на GitHub: InverseProblemService.cs

По обратной геодезической задаче на сфере — всё, нашли азимуты и расстояние. Теперь опишу их решение для эллипсоида.

Решение геодезических задач на эллипсоиде

Помните как мы искали широту промежуточной точки или координаты пересечения двух линий на эллипсоиде? Так у нас ещё были какие‑то итерационные алгоритмы. То есть ситуация следующая: есть эллипсоид — есть итерационный алгоритм, нет аналитического решения.

И здесь нам так же понадобится итерационный алгоритм. К счастью нам его не придётся выдумывать, его уже придумал Thaddeus Vincenty (Тадеуш Винсенти?)

На Вики, опять же, есть описание его решения: en:Vincenty's formulae. Пишут что имеет очень высокую точность, аж до 0,5 мм. Ну что же, начнём с обратной геодезической задачи.

Обратная геодезическая задача

Всё ещё пытаемся вычислить расстояние по известным координатам:

$\varphi_1$ и $\lambda_1$ — широта и долгота начальной точки отрезка

$\varphi_2$ и $\lambda_2$ — широта и долгота конечной точки отрезка

Начинается решение с пересчёта геодезических широт в приведённую широту, следующим образом:

$\begin{align} u_1=\arctan\left(\left(1-f\right)\tan\varphi_1\right)\\ u_2=\arctan\left(\left(1-f\right)\tan\varphi_2\right)\\ \end{align}$

Здесь — коэффициент полярного сжатия.

В коде мы сразу рассчитаем для них синусы и косинусы — пригодится.

double l = coord2.LonR - coord1.LonR; // Разность геодезических долгот

// Приведённая широта
double u1 = Math.Atan((1 - _ellipsoid.F) * Math.Tan(coord1.LatR));
double u2 = Math.Atan((1 - _ellipsoid.F) * Math.Tan(coord2.LatR));
double sinU1 = Math.Sin(u1), cosU1 = Math.Cos(u1);
double sinU2 = Math.Sin(u2), cosU2 = Math.Cos(u2);

Последовательным приближением рассчитывается разность долгот на эллипсоиде $\lambda$ и угловое расстояние $\sigma$ .

В первом приближении принимаем $\lambda=L$ , где $L=\lambda_2-\lambda_1$ и проводим итеративные вычисления, пока $\lambda$ не сойдётся.

$\begin{align} &\sin\sigma=\sqrt{\left(\cos{u_2}\sin\lambda\right)^2+\left(\cos{u_1}\sin{u_2}-\sin{u_1}cos{u_2}\cos\lambda\right)^2} \\ &\cos\sigma=\sin{u_1}\sin{u_2}+\cos{u_1}\cos{u_2}\cos\lambda \\ \end{align}$

Код:


var sinLambda = Math.Sin(lambda);
var cosLambda = Math.Cos(lambda);

sinSigma =
    Math.Sqrt(Math.Pow(cosU2 * sinLambda, 2) + Math.Pow(cosU1 * sinU2 - sinU1 * cosU2 * cosLambda, 2));

...

cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;

Из этих двух формул можем вывести $\sigma$ :

$\sigma=\arctan\left(\frac{\sin\sigma}{\cos\sigma}\right)$

$\sin\alpha=\frac{\cos{u_1}\cos{u_2}\sin\lambda}{\sin\sigma}$

В коде:

sigma = Math.Atan(sinSigma / cosSigma);
var sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;

По основному тригонометрическому тождеству $\cos^2\alpha=1-\sin^2\alpha$ :

$\cos{2\sigma_m}=\cos\sigma-\frac{2\sin{u_1}\sin{u_2}}{\cos^2\alpha}=\cos\sigma-\frac{2\sin{u_1}\sin{u_2}}{1-\sin^2\alpha}$ $C=\frac{f}{16}\cos^2\alpha\left(4+f\left(4-3\cos^2\alpha\right)\right)$ $\lambda=L+\left(1-C\right)f\sin\alpha\biggl(\sigma+C\sin\sigma\Bigl(\cos{2\sigma_m}+C\cos\sigma\bigl(-1+2\cos^2{2\sigma_m}\bigl)\Bigl)\biggl)$

Вычисления в программном коде:

cosSqAlpha = 1 - sinAlpha * sinAlpha;

if (Math.Abs(cosSqAlpha) > TOLERANCE)
    cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
else
    cos2SigmaM = 0;

double c = _ellipsoid.F / 16 * cosSqAlpha * (4 + _ellipsoid.F * (4 - 3 * cosSqAlpha));
lambdaP = lambda;
lambda = l +
         (1 - c) * _ellipsoid.F * sinAlpha *
         (sigma + c * sinSigma * (cos2SigmaM + c * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));

Приближение продолжается, пока изменения $\lambda$ в сравнении с предыдущей итераций не станут допустимо малы, например $10^{-12}$ , что соответствует примерно 0,06 мм.

do
{
  ...
} while (Math.Abs(lambda - lambdaP) > 1e-12 && --iterLimit > 0);

Теперь начинаем расчёт дистанции:

$u^2=\cos^2\alpha\left(\frac{a^2-b^2}{b^2}\right)$

где — экваториальный радиус, — полярный радиус.

Дальше идут магические числа и :

$\begin{align} &A=1+\frac{u^2}{16384}\biggl(4096+u^2\Bigl(-768+u^2\bigl(320-175u^2\bigl)\Bigl)\biggl) \\ &B=\frac{u^2}{1024}\biggl(256+u^2\Bigl(-128+u^2\bigl(74-47u^2\bigl)\Bigl)\biggl) \end{align}$

Код:

double uSq = cosSqAlpha * (Math.Pow(_ellipsoid.EquatorialRadius, 2) - Math.Pow(_ellipsoid.PolarRadius, 2)) /
             Math.Pow(_ellipsoid.PolarRadius, 2);
double a = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
double b = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));

Кстати, для вычисления и предложена модификация с более простыми формулами, давайте их тоже приведём, на всякий случай:

$\begin{align} &A=\frac{1+\frac{1}{4}{k_1}^2}{1-k_1} \\ &B=k_1\left(1-\frac{3}{8}{k_1}^2\right) \end{align}$

где вводится коэффициент:

$k_1=\frac{\sqrt{1+u^2}-1}{\sqrt{1+u^2}+1}$

В общем, выбирайте что больше нравится из этих двух вариантов, я выбрал первый.

Далее жуткая формула, надеюсь я ни кого не запутал своей расстановкой отступов:

$\begin{align} \Delta\sigma=B\sin\sigma\biggl(&\cos{2\sigma_m} \\ &+\frac{1}{4}B\Bigl(& \\ &&\cos\sigma\left(-1+2\cos^2{2\sigma_m}\right) \\ &&-\frac{B}{6}\left(\cos{2\sigma_m}\right)\left(-3+4\sin^2\sigma\right)\left(-3+4\cos^22\sigma_m\right) \\ &&\Bigl)\\ &\biggl) \end{align}$

Наконец, вычисляется расстояние между точками:

$s=bA\left(\sigma-\Delta\sigma\right)$

Снова немного кода:

double deltaSigma = b * sinSigma *
                    (cos2SigmaM +
                     b / 4 *
                     (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) -
                      b / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) *
                      (-3 + 4 * cos2SigmaM * cos2SigmaM)));
double s = _ellipsoid.PolarRadius * a * (sigma - deltaSigma);

Ну и какая обратная геодезическая задача без вычисления азимутов:

$\begin{align} &\tan\alpha_1=\frac{\cos{u_2}\sin\lambda}{\cos{u_1}\sin{u_2}-\sin{u_1}\cos{u_2}\cos\lambda} \\ &\tan\alpha_2=\frac{\cos{u_1}\sin\lambda}{\cos{u_1}\sin{u_2}\cos\lambda-\sin{u_1}\cos{u_2}} \end{align}$

Код для азимутов:

var a1 = Math.Atan(cosU2 * Math.Sin(lambda) / (cosU1 * sinU2 - sinU1 * cosU2 * Math.Cos(lambda))) * 180 /
         Math.PI;
var a2 = Math.Atan(cosU1 * Math.Sin(lambda) / (-sinU1 * cosU2 + cosU1 * sinU2 * Math.Cos(lambda))) * 180 /
         Math.PI;
a1 = Azimuth.AzimuthRecovery(coord1, coord2, a1);
a2 = Azimuth.AzimuthRecovery(coord2, coord1, a2);

Весь код вместе:

Скрытый текст

private InverseProblemAnswer OrthodromicEllipsoidDistance(Point coord1, Point coord2)
{
    double l = coord2.LonR - coord1.LonR; // Разность геодезических долгот

    double u1 = Math.Atan((1 - _ellipsoid.F) * Math.Tan(coord1.LatR)); // Приведённая широта
    double u2 = Math.Atan((1 - _ellipsoid.F) * Math.Tan(coord2.LatR));
    double sinU1 = Math.Sin(u1), cosU1 = Math.Cos(u1);
    double sinU2 = Math.Sin(u2), cosU2 = Math.Cos(u2);

    double lambda = l, lambdaP;
    double cosSqAlpha, sinSigma, cosSigma, cos2SigmaM, sigma;
    int iterLimit = 100;

    do
    {
        var sinLambda = Math.Sin(lambda);
        var cosLambda = Math.Cos(lambda);

        sinSigma =
            Math.Sqrt(Math.Pow(cosU2 * sinLambda, 2) + Math.Pow(cosU1 * sinU2 - sinU1 * cosU2 * cosLambda, 2));

        if (Math.Abs(sinSigma) < TOLERANCE)
            return new InverseProblemAnswer(0, 0, 0);

        cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
        sigma = Math.Atan(sinSigma / cosSigma);
        var sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
        cosSqAlpha = 1 - sinAlpha * sinAlpha;

        if (Math.Abs(cosSqAlpha) > TOLERANCE)
            cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;
        else
            cos2SigmaM = 0;

        double c = _ellipsoid.F / 16 * cosSqAlpha * (4 + _ellipsoid.F * (4 - 3 * cosSqAlpha));
        lambdaP = lambda;
        lambda = l +
                 (1 - c) * _ellipsoid.F * sinAlpha *
                 (sigma + c * sinSigma * (cos2SigmaM + c * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));
    } while (Math.Abs(lambda - lambdaP) > 1e-12 && --iterLimit > 0);

    double uSq = cosSqAlpha * (Math.Pow(_ellipsoid.EquatorialRadius, 2) - Math.Pow(_ellipsoid.PolarRadius, 2)) /
                 Math.Pow(_ellipsoid.PolarRadius, 2);
    double a = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
    double b = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));
    double deltaSigma = b * sinSigma *
                        (cos2SigmaM +
                         b / 4 *
                         (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) -
                          b / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) *
                          (-3 + 4 * cos2SigmaM * cos2SigmaM)));
    double s = _ellipsoid.PolarRadius * a * (sigma - deltaSigma);

    var a1 = Math.Atan(cosU2 * Math.Sin(lambda) / (cosU1 * sinU2 - sinU1 * cosU2 * Math.Cos(lambda))) * 180 /
             Math.PI;
    var a2 = Math.Atan(cosU1 * Math.Sin(lambda) / (-sinU1 * cosU2 + cosU1 * sinU2 * Math.Cos(lambda))) * 180 /
             Math.PI;
    a1 = Azimuth.AzimuthRecovery(coord1, coord2, a1);
    a2 = Azimuth.AzimuthRecovery(coord2, coord1, a2);

    return new InverseProblemAnswer(a1, a2, s);
}

Код метода на GitHub: InverseProblemService.cs

Ну, вроде всё, идём дальше.

Прямая геодезическая задача

Напоминаю, что задача состоит в том, чтобы по известной начальной координате $\{\varphi_1;\lambda_1\}$ , расстоянию , и азимуту $\alpha_1$ , найти вторую координату $\{\varphi_2;\lambda_2\}$ .

Начнём с пересчёта геодезической широты в приведённую широту:

$u_1=\arctan\left(\left(1-f\right)\tan\varphi_1\right)$

где, - всё ещё коэффициент сжатия.

Угловое расстояние между точкой и экватором:

$\sigma_1=\arctan\left(\frac{\tan{u_1}}{\cos{\alpha_1}}\right)$

Прямой азимут $\alpha$ геодезической линии на экваторе, если бы она была продолжена до него:

$\sin\alpha=\cos{u_1}\sin\alpha_1$

Код:

double u1 = Math.Atan((1 - _ellipsoid.F) * Math.Tan(coord.LatR));
double sigma1 = Math.Atan(Math.Tan(u1) / Math.Cos(a1));
var sinAlpha = Math.Cos(u1) * Math.Sin(a1);

Вычисляем вспомогательные величины:

$u^2=\cos^2\alpha\left(\frac{a^2-b^2}{b^2}\right)=\left(1-\sin^2\alpha\right)\left(\frac{a^2-b^2}{b^2}\right)$

, здесь - экваториальный радиус, - полярный радиус: $b=\left(1-f\right)a$

var cosSqAlpha = 1 - sinAlpha * sinAlpha;

double uSq = cosSqAlpha * (Math.Pow(_ellipsoid.EquatorialRadius, 2) - Math.Pow(_ellipsoid.PolarRadius, 2)) /
             Math.Pow(_ellipsoid.PolarRadius, 2);

И вновь магические числа, о которых мы уже писали, повторим ещё раз:

$\begin{align} &A=1+\frac{u^2}{16384}\biggl(4096+u^2\Bigl(-768+u^2\bigl(320-175u^2\bigl)\Bigl)\biggl) \\ &B=\frac{u^2}{1024}\biggl(256+u^2\Bigl(-128+u^2\bigl(74-47u^2\bigl)\Bigl)\biggl) \end{align}$

Зададим начальное значение $\sigma$ :

$\sigma=\frac{s}{bA}$

Код:

var si = s / _ellipsoid.PolarRadius / a;
double sigma = si;

Будем повторять следующие вычисления, пока изменения $\sigma$ (угловое расстояние между точками), не станут незначительными:

$2\sigma_m=2\sigma_1+\sigma$

$\sigma_m$ – угловое расстояние от экватора до середины линии.

Опять страшная формула:

Вся эта жуть в коде:

var sinSigma = Math.Sin(sigma);
var cosSigma = Math.Cos(sigma);

cos2SigmaM = Math.Cos(2 * sigma1 + sigma);
double deltaSigma = b * sinSigma *
                    (cos2SigmaM +
                     b / 4 *
                     (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) -
                      b / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) *
                      (-3 + 4 * cos2SigmaM * cos2SigmaM)));
sigmaP = sigma;
sigma = si + deltaSigma;

После того как $\sigma$ получена с достаточной точностью, вычислить широту второй точки:

$\varphi_2=\arctan\left(\frac{\sin{u_1}\cos\sigma+\cos{u_1}\sin\sigma\cos\alpha_1}{(1-f)\sqrt{\sin^2\alpha+(\sin{u_1}\sin\sigma-\cos{u_1}\cos\sigma\cos\alpha_1)^2}}\right)$

Можно упростить эту формулу, введя :

$\sin{u_2}=\sin{u_1}\cos\sigma+\cos{u_1}\sin\sigma\cos\alpha_1$

Тогда формула будет выглядеть:

$\varphi_2=\arctan\left(\frac{\sin{u_2}}{(1-f)\sqrt{1-\sin^2u_2}}\right)$

Код:

var sinU2 = Math.Sin(u1) * Math.Cos(sigma) + Math.Cos(u1) * Math.Sin(sigma) * Math.Cos(a1);
var lat2 = sinU2 /
           ((1 - _ellipsoid.F) *
            Math.Sqrt(sinAlpha * sinAlpha +
                      Math.Pow(
                          Math.Sin(u1) * Math.Sin(sigma) - Math.Cos(u1) * Math.Cos(sigma) * Math.Cos(a1), 2)));
lat2 = Math.Atan(lat2);

Теперь вычисляем долготу. Для этого рассчитывается разность долгот:

$\lambda=\arctan\left(\frac{\sin\sigma\sin\alpha_1}{\cos{u_1}\cos\sigma-\sin{u_1}\sin\sigma\cos\alpha_1}\right)$ $С=\frac{f}{16}\cos^2\alpha\left(4+f(4-3\cos^22\alpha)\right)$ $L=\lambda-(1-C)f\sin\alpha\biggl(\sigma+C\sin\sigma\Bigl(\cos{2\sigma_m}+C\cos\sigma\bigl(-1+2\cos^2{2\sigma_m}\bigl)\Bigl)\biggl)$ $\lambda_2=L+\lambda_1$

Код:

// Разность долгот
var lamda = Math.Atan(Math.Sin(sigma) * Math.Sin(a1) /
                      (Math.Cos(u1) * Math.Cos(sigma) - Math.Sin(u1) * Math.Sin(sigma) * Math.Cos(a1)));

double c = _ellipsoid.F / 16 * cosSqAlpha * (4 + _ellipsoid.F * (4 - 3 * cosSqAlpha));

var l = lamda - (1 - c) * _ellipsoid.F * sinAlpha *
        (sigma +
         c * Math.Sin(sigma) * (cos2SigmaM + c * Math.Cos(sigma) * (-1 + 2 * cos2SigmaM * cos2SigmaM)));
var lon2 = -l + coord.LonR;

Ну и обратный азимут:

$\alpha_2=\arctan\left(\frac{\sin\alpha}{-\sin{u_1}\sin\sigma+\cos{u_1}\cos\sigma\cos\alpha_1}\right)$

В программе:

var a2 = -Math.Atan(sinAlpha / (-Math.Sin(u1) * Math.Sin(sigma) + Math.Cos(u1) * Math.Cos(sigma) * Math.Cos(a1))) * 180 / Math.PI;
a2 = Azimuth.AzimuthRecovery(point1, coord, a2);

Под обратным азимутом, можно понимать, как азимут в направлении от первой точки, ко второй, так и наоборот, они отличаются на 180°, формула тангенса не фиксирует эти различия.

Снова код метода целиком:

Скрытый текст

private DirectProblemAnswer DirectProblemEllipsoid(Point coord, double a1, double s)
{
    var aDegree = a1;
    a1 = a1 * Math.PI / 180;

    double u1 = Math.Atan((1 - _ellipsoid.F) * Math.Tan(coord.LatR));
    double sigma1 = Math.Atan(Math.Tan(u1) / Math.Cos(a1));
    var sinAlpha = Math.Cos(u1) * Math.Sin(a1);

    var cosSqAlpha = 1 - sinAlpha * sinAlpha;

    double uSq = cosSqAlpha * (Math.Pow(_ellipsoid.EquatorialRadius, 2) - Math.Pow(_ellipsoid.PolarRadius, 2)) /
                 Math.Pow(_ellipsoid.PolarRadius, 2);
    double a = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)));
    double b = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)));

    var si = s / _ellipsoid.PolarRadius / a;
    double sigma = si;

    double sigmaP, cos2SigmaM;
    do
    {
        var sinSigma = Math.Sin(sigma);
        var cosSigma = Math.Cos(sigma);

        cos2SigmaM = Math.Cos(2 * sigma1 + sigma);
        double deltaSigma = b * sinSigma *
                            (cos2SigmaM +
                             b / 4 *
                             (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) -
                              b / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) *
                              (-3 + 4 * cos2SigmaM * cos2SigmaM)));
        sigmaP = sigma;
        sigma = si + deltaSigma;
    } while (Math.Abs(sigma - sigmaP) > 1e-12);

    var sinU2 = Math.Sin(u1) * Math.Cos(sigma) + Math.Cos(u1) * Math.Sin(sigma) * Math.Cos(a1);
    var lat2 = sinU2 /
               ((1 - _ellipsoid.F) *
                Math.Sqrt(sinAlpha * sinAlpha +
                          Math.Pow(
                              Math.Sin(u1) * Math.Sin(sigma) - Math.Cos(u1) * Math.Cos(sigma) * Math.Cos(a1), 2)));
    lat2 = Math.Atan(lat2);

    // Разность долгот
    var lamda = Math.Atan(Math.Sin(sigma) * Math.Sin(a1) /
                          (Math.Cos(u1) * Math.Cos(sigma) - Math.Sin(u1) * Math.Sin(sigma) * Math.Cos(a1)));

    double c = _ellipsoid.F / 16 * cosSqAlpha * (4 + _ellipsoid.F * (4 - 3 * cosSqAlpha));

    var l = lamda - (1 - c) * _ellipsoid.F * sinAlpha *
            (sigma +
             c * Math.Sin(sigma) * (cos2SigmaM + c * Math.Cos(sigma) * (-1 + 2 * cos2SigmaM * cos2SigmaM)));
    var lon2 = -l + coord.LonR;

    lon2 = lon2 * 180 / Math.PI;
    lat2 = lat2 * 180 / Math.PI;

    var point1 = IsPolisIntersect(coord, aDegree, s)
        ? new Point(lon2 + (lon2 < 0 ? 180 : -180), lat2)
        : new Point(lon2, lat2);

    var a2 = -Math.Atan(sinAlpha / (-Math.Sin(u1) * Math.Sin(sigma) + Math.Cos(u1) * Math.Cos(sigma) * Math.Cos(a1))) * 180 / Math.PI;
    a2 = Azimuth.AzimuthRecovery(point1, coord, a2);

    return new DirectProblemAnswer(point1, a2);
}

Код на GitHub: DirectProblemService.cs

Ну вот и всё, наконец-то, переходим к итогам.

Итоги

Голова идёт кругом от всей этой сферической геометрии. Давайте ещё раз сделаем краткий обзор того, что мы преодолели.

Для начала мы привели формулу вычисления широты по долготе, для сфероида. Потом мы вывели из неё формулу вычисления долготы по широте, опять же для сфероида, так, по фану.

Потом мы уже перешли непосредственно к задаче поиска точки пересечения, и воспользовавшись, всё теми же формулами вычисления широты по долготе, и долготы по широте, научились находить точку пересечения ортодромий на сфероиде.

Но нам то нужен эллипсоид, ибо Земля сплюснута. И мы проделали всё то же самое для эллипсоида: научились вычислять широту по долготе и долготу по широте. И снова вернулись к изначальной задаче - поиска точки пересечения двух линий, - определили алгоритм для эллипсоида.

Работа на эллипсоиде - это именно та часть, которую я и ставлю себе в заслугу. Magnum opus, так сказать. Остальное можно нагуглить, хоть и не без проблем, всё же вещь специфическая, либо, если не гуглится, можно аналитически вывести используя известные тригонометрические формулы.

Так или иначе при работе на эллипсоиде активно применяются главные геодезические задачи, их тоже нельзя было обойти стороной, поэтому описал и их, как для сфероида: прямая и обратная, так и для эллипсоида: обратная и прямая.

Решение для эллипсоида было предложено Тадеушем Винсенти и характеризуется высокой точностью. Однако высокая точность на эллипсоиде достигается за счёт итеративного приближения к ответу, таким образом, по сравнению с аналитическим решением, страдает производительность. Поэтому Вам остаётся только ответить на вопрос: а стоит ли игра свеч?

Комментарии (25)

aamonster
10.09.2024 12:49
#27272830
+6
Какая боль...

А ведь можно было свести задачу к пересечению плоскостей – и всё стало бы просто, линейно (пардон, почти: из-за сферы квадратное уравнение вылезает) и без всей этой жуткой тригонометриии.

Набросок:
1. Вычисляем декартовы координаты x, y, z концов "отрезка" (можно принять радиус сферы за 1, можно R – неважно)
2. "Отрезок" заменяем на плоскость, проходящую через его концы и центр сферы – она задаётся линейным уравнением вида a*x+b*y+c*z = 0
3. Два таких уравнения (для двух отрезков) задают прямую. Дополняем их уравнением сферы x²+y²+z²=R² и находим координаты двух точек пересечения.
4. Проверяем, попадает ли одна из двух точек пересечения на оба отрезка (она в плоскости отрезка, так что всё сведётся к проверке знаков векторных произведений, если надо – могу расписать подробнее)
5. Переводим декартовы координаты в эйлеровы или какие нам нужны.
1. anonymous
  10.09.2024 12:49
  #27272872
  НЛО прилетело и опубликовало эту надпись здесь
1. velon Автор
  10.09.2024 12:49
  #27272876
  Как-то у Вас всё однозначно. Вопрос же всегда в соотношении точности и скорости. Алгоритмы надо по этому принципу сравнивать, а не "просто"/"сложно".
  
  Любые преобразования крадут точность, да и на скорости положительно не скажутся.
  
  Если исходить из этих характеристик, то точность того что я описал (поска точки пересечения) будет близка к точности алгоритма Винсети, а скорость... ну ассимптотическая сложность логарифмическая
  1. aamonster
    10.09.2024 12:49
    #27273082
    +2
    Да, у меня всё однозначно: тригонометрические уравнения – всегда боль, если можно их избежать, всего лишь увеличив число переменных на 1 или 2 – следует это сделать. Посмотрите, к примеру, на 3d-графику: все представляют поворот матрицами или кватернионами, никто не использует углы в вычислениях.
    
    Правда, это всё про решение для шара. Для эллипсоида возникнут сложности – на нём, кажется, ортодромия может не лежать в одной плоскости. И, кстати, возможно, это сломает и ваш алгоритм (требует проверки) – не факт, что можно использовать прямую геодезическую задачу для вычисления координат промежуточной точки: не факт, что угол для промежуточной точки будет таким же, как угол для конечной.
    
    velon Автор
    10.09.2024 12:49
    #27273136
    Вы какие-то странные вещи пишете про ортодромию и геодезические задачи.
    
    Мне кажется факт что ортодромии в одной плоскости НЕ лежат. Из этого и исходим... но ничего не сломалось
    
    Прямую геодезическую задачу очень даже можно использовать для вычисления координат это факт - в этом смысл задачи.
    
    То что угол для промежуточной и конечной НЕ совпадёт, так это тоже факт, но что это меняет?!
    
    aamonster
    10.09.2024 12:49
    #27273172
    
    Используя обратную геодезическую задачу, определим расстояние между двумя известными точками и угол.
    
    Полученное расстояние разделим пополам, и, применив прямую геодезическую задачу, по найденному ранее углу, и полученному расстоянию, найдём координаты средней точки.
    
    Так вот, если угол для промежуточной точки не тот же, что для конечной – то нельзя в пункте 2 использовать угол из пункта 1. В этом сложность.
    (надо проверять, может, для эллипсоида угол и тот же – я не знаю).
    
    velon Автор
    10.09.2024 12:49
    #27273182
    Так здесь в обоих случаях же считаем от самой первой точки, только один раз мы берём полную дистанцию (т.е. вычисляем её), а во второй раз 1/2 прикладываем.
    
    Если формулировка вас запутала, предложите вариант - я исправлю
    
    П.С. Если на следующей итерации мы возьмём среднюю точку и конечную - то углы считаем заново
    
    aamonster
    10.09.2024 12:49
    #27273264
    Да нет, формулировка как раз не запутала. Но, кажется, тут действительно нет проблемы, "азимут" будет одним и тем же. Теперь надо проверять алгоритм решения прямой геодезической задачи :-)
    Ещё хотел написать: для эллипсоида можно использовать ваш же подход (проверив его) с алгоритмом половинного деления (или лучше заменить его на метод Ньютона), но всё, что только возможно, вычислять в декартовых координатах.
    
    Upd: И да, любопытно, какова будет погрешность вычисления координат точки пересечения, если её вычислять, пренебрегая "сплющенностью" Земли. Возможно, если точности не хватит – выгодней будет вначале найти точку пересечения "на сфере", а потом скорректировать.
    
    velon Автор
    10.09.2024 12:49
    #27273292
    Но ведь и так применяется метод половинного деления: он же, насколько я помню, заключается в том, чтобы искать тот подинтревал где функция меняет знак.
    
    Теперь попарно сравним $\Delta\varphi_f$ с $\Delta\varphi_m$ и $\Delta\varphi_m$ с $\Delta\varphi_s$
    
    ЕСЛИ у $\Delta\varphi_f$ и $\Delta\varphi_m$ разные знаки, ТО $\lambda_s=\lambda_m$
    
    ЕСЛИ у $\Delta\varphi_m$ и $\Delta\varphi_s$ разные знаки, ТО $\lambda_f=\lambda_m$
    
    Различие знаков говорит о том, что пересечение произошло на этом участке, между соответствующими значениями долготы. Алгоритм рекурсивно продолжается, переходим на ШАГ 2.
  1. mentin
    10.09.2024 12:49
    #27273986
    Вспомнился мультик: "Лучше день потерять… Запомни: лучше день потерять, потом за пять минут долететь."
    
    Так и тут, лучше преобразования сначала и в конце сделать, зато все основные вычисления без тригонометрии в векторах делать, точнее будет.
    
    velon Автор
    10.09.2024 12:49
    #27274006
    Что вы несёте? На чём вообще основано это утверждение что тригонометрия снижает как-то точность?
    
    mentin
    10.09.2024 12:49
    #27274018
    На опыте и оценках погрешностей вычислений.
    
    velon Автор
    10.09.2024 12:49
    #27274086
    Тогда от Филдсовской премии вас отделяет один шаг. Доказательство того, что тригонометрические функции врут - перевернёт мир вычислений
    
    mentin
    10.09.2024 12:49
    #27274290
    Врут не тригонометрические функции, а численные вычисления с ограниченной точностью.
    
    Как при вычислении корней квадратного уравнения, есть абсолютно правильная школьная формула. Но знакомые с численными вычислениями знают что для точных вычислений надо, вдобавок, использовать симметричную с дискриминантом в знаменателе.

AlexTmp8
10.09.2024 12:49
#27273180
А какая область применения?

Для какой задачи важна эллиптичность Земли? И какая ошибка будет если ей пренебречь?
1. velon Автор
  10.09.2024 12:49
  #27273268
  Мы использовали для триангуляции, применяли жадный алгоритм триангуляции, а он основан на поиске пересечений.
  
  Если грубо, то я хотел чтобы алгоритм разбивал на треугольники как в первом примере (слева), а он разбивал как во втором (справа).
  
  Это ещё красивый пример, на самом деле всё было гораздо хуже, получались треугольники с углами 170° , 7°, 3°, и высота у них была несколько метров при длине основания в несколько километров.
  
  Тогда решили всеми правдами и неправдами повысить точность.
  
  П.С. от проблемы полностью не избавились, но триангуляция стала строиться иначе
  
  П.П.С. мы в итоге остановились на упрощении до сферы, но наработки остались
  1. mentin
    10.09.2024 12:49
    #27273506
    Точность от всей этой тригонометрии сильно страдает. В статье есть оптимистическая фраза
    
    можно использовать координаты любой из двух пересекающихся линий - значение совпадёт
    
    Они не совпадут конечно, погрешность может быть существенная, особенно для отрезков около полюсов, почти-параллельных отрезков и так далее.
    
    Для вычисления на сфере есть гугловская библиотека S2 geometry, рекомендую глянуть как они это считают. Там оценивается точность, выбирается порядок вычислений дающий меньшую погрешность, и иногда все равно вынуждено использующую exact float для получения достаточно точного результата. Правда, она только для правильной сферы.
    
    velon Автор
    10.09.2024 12:49
    #27273688
    Совпадут на сфере, решение же аналитическое, если не совпали - значит ошибка в расчётах. С чего вообще в тригонометрии точность страдает?
    
    Вы конечно можете сказать что если 10 разделить на 3, а потом умножить на 3, то 10 мы уже не получим, - да, из-за подобных выкрутасов с плавающей точкой не совпадут, но я бы не называл эту погрешность существенной.
    
    П.С. я не понял что значит "правильная сфера".
    
    mentin
    10.09.2024 12:49
    #27274066
    да, из-за подобных выкрутасов с плавающей точкой не совпадут, но я бы не называл эту погрешность существенной.
    
    У меня требования когда эта погрешность бывает очень существенной. Например, из-за такой погрешности в вычислениях можно неправильно определить с какой стороны от отрезка лежит точка.
    
    velon Автор
    10.09.2024 12:49
    #27274096
    Тогда у меня для Вас плохие новости, вы только посмотрите сколько тригонометрии в преобразовании координат: ГОСТ 32453—2017
    
    mentin
    10.09.2024 12:49
    #27274962
    Можно очень аккуратно использовать тригонометрию, рассчитывая и указывая как в ГОСТ получаемую точность вычислений. Главное не считать это "аналитической" формулой, и понимать что не любая правильная с аналитической точки зрения формула даст такую же точность в любой точке глобуса.
    
    "Аналитическая" формула для пересечения на сфере у вас падает с делением на ноль если один из отрезков идёт строго вдоль меридиана, или если скажем результат попал на 90 градусов долготы. Уже это должно наводить на мысли что и рядом с этими точками точность вычислений сильно страдает.
    
    Вектора обычно лучше хотя бы тем, что избегают таких выделенных значений.

voldemar_d
10.09.2024 12:49
#27273500
Но нам то нужен эллипсоид, ибо Земля сплюснута.

Если точнее, то Земля - геоид, а не эллипсоид.
1. velon Автор
  10.09.2024 12:49
  #27273544
  Земля то конечно геоид, но математической формулой геоид не описать, поэтому эллипсоид принят как наиболее близкий к поверхности геоида
1. SadOcean
  10.09.2024 12:49
  #27275060
  Ну геоид - "подобный Земле", это по сути ни о чем не говорит.
  Если хочется просто и элегантно, нужно сводить к эллипсоиду или уж тогда просто фигячить модель на реальных данных с заданной точностью.

apro
10.09.2024 12:49
#27274578
И да, я его придумал потому что не смог нагуглить, может он где‑то и без меня описан был, может за эти 7 лет кто‑то написал что‑то похожее, а может я придумал какую‑то фигню, которую умные люди изобретать не станут...

Вот есть статья 2017 года про алгоритм для поиска пересечения двух ортодромий: https://www.researchgate.net/publication/321358300_Intersection_and_point-to-line_solutions_for_geodesics_on_the_ellipsoid

Изобретаю свой сложный способ поиска координат точки пересечения двух линий +3

Расчёты

Вывод формулы долготы промежуточной точки как функции широты

Поиск координат точки пересечения линий на сфероиде

Определение широты или долготы промежуточной точки на эллипсоиде

Поиск координат точки пересечения линий на эллипсоиде

Решение геодезических задач на сфероиде

Прямая геодезическая задача

Обратная геодезическая задача

Решение геодезических задач на эллипсоиде

Обратная геодезическая задача

Прямая геодезическая задача

Итоги

Комментарии (25)

velon Автор

velon Автор

velon Автор

velon Автор

velon Автор

velon Автор

velon Автор

velon Автор

velon Автор

velon Автор