def generate_n_vectors(N_max, dx=1, dy=1, dz=1, half_lattice=True):
    r"""
    Generate integer vectors, :math:`\boldsymbol{n}`, with
    :math:`|\boldsymbol{n}| < N_{\rm max}`.

    If ``half_lattice=True``, only return half of the three-dimensional
    lattice. If the set N = {(i,j,k)} defines the lattice, we restrict to
    the cases such that ``(k > 0)``, ``(k = 0, j > 0)``, and
    ``(k = 0, j = 0, i > 0)``.

    .. todo::

        Return shape should be (3,N) to be consistent.

    Parameters
    ----------
    N_max : int
        Maximum norm of the integer vector.
    dx : int
        Step size in x direction. Set to 1 for odd and even terms, set
        to 2 for just even terms.
    dy : int
        Step size in y direction. Set to 1 for odd and even terms, set
        to 2 for just even terms.
    dz : int
        Step size in z direction. Set to 1 for odd and even terms, set
        to 2 for just even terms.
    half_lattice : bool (optional)
        Only return half of the 3D lattice.

    Returns
    -------
    vecs : :class:`numpy.ndarray`
        A 2D array of integers with :math:`|\boldsymbol{n}| < N_{\rm max}`
        with shape (N,3).

    """
    vecs = np.meshgrid(np.arange(-N_max, N_max+1, dx),
                       np.arange(-N_max, N_max+1, dy),
                       np.arange(-N_max, N_max+1, dz))
    vecs = np.vstack(map(np.ravel, vecs)).T
    vecs = vecs[np.linalg.norm(vecs, axis=1) <= N_max]

    if half_lattice:
        ix = ((vecs[:, 2] > 0) |
              ((vecs[:, 2] == 0) &
               (vecs[:, 1] > 0)) |
              ((vecs[:, 2] == 0) &
               (vecs[:, 1] == 0) &
               (vecs[:, 0] > 0)))
        vecs = vecs[ix]

    vecs = np.array(sorted(vecs, key=lambda x: (x[0], x[1], x[2])))
    return vecs
r"""
    Generate integer vectors, :math:`\boldsymbol{n}`, with
    :math:`|\boldsymbol{n}| < N_{\rm max}`.

    If ``half_lattice=True``, only return half of the three-dimensional
    lattice. If the set N = {(i,j,k)} defines the lattice, we restrict to
    the cases such that ``(k > 0)``, ``(k = 0, j > 0)``, and
    ``(k = 0, j = 0, i > 0)``.

    .. todo::

        Return shape should be (3,N) to be consistent.

    Parameters
    ----------
    N_max : int
        Maximum norm of the integer vector.
    dx : int
        Step size in x direction. Set to 1 for odd and even terms, set
        to 2 for just even terms.
    dy : int
        Step size in y direction. Set to 1 for odd and even terms, set
        to 2 for just even terms.
    dz : int
        Step size in z direction. Set to 1 for odd and even terms, set
        to 2 for just even terms.
    half_lattice : bool (optional)
        Only return half of the 3D lattice.

    Returns
    -------
    vecs : :class:`numpy.ndarray`
        A 2D array of integers with :math:`|\boldsymbol{n}| < N_{\rm max}`
        with shape (N,3).