Calculate the distance between two locations on a sphere.

    lat1
      The latitude of the first location.
    lon1
      The longitude of the first location.
    lat2
      The latitude of the second location.
    lon2
      The longitude of the second location.

    Returns the separation in radians. All arguments are in radians as well.
    The arguments may be vectors.

    Note that the ordering of the arguments maps to the nonstandard ordering
    ``(Dec, RA)`` in equatorial coordinates. In a spherical projection it maps
    to ``(Y, X)`` which may also be unexpected.

    The distance is computed with the "specialized Vincenty formula". Faster
    but more error-prone formulae are possible; see Wikipedia on Great-circle
    Distance.

def sphdist (lat1, lon1, lat2, lon2):
    """Calculate the distance between two locations on a sphere.

    lat1
      The latitude of the first location.
    lon1
      The longitude of the first location.
    lat2
      The latitude of the second location.
    lon2
      The longitude of the second location.

    Returns the separation in radians. All arguments are in radians as well.
    The arguments may be vectors.

    Note that the ordering of the arguments maps to the nonstandard ordering
    ``(Dec, RA)`` in equatorial coordinates. In a spherical projection it maps
    to ``(Y, X)`` which may also be unexpected.

    The distance is computed with the "specialized Vincenty formula". Faster
    but more error-prone formulae are possible; see Wikipedia on Great-circle
    Distance.

    """
    cd = np.cos (lon2 - lon1)
    sd = np.sin (lon2 - lon1)
    c2 = np.cos (lat2)
    c1 = np.cos (lat1)
    s2 = np.sin (lat2)
    s1 = np.sin (lat1)
    a = np.sqrt ((c2 * sd)**2 + (c1 * s2 - s1 * c2 * cd)**2)
    b = s1 * s2 + c1 * c2 * cd
    return np.arctan2 (a, b)