// Способ первый
function GeoToReal(Grad, Min, Sec: Byte): Extended;
begin
  Result := Grad + Min/60 + Sec/3600;
end;

// Способ второй
function GeoToReal( Grad, Min, Sec: Byte ): Extended;
begin
  Result := Grad + (Min + Sec/60) / 60;
end;

// Способ первый
// Исходный градус как число	                        37,943547
//
// Непосредственно градус (берем целое)	                37
//
// Минуты и секунды как доли градуса (вычитаем градус)  0,943547
// они же в секундах (умножаем на 3600)	                3396,7692
//
// Округляем секунды	                                3397
//
// Вычисляем минуты (делим на 60)	                56,61666667
// берем целое	                                        56
//
// Секунды как доли минут	                        0,616666667
// переводим в секунды (умножаем на 60)	                37
//
// Таким образом получаем	                        37 гр. 56 мин. 37 сек.

procedure RealToGeo(Value: Extended; var Grad, Min, Sec: Byte);
begin
  Grad := Trunc(Value);
  Value := Round(Frac(Value) * 3600);
  Min := Trunc(Value / 60);
  Sec := Round(Frac(Value / 60) * 60);
end;

// Способ второй
// Минуты и секунды как доли градуса (вычитаем градус) 0,943547
//
// Умножаем остаток на 60.
// Целая часть равна минутам, дробная - десятым долям минуты.
//
// Умножаем дробную на 60.
// Целая часть равна секундам.

procedure RealToGeo(Value: Extended; var Grad, Min, Sec: Byte);
begin
  Grad := Trunc(Value);
  Value := Frac(Value) * 60;
  Min := Trunc(Value);
  Sec := Round(Frac(Value) * 60);
end;

// Способ первый
type
  TDistance = record
    Distance: Double; // расстояние в метрах
    Bearing: Double;  // смещение в градусах
  end;

implementation

uses
  {...,} Math;

/// Начальная широта – градусы и сотые доли
/// Начальная долгота – градусы и сотые доли
/// Конечная широта – градусы и сотые доли
/// Конечная долгота – градусы и сотые доли
function GetDistance( StartLat, StartLong, EndLat, EndLong: Double ): TDistance;
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;  // Выравнивание эллипсоида
var
  // Переменные, используемые для вычисления смещения и расстояния
  fPhimean: Double; // Средняя широта
  fdLambda: Double; // Разница между двумя значениями долготы
  fdPhi: Double;    // Разница между двумя значениями широты
  fAlpha: Double;   // Смещение
  fRho: Double;     // Меридианский радиус кривизны
  fNu: Double;      // Поперечный радиус кривизны
  fR: Double;       // Радиус сферы Земли
  fz: Double;       // Угловое расстояние от центра сфероида
  fTemp: Double;    // Временная переменная, используемая в вычислениях
  Distance: Double; // Вычисленное расстояния в метрах
  Bearing: Double;  // Вычисленное от и до смещение
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 := ArcSin(Cos(EndLat * D2R) * Sin(fdLambda) / Sin(fz));
  // Вычисляем радиус Земли
  // fR := fRho * fNu / ((fRho * Power(Sin(fAlpha), 2)) +
  //       (fNu * Power(Cos(fAlpha), 2)));
  fR := fRho * fNu / (fRho * Power(Sin(fAlpha), 2) + fNu *
        Power(Cos(fAlpha), 2));
  // Получаем смещение и расстояние
  Result.Distance := fz * fR;
  if ((StartLat < EndLat) and (StartLong < EndLong)) then
     Result.Bearing := Abs(fAlpha * R2D)
  else
  if ((StartLat < EndLat) and (StartLong > EndLong)) then
    Result.Bearing := 360 - Abs(fAlpha * R2D)
  else
  if ((StartLat > EndLat) and (StartLong > EndLong)) then
    Result.Bearing := 180 + Abs(fAlpha * R2D)
  else if ((StartLat > EndLat) and (StartLong < EndLong)) then
    Result.Bearing := 180 - Abs(fAlpha * R2D);
end;

procedure TForm1.Button1Click(Sender: TObject);
var
  Distance: TDistance;
begin
  Distance := GetDistance(55.756339, 37.623612, 55.758928, 37.626016);
  ShowMessage(FloatToStr(Distance.Distance));
end;

// Способ второй [XE8]
uses
  {...,} System.Math, FMX.Maps;

