distance[lat1_, long1_, lat2_, long2_] := (

phi1 = 90 - lat1;
theta1 = 360 - long1;
phi2 = 90 - lat2;
theta2 = 360 - long2;

rad[x_] := Pi x/180;
phi1 = rad[phi1];
theta1 = rad[theta1];
phi2 = rad[phi2];
theta2 = rad[theta2];

rho = 3960;
z1 = rho Cos[phi1];
z2 = rho Cos[phi2];
x1 = rho Sin[phi1] Cos[theta1];
x2 = rho Sin[phi2] Cos[theta2];
y1 = rho Sin[phi1] Sin[theta1];
y2 = rho Sin[phi2] Sin[theta2];

A = {x1, y1, z1};
B = {x2, y2, z2};

theta = ArcCos[(A . B)/(Norm[A] Norm[B])];
dist = rho theta
)

distance[34.06, 118.25, 45.50, 73.60]