haversin := theta->(1-cos(theta))/2;
distance := (theta1, phi1, theta2, phi2)->2*6378.14*arcsin( sqrt(haversin(theta2-theta1) + cos(theta1)*cos(theta2)*haversin(phi2-phi1)) );