Задачи на сфере: угловая засечка: различия между версиями
ErnieBoyd (обсуждение | вклад) м (→Алгоритм) |
ErnieBoyd (обсуждение | вклад) |
||
Строка 55: | Строка 55: | ||
== Пример программной реализации == | == Пример программной реализации == | ||
<syntaxhighlight lang="c"> | |||
/* | |||
* Решение угловой засечки | |||
* | |||
* Аргументы исходные: | |||
* pt1 - {широта, долгота} пункта Q1 | |||
* pt2 - {широта, долгота} пункта Q2 | |||
* | |||
* Аргументы определяемые: | |||
* pt3 - {широта, долгота} точки Q3 | |||
*/ | |||
int SphereAngular(double pt1[], double pt2[], double azi13, double azi23, | |||
double pt3[]) | |||
{ | |||
double azi12, dist12, azi21, dist13; | |||
double cos_beta1, sin_beta1, cos_beta2, sin_beta2, cos_dist12, sin_dist12; | |||
SphereInverse(pt2, pt1, &azi21, &dist12); | |||
SphereInverse(pt1, pt2, &azi12, &dist12); | |||
cos_beta1 = cos(azi13 - azi12); | |||
sin_beta1 = sin(azi13 - azi12); | |||
cos_beta2 = cos(azi21 - azi23); | |||
sin_beta2 = sin(azi21 - azi23); | |||
cos_dist12 = cos(dist12); | |||
sin_dist12 = sin(dist12); | |||
if (sin_beta1 == 0. && sin_beta2 == 0.) // Решение - любая точка | |||
return -1; // на большом круге Q1-Q2. | |||
else if (sin_beta1 == 0.) { | |||
pt3[0] = pt2[0]; // Решение - точка Q2. | |||
pt3[1] = pt2[1]; | |||
return 0; | |||
} else if (sin_beta2 == 0.) { // Решение - точка Q1. | |||
pt3[0] = pt1[0]; | |||
pt3[1] = pt1[1]; | |||
return 0; | |||
} else if (sin_beta1 * sin_beta2 < 0.) { // Лучи Q1-Q3 и Q2-Q3 направлены | |||
if (fabs(sin_beta1) >= fabs(sin_beta2)) { // в разные полусферы. | |||
cos_beta2 = -cos_beta2; // Выберем ближайшее решение: | |||
sin_beta2 = -sin_beta2; // развернём луч Q2-Q3 на 180°; | |||
} else { // иначе | |||
cos_beta1 = -cos_beta1; // развернём луч Q1-Q3 на 180°. | |||
sin_beta1 = -sin_beta1; | |||
} | |||
} | |||
dist13 = atan2(fabs(sin_beta2) * sin_dist12, | |||
cos_beta2 * fabs(sin_beta1) | |||
+ fabs(sin_beta2) * cos_beta1 * cos_dist12); | |||
SphereDirect(pt1, azi13, dist13, pt3); | |||
return 0; | |||
} | |||
</syntaxhighlight> | |||
== Ссылки == | == Ссылки == | ||
* [[Вычисление расстояния и начального азимута между двумя точками на сфере]] | * [[Вычисление расстояния и начального азимута между двумя точками на сфере]] | ||
* [http://gis-lab.info/qa/biangulation.html Нахождение точки пересечения двух линий по углам и двум известным точкам (биангуляция)] | * [http://gis-lab.info/qa/biangulation.html Нахождение точки пересечения двух линий по углам и двум известным точкам (биангуляция)] |
Версия от 12:17, 12 марта 2014
Линейная засечка — это нахождение положения точки по координатам двух исходных пунктов и значениям азимутов направлений с этих пунктов на определяемую точку.
Общие положения
В качестве модели Земли принимается сфера с радиусом R, равным среднему радиусу земного эллипсоида. Аналогом прямой линии на плоскости является геодезическая линия на поверхности. На сфере геодезическая линия — дуга большого круга.
Введём следующие обозначения:
- φ — географическая широта,
- λ — географическая долгота,
- α — азимут дуги большого круга,
- σ — сферическое расстояние (длина дуги большого круга, выраженная в долях радиуса шара).
Линейное расстояние по дуге большого круга s связано со сферическим расстоянием σ формулой s = R σ.
Постановка задачи
- Исходные данные
- координаты пунктов Q₁, Q₂ — φ₁, λ₁, φ₂, λ₂,
- начальные направления с пунктов Q₁, Q₂ на точку Q₃ — α₁₃, α₂₃.
- Определяемые величины
- координаты точки Q₃ — φ₃, λ₃.
Алгоритм
Решение любого вида засечек сводится к нахождению полярных координат искомой точки, т.е. начального направления и расстояния на неё с одного из исходных пунктов. На конечном этапе координаты находятся из решения прямой геодезической задачи. Поскольку в угловой засечке направления α₁₃ и α₂₃ уже заданы, остаётся определить расстояние σ₁₃ или σ₂₃.
На рисунке синим цветом выделены заданные элементы сферического треугольника, красным цветом неизвестные, зелёным — вспомогательные элементы. Очевидно, в треугольнике Q₁Q₂Q₃ нет ни одного известного элемента. Однако из решения обратной геодезической задачи для пунктов Q₁, Q₂ могут быть получены расстояние σ₁₂, а также азимуты α₁₂ и α₂₁, после чего углы β₁ и β₂ вычисляются как разности азимутов при соответствующих пунктах. Далее из решения треугольника Q₁Q₂Q₃ найдём сторону σ₁₃.
Последовательность действий:
- решить обратную геодезическую задачу для Q₁, Q₂: по φ₁, λ₁, φ₂, λ₂ получить α₁₂, α₂₁, σ₁₂;
- вычислить углы β₁, β₂;
- в треугольнике Q₁Q₂Q₃ по σ₁₂, β₁, β₂ вычислить σ₁₃;
- решить прямую геодезическую задачу для Q₁, Q₃: по φ₁, λ₁, α₁₃, σ₁₃ вычислить φ₃, λ₃.
Действия по первому и последнему пунктам рассмотрены в статьях Задачи на сфере: обратная геодезическая задача и Задачи на сфере: прямая геодезическая задача.
Углы β₁, β₂ и длина σ₁₃ вычисляются по формулам:
Правда, до вычисления длины σ₁₃ необходимо проанализировать полученные значения углов β₁ и β₂. Ниже в коде функции можно увидеть пример такого анализа:
- если линии Q₁Q₃ и Q₂Q₃ совпадают с Q₁Q₂, решение не определено, т.к. решением может быть любая точка геодезической линии Q₁Q₂;
- если одна из линий Q₁Q₃ и Q₂Q₃ совпадает с Q₁Q₂, а другая нет, решением является пункт, из которого выходит другая;
- если линии Q₁Q₃ и Q₂Q₃ уходят в разные полушария от Q₁Q₂, функция находит ближайшее к Q₁Q₂ «ложное пересечение» этих линий.
Здесь необходимо пояснить, что на сфере две несовпадающие геодезические линии всегда пересекаются в двух точках-антиподах. В традиционной постановке задачи направление на нужное пересечение задаётся явно. Если же прямое и обратное направления по условию равнозначны, возникает вопрос выбора одного из антиподов: φ₃, λ₃ или φ₃′ = −φ₃, λ₃′ = λ₃ ± 180°.
Пример программной реализации
/*
* Решение угловой засечки
*
* Аргументы исходные:
* pt1 - {широта, долгота} пункта Q1
* pt2 - {широта, долгота} пункта Q2
*
* Аргументы определяемые:
* pt3 - {широта, долгота} точки Q3
*/
int SphereAngular(double pt1[], double pt2[], double azi13, double azi23,
double pt3[])
{
double azi12, dist12, azi21, dist13;
double cos_beta1, sin_beta1, cos_beta2, sin_beta2, cos_dist12, sin_dist12;
SphereInverse(pt2, pt1, &azi21, &dist12);
SphereInverse(pt1, pt2, &azi12, &dist12);
cos_beta1 = cos(azi13 - azi12);
sin_beta1 = sin(azi13 - azi12);
cos_beta2 = cos(azi21 - azi23);
sin_beta2 = sin(azi21 - azi23);
cos_dist12 = cos(dist12);
sin_dist12 = sin(dist12);
if (sin_beta1 == 0. && sin_beta2 == 0.) // Решение - любая точка
return -1; // на большом круге Q1-Q2.
else if (sin_beta1 == 0.) {
pt3[0] = pt2[0]; // Решение - точка Q2.
pt3[1] = pt2[1];
return 0;
} else if (sin_beta2 == 0.) { // Решение - точка Q1.
pt3[0] = pt1[0];
pt3[1] = pt1[1];
return 0;
} else if (sin_beta1 * sin_beta2 < 0.) { // Лучи Q1-Q3 и Q2-Q3 направлены
if (fabs(sin_beta1) >= fabs(sin_beta2)) { // в разные полусферы.
cos_beta2 = -cos_beta2; // Выберем ближайшее решение:
sin_beta2 = -sin_beta2; // развернём луч Q2-Q3 на 180°;
} else { // иначе
cos_beta1 = -cos_beta1; // развернём луч Q1-Q3 на 180°.
sin_beta1 = -sin_beta1;
}
}
dist13 = atan2(fabs(sin_beta2) * sin_dist12,
cos_beta2 * fabs(sin_beta1)
+ fabs(sin_beta2) * cos_beta1 * cos_dist12);
SphereDirect(pt1, azi13, dist13, pt3);
return 0;
}