def as_unit(v, axis=1):
    """Return array of unit vectors parallel to vectors in `v`.

    Parameters
    ----------
    v : ndarray of float
    axis : int, optional
        Axis along which to normalize length.

    Returns
    -------
    ndarray of float : Unit vector of `v`, i.e. `v` divided by its
                       magnitude along `axis`.
    """
    u = np.array(v, dtype=np.float64, copy=True)
    if u.ndim == 1:
        sqmag = u.dot(u)
        if sqmag >= EPS:
            u /= sqmag ** 0.5
    else:
        if axis == 1:
            sqmag = np.einsum("...ij,...ij->...i", u, u)
        else:
            sqmag = np.einsum("...ij,...ij->...j", u, u)

        sqmag[sqmag < EPS] = 1.0
        u /= np.expand_dims(np.sqrt(sqmag), axis)
    return u