def METHOD_NAME(N_max, dx, dy, dz):
    from itertools import product
    NNx = range(-N_max, N_max+1, dx)
    NNy = range(-N_max, N_max+1, dy)
    NNz = range(-N_max, N_max+1, dz)
    n_vectors = np.array([[i, j, k] for (i, j, k) in product(NNx, NNy, NNz)
                          if (not (i == 0 and j == 0 and k == 0)        # exclude zero vector
                              and (k > 0                                # northern hemisphere
                                   or (k == 0 and j > 0)                # half of x-y plane
                                   or (k == 0 and j == 0 and i > 0))    # half of x axis
                              and np.sqrt(i*i + j*j + k*k) <= N_max)])  # inside sphere
    return n_vectors