Convert latitude and longitude values to alpha and beta

alphaConvert[{degrees_,minutes_}] := 
    Pi/2 - (degrees + If[degrees<0,-1,1]*minutes/60)*Pi/180;

betaConvert[{degrees_,minutes_}] := 
    (degrees + If[degrees<0,-1,1]*minutes/60)*Pi/180;

Convert Spherical coordinatres to rectangular

z[R_, alpha_, beta_] := R Cos[alpha];
r[R_, alpha_, beta_] := R Sin[alpha];
x[R_, alpha_, beta_] := r[R, alpha, beta] Cos[beta];
y[R_, alpha_, beta_] := r[R, alpha, beta] Sin[beta];

Tests of x, y, and z

{x[R,Pi/2,Pi], y[R,Pi/2,Pi], z[R,Pi/2,Pi]}

{x[R,0,Pi], y[R,0,Pi], z[R,0,Pi]}

{x[R,Pi/2,Pi/4], y[R,Pi/2,Pi/4], z[R,Pi/2,Pi/4]}

Surface Distance between two points

The following cell gives the cosine of the angle between the vectors

{x[R,alpha1,beta1],y[R,alpha1,beta1],z[R,alpha1,beta1]}.
         {x[R,alpha2,beta2],y[R,alpha2,beta2],z[R,alpha2,beta2]}/
         R^2

Now we do some simplification

Factor[%]

Simplify[%-%[[1]]]

angle = %+%%[[1]]

angleBetween[alpha1_,beta1_,alpha2_,beta2_] =
  ArcCos[angle]

distanceBetween [lat1_, lon1_, lat2_, lon2_] :=
3963 angleBetween[alphaConvert[lat1],betaConvert[lon1],
                  alphaConvert[lat2],betaConvert[lon2]]//N

Now some city data to run tests.

indyLong = {-86,9}; indyLat ={ 39,46};
chiLat = {41,53} ; chiLong={-87,38};
nyLat = {40,43}; nyLong = {-74,01};
parisLat = {48,42}; parisLong = {2,20};
laLat = {34,03}; laLong = {-118,15};
delhiLat = {28,40}; delhiLong={77,13};
rioLat= {-22,54}; rioLong = {-43,15};

distanceBetween[indyLat,indyLong,chiLat,chiLong]

distanceBetween[indyLat,indyLong,nyLat,nyLong]

distanceBetween[laLat,laLong,nyLat,nyLong]

distanceBetween[parisLat,parisLong,nyLat,nyLong]

distanceBetween[laLat,laLong,indyLat,indyLong]

distanceBetween[rioLat,rioLong,nyLat,nyLong]

distanceBetween[rioLat,rioLong,delhiLat,delhiLong]