def spherical_distance(pt0, pt1):
dtheta = pt1[<NUM_LIT:0>] - pt0[<NUM_LIT:0>]<EOL>dphi   = pt1[<NUM_LIT:1>] - pt0[<NUM_LIT:1>]<EOL>a = np.sin(dphi/<NUM_LIT:2>)**<NUM_LIT:2> + np.cos(pt0[<NUM_LIT:1>]) * np.cos(pt1[<NUM_LIT:1>]) * np.sin(dtheta/<NUM_LIT:2>)**<NUM_LIT:2><EOL>return <NUM_LIT:2> * np.arcsin(np.sqrt(a))<EOL>
spherical_distance(a, b) yields the angular distance between points a and b, both of which
  should be expressed in spherical coordinates as (longitude, latitude).
If a and/or b are (2 x n) matrices, then the calculation is performed over all columns.
The spherical_distance function uses the Haversine formula; accordingly it may suffer from
rounding errors in the case of nearly antipodal points.