Расчет расстояния между двумя точками на земной поверхности

Delphi , Синтаксис , Математика

Автор: Dimka Maslov
WEB-сайт: http://delphibase.endimus.com

{ **** UBPFD *********** by delphibase.endimus.com ****
>> Расчет расстояния между двумя точками на земной поверхности.

Расчет расстояния между 2мя точками на земной поверхности методом Винсенти.

Dimka Maslov:
Lat1, Lon1 - широта и долгота точки 1 в градусах
Lat2, Lon2 - широта и долгота точки 2 в градусах
Функция возвращает результат в метрах.

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

Зависимости: Math
Автор:       Вячеслав
Copyright:   Опубликован в Survey Review №175 за Апрель 1976г.
Дата:        19 марта 2003 г.
***************************************************** }

function Vincenty(Lat1, Lon1, Lat2, Lon2: Extended): Extended;
const // Параметры эллипсоида:
  a = 6378245.0;
  f = 1 / 298.3;
  b = (1 - f) * a;
  EPS = 0.5E-30;
var
  APARAM, BPARAM, CPARAM, OMEGA, TanU1, TanU2,
    Lambda, LambdaPrev, SinL, CosL, USQR, U1, U2,
    SinU1, CosU1, SinU2, CosU2, SinSQSigma, CosSigma,
    TanSigma, Sigma, SinAlpha, Cos2SigmaM, DSigma: Extended;
begin
  lon1 := lon1 * (PI / 180);
  lat1 := lat1 * (PI / 180);
  lon2 := lon2 * (PI / 180);
  lat2 := lat2 * (PI / 180); //Пересчет значений координат в радианы
  TanU1 := (1 - f) * Tan(lat1);
  TanU2 := (1 - f) * Tan(lat2);
  U1 := ArcTan(TanU1);
  U2 := ArcTan(TanU2);
  SinCos(U1, SinU1, CosU1);
  SinCos(U2, SinU2, CosU2);
  OMEGA := lon2 - lon1;
  lambda := OMEGA;
  repeat //Начало цикла итерации
    LambdaPrev := lambda;
    SinCos(lambda, SinL, CosL);
    SinSQSigma := (CosU2 * SinL * CosU2 * SinL) +
      (CosU1 * SinU2 - SinU1 * CosU2 * CosL) *
      (CosU1 * SinU2 - SinU1 * CosU2 * CosL);
    CosSigma := SinU1 * SinU2 + CosU1 * CosU2 * CosL;
    TanSigma := Sqrt(SinSQSigma) / CosSigma;
    if TanSigma > 0 then
      Sigma := ArcTan(TanSigma)
    else
      Sigma := ArcTan(TanSigma) + Pi;
    if SinSQSigma = 0 then
      SinAlpha := 0
    else
      SinAlpha := CosU1 * CosU2 * SinL / Sqrt(SinSQSigma);
    if (Cos(ArcSin(SinAlpha)) * Cos(ArcSin(SinAlpha))) = 0 then
      Cos2SigmaM := 0
    else
      Cos2SigmaM := CosSigma -
        (2 * SinU1 * SinU2 / (Cos(ArcSin(SinAlpha)) * Cos(ArcSin(SinAlpha))));
    CPARAM := (f / 16) * Cos(ArcSin(SinAlpha)) * Cos(ArcSin(SinAlpha)) *
      (4 + f * (4 - 3 * Cos(ArcSin(SinAlpha)) * Cos(ArcSin(SinAlpha))));
    lambda := OMEGA + (1 - CPARAM) * f * SinAlpha * (ArcCos(CosSigma) +
      CPARAM * Sin(ArcCos(CosSigma)) * (Cos2SigmaM + CPARAM * CosSigma *
      (-1 + 2 * Cos2SigmaM * Cos2SigmaM)));
  until Abs(lambda - LambdaPrev) < EPS; // Конец цикла итерации
  USQR := Cos(ArcSin(SinAlpha)) * Cos(ArcSin(SinAlpha)) *
    (a * a - b * b) / (b * b);
  APARAM := 1 + (USQR / 16384) *
    (4096 + USQR * (-768 + USQR * (320 - 175 * USQR)));
  BPARAM := (USQR / 1024) * (256 + USQR * (-128 + USQR * (74 - 47 * USQR)));
  DSigma := BPARAM * SQRT(SinSQSigma) * (Cos2SigmaM + BPARAM / 4 *
    (CosSigma * (-1 + 2 * Cos2SigmaM * Cos2SigmaM) - BPARAM / 6 * Cos2SigmaM *
    (-3 + 4 * SinSQSigma) * (-3 + 4 * Cos2SigmaM * Cos2SigmaM)));
  Result := b * APARAM * (Sigma - DSigma);
end;

Пример использования:

procedure TForm1.Button1Click(Sender: TObject);
var
  R: Extended;
begin
  R := Vicenty(59.8833, 30.2333, 55.7667, 37.5833);
  ShowMessageFmt('%g', [R]);
end;

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

Формула Винсента - это итеративный метод для расчета расстояния и азимута между двумя точками на эллипсоидальном модели Земли. Она более точна, чем другие методы, такие как формула Хаверсина, особенно при меньших расстояниях и вблизи полюсов.

Вот разбивка кода:

1. Функция Vincenty принимает четыре параметра: Lat1, Lon1, Lat2 и Lon2, которые представляют широту и долготу двух точек на поверхности Земли.
2. Функция конвертирует входные координаты из градусов в радианы с помощью коэффициента преобразования PI/180.
3. Затем функция рассчитывает различные промежуточные значения, такие как TanU1, TanU2, U1, U2, SinL, CosL, SinU1, CosU1 и т.д., используя тригонометрические функции и математические формулы.
4. Функция итерируется, пока разница между текущим значением lambda (промежуточной переменной) и ее предыдущим значением не будет меньше маленького порога (EPS).
5. В каждом шаге итерации функция обновляет значения SinL, CosL, SinU1, CosU1 и т.д., используя ранее рассчитанные значения.
6. Наконец, функция рассчитывает расстояние между двумя точками с помощью промежуточных значений и возвращает результат.

Пример использования показывает, как вызывать функцию Vincenty с четырьмя координатами (59.8833° N, 30.2333° E, 55.7667° N, 37.5833° E) и отображать результат в сообщении.

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

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

Расчет расстояния между двумя точками на земной поверхности методом Винсенти.

Комментарии и вопросы

Получайте свежие новости и обновления по Object Pascal, Delphi и Lazarus прямо в свой смартфон. Подпишитесь на наш Telegram-канал delphi_kansoftware и будьте в курсе последних тенденций в разработке под Linux, Windows, Android и iOS

:: Главная :: Математика ::