def METHOD_NAME(skycoord, rng):
    zhat = np.squeeze((skycoord.cartesian / skycoord.cartesian.norm()).xyz)
    # Random vector orthogonal to the pole:
    x = rng.uniform(size=3)
    x /= np.linalg.norm(x)
    xhat = x - (x @ zhat) * zhat
    xhat /= np.linalg.norm(xhat)
    return coord.SkyCoord(coord.CartesianRepresentation(xhat), frame=skycoord.frame)