#DIA=2*6372.8

Procedure.d Haversine(th1.d,ph1.d,th2.d,ph2.d)
  Define dx.d,
         dy.d,
         dz.d

  ph1=Radian(ph1-ph2)
  th1=Radian(th1)
  th2=Radian(th2)

  dz=Sin(th1)-Sin(th2)
  dx=Cos(ph1)*Cos(th1)-Cos(th2)
  dy=Sin(ph1)*Cos(th1)
  ProcedureReturn ASin(Sqr(Pow(dx,2)+Pow(dy,2)+Pow(dz,2))/2)*#DIA
EndProcedure

OpenConsole("Haversine distance")
Print("Haversine distance: ")
Print(StrD(Haversine(36.12,-86.67,33.94,-118.4),7)+" km.")
Input()