Как вычислить расстояние, имея широту и долготу

var
  // Передаваемые широта/долгота в градусах и сотых долях
  StartLat: double; // Начальная широта
  StartLong: double; // Начальная долгота
  EndLat: double; // Конечная широта
  EndLong: double; // Конечная долгота

  // Переменные, используемые для вычисления смещения и расстояния
  fPhimean: Double; // Средняя широта
  fdLambda: Double; // Разница между двумя значениями долготы
  fdPhi: Double; // Разница между двумя значениями широты
  fAlpha: Double; // Смещение
  fRho: Double; // Меридианский радиус кривизны
  fNu: Double; // Поперечный радиус кривизны
  fR: Double; // Радиус сферы Земли
  fz: Double; // Угловое расстояние от центра сфероида
  fTemp: Double; // Временная переменная, использующаяся в вычислениях
  Distance: Double; // Вычисленное расстояния в метрах
  Bearing: Double; // Вычисленное от и до смещение
end

const
  // Константы, используемые для вычисления смещения и расстояния
  D2R: Double = 0.017453; // Константа для преобразования градусов в радианы
  R2D: Double = 57.295781; // Константа для преобразования радиан в градусы
  a: Double = 6378137.0; // Основные полуоси
  b: Double = 6356752.314245; // Неосновные полуоси
  e2: Double = 0.006739496742337; // Квадрат эксцентричности эллипсоида
  f: Double = 0.003352810664747; // Выравнивание эллипсоида

begin
  // Вычисляем разницу между двумя долготами и широтами и получаем среднюю широту
  fdLambda := (StartLong - EndLong) * D2R;
  fdPhi := (StartLat - EndLat) * D2R;
  fPhimean := ((StartLat + EndLat) / 2.0) * D2R;

  // Вычисляем меридианные и поперечные радиусы кривизны средней широты
  fTemp := 1 - e2 * (Power(Sin(fPhimean), 2));
  fRho := (a * (1 - e2)) / Power(fTemp, 1.5);
  fNu := a / (Sqrt(1 - e2 * (Sin(fPhimean) * Sin(fPhimean))));

  // Вычисляем угловое расстояние
  fz :=
    Sqrt(Power(Sin(fdPhi / 2.0), 2) + Cos(EndLat * D2R) * Cos(StartLat * D2R) *
      Power(Sin(fdLambda / 2.0), 2));

  fz := 2 * ArcSin(fz);

  // Вычисляем смещение
  fAlpha := Cos(EndLat * D2R) * Sin(fdLambda) * 1 / Sin(fz);
  fAlpha := ArcSin(fAlpha);

  // Вычисляем радиус Земли
  fR := (fRho * fNu) / ((fRho * Power(Sin(fAlpha), 2)) + (fNu *
    Power(Cos(fAlpha), 2)));

  // Получаем смещение и расстояние
  Distance := (fz * fR);
  if ((StartLat < EndLat) and (StartLong < EndLong)) then
    Bearing := Abs(fAlpha * R2D)
  else if ((StartLat < EndLat) and (StartLong > EndLong)) then
    Bearing := 360 - Abs(fAlpha * R2D)
  else if ((StartLat > EndLat) and (StartLong > EndLong)) then
    Bearing := 180 + Abs(fAlpha * R2D)
  else if ((StartLat > EndLat) and (StartLong < EndLong)) then
    Bearing := 180 - Abs(fAlpha * R2D);
end;