diff --git a/kwant/continuum/__init__.py b/kwant/continuum/__init__.py
index fd08848f69cfae23f6e6612e74167427c675f8e3..bf6da67bdf8861443851148a4f3c25bf99546429 100644
--- a/kwant/continuum/__init__.py
+++ b/kwant/continuum/__init__.py
@@ -10,10 +10,12 @@ try:
     from .discretizer import discretize, discretize_symbolic, build_discretized
     from ._common import sympify, lambdify
     from ._common import momentum_operators, position_operators
+    from .landau_levels import to_landau_basis, discretize_landau, LandauLattice
 except ImportError as error:
     msg = ("'kwant.continuum' is not available because one or more of its "
            "dependencies is not installed.")
     raise ImportError(msg) from error
 
 __all__ = ['discretize', 'discretize_symbolic', 'build_discretized',
+           'to_landau_basis', 'discretize_landau', 'LandauLattice',
            'sympify', 'lambdify', 'momentum_operators', 'position_operators']
diff --git a/kwant/continuum/landau_levels.py b/kwant/continuum/landau_levels.py
new file mode 100644
index 0000000000000000000000000000000000000000..90544579724c1afa56cf061b6b94b96395e735a8
--- /dev/null
+++ b/kwant/continuum/landau_levels.py
@@ -0,0 +1,454 @@
+# Copyright 2011-2019 Kwant authors.
+#
+# This file is part of Kwant.  It is subject to the license terms in the file
+# LICENSE.rst found in the top-level directory of this distribution and at
+# http://kwant-project.org/license.  A list of Kwant authors can be found in
+# the file AUTHORS.rst at the top-level directory of this distribution and at
+# http://kwant-project.org/authors.
+
+from keyword import iskeyword
+import functools
+import operator
+import collections
+
+from math import sqrt
+import tinyarray as ta
+import numpy as np
+import sympy
+
+import kwant.lattice
+import kwant.builder
+
+import kwant.continuum
+import kwant.continuum._common
+import kwant.continuum.discretizer
+from kwant.continuum import momentum_operators, position_operators
+
+__all__ = ["to_landau_basis", "discretize_landau"]
+
+coordinate_vectors = dict(zip("xyz", np.eye(3)))
+
+ladder_lower, ladder_raise = sympy.symbols(r"a a^\dagger", commutative=False)
+
+
+def to_landau_basis(hamiltonian, momenta=None):
+    r"""Replace two momenta by Landau level ladder operators.
+
+    Replaces:
+
+        k_0 -> sqrt(B/2) * (a + a^\dagger)
+        k_1 -> 1j * sqrt(B/2) * (a - a^\dagger)
+
+    Parameters
+    ----------
+    hamiltonian : str or sympy expression
+        The Hamiltonian to which to apply the Landau level transformation.
+    momenta : sequence of str (optional)
+        The momenta to replace with Landau level ladder operators. If not
+        provided, 'k_x' and 'k_y' are used
+
+    Returns
+    -------
+    hamiltonian : sympy expression
+    momenta : sequence of sympy atoms
+        The momentum operators that have been replaced by ladder operators.
+    normal_coordinate : sympy atom
+        The remaining position coordinate. May or may not be present
+        in 'hamiltonian'.
+    """
+    hamiltonian = kwant.continuum.sympify(hamiltonian)
+    momenta = _normalize_momenta(momenta)
+    normal_coordinate = _find_normal_coordinate(hamiltonian, momenta)
+
+    # Substitute ladder operators for Landau level momenta
+    B = sympy.symbols("B")
+    hamiltonian = hamiltonian.subs(
+        {
+            momenta[0]: sympy.sqrt(abs(B) / 2) * (ladder_raise + ladder_lower),
+            momenta[1]: sympy.I
+            * sympy.sqrt(abs(B) / 2)
+            * (ladder_lower - ladder_raise),
+        }
+    )
+
+    return hamiltonian, momenta, normal_coordinate
+
+
+def discretize_landau(hamiltonian, N, momenta=None, grid_spacing=1):
+    """Discretize a Hamiltonian in a basis of Landau levels.
+
+    Parameters
+    ----------
+    hamiltonian : str or sympy expression
+    N : positive integer
+        The number of Landau levels in the basis.
+    momenta : sequence of str (optional)
+        The momenta defining the plane perpendicular to the magnetic field.
+        If not provided, "k_x" and "k_y" are used.
+    grid_spacing : float, default: 1
+        The grid spacing to use when discretizing the normal direction
+        (parallel to the magnetic field).
+
+    Returns
+    -------
+    builder : `~kwant.builder.Builder`
+        'hamiltonian' discretized in a basis of Landau levels in the plane
+        defined by 'momenta'. If a third coordinate is present in 'hamiltonian',
+        then this also has a translational symmetry in that coordinate direction.
+    """
+
+    if N <= 0:
+        raise ValueError("N must be positive")
+
+    hamiltonian, momenta, normal_coordinate = to_landau_basis(hamiltonian, momenta)
+
+    # Discretize in normal direction and split terms for onsites/hoppings into
+    # monomials in ladder operators.
+    tb_hamiltonian, _ = kwant.continuum.discretize_symbolic(
+        hamiltonian, coords=[normal_coordinate.name]
+    )
+    tb_hamiltonian = {
+        key: kwant.continuum._common.monomials(value, gens=(ladder_lower, ladder_raise))
+        for key, value in tb_hamiltonian.items()
+    }
+
+    # Replace ladder operator monomials by tuple of integers:
+    # e.g. a^\dagger a^2 -> (+1, -2).
+    tb_hamiltonian = {
+        outer_key: {
+            _ladder_term(inner_key): inner_value
+            for inner_key, inner_value in outer_value.items()
+        }
+        for outer_key, outer_value in tb_hamiltonian.items()
+    }
+
+    # Construct map from LandauLattice HoppingKinds to a sequence of pairs
+    # that encode the ladder operators and multiplying expression.
+    tb_hoppings = collections.defaultdict(list)
+    for outer_key, outer_value in tb_hamiltonian.items():
+        for inner_key, inner_value in outer_value.items():
+            tb_hoppings[(*outer_key, sum(inner_key))] += [(inner_key, inner_value)]
+    # Extract the number of orbitals on each site/hopping
+    random_element = next(iter(tb_hoppings.values()))[0][1]
+    norbs = 1 if isinstance(random_element, sympy.Expr) else random_element.shape[0]
+    tb_onsite = tb_hoppings.pop((0, 0), None)
+
+    # Construct Builder
+    if _has_coordinate(normal_coordinate, hamiltonian):
+        sym = kwant.lattice.TranslationalSymmetry([grid_spacing, 0])
+    else:
+        sym = kwant.builder.NoSymmetry()
+    lat = LandauLattice(grid_spacing, norbs=norbs)
+    syst = kwant.Builder(sym)
+
+    # Add onsites
+    landau_sites = (lat(0, j) for j in range(N))
+    if tb_onsite is None:
+        syst[landau_sites] = ta.zeros((norbs, norbs))
+    else:
+        syst[landau_sites] = _builder_value(
+            tb_onsite, normal_coordinate.name, grid_spacing, is_onsite=True
+        )
+
+    # Add zero hoppings between adjacent Landau levels.
+    # Necessary to be able to use the Landau level builder
+    # to populate another builder using builder.fill().
+    syst[kwant.builder.HoppingKind((0, 1), lat)] = ta.zeros((norbs, norbs))
+
+    # Add the hoppings from the Hamiltonian
+    for hopping, parts in tb_hoppings.items():
+        syst[kwant.builder.HoppingKind(hopping, lat)] = _builder_value(
+            parts, normal_coordinate.name, grid_spacing, is_onsite=False
+        )
+
+    return syst
+
+
+# This has to subclass lattice so that it will work with TranslationalSymmetry.
+class LandauLattice(kwant.lattice.Monatomic):
+    """
+    A `~kwant.lattice.Monatomic` lattice with a Landau level index per site.
+
+    Site tags (see `~kwant.builder.SiteFamily`) are pairs of integers, where
+    the first integer describes the real space position and the second the
+    Landau level index.
+
+    The real space Bravais lattice is one dimensional, oriented parallel
+    to the magnetic field. The Landau level index represents the harmonic
+    oscillator basis states used for the Landau quantization in the plane
+    normal to the field.
+
+    Parameters
+    ----------
+    grid_spacing : float
+        Real space lattice spacing (parallel to the magnetic field).
+    offset : float (optional)
+        Displacement of the lattice origin from the real space
+        coordinates origin.
+    """
+
+    def __init__(self, grid_spacing, offset=None, name="", norbs=None):
+        if offset is not None:
+            offset = [offset, 0]
+        # The second vector and second coordinate do not matter (they are
+        # not used in pos())
+        super().__init__([[grid_spacing, 0], [0, 1]], offset, name, norbs)
+
+    def pos(self, tag):
+        return ta.array((self.prim_vecs[0, 0] * tag[0] + self.offset[0],))
+
+    def landau_index(self, tag):
+        return tag[-1]
+
+
+def _builder_value(terms, normal_coordinate, grid_spacing, is_onsite):
+    """Construct an onsite/hopping value function from a list of terms
+
+    Parameters
+    ----------
+    terms : list
+        Each element is a pair (ladder_term, sympy expression/matrix).
+        ladder_term is a tuple of integers that encodes a string of
+        Landau raising/lowering operators and the sympy expression
+        is the rest
+    normal_coordinate : str
+    grid_spacing : float
+        The grid spacing in the normal direction
+    is_onsite : bool
+        True if we are constructing an onsite value function
+    """
+    ladder_term_symbols = [sympy.Symbol(_ladder_term_name(lt)) for lt, _ in terms]
+    # Construct a single expression from the terms, multiplying each part
+    # by the placeholder that represents the prefactor from the ladder operator term.
+    (ladder, (_, part)), *rest = zip(ladder_term_symbols, terms)
+    expr = ladder * part
+    for ladder, (_, part) in rest:
+        expr += ladder * part
+    expr = expr.subs(
+        {kwant.continuum.discretizer._displacements[normal_coordinate]: grid_spacing}
+    )
+    # Construct the return string and temporary variable names
+    # for function calls.
+    return_string, map_func_calls, const_symbols, _cache = kwant.continuum.discretizer._return_string(
+        expr, coords=[normal_coordinate]
+    )
+
+    # Remove the ladder term placeholders, as these are not parameters
+    const_symbols = set(const_symbols)
+    for ladder_term in ladder_term_symbols:
+        const_symbols.discard(ladder_term)
+
+    # Construct the argument part of signature. Arguments
+    # consist of any constants and functions in the return string.
+    arg_names = set.union(
+        {s.name for s in const_symbols}, {str(k.func) for k in map_func_calls}
+    )
+    arg_names.discard("Abs")  # Abs function is not a parameter
+    for arg_name in arg_names:
+        if not (arg_name.isidentifier() and not iskeyword(arg_name)):
+            raise ValueError(
+                "Invalid name in used symbols: {}\n"
+                "Names of symbols used in Hamiltonian "
+                "must be valid Python identifiers and "
+                "may not be keywords".format(arg_name)
+            )
+    arg_names = ", ".join(sorted(arg_names))
+    # Construct site part of the function signature
+    site_string = "from_site" if is_onsite else "to_site, from_site"
+    # Construct function signature
+    if arg_names:
+        function_header = "def _({}, {}):".format(site_string, arg_names)
+    else:
+        function_header = "def _({}):".format(site_string)
+    # Construct function body
+    function_body = []
+    if "B" not in arg_names:
+        # B is not a parameter for terms with no ladder operators but we still
+        # need something to pass to _evaluate_ladder_term
+        function_body.append("B = +1")
+    function_body.extend(
+        [
+            "{}, = from_site.pos".format(normal_coordinate),
+            "_ll_index = from_site.family.landau_index(from_site.tag)",
+        ]
+    )
+    # To get the correct hopping if B < 0, we need to set the Hermitian
+    # conjugate of the ladder operator matrix element, which swaps the
+    # from_site and to_site Landau level indices.
+    if not is_onsite:
+        function_body.extend(
+            ["if B < 0:",
+             "    _ll_index = to_site.family.landau_index(to_site.tag)"
+            ]
+        )
+    function_body.extend(
+        "{} = _evaluate_ladder_term({}, _ll_index, B)".format(symb.name, lt)
+        for symb, (lt, _) in zip(ladder_term_symbols, terms)
+    )
+    function_body.extend(
+        "{} = {}".format(v, kwant.continuum.discretizer._print_sympy(k))
+        for k, v in map_func_calls.items()
+    )
+    function_body.append(return_string)
+    func_code = "\n    ".join([function_header] + function_body)
+    # Add "I" to namespace just in case sympy again would miss to replace it
+    # with Python's 1j as it was the case with SymPy 1.2 when "I" was argument
+    # of some function.
+    namespace = {
+        "pi": np.pi,
+        "I": 1j,
+        "_evaluate_ladder_term": _evaluate_ladder_term,
+        "Abs": abs,
+    }
+    namespace.update(_cache)
+    # Construct full source, including cached arrays
+    source = []
+    for k, v in _cache.items():
+        source.append("{} = (\n{})\n".format(k, repr(np.array(v))))
+    source.append(func_code)
+    exec(func_code, namespace)
+    f = namespace["_"]
+    f._source = "".join(source)
+
+    return f
+
+
+def _ladder_term(operator_string):
+    r"""Return a tuple of integers representing a string of ladder operators
+
+    Parameters
+    ----------
+    operator_string : Sympy expression
+        Monomial in ladder operators, e.g. a^\dagger a^2 a^\dagger.
+
+    Returns
+    -------
+    ladder_term : tuple of int
+        e.g. a^\dagger a^2 -> (+1, -2)
+    """
+    ret = []
+    for factor in operator_string.as_ordered_factors():
+        ladder_op, exponent = factor.as_base_exp()
+        if ladder_op == ladder_lower:
+            sign = -1
+        elif ladder_op == ladder_raise:
+            sign = +1
+        else:
+            sign = 0
+        ret.append(sign * int(exponent))
+    return tuple(ret)
+
+
+def _ladder_term_name(ladder_term):
+    """
+    Parameters
+    ----------
+    ladder_term : tuple of int
+
+    Returns
+    -------
+    ladder_term_name : str
+    """
+
+    def ns(i):
+        if i >= 0:
+            return str(i)
+        else:
+            return "_" + str(-i)
+
+    return "_ladder_{}".format("_".join(map(ns, ladder_term)))
+
+
+def _evaluate_ladder_term(ladder_term, n, B):
+    r"""Evaluates the prefactor for a ladder operator on a landau level.
+
+    Example: a^\dagger a^2 -> (n - 1) * sqrt(n)
+
+    Parameters
+    ----------
+    ladder_term : tuple of int
+        Represents a string of ladder operators. Positive
+        integers represent powers of the raising operator,
+        negative integers powers of the lowering operator.
+    n : non-negative int
+        Landau level index on which to act with ladder_term.
+    B : float
+        Magnetic field with sign
+
+    Returns
+    -------
+    ladder_term_prefactor : float
+    """
+    assert n >= 0
+    # For negative B we swap a and a^\dagger.
+    if B < 0:
+        ladder_term = tuple(-i for i in ladder_term)
+    ret = 1
+    for m in reversed(ladder_term):
+        if m > 0:
+            factors = range(n + 1, n + m + 1)
+        elif m < 0:
+            factors = range(n + m + 1, n + 1)
+            if n == 0:
+                return 0  # a|0> = 0
+        else:
+            factors = (1,)
+        ret *= sqrt(functools.reduce(operator.mul, factors))
+        n += m
+    return ret
+
+
+def _normalize_momenta(momenta=None):
+    """Return Landau level momenta as Sympy atoms
+
+    Parameters
+    ----------
+    momenta : None or pair of int or pair of str
+        The momenta to choose. If None then 'k_x' and 'k_y'
+        are chosen. If integers, then these are the indices
+        of the momenta: 0 → k_x, 1 → k_y, 2 → k_z. If strings,
+        then these name the momenta.
+
+    Returns
+    -------
+    The specified momenta as sympy atoms.
+    """
+
+    # Specify which momenta to substitute for the Landau level basis.
+    if momenta is None:
+        # Use k_x and k_y by default
+        momenta = momentum_operators[:2]
+    else:
+        if len(momenta) != 2:
+            raise ValueError("Two momenta must be specified.")
+
+        k_names = [k.name for k in momentum_operators]
+
+        if all([type(i) is int for i in momenta]) and all(
+            [i >= 0 and i < 3 for i in momenta]
+        ):
+            momenta = [momentum_operators[i] for i in momenta]
+        elif all([isinstance(momentum, str) for momentum in momenta]) and all(
+            [momentum in k_names for momentum in momenta]
+        ):
+            momenta = [
+                momentum_operators[k_names.index(momentum)] for momentum in momenta
+            ]
+        else:
+            raise ValueError("Momenta must all be integers or strings.")
+
+    return tuple(momenta)
+
+
+def _find_normal_coordinate(hamiltonian, momenta):
+    discrete_momentum = next(
+        momentum for momentum in momentum_operators if momentum not in momenta
+    )
+    normal_coordinate = position_operators[momentum_operators.index(discrete_momentum)]
+    return normal_coordinate
+
+
+def _has_coordinate(coord, expr):
+    momentum = momentum_operators[position_operators.index(coord)]
+    atoms = set(expr.atoms())
+    return coord in atoms or momentum in atoms
diff --git a/kwant/continuum/tests/test_landau_levels.py b/kwant/continuum/tests/test_landau_levels.py
new file mode 100644
index 0000000000000000000000000000000000000000..8a568e6f4a4a716f1daa2c82542f47ea237a0f9d
--- /dev/null
+++ b/kwant/continuum/tests/test_landau_levels.py
@@ -0,0 +1,252 @@
+# Copyright 2011-2019 Kwant authors.
+#
+# This file is part of Kwant.  It is subject to the license terms in the file
+# LICENSE.rst found in the top-level directory of this distribution and at
+# http://kwant-project.org/license.  A list of Kwant authors can be found in
+# the file AUTHORS.rst at the top-level directory of this distribution and at
+# http://kwant-project.org/authors.
+
+from math import sqrt
+
+import numpy as np
+import sympy
+import pytest
+import itertools
+
+import kwant.builder
+import kwant.lattice
+
+from .._common import position_operators, momentum_operators, sympify
+from ..landau_levels import (
+    to_landau_basis,
+    discretize_landau,
+    _ladder_term,
+    _evaluate_ladder_term,
+    ladder_lower,
+    ladder_raise,
+    LandauLattice,
+)
+
+x, y, z = position_operators
+k_x, k_y, k_z = momentum_operators
+B = sympy.symbols("B")
+V = sympy.symbols("V", cls=sympy.Function)
+a, a_dag = ladder_lower, ladder_raise
+
+
+def test_ladder_term():
+    assert _ladder_term(a ** 2 * a_dag) == (-2, 1)
+    assert _ladder_term(a_dag ** 5 * a ** 3 * a_dag) == (5, -3, 1)
+    # non-ladder operators give a shift of 0
+    assert _ladder_term(B * a ** 2) == (0, -2)
+
+
+def test_evaluate_ladder_term():
+    assert np.isclose(_evaluate_ladder_term((+1, -1), 1, +1), 1)
+    assert np.isclose(_evaluate_ladder_term((+1, -1), 2, +1), 2)
+    assert np.isclose(_evaluate_ladder_term((-2, +3, -2), 5, +1), 4 * 5 * 6 * sqrt(5))
+    # annihilating |0> is always 0
+    assert _evaluate_ladder_term((-1,), 0, +1) == 0
+    assert _evaluate_ladder_term((-2,), 1, +1) == 0
+    assert _evaluate_ladder_term((-3,), 1, +1) == 0
+    assert _evaluate_ladder_term((+3, -2), 1, +1) == 0
+    assert _evaluate_ladder_term((-3, -2, +3), 1, +1) == 0
+    # negative B swaps creation and annihilation operators
+    assert _evaluate_ladder_term((+1, -1), 2, +1) == _evaluate_ladder_term(
+        (-1, +1), 2, -1
+    )
+    assert _evaluate_ladder_term((-2, +3, -2), 5, +1) == _evaluate_ladder_term(
+        (+2, -3, +2), 5, -1
+    )
+
+
+def test_to_landau_basis():
+    # test basic usage
+    ham, momenta, normal_coord = to_landau_basis("k_x**2 + k_y**2")
+    assert sympy.expand(ham) == abs(B) * a * a_dag + abs(B) * a_dag * a
+    assert momenta == (k_x, k_y)
+    assert normal_coord == z
+
+    # test that hamiltonian can be specified as a sympy expression
+    ham, momenta, normal_coord = to_landau_basis(sympify("k_x**2 + k_y**2"))
+    assert sympy.expand(ham) == abs(B) * a * a_dag + abs(B) * a_dag * a
+    assert momenta == (k_x, k_y)
+    assert normal_coord == z
+
+    # test that
+    ham, momenta, normal_coord = to_landau_basis("k_x**2 + k_y**2 + k_z**2 + V(z)")
+    assert sympy.expand(ham) == (
+        abs(B) * a * a_dag + abs(B) * a_dag * a + k_z ** 2 + V(z)
+    )
+    assert momenta == (k_x, k_y)
+    assert normal_coord == z
+
+    # test for momenta explicitly specified
+    ham, momenta, normal_coord = to_landau_basis(
+        "k_x**2 + k_y**2 + k_z**2 + k_x*k_y", momenta=("k_z", "k_y")
+    )
+    assert sympy.expand(ham) == (
+        abs(B) * a * a_dag
+        + abs(B) * a_dag * a
+        + k_x ** 2
+        + sympy.I * sympy.sqrt(abs(B) / 2) * k_x * a
+        - sympy.I * sympy.sqrt(abs(B) / 2) * k_x * a_dag
+    )
+    assert momenta == (k_z, k_y)
+    assert normal_coord == x
+
+
+def test_discretize_landau():
+    n_levels = 10
+    magnetic_field = 1 / 3  # a suitably arbitrary value
+    # Ensure that we can handle products of ladder operators by truncating
+    # several levels higher than the number of levels we actually want.
+    a = np.diag(np.sqrt(np.arange(1, n_levels + 5)), k=1)
+    a_dag = a.conjugate().transpose()
+    k_x = sqrt(magnetic_field / 2) * (a + a_dag)
+    k_y = 1j * sqrt(magnetic_field / 2) * (a - a_dag)
+    sigma_0 = np.eye(2)
+    sigma_x = np.array([[0, 1], [1, 0]])
+    sigma_y = np.array([[0, -1j], [1j, 0]])
+
+    # test that errors are raised on invalid input
+    with pytest.raises(ValueError):
+        discretize_landau("k_x**2 + k_y**2", N=0)
+    with pytest.raises(ValueError):
+        discretize_landau("k_x**2 + k_y**2", N=-1)
+
+    # test a basic Hamiltonian with no normal coordinate dependence
+    syst = discretize_landau("k_x**2 + k_y**2", N=n_levels)
+    lat = LandauLattice(1, norbs=1)
+    assert isinstance(syst.symmetry, kwant.builder.NoSymmetry)
+    syst = syst.finalized()
+    assert set(syst.sites) == {lat(0, j) for j in range(n_levels)}
+    assert np.allclose(
+        syst.hamiltonian_submatrix(params=dict(B=0)), np.zeros((n_levels, n_levels))
+    )
+    should_be = k_x @ k_x + k_y @ k_y
+    assert np.allclose(
+        syst.hamiltonian_submatrix(params=dict(B=magnetic_field)),
+        should_be[:n_levels, :n_levels],
+    )
+
+    # test negative magnetic field
+    assert np.allclose(
+        syst.hamiltonian_submatrix(params=dict(B=-magnetic_field)),
+        should_be[:n_levels, :n_levels],
+    )
+
+    # test hamiltonian with no onsite elements
+    syst = discretize_landau("k_x", N=n_levels)
+    syst = syst.finalized()
+    assert np.allclose(
+        syst.hamiltonian_submatrix(params=dict(B=magnetic_field)),
+        k_x[:n_levels, :n_levels],
+    )
+
+    # test a basic Hamiltonian with normal coordinate dependence
+    grid = 1 / 7  # a suitably arbitrary value
+    syst = discretize_landau(
+        "k_x**2 + k_y**2 + k_z**2 + V(z)", N=n_levels, grid_spacing=grid
+    )
+    assert isinstance(syst.symmetry, kwant.lattice.TranslationalSymmetry)
+    syst = syst.finalized()
+    zero_potential = syst.cell_hamiltonian(params=dict(B=magnetic_field, V=lambda z: 0))
+    with_potential = syst.cell_hamiltonian(params=dict(B=magnetic_field, V=lambda z: 1))
+    # extra +2/grid**2 from onsite kinetic energy
+    assert np.allclose(
+        zero_potential,
+        should_be[:n_levels, :n_levels] + (2 / grid ** 2) * np.eye(n_levels),
+    )
+    # V(z) just adds the same energy to each onsite
+    assert np.allclose(with_potential - zero_potential, np.eye(n_levels))
+    # hopping matrix does not exchange landau levels
+    assert np.allclose(
+        syst.inter_cell_hopping(params=dict(B=magnetic_field, V=lambda z: 0)),
+        -np.eye(n_levels) / grid ** 2,
+    )
+
+    # test a Hamiltonian with mixing between Landau levels
+    # and spatial degrees of freedom.
+    syst = discretize_landau("k_x**2 + k_y**2 + k_x*k_z", N=n_levels)
+    syst = syst.finalized()
+    assert np.allclose(
+        syst.inter_cell_hopping(params=dict(B=magnetic_field)),
+        -1j * k_x[:n_levels, :n_levels] / 2,
+    )
+
+    # test a Hamiltonian with extra degrees of freedom
+    syst = discretize_landau("sigma_0 * k_x**2 + sigma_x * k_y**2", N=n_levels)
+    syst = syst.finalized()
+    assert syst.sites[0].family.norbs == 2
+    should_be = np.kron(k_x @ k_x, sigma_0) + np.kron(k_y @ k_y, sigma_x)
+    assert np.allclose(
+        syst.hamiltonian_submatrix(params=dict(B=magnetic_field)),
+        should_be[: 2 * n_levels, : 2 * n_levels],
+    )
+
+    # test a linear Hamiltonian
+    syst = discretize_landau("sigma_y * k_x - sigma_x * k_y", N=n_levels)
+    syst = syst.finalized()
+    should_be = np.kron(k_x, sigma_y) - np.kron(k_y, sigma_x)
+    assert np.allclose(
+        syst.hamiltonian_submatrix(params=dict(B=magnetic_field)),
+        should_be[: 2 * n_levels, : 2 * n_levels],
+    )
+    assert np.allclose(
+        syst.hamiltonian_submatrix(params=dict(B=magnetic_field)),
+        syst.hamiltonian_submatrix(params=dict(B=-magnetic_field)),
+    )
+
+
+def test_analytical_spectrum():
+    hamiltonian = """(k_x**2 + k_y**2) * sigma_0 +
+                    alpha * (k_x * sigma_y - k_y * sigma_x) +
+                    g * B * sigma_z"""
+
+    def exact_Es(n, B, alpha, g):
+        # See e.g. R. Winkler (2003), section 8.4.1
+        sign_B = np.sign(B)
+        B = np.abs(B)
+        Ep = 2*B*(n+1) - 0.5*np.sqrt((2*B - sign_B*2*g*B)**2 + 8*B*alpha**2*(n+1))
+        Em = 2*B*n + 0.5*np.sqrt((2*B - sign_B*2*g*B)**2 + 8*B*alpha**2*n)
+        return Ep, Em
+
+    N = 20
+    syst = discretize_landau(hamiltonian, N)
+    syst = syst.finalized()
+    params = dict(alpha=0.07, g=0.04)
+    for _ in range(5):
+        B = 0.01 + np.random.rand()*3
+        params['B'] = B
+        exact = [exact_Es(n, B, params['alpha'], params['g']) for n in range(N)]
+        # We only check the first N levels - the SOI couples adjacent levels,
+        # so the higher ones are not necessarily determined accurately in the
+        # discretization
+        exact = np.sort([energy for energies in exact for energy in energies])[:N]
+        ll_spect = np.linalg.eigvalsh(syst.hamiltonian_submatrix(params=params))[:len(exact)]
+        assert np.allclose(ll_spect, exact)
+
+
+
+def test_fill():
+
+    def shape(site, lower, upper):
+        (z, )= site.pos
+        return lower <= z < upper
+
+    hamiltonian = "k_x**2 + k_y**2 + k_z**2"
+    N = 6
+    template = discretize_landau(hamiltonian, N)
+
+    syst = kwant.Builder()
+    width = 4
+    syst.fill(template, lambda site: shape(site, 0, width), (0, ));
+
+    correct_tags = [(coordinate, ll_index) for coordinate, ll_index
+                    in itertools.product(range(width), range(N))]
+
+    syst_tags = [site.tag for site in syst.sites()]
+
+    assert len(syst_tags) == len(correct_tags)
+    assert all(tag in correct_tags for tag in syst_tags)
diff --git a/ll.ipynb b/ll.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..b6f7dda7634d86d787ed7908ae7e8ec99d9b2636
--- /dev/null
+++ b/ll.ipynb
@@ -0,0 +1,1080 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stderr",
+     "output_type": "stream",
+     "text": [
+      "/home/tinkerer/kwant/kwant/solvers/default.py:18: RuntimeWarning: MUMPS is not available, SciPy built-in solver will be used as a fallback. Performance can be very poor in this case.\n",
+      "  \"Performance can be very poor in this case.\", RuntimeWarning)\n"
+     ]
+    }
+   ],
+   "source": [
+    "import numpy as np\n",
+    "import sympy\n",
+    "import operator\n",
+    "import inspect\n",
+    "import itertools\n",
+    "import functools\n",
+    "from copy import deepcopy\n",
+    "\n",
+    "import kwant\n",
+    "import kwant.continuum"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "import sympy\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "sympy.init_printing(print_builtin=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from kwant.continuum import landau_levels"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "builder = kwant.continuum.discretize('k_x**2 + k_y**2 + V(x, y) + sqrt(B)')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "# Discrete coordinates: x y\n",
+      "\n",
+      "# Onsite element:\n",
+      "def onsite(site, B, V):\n",
+      "    (x, y, ) = site.pos\n",
+      "    _const_0 = (V(x, y))\n",
+      "    return (B**(1/2) + _const_0 + 4.0)\n",
+      "\n",
+      "# Hopping from (1, 0):\n",
+      "(-1+0j)\n",
+      "\n",
+      "# Hopping from (0, 1):\n",
+      "(-1+0j)\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(builder)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 35,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "a_dag = sympy.symbols('a^\\dagger')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 36,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAABIAAAATBAMAAABvvEDBAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIpmJdu8QRM1mu90yVKvMIHo8AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAdUlEQVQIHWNggIDnUJqBAcEyg4mtWBQLY7olwFhvYAyGtVhYb+FiIPMYlVXCGBh0gKxSAY4PDAyGDAxskQxMCxgY3BgYmL4xcDswMAgwMLAuYMg/wAAC+QYM98EMhvwLDC94gXIMDFwOnPuYwYJsT3zk3EEsAMsTFcMrnqiMAAAAAElFTkSuQmCC\n",
+      "text/latex": [
+       "$$a^\\dagger$$"
+      ],
+      "text/plain": [
+       "a__\\dagger"
+      ]
+     },
+     "execution_count": 36,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "a_dag"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 24,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "chain = kwant.lattice.chain()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 25,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "kwant.builder.Site"
+      ]
+     },
+     "execution_count": 25,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "type(chain(1))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 26,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "Site(kwant.lattice.Monatomic([[1.0]], [0.0], '', None), array([2]))"
+      ]
+     },
+     "execution_count": 26,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "chain(2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "Site(kwant.lattice.Monatomic([[1.0]], [0.0], '', None), array([1]))"
+      ]
+     },
+     "execution_count": 17,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "kwant.builder.Site(chain, (1,))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 28,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "s = kwant.builder.Site(chain, (40,))  # Second argument is the tag, for the chain its the position."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 30,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([40.0])"
+      ]
+     },
+     "execution_count": 30,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "s.pos"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "\u001b[0;31mInit signature:\u001b[0m \u001b[0mkwant\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuilder\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mSite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfamily\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtag\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_i_know_what_i_do\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+       "\u001b[0;31mDocstring:\u001b[0m     \n",
+       "A site, member of a `SiteFamily`.\n",
+       "\n",
+       "Sites are the vertices of the graph which describes the tight binding\n",
+       "system in a `Builder`.\n",
+       "\n",
+       "A site is uniquely identified by its family and its tag.\n",
+       "\n",
+       "Parameters\n",
+       "----------\n",
+       "family : an instance of `SiteFamily`\n",
+       "    The 'type' of the site.\n",
+       "tag : a hashable python object\n",
+       "    The unique identifier of the site within the site family, typically a\n",
+       "    vector of integers.\n",
+       "\n",
+       "Raises\n",
+       "------\n",
+       "ValueError\n",
+       "    If `tag` is not a proper tag for `family`.\n",
+       "\n",
+       "Notes\n",
+       "-----\n",
+       "For convenience, ``family(*tag)`` can be used instead of ``Site(family,\n",
+       "tag)`` to create a site.\n",
+       "\n",
+       "The parameters of the constructor (see above) are stored as instance\n",
+       "variables under the same names.  Given a site ``site``, common things to\n",
+       "query are thus ``site.family``, ``site.tag``, and ``site.pos``.\n",
+       "\u001b[0;31mFile:\u001b[0m           ~/kwant/kwant/builder.py\n",
+       "\u001b[0;31mType:\u001b[0m           type\n",
+       "\u001b[0;31mSubclasses:\u001b[0m     \n"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "kwant.builder.Site?"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "\u001b[0;31mSignature:\u001b[0m      \u001b[0mchain\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mtag\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
+       "\u001b[0;31mType:\u001b[0m           Monatomic\n",
+       "\u001b[0;31mString form:\u001b[0m    <unnamed Monatomic lattice, vectors [1.0], origin [0.0]>\n",
+       "\u001b[0;31mFile:\u001b[0m           ~/kwant/kwant/lattice.py\n",
+       "\u001b[0;31mDocstring:\u001b[0m     \n",
+       "A Bravais lattice with a single site in the basis.\n",
+       "\n",
+       "Instances of this class provide the `~kwant.builder.SiteFamily` interface.\n",
+       "Site tags (see `~kwant.builder.SiteFamily`) are sequences of integers and\n",
+       "describe the lattice coordinates of a site.\n",
+       "\n",
+       "``Monatomic`` instances are used as site families on their own or as\n",
+       "sublattices of `Polyatomic` lattices.\n",
+       "\n",
+       "Parameters\n",
+       "----------\n",
+       "prim_vecs : 2d array-like of floats\n",
+       "    Primitive vectors of the Bravais lattice.\n",
+       "\n",
+       "offset : vector of floats\n",
+       "    Displacement of the lattice origin from the real space\n",
+       "    coordinates origin.\n",
+       "\n",
+       "Attributes\n",
+       "----------\n",
+       "``offset`` : vector\n",
+       "    Displacement of the lattice origin from the real space coordinates origin\n",
+       "\u001b[0;31mCall docstring:\u001b[0m\n",
+       "A convenience function.\n",
+       "\n",
+       "This function allows to write fam(1, 2) instead of Site(fam, (1, 2)).\n"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "chain?"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "H1 = \"k_x**2 + 10*k_y**2\"\n",
+    "H2 = \"t_1*k_x**2 + t_2*k_y**2 + t_3*k_z**2\"\n",
+    "H3 = \"sigma_z*k_x**2 + 10*sigma_x*k_y**2\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ll = landau_levels.ll_hamiltonian(H3)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS0AAAAyBAMAAAD7OpAnAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEImZRO/dMlQiu6vNZnZmcXX2AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAFhElEQVRYCcWZTWhcVRTHz3uZJDOZTOZRaRZ+zduoiKWOQhcq6uDCjUJalPpBIRWCtRRkNroSMy4qKkhSEDEiJBtBLEgUpCDaDmTjwkVocZ2JXZlSM7UFU6uO99x73v18b2buE5m7eHPuuff+z++dd+fOOwwAa/vqeHXbeK8Xu17lCXu9w6rX35rv/d5/ApDap61kXuFoMzHNz/H3VyLTY/aClU+G5npx5Yy52OmRWuXlJOZS3ZkjHOPNjAHpnhya6wW5JtsQah+/ImaU/syamcl18CIRe3OFV+9LiybcQq10S8woX0ubiT7OVTziDE+3gpPC6c31I8wuW3qnIwDhJrUnxU1X29ZE2UWuOx5wsX8GOCAm+XJVDkNohgu+uRABuUltscPFq+sihnvl+aq4XGcBlthNsubLNbUORVtvKwJyk1qtzrVr3lxsA8y1+FpfLpaDwh5fqS6Mi9wJlwDy5gquM65OLi6WisI/ColbjIvc/5GrxE5I2gO++VqNYTKFi9yZXOGjD+k3k7G/SixfizGf2IfL0hLn12odJu1zn+WL3Jlc98P3Q3AN9xwtLcGV8zmehMVlDSwjX8D2/VJrUL4sLcHFNngxfd+jOzNfC7DaGILrSwCWfGx9nqOlJbim1sA5d5gUuTO55Lkk4LLyNeS5Smec0BJcFedc5bdI7j5crwsVjWvGPgcBJlqBCNMvXwCGFi14DfY39RDM3mWpF+40rtuP4vSwri/CfI19+9fjug/t4OAVEk99jsHnx3GWqUVc+347h2Nae2x+YxOEO40LNnDuJW0B/W4bHruTygXT/O5MLeKyBYx+Khd+TQpxgRLB5/P9Zay0O+lc1RabZ2l5cNEPNjtsUQa30VvvfRBhh1perjkUsbR8uNZ4/K0GfoTH2WW+10M7aXm57kUBS8uDa7rLAV7l1+l68dIyt9hlYp1bebk24MNHEik4JSwPruLfuCS8yRfWWr/ObQoJgP1tbuXkCvZmPuLSKBIc41LgwQUXMEH3NPjCpbuju4SCuubkKl7bDuQtJmo+XOXnAe78Wqw8VE8U1GdOrrHnlIS0fLjg3Q7M0pvNqZ2O1EiMnFwT579KFNSnF9ebHYBDTVwc3CjdUCJk5eSqNVYbjpYPV8ieI4y9hBqVbrAX2mI5ubag1optLR+uB3mq+DlRXoM/Qggvv30OBfXvo6xh7UjYF+d9rW2MfQfVZly4AhMRczvfx5TKNglBvx5UAew22PKpFlzehu1yewFDaOeXqmFxQGunMapQmqDKnUa/gJknorAavcH7+vmFa9zKVoYgrqkuXyh+h7jZnN4kbf6Bz1G9a+kjwMvRhAveMcZ458zTsKN72XPUS1h9SIYgLiq02Tu3bNVlaTIDuc5a74pynL+xklIsvcp4GM6rDohzFV9N3cpWhiA1Khz1+nERGIpsyKVqWOkWhsYVdKwx7C6YRzzf92wNe1mwK1sZIour1D0AsRaCcWm1jzaApsZVsYZ490Shq7sTLrciUiGyuIIfnvol0sQYl1bDagNoalzWiOg+s9PW/QkX281WZatCZHHpOmgj13VZw1qjA7iC5alNfYXkcipbFcKDSyVZBAl+4o1tvAFc5fhZHUvu+zzP0RBiHdr3Sy17APsDuAoXY2NVki+27+3KVpbJHvkCVcMaYQZyWbNlvlIqWxnCh0seenagAfmypyf5ohJWH5YhfLhUDatLMRvL0f51rbGAc+Eat7KVIXy4VA1rhAFejnpyaSWsLiZD+HDpAq5NSu6A4+H5crymw+Ri35CMht/H/u3/4arzoOPtrNgj4lqMOVBZFGopcCPi2hWPKfnfwwUbDVeQ/C+023GRuGc0XLNHCKdyK0oHG8n/fMVjkua29QyuUfwv+tkmwvwLRnzCe2DCMugAAAAASUVORK5CYII=\n",
+      "text/latex": [
+       "$$\\left \\{ k_{x}^{2} : \\left[\\begin{matrix}1 & 0\\\\0 & -1\\end{matrix}\\right], \\quad k_{y}^{2} : \\left[\\begin{matrix}0 & 10\\\\10 & 0\\end{matrix}\\right]\\right \\}$$"
+      ],
+      "text/plain": [
+       "⎧  2  ⎡1  0 ⎤     2  ⎡0   10⎤⎫\n",
+       "⎨kₓ : ⎢     ⎥, k_y : ⎢      ⎥⎬\n",
+       "⎩     ⎣0  -1⎦        ⎣10  0 ⎦⎭"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "ll.monomials"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "True"
+      ]
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "ll.normal_coordinate == None"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ll.set_ll_momenta([1, 2])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAE8AAAAVBAMAAAD1D64kAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAds0yiUTdELvvmVRmIqtiKJjQAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABQklEQVQoFYXTMUjDQBQG4D9tktM2aYODFUEoCAqdsrlmUOe4C8nmJnEQBJVEN9HBScQpjiJCNpeCHZxcFIquWujgKEjJEAdDLgc5k8NbXu5/H/cuw2FmEf8vbcnGGrB3IaDkxM07cgbRE0DUQh4mItiOOKjEIuh4HJRdEVxnDXrHWth8emERV3vo32UBhX604Aw4kG9IIr3+FODDjndQ5dCMO2RQgDdhJQMaV6xBR2/PGSzgqzo6yoMMkok+ATY+sMUzwLffbeBteYgMal8kkWHF9fO/8Bl+ZOH20aCwFeBb1nddKYDf5ewn2qYF7IPCqQjzHTSslgGVn34K6doDScfTn8lOUe1pE7jnTqSbIVY4aDppbpWhPu4HRaiMjgFilKFWv/SK8BBngFZ2NEnvOJs/hbHSFSlASp8CW5ur7Ku6/gLiIEbRKeWmjQAAAABJRU5ErkJggg==\n",
+      "text/latex": [
+       "$$\\left [ k_{y}, \\quad k_{z}\\right ]$$"
+      ],
+      "text/plain": [
+       "[k_y, k_z]"
+      ]
+     },
+     "execution_count": 16,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "ll.ll_momenta"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATcAAAAzBAMAAAAA6lKjAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEImZRO/dMlQiu6vNZnZmcXX2AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAFwUlEQVRYCbWZTYgcRRTHX/fO7szs7GSHqAt+ZQeDiBji5JCLIg65eNJNQOIHgU0gBENA5uRJ3Lko0cuOF3FF2BVvfsAaEC9qBnLJwcNg8JyJOanoTmLAaKJj1Xuvvrqre6t3mT5Mv3rV71+/qeruV28GYJsjfurQNldMuHtvKzXA8vgm+R6Hb1OdjiMej486jpzG9HjczOkWXaz2UVddVjreUaY+v7z2PtlnYaWnvT4jWvswHO7dtYZPw/hYrf6qum61ZTqV9ZIyTsN6W9kZ53I4XCdDwnKT2genyFX9x+pSpoR74Di2VtWXUH3yHP/+mG4WhrODtYrSJLXqHeqo3bAuUCbO3CVsva589vkyLPRUuzCcHcwib4sZIDerPUuTPN9Xo1hnhLstHXHLciuzfhTivmoUhXOCUSS6cLEB7Ga1lSH2zG+qUayzhCvhlF6xvNqc3YSKnvCicE4wS15tALtZbZHmZDELLj4pIkvNUkczaUN8oRLOq/QUhXOCWVLAsVvBEVUm3FyrcqX35vn3GppJG+J7lf5TraJwTrCBY3cg3GL3l6XB8nisIKzzehPKO4Zzgg0cuwPhVvc1HraAbHO9BWVOIbnLmkh+03iHOMEWHGkGwh1u2TyO7axMzrImkh/BOcEGjt2BcOd+HTpEVkPcvJWQByKR/AjOCTZw7A6Di25Vb1k8jjm7AfWQV0ki+RGcE2zg2B0GVx9Ft2MHyTT4hUmOnGUFWLWfdYJzgg0cu8PgahvwVwzx9be+MVDaeg3uw5tbOnLhnORHcGAHs+KW+A7kDoOb7cL1a3Ct1j+tkYyx9w+D7IejbYOb/BjODibJp5cvDYDcPrjKMTOysDC3Sk9nbuB0eBp+OMBtg5v8GM6jYVweuAef0Pc3XqfhYL5nAv1WBpx8mhPJrwCceHblId7M4jAPn2yZmYMV6KAn+8MPh9uGRPIrAreBA15ty1MGXHV0AJqyP+fww+G2IZH8CsDNjXDEM/ip4M4RhVrW6LsjPzdywGSXHw63DSqSVQvAVf6VsfHfqMBw0QnSU3DUyv/0w+G2gQOVagE4uNgTwY+2UUHNHMvtHs63bSgCV3sR4KGvCGdncPs/FWvunznftqEIHLwzhIVDu4CLe7V+Fpxv21AI7o0hwOEO0u1o5mZ6ZXHHemfOu20oAheLZYWpV3YB15YbFy+cd9tQBO5JnDR6leyhDKGeK34gDv5I84r8zgdX3LIGI7jortNP24YXNuGUcCtVA5cWVpU2f1UuUbbaInzq67vPoLjznpvrRmfR6/mgintmqODgjHg2zEHbhuYNGst+z1WOAXiEVaXNcLMj1KL0ZWTRopn7CeBAokM3aSO8X7RZruLA4XXVe0/u2dARAHLmMIunhfUmj9W41Bd79/RBcJ8nNovWdVRx12Usy/k2plPN2tAKwmWVj15aeFaV6azGBWtm3Qogfk1Z6lrilompE47Acxpun9WrzJm2uc2ET8OlhcUuhMr0ULjoTwE3VAO5Z0yd5cvnP9FwA/cCbM10lmyvgvMI65IsFK4qStOVoa1ubEyd4pfKkYYzfcYqff+FaZiZ8whzSa3Vtl3Wqpi5laatbmwrdfJ3NX3Guh++NA0LLi28rsr00JnzzH70Ax4dsFJnDtyFUt8H5xEuvKzygVjt2urGtlJnDtzzj5gIYal7ziPMJXX4ssJnAOK3Kd9hp84cuESohksL60o7dFkh/a5Uo9mpcwdwaeHkS3jbBwJmulHGrpNSJ1XcBeFkFvcIq0o7eOaig7911Fy5Z7viLgaHWdwjrCrtYDgXKN2iirsYXFrF9bhwXL66l2SspnsRUMU9GbgWjjXdTw4p2oFwWHFPBI5f/jWqDV3CIDiuuCcCt0U3u/onx6ELguOKexJwkfrPa2vocGEjCI7DJgG3ILbLeNTvpLOA/kuTr8k+TeQvzcoJjXTPZmps/WdwqifpmMifwR8P5DD/AzdT3d3vz5R8AAAAAElFTkSuQmCC\n",
+      "text/latex": [
+       "$$\\left \\{ 1 : \\left[\\begin{matrix}k_{x}^{2} & 0\\\\0 & - k_{x}^{2}\\end{matrix}\\right], \\quad k_{y}^{2} : \\left[\\begin{matrix}0 & 10\\\\10 & 0\\end{matrix}\\right]\\right \\}$$"
+      ],
+      "text/plain": [
+       "⎧   ⎡  2      ⎤                ⎫\n",
+       "⎪   ⎢kₓ    0  ⎥     2  ⎡0   10⎤⎪\n",
+       "⎨1: ⎢         ⎥, k_y : ⎢      ⎥⎬\n",
+       "⎪   ⎢        2⎥        ⎣10  0 ⎦⎪\n",
+       "⎩   ⎣ 0   -kₓ ⎦                ⎭"
+      ]
+     },
+     "execution_count": 17,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "ll.monomials"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 18,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAAsAAAAJBAMAAAAWSsseAAAALVBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADAOrOgAAAADnRSTlMAEHarIkSJZt3NuzJUmW693xMAAAAJcEhZcwAADsQAAA7EAZUrDhsAAABASURBVAgdY2AQUnZVU2BgTGBv4pjAwCbA9pDVgYGRgWsBAwjwKYCpfRuAFI+AHgOTAEPcgXUM7gwMwkpC1wsYABfXCcn8wW65AAAAAElFTkSuQmCC\n",
+      "text/latex": [
+       "$$x$$"
+      ],
+      "text/plain": [
+       "x"
+      ]
+     },
+     "execution_count": 18,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "ll.normal_coordinate"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 19,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# ll.set_ll_momenta(['k_x', 'k_y'])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 20,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAAsAAAAJBAMAAAAWSsseAAAALVBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADAOrOgAAAADnRSTlMAEHarIkSJZt3NuzJUmW693xMAAAAJcEhZcwAADsQAAA7EAZUrDhsAAABASURBVAgdY2AQUnZVU2BgTGBv4pjAwCbA9pDVgYGRgWsBAwjwKYCpfRuAFI+AHgOTAEPcgXUM7gwMwkpC1wsYABfXCcn8wW65AAAAAElFTkSuQmCC\n",
+      "text/latex": [
+       "$$x$$"
+      ],
+      "text/plain": [
+       "x"
+      ]
+     },
+     "execution_count": 20,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "ll.normal_coordinate"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 21,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATcAAAAzBAMAAAAA6lKjAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEImZRO/dMlQiu6vNZnZmcXX2AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAFwUlEQVRYCbWZTYgcRRTHX/fO7szs7GSHqAt+ZQeDiBji5JCLIg65eNJNQOIHgU0gBENA5uRJ3Lko0cuOF3FF2BVvfsAaEC9qBnLJwcNg8JyJOanoTmLAaKJj1Xuvvrqre6t3mT5Mv3rV71+/qeruV28GYJsjfurQNldMuHtvKzXA8vgm+R6Hb1OdjiMej486jpzG9HjczOkWXaz2UVddVjreUaY+v7z2PtlnYaWnvT4jWvswHO7dtYZPw/hYrf6qum61ZTqV9ZIyTsN6W9kZ53I4XCdDwnKT2genyFX9x+pSpoR74Di2VtWXUH3yHP/+mG4WhrODtYrSJLXqHeqo3bAuUCbO3CVsva589vkyLPRUuzCcHcwib4sZIDerPUuTPN9Xo1hnhLstHXHLciuzfhTivmoUhXOCUSS6cLEB7Ga1lSH2zG+qUayzhCvhlF6xvNqc3YSKnvCicE4wS15tALtZbZHmZDELLj4pIkvNUkczaUN8oRLOq/QUhXOCWVLAsVvBEVUm3FyrcqX35vn3GppJG+J7lf5TraJwTrCBY3cg3GL3l6XB8nisIKzzehPKO4Zzgg0cuwPhVvc1HraAbHO9BWVOIbnLmkh+03iHOMEWHGkGwh1u2TyO7axMzrImkh/BOcEGjt2BcOd+HTpEVkPcvJWQByKR/AjOCTZw7A6Di25Vb1k8jjm7AfWQV0ki+RGcE2zg2B0GVx9Ft2MHyTT4hUmOnGUFWLWfdYJzgg0cu8PgahvwVwzx9be+MVDaeg3uw5tbOnLhnORHcGAHs+KW+A7kDoOb7cL1a3Ct1j+tkYyx9w+D7IejbYOb/BjODibJp5cvDYDcPrjKMTOysDC3Sk9nbuB0eBp+OMBtg5v8GM6jYVweuAef0Pc3XqfhYL5nAv1WBpx8mhPJrwCceHblId7M4jAPn2yZmYMV6KAn+8MPh9uGRPIrAreBA15ty1MGXHV0AJqyP+fww+G2IZH8CsDNjXDEM/ip4M4RhVrW6LsjPzdywGSXHw63DSqSVQvAVf6VsfHfqMBw0QnSU3DUyv/0w+G2gQOVagE4uNgTwY+2UUHNHMvtHs63bSgCV3sR4KGvCGdncPs/FWvunznftqEIHLwzhIVDu4CLe7V+Fpxv21AI7o0hwOEO0u1o5mZ6ZXHHemfOu20oAheLZYWpV3YB15YbFy+cd9tQBO5JnDR6leyhDKGeK34gDv5I84r8zgdX3LIGI7jortNP24YXNuGUcCtVA5cWVpU2f1UuUbbaInzq67vPoLjznpvrRmfR6/mgintmqODgjHg2zEHbhuYNGst+z1WOAXiEVaXNcLMj1KL0ZWTRopn7CeBAokM3aSO8X7RZruLA4XXVe0/u2dARAHLmMIunhfUmj9W41Bd79/RBcJ8nNovWdVRx12Usy/k2plPN2tAKwmWVj15aeFaV6azGBWtm3Qogfk1Z6lrilompE47Acxpun9WrzJm2uc2ET8OlhcUuhMr0ULjoTwE3VAO5Z0yd5cvnP9FwA/cCbM10lmyvgvMI65IsFK4qStOVoa1ubEyd4pfKkYYzfcYqff+FaZiZ8whzSa3Vtl3Wqpi5laatbmwrdfJ3NX3Guh++NA0LLi28rsr00JnzzH70Ax4dsFJnDtyFUt8H5xEuvKzygVjt2urGtlJnDtzzj5gIYal7ziPMJXX4ssJnAOK3Kd9hp84cuESohksL60o7dFkh/a5Uo9mpcwdwaeHkS3jbBwJmulHGrpNSJ1XcBeFkFvcIq0o7eOaig7911Fy5Z7viLgaHWdwjrCrtYDgXKN2iirsYXFrF9bhwXL66l2SspnsRUMU9GbgWjjXdTw4p2oFwWHFPBI5f/jWqDV3CIDiuuCcCt0U3u/onx6ELguOKexJwkfrPa2vocGEjCI7DJgG3ILbLeNTvpLOA/kuTr8k+TeQvzcoJjXTPZmps/WdwqifpmMifwR8P5DD/AzdT3d3vz5R8AAAAAElFTkSuQmCC\n",
+      "text/latex": [
+       "$$\\left \\{ 1 : \\left[\\begin{matrix}k_{x}^{2} & 0\\\\0 & - k_{x}^{2}\\end{matrix}\\right], \\quad k_{y}^{2} : \\left[\\begin{matrix}0 & 10\\\\10 & 0\\end{matrix}\\right]\\right \\}$$"
+      ],
+      "text/plain": [
+       "⎧   ⎡  2      ⎤                ⎫\n",
+       "⎪   ⎢kₓ    0  ⎥     2  ⎡0   10⎤⎪\n",
+       "⎨1: ⎢         ⎥, k_y : ⎢      ⎥⎬\n",
+       "⎪   ⎢        2⎥        ⎣10  0 ⎦⎪\n",
+       "⎩   ⎣ 0   -kₓ ⎦                ⎭"
+      ]
+     },
+     "execution_count": 21,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "ll.monomials"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 22,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ll.shape_func"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 23,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "N = 3"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 24,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ham = ll.hamiltonian_matrix(N)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 25,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATcAAAAzBAMAAAAA6lKjAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEImZRO/dMlQiu6vNZnZmcXX2AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAFwUlEQVRYCbWZTYgcRRTHX/fO7szs7GSHqAt+ZQeDiBji5JCLIg65eNJNQOIHgU0gBENA5uRJ3Lko0cuOF3FF2BVvfsAaEC9qBnLJwcNg8JyJOanoTmLAaKJj1Xuvvrqre6t3mT5Mv3rV71+/qeruV28GYJsjfurQNldMuHtvKzXA8vgm+R6Hb1OdjiMej486jpzG9HjczOkWXaz2UVddVjreUaY+v7z2PtlnYaWnvT4jWvswHO7dtYZPw/hYrf6qum61ZTqV9ZIyTsN6W9kZ53I4XCdDwnKT2genyFX9x+pSpoR74Di2VtWXUH3yHP/+mG4WhrODtYrSJLXqHeqo3bAuUCbO3CVsva589vkyLPRUuzCcHcwib4sZIDerPUuTPN9Xo1hnhLstHXHLciuzfhTivmoUhXOCUSS6cLEB7Ga1lSH2zG+qUayzhCvhlF6xvNqc3YSKnvCicE4wS15tALtZbZHmZDELLj4pIkvNUkczaUN8oRLOq/QUhXOCWVLAsVvBEVUm3FyrcqX35vn3GppJG+J7lf5TraJwTrCBY3cg3GL3l6XB8nisIKzzehPKO4Zzgg0cuwPhVvc1HraAbHO9BWVOIbnLmkh+03iHOMEWHGkGwh1u2TyO7axMzrImkh/BOcEGjt2BcOd+HTpEVkPcvJWQByKR/AjOCTZw7A6Di25Vb1k8jjm7AfWQV0ki+RGcE2zg2B0GVx9Ft2MHyTT4hUmOnGUFWLWfdYJzgg0cu8PgahvwVwzx9be+MVDaeg3uw5tbOnLhnORHcGAHs+KW+A7kDoOb7cL1a3Ct1j+tkYyx9w+D7IejbYOb/BjODibJp5cvDYDcPrjKMTOysDC3Sk9nbuB0eBp+OMBtg5v8GM6jYVweuAef0Pc3XqfhYL5nAv1WBpx8mhPJrwCceHblId7M4jAPn2yZmYMV6KAn+8MPh9uGRPIrAreBA15ty1MGXHV0AJqyP+fww+G2IZH8CsDNjXDEM/ip4M4RhVrW6LsjPzdywGSXHw63DSqSVQvAVf6VsfHfqMBw0QnSU3DUyv/0w+G2gQOVagE4uNgTwY+2UUHNHMvtHs63bSgCV3sR4KGvCGdncPs/FWvunznftqEIHLwzhIVDu4CLe7V+Fpxv21AI7o0hwOEO0u1o5mZ6ZXHHemfOu20oAheLZYWpV3YB15YbFy+cd9tQBO5JnDR6leyhDKGeK34gDv5I84r8zgdX3LIGI7jortNP24YXNuGUcCtVA5cWVpU2f1UuUbbaInzq67vPoLjznpvrRmfR6/mgintmqODgjHg2zEHbhuYNGst+z1WOAXiEVaXNcLMj1KL0ZWTRopn7CeBAokM3aSO8X7RZruLA4XXVe0/u2dARAHLmMIunhfUmj9W41Bd79/RBcJ8nNovWdVRx12Usy/k2plPN2tAKwmWVj15aeFaV6azGBWtm3Qogfk1Z6lrilompE47Acxpun9WrzJm2uc2ET8OlhcUuhMr0ULjoTwE3VAO5Z0yd5cvnP9FwA/cCbM10lmyvgvMI65IsFK4qStOVoa1ubEyd4pfKkYYzfcYqff+FaZiZ8whzSa3Vtl3Wqpi5laatbmwrdfJ3NX3Guh++NA0LLi28rsr00JnzzH70Ax4dsFJnDtyFUt8H5xEuvKzygVjt2urGtlJnDtzzj5gIYal7ziPMJXX4ssJnAOK3Kd9hp84cuESohksL60o7dFkh/a5Uo9mpcwdwaeHkS3jbBwJmulHGrpNSJ1XcBeFkFvcIq0o7eOaig7911Fy5Z7viLgaHWdwjrCrtYDgXKN2iirsYXFrF9bhwXL66l2SspnsRUMU9GbgWjjXdTw4p2oFwWHFPBI5f/jWqDV3CIDiuuCcCt0U3u/onx6ELguOKexJwkfrPa2vocGEjCI7DJgG3ILbLeNTvpLOA/kuTr8k+TeQvzcoJjXTPZmps/WdwqifpmMifwR8P5DD/AzdT3d3vz5R8AAAAAElFTkSuQmCC\n",
+      "text/latex": [
+       "$$\\left \\{ 1 : \\left[\\begin{matrix}k_{x}^{2} & 0\\\\0 & - k_{x}^{2}\\end{matrix}\\right], \\quad k_{y}^{2} : \\left[\\begin{matrix}0 & 10\\\\10 & 0\\end{matrix}\\right]\\right \\}$$"
+      ],
+      "text/plain": [
+       "⎧   ⎡  2      ⎤                ⎫\n",
+       "⎪   ⎢kₓ    0  ⎥     2  ⎡0   10⎤⎪\n",
+       "⎨1: ⎢         ⎥, k_y : ⎢      ⎥⎬\n",
+       "⎪   ⎢        2⎥        ⎣10  0 ⎦⎪\n",
+       "⎩   ⎣ 0   -kₓ ⎦                ⎭"
+      ]
+     },
+     "execution_count": 25,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "ll.monomials"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 28,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "params = dict(t_1 = 1, t_2 = 2, t_3 = 3, k_x=0)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 29,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "matrix([[ 0.        +0.j,  5.        +0.j,  0.        +0.j,\n",
+       "          0.        +0.j,  0.        +0.j,  7.07106781+0.j],\n",
+       "        [ 5.        +0.j,  0.        +0.j,  0.        +0.j,\n",
+       "          0.        +0.j,  7.07106781+0.j,  0.        +0.j],\n",
+       "        [ 0.        +0.j,  0.        +0.j,  0.        +0.j,\n",
+       "         15.        +0.j,  0.        +0.j,  0.        +0.j],\n",
+       "        [ 0.        +0.j,  0.        +0.j, 15.        +0.j,\n",
+       "          0.        +0.j,  0.        +0.j,  0.        +0.j],\n",
+       "        [ 0.        +0.j,  7.07106781+0.j,  0.        +0.j,\n",
+       "          0.        +0.j,  0.        +0.j, 25.        +0.j],\n",
+       "        [ 7.07106781+0.j,  0.        +0.j,  0.        +0.j,\n",
+       "          0.        +0.j, 25.        +0.j,  0.        +0.j]])"
+      ]
+     },
+     "execution_count": 29,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "ham(B=1, params=params).todense()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "H1 = \"k_x**2 + 10*k_y**2\"\n",
+    "H2 = \"t_1*k_x**2 + t_2*k_y**2 + t_3*k_z**2\"\n",
+    "H3 = \"sigma_z*k_x**2 + 10*sigma_x*k_y**2\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "N = 2\n",
+    "ham, monomials = ll_hamiltonian(H3, N, momenta=None)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<function kwant.continuum.landau_levels.combine.<locals>.final_wrapper(*, B, params=None)>"
+      ]
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "ham"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "matrix([[ 0.5+0.j,  5. +0.j,  0. +0.j,  0. +0.j],\n",
+       "        [ 5. +0.j, -0.5+0.j,  0. +0.j,  0. +0.j],\n",
+       "        [ 0. +0.j,  0. +0.j,  1.5+0.j, 15. +0.j],\n",
+       "        [ 0. +0.j,  0. +0.j, 15. +0.j, -1.5+0.j]])"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "ham(B=1, params=dict(t_1=1, t_2=1, t_3=1)).todense()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{k_x**2: <function __main__.ll_hamiltonian.<locals>.<lambda>.<locals>.<lambda>(params=None)>,\n",
+       " k_y**2: <function __main__.ll_hamiltonian.<locals>.<lambda>.<locals>.<lambda>(params=None)>}"
+      ]
+     },
+     "execution_count": 9,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "monomials"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "funcs = list(monomials.values())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<function __main__.ll_hamiltonian.<locals>.<lambda>.<locals>.<lambda>(params=None)>"
+      ]
+     },
+     "execution_count": 11,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "funcs[0]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# funcs[0]().todense()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "k_x**2 [[ 1.+0.j  0.+0.j]\n",
+      " [ 0.+0.j -1.+0.j]]\n",
+      "k_y**2 [[ 0.+0.j 10.+0.j]\n",
+      " [10.+0.j  0.+0.j]]\n"
+     ]
+    }
+   ],
+   "source": [
+    "params = dict(t_1 = 1, t_2=1, t_3=1)\n",
+    "for key, value in monomials.items():\n",
+    "    print(key, value(params=params).todense())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 60,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "f1 = lambda x, y=1, z=2: x*y + z\n",
+    "f2 = lambda l, y=1, z=2: l*y + z"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 61,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "funcs = [f1, f2]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 62,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<function __main__.combine.<locals>.final_wrapper(*, l, x, y=1, z=2)>"
+      ]
+     },
+     "execution_count": 62,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "combine(operator.add, funcs, args_names=None)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "False"
+      ]
+     },
+     "execution_count": 15,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "('N_1', None) == ('N_2', None)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 81,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "H = \"k_x**2 + k_y**2 + k_z**2\""
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 87,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "disc = kwant.continuum.discretize(H, coords=['y', 'z'])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 88,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "# Discrete coordinates: y z\n",
+       "\n",
+       "# Onsite element:\n",
+       "def onsite(site, k_x):\n",
+       "    return (4.0 + k_x**2)\n",
+       "\n",
+       "# Hopping from (1, 0):\n",
+       "(-1+0j)\n",
+       "\n",
+       "# Hopping from (0, 1):\n",
+       "(-1+0j)"
+      ],
+      "text/plain": [
+       "<kwant.continuum.discretizer._DiscretizedBuilder at 0x7f8e4a178630>"
+      ]
+     },
+     "execution_count": 88,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "disc"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 89,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[0.0 0.0]\n"
+     ]
+    }
+   ],
+   "source": [
+    "for site in disc.sites():\n",
+    "    print(site.pos)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 71,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fsyst = disc.finalized()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 74,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[ 2.+0.j, -1.-0.j],\n",
+       "       [-1.+0.j,  2.+0.j]])"
+      ]
+     },
+     "execution_count": 74,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "fsyst.hamiltonian_submatrix()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "l = ['a', 'b', 'c', 'd']\n",
+    "inds = [0, 3]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 22,
+   "metadata": {},
+   "outputs": [
+    {
+     "ename": "TypeError",
+     "evalue": "zip argument #1 must support iteration",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
+      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
+      "\u001b[0;32m<ipython-input-22-39657fc954cd>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0ml\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mzip\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0minds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
+      "\u001b[0;31mTypeError\u001b[0m: zip argument #1 must support iteration"
+     ]
+    }
+   ],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "work_env",
+   "language": "python",
+   "name": "work_env"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.1"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}