function GetDistance(const aStart, aEnd: TMapCoordinate): Real;
const
  EarthRadius = 6372795;
var
  CosLatStart, SinLatStart, CosLatEnd, SinLatEnd,
  Delta, CosDelta, SinDelta, X, Y: Real;
begin
  try
    CosLatStart := Cos(DegToRad(aStart.Latitude));
    CosLatEnd := Cos(DegToRad(aEnd.Latitude));
    SinLatStart := Sin(DegToRad(aStart.Latitude));
    SinLatEnd := Sin(DegToRad(aEnd.Latitude));
    Delta := (DegToRad(aEnd.Longitude)) - (DegToRad(aStart.Longitude));
    CosDelta := Cos(Delta);
    SinDelta := Sin(Delta);
    Y := Sqrt(((CosLatEnd * SinDelta) * (CosLatEnd * SinDelta)) +
        ((CosLatStart * SinLatEnd - SinLatStart * CosLatEnd * CosDelta) *
        (CosLatStart * SinLatEnd - SinLatStart * CosLatEnd * CosDelta)));
    X := SinLatStart * SinLatEnd + CosLatStart * CosLatEnd * CosDelta;
    Result := Round(ArcTan2(Y, X) * EarthRadius);
  except
    Result := -1;
  end;
end;

procedure TForm1.Button1Click(Sender: TObject);
var
  aStart, aEnd: TMapCoordinate;
begin
  aStart.Latitude := 55.756339;
  aStart.Longitude := 37.623612;
  aEnd.Latitude := 55.758928;
  aEnd.Longitude := 37.626016;
  ShowMessage(FloatToStr(GetDistance(aStart, aEnd)));
end;

// Способ третий
// Для модели эллипсоида WGS84
uses
  {...,} System.Math;

function GetDistance(Lat1, Lon1, Lat2, Lon2: Double): Double;
var
  fdLambda: Double;
  fdPhi: Double;
  fPhimean: Double;
  fTemp: Double;
  fRho: Double;
  fNu: Double;
  fz: Double;
  fAlpha: Double;
  fR: Double;
const
  e2 = 0.0067394967423;
  a = 6378137;
begin
  Result := 0;

  if (Abs(Lon1-Lon2) <= 0.0000014) and (Abs(Lat1-Lat2) <= 0.0000008) then
    Exit
  else
  begin
    Lat1 := Lat1 * Pi / 180;
    Lon1 := Lon1 * Pi / 180;
    Lat2 := Lat2 * Pi / 180;
    Lon2 := Lon2 * Pi / 180;

    fdLambda := Lon1 - Lon2;
    fdPhi := Lat1 - Lat2;
    fPhimean := (Lat1 + Lat2) / 2;
    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 := 2 * ArcSin(Sqrt(Power(Sin(fdPhi / 2.0), 2) + Cos(Lat2) *
       Cos(Lat1) * Power(Sin(fdLambda / 2.0), 2)));
    fAlpha := ArcSin(Cos(Lat2) * Sin(fdLambda) / Sin(fz));
    fR := fRho * fNu / (fRho * Power(Sin(fAlpha), 2) + fNu * Power(Cos(fAlpha), 2));
    Result := fz * fR;
  end;
end;

procedure TForm1.Button3Click(Sender: TObject);
begin
  ShowMessage(FloatToStr(GetDistance(55.756339, 37.623612, 55.758928, 37.626016)));
end;

// Способ четвертый
uses
  {...,} System.Math;

type
  TPointF = record
    X: Double;
    Y: Double;
  end;

function GetGeoDistance(P1, P2: TPointF): Double;
var
  x1, x2, y1, y2: Double;
begin
  try
    x1 := P1.X * pi / 180;
    y1 := P1.Y * pi / 180;
    x2 := P2.X * pi / 180;
    y2 := P2.Y * pi / 180;

    Result := 6378137 * 2 * ArcSin(Sqrt(Sqr(Sin((y1 - y2) / 2)) +
       Cos(y1) * Cos(y2) * Sqr(Sin((x1 - x2) / 2)))) ;
  except
    Result := 0;
  end;
end;

procedure TForm1.Button4Click(Sender: TObject);
var
  p1, p2: TPointF;
begin
  p1.X := 37.623612;
  p1.Y := 55.756339;
  p2.X := 37.626016;
  p2.Y := 55.758928;

  ShowMessage(FloatToStr(GetGeoDistance(p1, p2)));
end;