distance[{theta1_, phi1_}, {theta2_, phi2_}] :=
 2*6378.14 ArcSin@
   Sqrt[Haversine[(theta2 - theta1) Degree] +
     Cos[theta1*Degree] Cos[theta2*Degree] Haversine[(phi2 - phi1) Degree]]