Geodesic distance between two points

Geodesic distance

When we talk about distance between two points, we usually refer to the straight line distance between them. However when we want to calculate the distance between two points which are located on the surface of a sphere (here in the surface of the Earth) we are talking about the geodesic distance.

I'm just putting this function here because I needed it at work and I was using a COM library which was extremely slow (about 4 miliseconds by 1 microsecond of this solution). I don't know how exactly the code works but basically ti applies the arc length formula applying also a deviation ff according to the flattened of the Earth) given the longitude and latitude of both points. The function is adapted from a Visual Basic code from Jay Tanner to Delphi.

function GetDistanceBetween(long1, lat1, long2, lat2 : Double) : Double;
var
  F,G,L : Double;
  SF, SG, SL,
  CF, CG, CL : Double;
  W1, W2 : Double;
  S, C : Double;
  O,R,D : Double;
  H1, H2 : Double;
  ff : Double;
  temp : Double;
begin
  ff := 1 / 298.257;
  F := (lat1 + lat2) / 2;
  G := (lat1 - lat2) / 2;
  L := (long1 - long2) / 2;

  SF := Sin(F*Pi/180);
  SG := Sin(G*Pi/180);
  SL := Sin(L*Pi/180);
  CF := Cos(F*Pi/180);
  CG := Cos(G*Pi/180);
  CL := Cos(L*Pi/180);

  W1 := sqr(SG * CL);
  W2 := sqr(CF * SL);
  S := W1 + W2;

  W1 := sqr(CG * CL);
  W2 := sqr(SF * SL);
  C := W1 + W2;
  O := ArcTan(Sqrt(S/C));
  R := Sqrt(S*C) / O;
  D := 2 * O * 6378.14;

  H1 := (3*R-1) / (2*C);
  H2 := (3*R+1) / (2*S);

  W1 := sqr(SF * CG) * H1 * ff + 1;
  W2 := sqr(CF * SG) * H2 * ff;

  result := D * (W1 - W2);
end;

10
Average: 10 (2 votes)
Your rating: None