def get_indices(Id, N0, N1=None, N2=None):
    """Gets 3D indices from 1D index. Assumes row-major ordering.

    Parameters
    ----------
    Id : int
        1D index
    N0 : int
        size of the array x
    N1 : int, optional, default=N0
        size of the array y
    N2 : int, optional, default=N0
        size of the array z

    Returns
    (i,j,k) : array, int, dimension=3
        corresponding 3D indices

    """
    N1 = N1 or N0
    N2 = N2 or N0
    i = int(Id / (N1 * N2) % N0)
    j = int((Id - N1 * N2 * i) / N2 % N1)
    k = int((Id - N2 * j - N2 * N1 * i) % N2)
    return i, j, k
