def calculate_angles(vec_arr, ref, ref_norm=None):
    """Calculate angles between vectors in `vec_arr` and `ref` vector.

    If `ref_norm` is not provided, angle ranges between 0 and pi. If it is
    provided, angle ranges between 0 and 2pi. Note that if `ref_norm` is
    orthogonal to `vec_arr` and `ref`, then the angle is rotation around the
    axis, but if a non-orthogonal axis is provided, this may not be the case.

    Parameters
    ----------
    vec_arr : Nx3 array of float
        Array of N 3D vectors.
    ref : 1x3 array of float
        Reference vector
    ref_norm : 1x3 array of float
        Normal vector.

    Returns
    -------
    1-D array
        Array of N angles
    """
    unit_vec_arr = as_unit(vec_arr)
    unit_ref = as_unit(ref).flatten()
    ang = np.arccos(np.clip(np.dot(unit_vec_arr, unit_ref), -1.0, 1.0))
    # handle cases where a vector is the origin
    ang[np.all(unit_vec_arr == np.zeros(3), axis=1)] = 0.0
    if ref_norm is not None:
        sign = np.sign(
            np.dot(ref_norm, np.cross(unit_vec_arr, unit_ref).T)
        ).flatten()
        sign[sign == 0] = 1
        ang = rotate_angles(sign * ang, 2 * np.pi)
    return ang