Задачи на сфере: линейная засечка: различия между версиями

Материал из GIS-Lab
Перейти к навигации Перейти к поиску
Строка 48: Строка 48:
Если величина косинуса превышает единицу, задача поставлена некорректно, не выполняется закон «Длина стороны не может превышать сумму длин других сторон».
Если величина косинуса превышает единицу, задача поставлена некорректно, не выполняется закон «Длина стороны не может превышать сумму длин других сторон».


В общем случае имеется два решения, расположенных симметрично относительно большого круга ''Q''₁''Q''₂. Следует явно определить, с какой стороны от направления ''Q''₁''Q''₂ находится точка ''Q''₃: если слева, как на рисунке, в последней формуле ставим знак минус, если же справа — знак плюс.
В общем случае имеется два решения, расположенных симметрично относительно большого круга ''Q''₁''Q''₂. Следует явно определить, с какой стороны от направления ''Q''₁''Q''₂ находится точка ''Q''₃: если слева, как на рисунке, то в последней формуле ставим знак минус, если же справа — знак плюс.


== Пример программной реализации ==
== Пример программной реализации ==

Версия от 19:47, 12 марта 2014

Эта страница является черновиком статьи.


Линейная засечка — это нахождение положения точки по координатам двух исходных пунктов и расстояниям от этих пунктов до определяемой точки.

Общие положения

В качестве модели Земли принимается сфера с радиусом R, равным среднему радиусу земного эллипсоида. Аналогом прямой линии на плоскости является геодезическая линия на поверхности. На сфере геодезическая линия — дуга большого круга.

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

  • φ — географическая широта,
  • λ — географическая долгота,
  • α — азимут дуги большого круга,
  • σ — сферическое расстояние (длина дуги большого круга, выраженная в долях радиуса шара).

Линейное расстояние по дуге большого круга s связано со сферическим расстоянием σ формулой s = R σ.

Постановка задачи

Исходные данные
координаты пунктов Q₁, Q₂ — φ₁, λ₁, φ₂, λ₂,
расстояния от пунктов Q₁, Q₂ до точки Q₃ — σ₁₃, σ₂₃.
Определяемые величины
координаты точки Q₃ — φ₃, λ₃.

Алгоритм

Файл:Sph lin.png
Угловая засечка

Решение любого вида засечек сводится к нахождению полярных координат искомой точки, т.е. начального направления и расстояния на неё с одного из исходных пунктов. На конечном этапе координаты находятся из решения прямой геодезической задачи. Поскольку в линейной засечке расстояния σ₁₃ и σ₂₃ уже заданы, остаётся определить направление α₁₃ или α₂₃.

На рисунке синим цветом выделены заданные элементы сферического треугольника, красным цветом неизвестные, зелёным — вспомогательные элементы. Итак, в треугольнике QQQ₃ известны только два элемента — стороны σ₁₃ и σ₂₃. Из решения обратной геодезической задачи для пунктов Q₁, Q₂ можно получить недостающий третий элемент — расстояние σ₁₂, а также азимут α₁₂.

Последовательность действий:

  1. решить обратную геодезическую задачу для Q₁, Q₂: по φ₁, λ₁, φ₂, λ₂ получить α₁₂, σ₁₂;
  2. в треугольнике QQQ₃ по σ₁₂, σ₁₃, σ₂₃ вычислить угол β₁;
  3. вычислить азимут α₁₃;
  4. решить прямую геодезическую задачу для Q₁, Q₃: по φ₁, λ₁, α₁₃, σ₁₃ вычислить φ₃, λ₃.

Действия по первому и последнему пунктам рассмотрены в статьях Задачи на сфере: обратная геодезическая задача и Задачи на сфере: прямая геодезическая задача.

Угол β₁ и азимут α₁₃ вычисляются по формулам:

Если величина косинуса превышает единицу, задача поставлена некорректно, не выполняется закон «Длина стороны не может превышать сумму длин других сторон».

В общем случае имеется два решения, расположенных симметрично относительно большого круга QQ₂. Следует явно определить, с какой стороны от направления QQ₂ находится точка Q₃: если слева, как на рисунке, то в последней формуле ставим знак минус, если же справа — знак плюс.

Пример программной реализации

Пример функции SphereAngular на языке Си, реализующей вышеизоложенный алгоритм:

/*
 * Решение угловой засечки
 *
 * Аргументы исходные:
 *     pt1   - {широта, долгота} пункта Q1
 *     pt2   - {широта, долгота} пункта Q2
 *     azi13 - азимут направления Q1-Q3
 *     azi23 - азимут направления Q2-Q3
 *
 * Аргументы определяемые:
 *     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;
}

Этот код находится в архиве Sph.tar.gz в файле sph.c. Кроме того, в файл sph.h включены следующие определения:

#define A_E 6371.0				// радиус Земли в километрах
#define Degrees(x) (x * 57.29577951308232)	// радианы -> градусы
#define Radians(x) (x / 57.29577951308232)	// градусы -> радианы

Теперь напишем программу, которая обращается к функции SphereAngular для решения угловой засечки:

#include <stdio.h>
#include <stdlib.h>
#include "sph.h"

int main(int argc, char *argv[])
{
  char buf[1024];
  double pt1[2], pt2[2], pt3[2];
  double lat1, lon1, lat2, lon2, azi13, azi23;

  while (fgets(buf, 1024, stdin) != NULL) {
    sscanf(buf, "%lf %lf %lf %lf %lf %lf",
	   &lat1, &lon1, &lat2, &lon2, &azi13, &azi23);
    pt1[0] = Radians(lat1);
    pt1[1] = Radians(lon1);
    pt2[0] = Radians(lat2);
    pt2[1] = Radians(lon2);
    if (SphereAngular(pt1, pt2, Radians(azi13), Radians(azi23), pt3))
      puts("\t"); /* Бесконечно много решений на большом круге Q1-Q2 */
    else
      printf("%f\t%f\n", Degrees(pt3[0]), Degrees(pt3[1]));
  }
  return 0;
}

В архиве Sph.tar.gz этот код находится в файле ang.c. Создадим исполняемый модуль ang компилятором gcc:

$ gcc -o ang ang.c sph.c -lm

Впрочем, в архиве есть Makefile. Для MS Windows готовую программу ang.exe можно найти в архиве Sph-win32.zip.

Программа читает данные из стандартного ввода консоли и отправляет результаты на стандартный вывод. Для чтения и записи файлов используются символы перенаправления потока «>» и «<» соответственно. Из каждой строки ввода программа считывает координаты первого и второго пунктов φ₁, λ₁, φ₂, λ₂, начальные азимуты α₁₃ и α₂₃ в градусах; решает угловую засечку; выводит в строку вывода координаты третьей точки φ₃, λ₃ в градусах.

Создадим файл ang.dat, содержащий одну строку данных:

30 0 60 30 44.80406 110.389945

После запуска программы

$ ang < ang.dat

получим φ₃, λ₃:

52.000000 54.000000

Ссылки