def angdist(v1, v2, zenith=True):
    if (v1.shape[0] == 2):
        v1 = sph2vec(v1, zenith=zenith)
    if (v2.shape[0] == 2):
        v2 = sph2vec(v2, zenith=zenith)
    v1 /= np.linalg.norm(v1, axis=0)
    v2 /= np.linalg.norm(v2, axis=0)
    if ((v1.ndim > 1) or (v2.ndim > 1)):
        d = np.einsum('ij,ij->j', v1, v2)
    else:
        d = np.dot(v1.T, v2)
    return np.absolute(np.arccos(d))
def angdist(v1, v2, zenith=True):
    if (v1.shape[0] == 2):
        v1 = sph2vec(v1, zenith=zenith)
    if (v2.shape[0] == 2):
        v2 = sph2vec(v2, zenith=zenith)
    v1 /= np.linalg.norm(v1, axis=0)
    v2 /= np.linalg.norm(v2, axis=0)
    if ((v1.ndim > 1) or (v2.ndim > 1)):
        d = np.einsum('ij,ij->j', v1, v2)
    else:
        d = np.dot(v1.T, v2)
    return np.absolute(np.arccos(d))