doc/source/code/figure/quantum_well.py.diff

-@@ -1,88 +1,95 @@
- # Tutorial 2.3.2. Spatially dependent values through functions
- # ============================================================
- #
- # Physics background
- # ------------------
- #  transmission through a quantum well
- #
- # Kwant features highlighted
- # --------------------------
- #  - Functions as values in Builder
- 
-+import _defs
- import kwant
- 
- # For plotting
- from matplotlib import pyplot
- 
- 
- #HIDDEN_BEGIN_ehso
- def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
-     # Start with an empty tight-binding system and a single square lattice.
-     # `a` is the lattice constant (by default set to 1 for simplicity.
-     lat = kwant.lattice.square(a)
- 
-     syst = kwant.Builder()
- 
-     #### Define the scattering region. ####
-     # Potential profile
-     def potential(site, pot):
-         (x, y) = site.pos
-         if (L - L_well) / 2 < x < (L + L_well) / 2:
-             return pot
-         else:
-             return 0
- #HIDDEN_END_ehso
- 
- #HIDDEN_BEGIN_coid
-     def onsite(site, pot):
-         return 4 * t + potential(site, pot)
- 
-     syst[(lat(x, y) for x in range(L) for y in range(W))] = onsite
-     syst[lat.neighbors()] = -t
- #HIDDEN_END_coid
- 
-     #### Define and attach the leads. ####
-     lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
-     lead[(lat(0, j) for j in range(W))] = 4 * t
-     lead[lat.neighbors()] = -t
-     syst.attach_lead(lead)
-     syst.attach_lead(lead.reversed())
- 
-     return syst
- 
- 
- def plot_conductance(syst, energy, welldepths):
- #HIDDEN_BEGIN_sqvr
- 
-     # Compute conductance
-     data = []
-     for welldepth in welldepths:
-         smatrix = kwant.smatrix(syst, energy, params=dict(pot=-welldepth))
-         data.append(smatrix.transmission(1, 0))
- 
--    pyplot.figure()
-+    fig = pyplot.figure()
-     pyplot.plot(welldepths, data)
--    pyplot.xlabel("well depth [t]")
--    pyplot.ylabel("conductance [e^2/h]")
--    pyplot.show()
-+    pyplot.xlabel("well depth [t]",
-+                 fontsize=_defs.mpl_label_size)
-+    pyplot.ylabel("conductance [e^2/h]",
-+                 fontsize=_defs.mpl_label_size)
-+    pyplot.setp(fig.get_axes()[0].get_xticklabels(),
-+               fontsize=_defs.mpl_tick_size)
-+    pyplot.setp(fig.get_axes()[0].get_yticklabels(),
-+               fontsize=_defs.mpl_tick_size)
-+    fig.set_size_inches(_defs.mpl_width_in, _defs.mpl_width_in * 3. / 4.)
-+    fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-+    for extension in ('pdf', 'png'):
-+        fig.savefig("quantum_well_result." + extension, dpi=_defs.dpi)
- #HIDDEN_END_sqvr
- 
- 
- def main():
-     syst = make_system()
- 
--    # Check that the system looks as intended.
--    kwant.plot(syst)
--
-     # Finalize the system.
-     syst = syst.finalized()
- 
-     # We should see conductance steps.
-     plot_conductance(syst, energy=0.2,
-                      welldepths=[0.01 * i for i in range(100)])
- 
- 
- # Call the main function if the script gets executed (as opposed to imported).
- # See <http://docs.python.org/library/__main__.html>.
- if __name__ == '__main__':
-     main()

doc/source/code/figure/quantum_wire.py.diff

-@@ -1,121 +1,130 @@
- # Tutorial 2.2.2. Transport through a quantum wire
- # ================================================
- #
- # Physics background
- # ------------------
- #  Conductance of a quantum wire; subbands
- #
- # Kwant features highlighted
- # --------------------------
- #  - Builder for setting up transport systems easily
- #  - Making scattering region and leads
- #  - Using the simple sparse solver for computing Landauer conductance
- 
-+import _defs
- from matplotlib import pyplot
- #HIDDEN_BEGIN_dwhx
- import kwant
- #HIDDEN_END_dwhx
- 
- # First, define the tight-binding system
- 
- #HIDDEN_BEGIN_goiq
- syst = kwant.Builder()
- #HIDDEN_END_goiq
- 
- # Here, we are only working with square lattices
- #HIDDEN_BEGIN_suwo
- a = 1
- lat = kwant.lattice.square(a)
- #HIDDEN_END_suwo
- 
- #HIDDEN_BEGIN_zfvr
- t = 1.0
- W = 10
- L = 30
- 
- # Define the scattering region
- 
- for i in range(L):
-     for j in range(W):
-         # On-site Hamiltonian
-         syst[lat(i, j)] = 4 * t
- 
-         # Hopping in y-direction
-         if j > 0:
-             syst[lat(i, j), lat(i, j - 1)] = -t
- 
-         # Hopping in x-direction
-         if i > 0:
-             syst[lat(i, j), lat(i - 1, j)] = -t
- #HIDDEN_END_zfvr
- 
- # Then, define and attach the leads:
- 
- # First the lead to the left
- # (Note: TranslationalSymmetry takes a real-space vector)
- #HIDDEN_BEGIN_xcmc
- sym_left_lead = kwant.TranslationalSymmetry((-a, 0))
- left_lead = kwant.Builder(sym_left_lead)
- #HIDDEN_END_xcmc
- 
- #HIDDEN_BEGIN_ndez
- for j in range(W):
-     left_lead[lat(0, j)] = 4 * t
-     if j > 0:
-         left_lead[lat(0, j), lat(0, j - 1)] = -t
-     left_lead[lat(1, j), lat(0, j)] = -t
- #HIDDEN_END_ndez
- 
- #HIDDEN_BEGIN_fskr
- syst.attach_lead(left_lead)
- #HIDDEN_END_fskr
- 
- # Then the lead to the right
- #HIDDEN_BEGIN_xhqc
- sym_right_lead = kwant.TranslationalSymmetry((a, 0))
- right_lead = kwant.Builder(sym_right_lead)
- 
- for j in range(W):
-     right_lead[lat(0, j)] = 4 * t
-     if j > 0:
-         right_lead[lat(0, j), lat(0, j - 1)] = -t
-     right_lead[lat(1, j), lat(0, j)] = -t
- 
- syst.attach_lead(right_lead)
- #HIDDEN_END_xhqc
- 
- # Plot it, to make sure it's OK
- #HIDDEN_BEGIN_wsgh
--kwant.plot(syst)
-+size = (_defs.figwidth_in, 0.3 * _defs.figwidth_in)
-+for extension in ('pdf', 'png'):
-+    kwant.plot(syst, file="quantum_wire_syst." + extension,
-+               fig_size=size, dpi=_defs.dpi)
- #HIDDEN_END_wsgh
- 
- # Finalize the system
- #HIDDEN_BEGIN_dngj
- syst = syst.finalized()
- #HIDDEN_END_dngj
- 
- # Now that we have the system, we can compute conductance
- #HIDDEN_BEGIN_buzn
- energies = []
- data = []
- for ie in range(100):
-     energy = ie * 0.01
- 
-     # compute the scattering matrix at a given energy
-     smatrix = kwant.smatrix(syst, energy)
- 
-     # compute the transmission probability from lead 0 to
-     # lead 1
-     energies.append(energy)
-     data.append(smatrix.transmission(1, 0))
- #HIDDEN_END_buzn
- 
- # Use matplotlib to write output
- # We should see conductance steps
- #HIDDEN_BEGIN_lliv
--pyplot.figure()
-+fig = pyplot.figure()
- pyplot.plot(energies, data)
--pyplot.xlabel("energy [t]")
--pyplot.ylabel("conductance [e^2/h]")
--pyplot.show()
-+pyplot.xlabel("energy [t]", fontsize=_defs.mpl_label_size)
-+pyplot.ylabel("conductance [e^2/h]", fontsize=_defs.mpl_label_size)
-+pyplot.setp(fig.get_axes()[0].get_xticklabels(), fontsize=_defs.mpl_tick_size)
-+pyplot.setp(fig.get_axes()[0].get_yticklabels(), fontsize=_defs.mpl_tick_size)
-+fig.set_size_inches(_defs.mpl_width_in, _defs.mpl_width_in * 3. / 4.)
-+fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-+for extension in ('pdf', 'png'):
-+    fig.savefig("quantum_wire_result." + extension, dpi=_defs.dpi)
- #HIDDEN_END_lliv

doc/source/code/figure/spin_orbit.py.diff

-@@ -1,104 +1,110 @@
- # Tutorial 2.3.1. Matrix structure of on-site and hopping elements
- # ================================================================
- #
- # Physics background
- # ------------------
- #  Gaps in quantum wires with spin-orbit coupling and Zeeman splititng,
- #  as theoretically predicted in
- #   http://prl.aps.org/abstract/PRL/v90/i25/e256601
- #  and (supposedly) experimentally oberved in
- #   http://www.nature.com/nphys/journal/v6/n5/abs/nphys1626.html
- #
- # Kwant features highlighted
- # --------------------------
- #  - Numpy matrices as values in Builder
- 
-+import _defs
- import kwant
- 
- # For plotting
- from matplotlib import pyplot
- 
- # For matrix support
- #HIDDEN_BEGIN_xumz
- import tinyarray
- #HIDDEN_END_xumz
- 
- # define Pauli-matrices for convenience
- #HIDDEN_BEGIN_hwbt
- sigma_0 = tinyarray.array([[1, 0], [0, 1]])
- sigma_x = tinyarray.array([[0, 1], [1, 0]])
- sigma_y = tinyarray.array([[0, -1j], [1j, 0]])
- sigma_z = tinyarray.array([[1, 0], [0, -1]])
- #HIDDEN_END_hwbt
- 
- 
- def make_system(t=1.0, alpha=0.5, e_z=0.08, W=10, L=30):
-     # Start with an empty tight-binding system and a single square lattice.
-     # `a` is the lattice constant (by default set to 1 for simplicity).
-     lat = kwant.lattice.square()
- 
-     syst = kwant.Builder()
- 
-     #### Define the scattering region. ####
- #HIDDEN_BEGIN_uxrm
-     syst[(lat(x, y) for x in range(L) for y in range(W))] = \
-         4 * t * sigma_0 + e_z * sigma_z
-     # hoppings in x-direction
-     syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = \
-         -t * sigma_0 + 1j * alpha * sigma_y / 2
-     # hoppings in y-directions
-     syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
-         -t * sigma_0 - 1j * alpha * sigma_x / 2
- #HIDDEN_END_uxrm
- 
-     #### Define the left lead. ####
-     lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
- 
- #HIDDEN_BEGIN_yliu
-     lead[(lat(0, j) for j in range(W))] = 4 * t * sigma_0 + e_z * sigma_z
-     # hoppings in x-direction
-     lead[kwant.builder.HoppingKind((1, 0), lat, lat)] = \
-         -t * sigma_0 + 1j * alpha * sigma_y / 2
-     # hoppings in y-directions
-     lead[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
-         -t * sigma_0 - 1j * alpha * sigma_x / 2
- #HIDDEN_END_yliu
- 
-     #### Attach the leads and return the finalized system. ####
-     syst.attach_lead(lead)
-     syst.attach_lead(lead.reversed())
- 
-     return syst
- 
- 
- def plot_conductance(syst, energies):
-     # Compute conductance
-     data = []
-     for energy in energies:
-         smatrix = kwant.smatrix(syst, energy)
-         data.append(smatrix.transmission(1, 0))
- 
--    pyplot.figure()
-+    fig = pyplot.figure()
-     pyplot.plot(energies, data)
--    pyplot.xlabel("energy [t]")
--    pyplot.ylabel("conductance [e^2/h]")
--    pyplot.show()
-+    pyplot.xlabel("energy [t]", fontsize=_defs.mpl_label_size)
-+    pyplot.ylabel("conductance [e^2/h]",
-+                 fontsize=_defs.mpl_label_size)
-+    pyplot.setp(fig.get_axes()[0].get_xticklabels(),
-+               fontsize=_defs.mpl_tick_size)
-+    pyplot.setp(fig.get_axes()[0].get_yticklabels(),
-+               fontsize=_defs.mpl_tick_size)
-+    fig.set_size_inches(_defs.mpl_width_in, _defs.mpl_width_in * 3. / 4.)
-+    fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-+    for extension in ('pdf', 'png'):
-+        fig.savefig("spin_orbit_result." + extension, dpi=_defs.dpi)
- 
- 
- def main():
-     syst = make_system()
- 
--    # Check that the system looks as intended.
--    kwant.plot(syst)
--
-     # Finalize the system.
-     syst = syst.finalized()
- 
-     # We should see non-monotonic conductance steps.
-     plot_conductance(syst, energies=[0.01 * i - 0.3 for i in range(100)])
- 
- 
- # Call the main function if the script gets executed (as opposed to imported).
- # See <http://docs.python.org/library/__main__.html>.
- if __name__ == '__main__':
-     main()

doc/source/code/figure/superconductor.py.diff

-@@ -1,132 +1,141 @@
- # Tutorial 2.6. "Superconductors": orbitals, conservation laws and symmetries
- # ===========================================================================
- #
- # Physics background
- # ------------------
- # - conductance of a NS-junction (Andreev reflection, superconducting gap)
- #
- # Kwant features highlighted
- # --------------------------
- # - Implementing electron and hole ("orbital") degrees of freedom
- #   using conservation laws.
- # - Use of discrete symmetries to relate scattering states.
- 
-+import _defs
- import kwant
- 
- import tinyarray
- import numpy as np
- 
- # For plotting
- from matplotlib import pyplot
-+from contextlib import redirect_stdout
- 
- tau_x = tinyarray.array([[0, 1], [1, 0]])
- tau_y = tinyarray.array([[0, -1j], [1j, 0]])
- tau_z = tinyarray.array([[1, 0], [0, -1]])
- 
- #HIDDEN_BEGIN_nbvn
- def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
-                 mu=0.4, Delta=0.1, Deltapos=4, t=1.0):
-     # Start with an empty tight-binding system. On each site, there
-     # are now electron and hole orbitals, so we must specify the
-     # number of orbitals per site. The orbital structure is the same
-     # as in the Hamiltonian.
-     lat = kwant.lattice.square(norbs=2)
-     syst = kwant.Builder()
- 
-     #### Define the scattering region. ####
-     # The superconducting order parameter couples electron and hole orbitals
-     # on each site, and hence enters as an onsite potential.
-     # The pairing is only included beyond the point 'Deltapos' in the scattering region.
-     syst[(lat(x, y) for x in range(Deltapos) for y in range(W))] = (4 * t - mu) * tau_z
-     syst[(lat(x, y) for x in range(Deltapos, L) for y in range(W))] = (4 * t - mu) * tau_z + Delta * tau_x
- 
-     # The tunnel barrier
-     syst[(lat(x, y) for x in range(barrierpos[0], barrierpos[1])
-          for y in range(W))] = (4 * t + barrier - mu) * tau_z
- 
-     # Hoppings
-     syst[lat.neighbors()] = -t * tau_z
- #HIDDEN_END_nbvn
- #HIDDEN_BEGIN_ttth
-     #### Define the leads. ####
-     # Left lead - normal, so the order parameter is zero.
-     sym_left = kwant.TranslationalSymmetry((-a, 0))
-     # Specify the conservation law used to treat electrons and holes separately.
-     # We only do this in the left lead, where the pairing is zero.
-     lead0 = kwant.Builder(sym_left, conservation_law=-tau_z, particle_hole=tau_y)
-     lead0[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z
-     lead0[lat.neighbors()] = -t * tau_z
- #HIDDEN_END_ttth
- #HIDDEN_BEGIN_zuuw
-     # Right lead - superconducting, so the order parameter is included.
-     sym_right = kwant.TranslationalSymmetry((a, 0))
-     lead1 = kwant.Builder(sym_right)
-     lead1[(lat(0, j) for j in range(W))] = (4 * t - mu) * tau_z + Delta * tau_x
-     lead1[lat.neighbors()] = -t * tau_z
- 
-     #### Attach the leads and return the system. ####
-     syst.attach_lead(lead0)
-     syst.attach_lead(lead1)
- 
-     return syst
- #HIDDEN_END_zuuw
- 
- #HIDDEN_BEGIN_jbjt
- def plot_conductance(syst, energies):
-     # Compute conductance
-     data = []
-     for energy in energies:
-         smatrix = kwant.smatrix(syst, energy)
-         # Conductance is N - R_ee + R_he
-         data.append(smatrix.submatrix((0, 0), (0, 0)).shape[0] -
-                     smatrix.transmission((0, 0), (0, 0)) +
-                     smatrix.transmission((0, 1), (0, 0)))
- #HIDDEN_END_jbjt
--    pyplot.figure()
-+    fig = pyplot.figure()
-     pyplot.plot(energies, data)
-     pyplot.xlabel("energy [t]")
-     pyplot.ylabel("conductance [e^2/h]")
--    pyplot.show()
-+    pyplot.setp(fig.get_axes()[0].get_xticklabels(),
-+                fontsize=_defs.mpl_tick_size)
-+    pyplot.setp(fig.get_axes()[0].get_yticklabels(),
-+                fontsize=_defs.mpl_tick_size)
-+    fig.set_size_inches(_defs.mpl_width_in, _defs.mpl_width_in * 3. / 4.)
-+    fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-+    for extension in ('pdf', 'png'):
-+        fig.savefig("superconductor_transport_result." + extension,
-+                    dpi=_defs.dpi)
- 
- #HIDDEN_BEGIN_pqmp
- def check_PHS(syst):
-     # Scattering matrix
-     s = kwant.smatrix(syst, energy=0)
-     # Electron to electron block
-     s_ee = s.submatrix((0,0), (0,0))
-     # Hole to hole block
-     s_hh = s.submatrix((0,1), (0,1))
-     print('s_ee: \n', np.round(s_ee, 3))
-     print('s_hh: \n', np.round(s_hh[::-1, ::-1], 3))
-     print('s_ee - s_hh^*: \n',
-           np.round(s_ee - s_hh[::-1, ::-1].conj(), 3), '\n')
-     # Electron to hole block
-     s_he = s.submatrix((0,1), (0,0))
-     # Hole to electron block
-     s_eh = s.submatrix((0,0), (0,1))
-     print('s_he: \n', np.round(s_he, 3))
-     print('s_eh: \n', np.round(s_eh[::-1, ::-1], 3))
-     print('s_he + s_eh^*: \n',
-           np.round(s_he + s_eh[::-1, ::-1].conj(), 3))
- #HIDDEN_END_pqmp
- 
- def main():
-     syst = make_system(W=10)
- 
--    # Check that the system looks as intended.
--    kwant.plot(syst)
--
-     # Finalize the system.
-     syst = syst.finalized()
- 
-     # Check particle-hole symmetry of the scattering matrix
--    check_PHS(syst)
-+    with open('check_PHS_out.txt', 'w') as f:
-+        with redirect_stdout(f):
-+            check_PHS(syst)
- 
-     # Compute and plot the conductance
-     plot_conductance(syst, energies=[0.002 * i for i in range(-10, 100)])
- 
- 
- # Call the main function if the script gets executed (as opposed to imported).
- # See <http://docs.python.org/library/__main__.html>.
- if __name__ == '__main__':
-     main()

doc/source/code/include/quantum_wire_revisited.py

-# Tutorial 2.2.3. Building the same system with less code
-# =======================================================
-#
-# Physics background
-# ------------------
-#  Conductance of a quantum wire; subbands
-#
-# Kwant features highlighted
-# --------------------------
-#  - Using iterables and builder.HoppingKind for making systems
-#  - introducing `reversed()` for the leads
-#
-# Note: Does the same as tutorial1a.py, but using other features of Kwant.
-
-#HIDDEN_BEGIN_xkzy
-import kwant
-
-# For plotting
-from matplotlib import pyplot
-
-
-def make_system(a=1, t=1.0, W=10, L=30):
-    # Start with an empty tight-binding system and a single square lattice.
-    # `a` is the lattice constant (by default set to 1 for simplicity.
-    lat = kwant.lattice.square(a)
-
-    syst = kwant.Builder()
-#HIDDEN_END_xkzy
-
-    #### Define the scattering region. ####
-#HIDDEN_BEGIN_vvjt
-    syst[(lat(x, y) for x in range(L) for y in range(W))] = 4 * t
-#HIDDEN_END_vvjt
-#HIDDEN_BEGIN_nooi
-    syst[lat.neighbors()] = -t
-#HIDDEN_END_nooi
-
-    #### Define and attach the leads. ####
-    # Construct the left lead.
-#HIDDEN_BEGIN_iepx
-    lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 0)))
-    lead[(lat(0, j) for j in range(W))] = 4 * t
-    lead[lat.neighbors()] = -t
-#HIDDEN_END_iepx
-
-    # Attach the left lead and its reversed copy.
-#HIDDEN_BEGIN_yxot
-    syst.attach_lead(lead)
-    syst.attach_lead(lead.reversed())
-
-    return syst
-#HIDDEN_END_yxot
-
-
-#HIDDEN_BEGIN_ayuk
-def plot_conductance(syst, energies):
-    # Compute conductance
-    data = []
-    for energy in energies:
-        smatrix = kwant.smatrix(syst, energy)
-        data.append(smatrix.transmission(1, 0))
-
-    pyplot.figure()
-    pyplot.plot(energies, data)
-    pyplot.xlabel("energy [t]")
-    pyplot.ylabel("conductance [e^2/h]")
-    pyplot.show()
-#HIDDEN_END_ayuk
-
-
-#HIDDEN_BEGIN_cjel
-def main():
-    syst = make_system()
-
-    # Check that the system looks as intended.
-    kwant.plot(syst)
-
-    # Finalize the system.
-    syst = syst.finalized()
-
-    # We should see conductance steps.
-    plot_conductance(syst, energies=[0.01 * i for i in range(100)])
-#HIDDEN_END_cjel
-
-
-# Call the main function if the script gets executed (as opposed to imported).
-# See <http://docs.python.org/library/__main__.html>.
-#HIDDEN_BEGIN_ypbj
-if __name__ == '__main__':
-    main()
-#HIDDEN_END_ypbj

doc/source/conf.py

@@ -15,8 +15,15 @@ import sys
 import os
 import string
 from distutils.util import get_platform
-sys.path.insert(0, "../../build/lib.{0}-{1}.{2}".format(
-        get_platform(), *sys.version_info[:2]))
+
+package_path = os.path.abspath(
+    "../../build/lib.{0}-{1}.{2}"
+    .format(get_platform(), *sys.version_info[:2]))
+
+# Insert into sys.path so that we can import kwant here
+sys.path.insert(0, package_path)
+# Insert into PYTHONPATH so that jupyter-sphinx will pick it up
+os.environ['PYTHONPATH'] = ':'.join((package_path, os.environ.get('PYTHONPATH','')))
 
 import kwant
 import kwant.qsymm
@@ -31,7 +38,8 @@ sys.path.insert(0, os.path.abspath('../sphinxext'))
 
 extensions = ['sphinx.ext.autodoc', 'sphinx.ext.autosummary',
               'sphinx.ext.todo', 'sphinx.ext.mathjax', 'numpydoc',
-              'kwantdoc', 'sphinx.ext.linkcode']
+              'kwantdoc', 'sphinx.ext.linkcode', 'jupyter_sphinx.execute',
+              'sphinxcontrib.rsvgconverter']
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['../templates']
@@ -99,6 +107,15 @@ pygments_style = 'sphinx'
 # Do not show all class members automatically in the class documentation
 numpydoc_show_class_members = False
 
+# Jupyter Sphinx config
+
+jupyter_sphinx_thebelab_config = {
+    "binderOptions": {
+        "repo": "kwant-project/binder",
+        "ref": "master",
+    }
+}
+
 # -- Options for HTML output ---------------------------------------------------
 
 # The theme to use for HTML and HTML Help pages.  Major themes that come with
@@ -218,7 +235,7 @@ latex_documents = [
 
 latex_engine = 'xelatex'
 
-latex_use_xindy = False
+latex_use_xindy = False  # Xindy not installable in CI environment
 
 # The name of an image file (relative to this directory) to place at the top of
 # the title page.
@@ -239,8 +256,9 @@ latex_domain_indices = False
 autosummary_generate = True
 
 autoclass_content = "both"
-autodoc_default_flags = ['show-inheritance']
-
+autodoc_default_options = {
+    'show-inheritance': True,
+}
 
 # -- Teach Sphinx to document bound methods like functions ---------------------
 import types

doc/source/pre/whatsnew/0.2.rst

@@ -42,7 +42,7 @@ description check `~kwant.plotter.plot` documentation.
 
 The behavior of `~kwant.plotter.plot` with low level systems has changed.
 Arguments of plot which are functions are given site numbers in place of
-`~kwant.builder.Site` objects when plotting a low level system.  This
+`~kwant.system.Site` objects when plotting a low level system.  This
 provides an easy way to make the appearance of lines and symbols depend on
 computation results.
 

doc/source/pre/whatsnew/1.0.rst

@@ -39,7 +39,7 @@ now it suffices to write ::
     syst[lat.neighbors()] = t
 
 This is possible because `~kwant.builder.Builder` now accepts *functions* as
-keys in addition to `~kwant.builder.Site` objects and tuples of them
+keys in addition to `~kwant.system.Site` objects and tuples of them
 (hoppings).  These functions are expected to yield either sites or hoppings,
 when given a builder instance as the sole argument. The use of such keys is to
 implement sets of sites or hoppings that depend on what is already present in

doc/source/pre/whatsnew/1.3.rst

@@ -28,11 +28,43 @@ tight-binding band structures or construct systems with different/lower
 symmetry. For example in just a few lines we can construct a two-band model that exhibits
 a quantum anomalous spin Hall phase:
 
-.. literalinclude:: ../../code/include/plot_qahe.py
-    :start-after: HIDDEN_BEGIN_model
-    :end-before: HIDDEN_END_model
+.. jupyter-kernel::
+    :id: plot_qahe
+
+.. jupyter-execute::
+    :hide-code:
+
+    # Comprehensive example: quantum anomalous Hall effect
+    # ====================================================
+    #
+    # Physics background
+    # ------------------
+    # + Quantum anomalous Hall effect
+    #
+    # Features highlighted
+    # --------------------
+    # + Use of `kwant.continuum` to discretize a continuum Hamiltonian
+    # + Use of `kwant.operator` to compute local current
+    # + Use of `kwant.plotter.current` to plot local current
+
+    import math
+    import matplotlib.pyplot
+    import kwant
+    import kwant.continuum
+
+.. jupyter-execute:: ../../tutorial/boilerplate.py
+    :hide-code:
+
+.. jupyter-execute::
 
-From: :download:`QAHE example script <../../code/download/plot_qahe.py>`
+    def make_model(a):
+        ham = ("alpha * (k_x * sigma_x - k_y * sigma_y)"
+               "+ (m + beta * kk) * sigma_z"
+               "+ (gamma * kk + U) * sigma_0")
+        subs = {"kk": "k_x**2 + k_y**2"}
+        return kwant.continuum.discretize(ham, locals=subs, grid=a)
+
+From: :jupyter-download:script:`plot_qahe`
 
 See the tutorial: :doc:`../../tutorial/discretize`
 
@@ -71,13 +103,47 @@ The example below shows edge states of a quantum anomalous Hall phase
 of the two-band model shown in the `above section
 <#tools-for-continuum-hamiltonians>`_:
 
-.. literalinclude:: ../../code/include/plot_qahe.py
-    :start-after: HIDDEN_BEGIN_current
-    :end-before: HIDDEN_END_current
+.. jupyter-execute::
+    :hide-code:
+
+    def make_system(model, L):
+        def lead_shape(site):
+            x, y = site.pos / L
+            return abs(y) < 0.5
+
+        # QPC shape: a rectangle with 2 gaussians
+        # etched out of the top and bottom edge.
+        def central_shape(site):
+            x, y = site.pos / L
+            return abs(x) < 3/5 and abs(y) < 0.5 - 0.4 * math.exp(-40 * x**2)
+
+        lead = kwant.Builder(kwant.TranslationalSymmetry(
+            model.lattice.vec((-1, 0))))
+        lead.fill(model, lead_shape, (0, 0))
+
+        syst = kwant.Builder()
+        syst.fill(model, central_shape, (0, 0))
+        syst.attach_lead(lead)
+        syst.attach_lead(lead.reversed())
+
+        return syst.finalized()
+
+    # Set up our model and system, and define the model parameters.
+    params = dict(alpha=0.365, beta=0.686, gamma=0.512, m=-0.01, U=0)
+    model = make_model(1)
+    syst = make_system(model, 70)
+
+    # Calculate the scattering states at energy 'm' coming from the left
+    # lead, and the associated particle current.
+    psi = kwant.wave_function(syst, energy=params['m'], params=params)(0)
+
+.. jupyter-execute::
 
-.. image:: ../../code/figure/plot_qahe_current.*
+    J = kwant.operator.Current(syst).bind(params=params)
+    current = sum(J(p) for p in psi)
+    kwant.plotter.current(syst, current);
 
-From: :download:`QAHE example script <../../code/download/plot_qahe.py>`
+From: :jupyter-download:script:`plot_qahe`
 
 Scattering states with discrete symmetries and conservation laws
 ----------------------------------------------------------------

doc/source/pre/whatsnew/1.5.rst

+What's new in Kwant 1.5
+=======================
+
+This article explains the user-visible changes in Kwant 1.5.0.
+
+
+Value functions may now be vectorized
+-------------------------------------
+It is now possible to define value functions that act on whole
+arrays of sites at a time.
+::
+
+    def onsite(sites):
+        r = sites.positions()
+        return np.exp(-np.linalg.norm(r, axis=1)**2)
+
+    def hopping(to_sites, from_sites):
+        dr = to_sites.positions() - from_sites.positions()
+        return 1 / np.linalg.norm(dr, axis=1)**2
+
+
+    lat = kwant.lattice.square(norbs=1)
+    syst = kwant.Builder(vectorize=True)
+    syst[(lat(i, j) for i in range(4) for j in range(10)] = onsite
+    syst[lat.neighbors()] = hopping
+    syst = syst.finalized()
+
+    syst.hamiltonian_submatrix()
+
+Notice that when creating the Builder we provide the ``vectorize=True`` flag,
+and that the value functions now take `~kwant.system.SiteArray` objects
+(which have a ``positions`` method, rather than a ``pos`` property as for
+`~kwant.system.Site`). We can now operate on the array of site positions
+using ``numpy``. Using vectorized value functions in this way can produce an
+order of magnitude speedup when evaluating system Hamiltonians (though this
+speedup may be masked by other parts of your computation e.g. calculating
+the scattering matrix).
+
+
+Deprecation of leaving 'norbs' unset for site families
+------------------------------------------------------
+When constructing site families (e.g. lattices) it is now deprecated to
+leave the 'norbs' parameter unset. This will now raise a
+KwantDeprecationWarning and will be disallowed in a future version of
+Kwant. For example, when constructing a square lattice with 1 orbital
+per site, use::
+
+    kwant.lattice.square(norbs=1)
+
+rather than::
+
+    kwant.lattice.square()
+
+
+Automatic addition of Peierls phase terms to Builders
+-----------------------------------------------------
+Kwant 1.4 introduced `kwant.physics.magnetic_gauge` for computing Peierls
+phases for arbitrary geometries and for systems with leads. Using this
+functionality requires that the system value functions are equipped to
+take the required Peierls phase parameters, which is not possible when
+you are not in direct control of the value functions (e.g. discretized
+systems). In Kwant 1.5 we have added the missing piece of functionality,
+`kwant.builder.add_peierls_phase`, which properly adds the Peierls phase
+parameters to a Builder's value functions::
+
+    syst = make_system()
+    lead = make_lead()
+    syst.attach_lead(lead)
+    syst.attach_lead(lead.reversed())
+
+    fsyst, phase = kwant.builder.add_peierls_phase(syst)
+
+    def B_syst(pos):
+        return np.exp(-np.sum(pos * pos))
+
+    kwant.smatrix(fsyst, parameters=dict(t=-1, **phase(B_syst, 0, 0)))
+
+
+Improved tutorial building
+--------------------------
+Improving or adding to Kwant's tutorial is now much simpler. Now
+the text and code for each tutorial is kept in the same file, making
+it easy to see where changes need to be made, and images generated by
+the code are inserted directly into the document thanks to the magic of
+`jupyter-sphinx <https://github.com/jupyter-widgets/jupyter-sphinx/>`_.
+It has never been easier to get started contributing to Kwant by
+helping us improve our documentation.
+
+Discretization onto a Landau level basis
+----------------------------------------
+The Hamiltonian for a system infinite in at least two dimensions and with
+a constant applied magnetic field may be expressed in a basis of Landau levels.
+The momenta in the plane perpendicular to the magnetic field direction are
+written in terms of the Landau level ladder operators:
+
+.. math::
+    k_x = \sqrt{\frac{B}{2}} (a + a^\dagger) \quad\quad
+    k_y = i \sqrt{\frac{B}{2}} (a - a^\dagger)
+
+The Hamiltonian is then expressed in terms of these ladder operators, which
+allows for a straight-forward discretization in the basis of Landau levels,
+provided that the basis is truncated after $N$ levels.
+`kwant.continuum.discretize_landau` makes this procedure simple::
+
+    syst = kwant.continuum.discretize_landau("k_x**2 + k_y**2", N)
+    syst.finalized().hamiltonian_submatrix(params=dict(B=0.5))
+
+`~kwant.continuum.discretize_landau` produces a regular Kwant Builder that
+can be inspected or finalized as usual. 3D Hamiltonians for systems that
+extend into the direction perpendicular to the magnetic field are also
+possible::
+
+    template = kwant.continuum.discretize_landau("k_x**2 + k_y**2 + k_z**2 + V(z)", N)
+
+This will produce a Builder with a single translational symmetry, which can be
+directly finalized, or can be used as a template for e.g. a heterostructure stack
+in the direction of the magnetic field::
+
+    def stack(site):
+        z, = site.pos
+        return 0 <= z < 10
+
+    template = kwant.continuum.discretize_landau("k_x**2 + k_y**2 + k_z**2 + V(z)", N)
+    syst = kwant.Builder()
+    syst.fill(template, stack, (0,))

doc/source/pre/whatsnew/index.rst

@@ -2,6 +2,7 @@ What's new in Kwant
 ===================
 
 .. toctree::
+    1.5
     1.4
     1.3
     1.2

doc/source/reference/kwant.builder.rst

@@ -9,7 +9,6 @@ Types
    :toctree: generated/
 
    Builder
-   Site
    HoppingKind
    SimpleSiteFamily
    BuilderLead
@@ -17,12 +16,28 @@ Types
    ModesLead
    FiniteSystem
    InfiniteSystem
+   FiniteVectorizedSystem
+   InfiniteVectorizedSystem
 
 Abstract base classes
 ---------------------
 .. autosummary::
    :toctree: generated/
 
-   SiteFamily
    Symmetry
    Lead
+
+Functions
+---------
+.. autosummary::
+   :toctree: generated/
+
+   add_peierls_phase
+
+Mixin Classes
+-------------
+.. autosummary::
+   :toctree: generated/
+
+   _FinalizedBuilderMixin
+   _VectorizedFinalizedBuilderMixin

doc/source/reference/kwant.continuum.rst

@@ -11,6 +11,7 @@ Discretizer
     discretize
     discretize_symbolic
     build_discretized
+    discretize_landau
 
 Symbolic helpers
 ----------------
@@ -19,3 +20,11 @@ Symbolic helpers
 
     sympify
     lambdify
+    to_landau_basis
+
+Other
+-----
+.. autosummary::
+    :toctree: generated/
+
+    LandauLattice

doc/source/reference/kwant.system.rst

@@ -19,10 +19,32 @@ In practice, very often Kwant systems are finalized
 `kwant.builder.FiniteSystem` or `kwant.builder.InfiniteSystem`) and offer some
 additional attributes.
 
+Sites
+-----
+.. autosummary::
+   :toctree: generated/
+
+   Site
+   SiteArray
+   SiteFamily
+
+Systems
+-------
 .. autosummary::
    :toctree: generated/
 
    System
+   VectorizedSystem
    InfiniteSystem
+   InfiniteVectorizedSystem
    FiniteSystem
+   FiniteVectorizedSystem
    PrecalculatedLead
+
+Mixin Classes
+-------------
+.. autosummary::
+   :toctree: generated/
+
+   FiniteSystemMixin
+   InfiniteSystemMixin

doc/source/tutorial/boilerplate.py

+import matplotlib
+import matplotlib.pyplot
+from IPython.display import set_matplotlib_formats
+
+matplotlib.rcParams['figure.figsize'] = matplotlib.pyplot.figaspect(1) * 2
+set_matplotlib_formats('svg')

doc/source/tutorial/discretize.rst

@@ -16,10 +16,40 @@ The `~kwant.continuum` sub-package aims to be a solution to this problem.
 It is a collection of tools for working with
 continuum models and for discretizing them into tight-binding models.
 
+.. seealso::
+    You can execute the code examples live in your browser by
+    activating thebelab:
+
+    .. thebe-button:: Activate Thebelab
+
 .. seealso::
     The complete source code of this tutorial can be found in
-    :download:`discretize.py </code/download/discretize.py>`
+    :jupyter-download:script:`discretize`
+
+.. jupyter-kernel::
+    :id: discretize
+
+.. jupyter-execute::
+    :hide-code:
+
+    # Tutorial 2.9. Processing continuum Hamiltonians with discretize
+    # ===============================================================
+    #
+    # Physics background
+    # ------------------
+    #  - tight-binding approximation of continuous Hamiltonians
+    #
+    # Kwant features highlighted
+    # --------------------------
+    #  - kwant.continuum.discretize
+
+    import matplotlib as mpl
+    from matplotlib import pyplot
 
+    import kwant
+
+.. jupyter-execute:: boilerplate.py
+    :hide-code:
 
 .. _tutorial_discretizer_introduction:
 
@@ -55,9 +85,16 @@ Using `~kwant.continuum.discretize` to obtain a template
 ........................................................
 First we must explicitly import the `kwant.continuum` package:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_import
-    :end-before: #HIDDEN_END_import
+.. jupyter-execute::
+
+    import kwant.continuum
+
+.. jupyter-execute::
+    :hide-code:
+
+    import scipy.sparse.linalg
+    import scipy.linalg
+    import numpy as np
 
 The function `kwant.continuum.discretize` takes a symbolic Hamiltonian and
 turns it into a `~kwant.builder.Builder` instance with appropriate spatial
@@ -65,17 +102,16 @@ symmetry that serves as a template.
 (We will see how to use the template to build systems with a particular
 shape later).
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_symbolic_discretization
-    :end-before: #HIDDEN_END_symbolic_discretization
+.. jupyter-execute::
+
+    template = kwant.continuum.discretize('k_x * A(x) * k_x')
+    print(template)
 
 It is worth noting that ``discretize`` treats ``k_x`` and ``x`` as
 non-commuting operators, and so their order is preserved during the
 discretization process.
 
-The builder produced by ``discretize`` may be printed to show the source code of its onsite and hopping functions (this is a special feature of builders returned by ``discretize``):
-
-.. literalinclude:: /code/figure/discretizer_intro_verbose.txt
+Printing the Builder produced by ``discretize`` shows the source code of its onsite and hopping functions (this is a special feature of builders returned by ``discretize``).
 
 .. specialnote:: Technical details
 
@@ -139,26 +175,45 @@ where :math:`V(x, y)` is some arbitrary potential.
 First, use ``discretize`` to obtain a
 builder that we will use as a template:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_template
-    :end-before: #HIDDEN_END_template
+.. jupyter-execute::
+
+    hamiltonian = "k_x**2 + k_y**2 + V(x, y)"
+    template = kwant.continuum.discretize(hamiltonian)
+    print(template)
 
 We now use this system with the `~kwant.builder.Builder.fill`
 method of `~kwant.builder.Builder` to construct the system we
 want to investigate:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_fill
-    :end-before: #HIDDEN_END_fill
+.. jupyter-execute::
+
+    def stadium(site):
+        (x, y) = site.pos
+        x = max(abs(x) - 20, 0)
+        return x**2 + y**2 < 30**2
+
+    syst = kwant.Builder()
+    syst.fill(template, stadium, (0, 0));
+    syst = syst.finalized()
 
 After finalizing this system, we can plot one of the system's
 energy eigenstates:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_plot_eigenstate
-    :end-before: #HIDDEN_END_plot_eigenstate
+.. jupyter-execute::
+
+    def plot_eigenstate(syst, n=2, Vx=.0003, Vy=.0005):
+
+        def potential(x, y):
+            return Vx * x + Vy * y
+
+        ham = syst.hamiltonian_submatrix(params=dict(V=potential), sparse=True)
+        evecs = scipy.sparse.linalg.eigsh(ham, k=10, which='SM')[1]
+        kwant.plotter.density(syst, abs(evecs[:, n])**2, show=False)
+
+.. jupyter-execute::
+    :hide-code:
 
-.. image:: /code/figure/discretizer_gs.*
+    plot_eigenstate(syst)
 
 Note in the above that we pass the spatially varying potential *function*
 to our system via a parameter called ``V``, because the symbol :math:`V`
@@ -168,6 +223,8 @@ In addition, the function passed as ``V`` expects two input parameters ``x``
 and ``y``, the same as in the initial continuum Hamiltonian.
 
 
+.. _discretize-bhz-model:
+
 Models with more structure: Bernevig-Hughes-Zhang
 .................................................
 
@@ -176,32 +233,95 @@ model [1]_ [2]_, one can provide matrix input to `~kwant.continuum.discretize`
 using ``identity`` and ``kron``. For example, the definition of the BHZ model can be
 written succinctly as:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_define_qsh
-    :end-before: #HIDDEN_END_define_qsh
+.. jupyter-execute::
+
+    hamiltonian = """
+       + C * identity(4) + M * kron(sigma_0, sigma_z)
+       - B * (k_x**2 + k_y**2) * kron(sigma_0, sigma_z)
+       - D * (k_x**2 + k_y**2) * kron(sigma_0, sigma_0)
+       + A * k_x * kron(sigma_z, sigma_x)
+       - A * k_y * kron(sigma_0, sigma_y)
+    """
+
+    a = 20
+
+    template = kwant.continuum.discretize(hamiltonian, grid=a)
 
 We can then make a ribbon out of this template system:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_define_qsh_build
-    :end-before: #HIDDEN_END_define_qsh_build
+.. jupyter-execute::
+
+    L, W = 2000, 1000
+
+    def shape(site):
+        (x, y) = site.pos
+        return (0 <= y < W and 0 <= x < L)
+
+    def lead_shape(site):
+        (x, y) = site.pos
+        return (0 <= y < W)
+
+    syst = kwant.Builder()
+    syst.fill(template, shape, (0, 0))
+
+    lead = kwant.Builder(kwant.TranslationalSymmetry([-a, 0]))
+    lead.fill(template, lead_shape, (0, 0))
+
+    syst.attach_lead(lead)
+    syst.attach_lead(lead.reversed())
+
+    syst = syst.finalized()
 
 and plot its dispersion using `kwant.plotter.bands`:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_plot_qsh_band
-    :end-before: #HIDDEN_END_plot_qsh_band
+.. jupyter-execute::
+
+    params = dict(A=3.65, B=-68.6, D=-51.1, M=-0.01, C=0)
+
+    kwant.plotter.bands(syst.leads[0], params=params,
+                        momenta=np.linspace(-0.3, 0.3, 201), show=False)
 
-.. image:: /code/figure/discretizer_qsh_band.*
+    pyplot.grid()
+    pyplot.xlim(-.3, 0.3)
+    pyplot.ylim(-0.05, 0.05)
+    pyplot.xlabel('momentum [1/A]')
+    pyplot.ylabel('energy [eV]')
+    pyplot.show()
 
 In the above we see the edge states of the quantum spin Hall effect, which
-we can visualize using `kwant.plotter.map`:
+we can visualize using `kwant.plotter.density`:
+
+.. jupyter-execute::
+
+    # get scattering wave functions at E=0
+    wf = kwant.wave_function(syst, energy=0, params=params)
+
+    # prepare density operators
+    sigma_z = np.array([[1, 0], [0, -1]])
+    prob_density = kwant.operator.Density(syst, np.eye(4))
+    spin_density = kwant.operator.Density(syst, np.kron(sigma_z, np.eye(2)))
+
+    # calculate expectation values and plot them
+    wf_sqr = sum(prob_density(psi) for psi in wf(0))  # states from left lead
+    rho_sz = sum(spin_density(psi) for psi in wf(0))  # states from left lead
+
+    fig, (ax1, ax2) = pyplot.subplots(1, 2, sharey=True, figsize=(16, 4))
+    kwant.plotter.density(syst, wf_sqr, ax=ax1)
+    kwant.plotter.density(syst, rho_sz, ax=ax2)
+
+    ax = ax1
+    im = [obj for obj in ax.get_children()
+          if isinstance(obj, mpl.image.AxesImage)][0]
+    fig.colorbar(im, ax=ax)
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_plot_qsh_wf
-    :end-before: #HIDDEN_END_plot_qsh_wf
+    ax = ax2
+    im = [obj for obj in ax.get_children()
+          if isinstance(obj, mpl.image.AxesImage)][0]
+    fig.colorbar(im, ax=ax)
 
-.. image:: /code/figure/discretizer_qsh_wf.*
+    ax1.set_title('Probability density')
+    ax2.set_title('Spin density')
+    pyplot.show()
 
 
 Limitations of discretization
@@ -238,29 +358,62 @@ example. Let us start from the continuum Hamiltonian
 We start by defining this model as a string and setting the value of the
 :math:`α` parameter:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_ls_def
-    :end-before: #HIDDEN_END_ls_def
+.. jupyter-execute::
+
+    hamiltonian = "k_x**2 * identity(2) + alpha * k_x * sigma_y"
+    params = dict(alpha=.5)
 
 Now we can use `kwant.continuum.lambdify` to obtain a function that computes
 :math:`H(k)`:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_ls_hk_cont
-    :end-before: #HIDDEN_END_ls_hk_cont
+.. jupyter-execute::
+
+    h_k = kwant.continuum.lambdify(hamiltonian, locals=params)
+    k_cont = np.linspace(-4, 4, 201)
+    e_cont = [scipy.linalg.eigvalsh(h_k(k_x=ki)) for ki in k_cont]
 
 We can also construct a discretized approximation using
 `kwant.continuum.discretize`, in a similar manner to previous examples:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_ls_hk_tb
-    :end-before: #HIDDEN_END_ls_hk_tb
+.. jupyter-execute::
+
+    def plot(ax, a=1):
+        template = kwant.continuum.discretize(hamiltonian, grid=a)
+        syst = kwant.wraparound.wraparound(template).finalized()
+
+        def h_k(k_x):
+            p = dict(k_x=k_x, **params)
+            return syst.hamiltonian_submatrix(params=p)
+
+        k_tb = np.linspace(-np.pi/a, np.pi/a, 201)
+        e_tb = [scipy.linalg.eigvalsh(h_k(k_x=a*ki)) for ki in k_tb]
+
+        ax.plot(k_cont, e_cont, 'r-')
+        ax.plot(k_tb, e_tb, 'k-')
+
+        ax.plot([], [], 'r-', label='continuum')
+        ax.plot([], [], 'k-', label='tight-binding')
+
+        ax.set_xlim(-4, 4)
+        ax.set_ylim(-1, 14)
+        ax.set_title('a={}'.format(a))
+
+        ax.set_xlabel('momentum [a.u.]')
+        ax.set_ylabel('energy [a.u.]')
+        ax.grid()
+        ax.legend()
 
 Below we can see the continuum and tight-binding dispersions for two
 different values of the discretization grid spacing :math:`a`:
 
-.. image:: /code/figure/discretizer_lattice_spacing.*
+.. jupyter-execute::
+    :hide-code:
 
+    _, (ax1, ax2) = pyplot.subplots(1, 2, sharey=True, figsize=(12, 4))
+
+    plot(ax1, a=1)
+    plot(ax2, a=.25)
+    pyplot.show()
 
 We clearly see that the smaller grid spacing is, the better we approximate
 the original continuous dispersion. It is also worth remembering that the
@@ -284,20 +437,26 @@ It is possible to use ``identity`` (for identity matrix), ``kron`` (for Kronecke
 expressions involving matrices. Matrices can also be provided explicitly using
 square ``[]`` brackets. For example, all following expressions are equivalent:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_subs_1
-    :end-before: #HIDDEN_END_subs_1
+.. jupyter-execute::
+
+    sympify = kwant.continuum.sympify
+    subs = {'sx': [[0, 1], [1, 0]], 'sz': [[1, 0], [0, -1]]}
+
+    e = (
+        sympify('[[k_x**2, alpha * k_x], [k_x * alpha, -k_x**2]]'),
+        sympify('k_x**2 * sigma_z + alpha * k_x * sigma_x'),
+        sympify('k_x**2 * sz + alpha * k_x * sx', locals=subs),
+    )
 
-.. literalinclude:: /code/figure/discretizer_subs_1.txt
+    print(e[0] == e[1] == e[2])
 
 We can use the ``locals`` keyword parameter to substitute expressions
 and numerical values:
 
-.. literalinclude:: /code/include/discretize.py
-    :start-after: #HIDDEN_BEGIN_subs_2
-    :end-before: #HIDDEN_END_subs_2
+.. jupyter-execute::
 
-.. literalinclude:: /code/figure/discretizer_subs_2.txt
+    subs = {'A': 'A(x) + B', 'V': 'V(x) + V_0', 'C': 5}
+    print(sympify('k_x * A * k_x + V + C', locals=subs))
 
 Symbolic expressions obtained in this way can be directly passed to all
 ``discretizer`` functions.

doc/source/tutorial/faq.rst

@@ -5,6 +5,23 @@ questions that are discussed there.  The `Kwant paper
 <https://downloads.kwant-project.org/doc/kwant-paper.pdf>`_ also digs deeper
 into Kwant's structure.
 
+.. seealso::
+    You can execute the code examples live in your browser by
+    activating thebelab:
+
+    .. thebe-button:: Activate Thebelab
+
+.. jupyter-execute::
+    :hide-code:
+
+    import kwant
+    import numpy as np
+    import tinyarray
+    import matplotlib
+    from matplotlib import pyplot as plt
+
+.. jupyter-execute:: boilerplate.py
+    :hide-code:
 
 What is a system, and what is a builder?
 ========================================
@@ -48,11 +65,16 @@ combination of family and tag uniquely defines a site.
 
 For example let us create an empty tight binding system and add two sites:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_site
-    :end-before: #HIDDEN_END_site
+.. jupyter-execute::
+
+    a = 1
+    lat = kwant.lattice.square(a, norbs=1)
+    syst = kwant.Builder()
+
+    syst[lat(1, 0)] = 4
+    syst[lat(1, 1)] = 4
 
-.. image:: /code/figure/faq_site.*
+    kwant.plot(syst);
 
 In the above snippet we added 2 sites: ``lat(1, 0)`` and ``lat(0, 1)``. Both
 of these sites belong to the same family, ``lat``, but have different tags:
@@ -79,11 +101,11 @@ tuples, for example lists, are not treated as hoppings.
 Starting from the example code from `What is a site?`_, we can add a hopping
 to our system in the following way:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_hopping
-    :end-before: #HIDDEN_END_hopping
+.. jupyter-execute::
 
-.. image:: /code/figure/faq_hopping.*
+    syst[(lat(1, 0), lat(1, 1))] = 1j
+
+    kwant.plot(syst);
 
 Visually, a hopping is represented as a line that joins two sites.
 
@@ -139,18 +161,29 @@ Let's create two monatomic lattices (``lat_a`` and ``lat_b``).  ``(1, 0)``
 and ``(0, 1)`` will be the primitive vectors and ``(0, 0)`` and ``(0.5, 0.5)``
 the origins of the two lattices:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_lattice_monatomic
-    :end-before: #HIDDEN_END_lattice_monatomic
+.. jupyter-execute::
+
+    # Two monatomic lattices
+    primitive_vectors = [(1, 0), (0, 1)]
+    lat_a = kwant.lattice.Monatomic(primitive_vectors, offset=(0, 0), norbs=1)
+    lat_b = kwant.lattice.Monatomic(primitive_vectors, offset=(0.5, 0.5), norbs=1)
+    # lat1 is equivalent to kwant.lattice.square(norbs=1)
+
+    syst = kwant.Builder()
 
-.. image:: /code/figure/faq_lattice.*
+    syst[lat_a(0, 0)] = 4
+    syst[lat_b(0, 0)] = 4
+
+    kwant.plot(syst);
 
 We can also create a ``Polyatomic`` lattice with the same primitive vectors and
 two sites in the basis:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_lattice_polyatomic
-    :end-before: #HIDDEN_END_lattice_polyatomic
+.. jupyter-execute::
+
+    # One polyatomic lattice containing two sublattices
+    lat = kwant.lattice.Polyatomic([(1, 0), (0, 1)], [(0, 0), (0.5, 0.5)], norbs=1)
+    sub_a, sub_b = lat.sublattices
 
 The two sublattices ``sub_a`` and ``sub_b`` are nothing else than ``Monatomic``
 instances, and are equivalent to ``lat_a`` and ``lat_b`` that we created
@@ -169,17 +202,40 @@ When plotting, how to color the different sublattices differently?
 ==================================================================
 In the following example we shall use a kagome lattice, which has three sublattices.
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_colors1
-    :end-before: #HIDDEN_END_colors1
+.. jupyter-execute::
+
+    lat = kwant.lattice.kagome(norbs=1)
+    syst = kwant.Builder()
+
+    a, b, c = lat.sublattices  # The kagome lattice has 3 sublattices
 
 As we can see below, we create a new plotting function that assigns a color for each family, and a different size for the hoppings depending on the family of the two sites. Finally we add sites and hoppings to our system and plot it with the new function.
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_colors2
-    :end-before: #HIDDEN_END_colors2
+.. jupyter-execute::
+
+    # Plot sites from different families in different colors
+    def family_color(site):
+        if site.family == a:
+            return 'red'
+        if site.family == b:
+            return 'green'
+        else:
+            return 'blue'
+
+    def plot_system(syst):
+        kwant.plot(syst, site_lw=0.1, site_color=family_color)
+
+    ## Add sites and hoppings.
+    for i in range(4):
+        for j in range (4):
+            syst[a(i, j)] = 4
+            syst[b(i, j)] = 4
+            syst[c(i, j)] = 4
+
+    syst[lat.neighbors()] = -1
 
-.. image:: /code/figure/faq_colors.*
+    ## Plot the system.
+    plot_system(syst)
 
 
 How to create many similar hoppings in one go?
@@ -197,31 +253,64 @@ integers).
 
 The following example shows how this can be used:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_direction1
-    :end-before: #HIDDEN_END_direction1
+.. jupyter-execute::
 
-.. image:: /code/figure/faq_direction1.*
+    # Create hopping between neighbors with HoppingKind
+    a = 1
+    syst = kwant.Builder()
+    lat = kwant.lattice.square(a, norbs=1)
+    syst[ (lat(i, j) for i in range(5) for j in range(5)) ] = 4
+
+    syst[kwant.builder.HoppingKind((1, 0), lat)] = -1
+    kwant.plot(syst);
 
 Note that ``HoppingKind`` only works with site families so you cannot use
 them directly with ``Polyatomic`` lattices; you have to explicitly specify
 the sublattices when creating a ``HoppingKind``:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_direction2
-    :end-before: #HIDDEN_END_direction2
+.. jupyter-execute::
+    :hide-code:
+
+    lat = kwant.lattice.kagome(norbs=1)
+    syst = kwant.Builder()
+
+    a, b, c = lat.sublattices  # The kagome lattice has 3 sublattices
+
+    for i in range(4):
+        for j in range (4):
+            syst[a(i, j)] = 4
+            syst[b(i, j)] = 4
+            syst[c(i, j)] = 4
+
+.. jupyter-execute::
+
+    # equivalent to syst[kwant.builder.HoppingKind((0, 1), b)] = -1
+    syst[kwant.builder.HoppingKind((0, 1), b, b)] = -1
 
 Here, we want the hoppings between the sites from sublattice b with a direction of (0,1) in the lattice coordinates.
 
-.. image:: /code/figure/faq_direction2.*
+.. jupyter-execute::
+    :hide-code:
+
+    plot_system(syst)
+
+.. jupyter-execute::
+    :hide-code:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_direction3
-    :end-before: #HIDDEN_END_direction3
+    # Delete the hoppings previously created
+    del syst[kwant.builder.HoppingKind((0, 1), b, b)]
+
+.. jupyter-execute::
+
+    syst[kwant.builder.HoppingKind((0, 0), a, b)] = -1
+    syst[kwant.builder.HoppingKind((0, 0), a, c)] = -1
+    syst[kwant.builder.HoppingKind((0, 0), c, b)] = -1
 
 Here, we create hoppings between the sites of the same lattice coordinates but from different families.
 
-.. image:: /code/figure/faq_direction3.*
+.. jupyter-execute::
+
+    plot_system(syst)
 
 
 How to set the hoppings between adjacent sites?
@@ -230,32 +319,53 @@ How to set the hoppings between adjacent sites?
 that returns a list of ``HoppingKind`` instances that connect sites with their
 (n-nearest) neighors:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_adjacent1
-    :end-before: #HIDDEN_END_adjacent1
+.. jupyter-execute::
 
-.. image:: /code/figure/faq_adjacent1.*
-.. image:: /code/figure/faq_adjacent2.*
+    # Create hoppings with lat.neighbors()
+    syst = kwant.Builder()
+    lat = kwant.lattice.square(norbs=1)
+    syst[(lat(i, j) for i in range(3) for j in range(3))] = 4
 
-As we can see in the figure above, ``lat.neighbors()`` (on the left) returns the hoppings between the first nearest neighbors and ``lat.neighbors(2)`` (on the right) returns the hoppings between the second nearest neighbors.
+    syst[lat.neighbors()] = -1  # Equivalent to lat.neighbors(1)
+    kwant.plot(syst);
+
+.. jupyter-execute::
+
+    del syst[lat.neighbors()]  # Delete all nearest-neighbor hoppings
+    syst[lat.neighbors(2)] = -1
+
+    kwant.plot(syst);
+
+As we can see in the figures above, ``lat.neighbors()`` returns the hoppings between the first nearest neighbors and ``lat.neighbors(2)`` returns the hoppings between the second nearest neighbors.
 
 When using a ``Polyatomic`` lattice ``neighbors()`` knows about the different
 sublattices:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_adjacent2
-    :end-before: #HIDDEN_END_adjacent2
+.. jupyter-execute::
+
+    # Create the system
+    lat = kwant.lattice.kagome(norbs=1)
+    syst = kwant.Builder()
+    a, b, c = lat.sublattices  # The kagome lattice has 3 sublattices
 
-.. image:: /code/figure/faq_adjacent3.*
+    for i in range(4):
+        for j in range (4):
+            syst[a(i, j)] = 4  # red
+            syst[b(i, j)] = 4  # green
+            syst[c(i, j)] = 4  # blue
+
+    syst[lat.neighbors()] = -1
+
+    plot_system(syst)
 
 However, if we use the ``neighbors()`` method of a single sublattice, we will
 only get the neighbors *on that sublattice*:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_adjacent3
-    :end-before: #HIDDEN_END_adjacent3
+.. jupyter-execute::
 
-.. image:: /code/figure/faq_adjacent4.*
+    del syst[lat.neighbors()]  # Delete the hoppings previously created
+    syst[a.neighbors()] = -1
+    plot_system(syst)
 
 Note in the above that there are *only* hoppings between the red sites. This
 is an artifact of the visualisation: the blue and green sites just happen to lie
@@ -268,12 +378,35 @@ To make a hole in the system, use ``del syst[site]``, just like with any other
 mapping. In the following example we remove all sites inside some "hole"
 region:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_hole
-    :end-before: #HIDDEN_END_hole
+.. jupyter-execute::
+
+    # Define the lattice and the (empty) system
+    a = 2
+    lat = kwant.lattice.cubic(a, norbs=1)
+    syst = kwant.Builder()
+
+    L = 10
+    W = 10
+    H = 2
+
+    # Add sites to the system in a cuboid
+
+    syst[(lat(i, j, k) for i in range(L) for j in range(W) for k in range(H))] = 4
+    kwant.plot(syst);
+
+.. jupyter-execute::
+
+    # Delete sites to create a hole
+
+    def in_hole(site):
+        x, y, z = site.pos / a - (L/2, W/2, H/2)  # position relative to centre
+        return abs(x) < L / 4 and abs(y) < W / 4
+
+    for site in filter(in_hole, list(syst.sites())):
+        del syst[site]
+
+    kwant.plot(syst);
 
-.. image:: /code/figure/faq_hole1.*
-.. image:: /code/figure/faq_hole2.*
 
 ``del syst[site]`` also works after hoppings have been added to the system.
 If a site is deleted, then all the hoppings to/from that site are also deleted.
@@ -287,9 +420,18 @@ what is a builder?`_ if you do not know the difference).
 
 We can access the sites of a ``Builder`` by using its `~kwant.builder.Builder.sites` method:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_sites1
-    :end-before: #HIDDEN_END_sites1
+.. jupyter-execute::
+    :hide-code:
+
+    builder = kwant.Builder()
+    lat = kwant.lattice.square(norbs=1)
+    builder[(lat(i, j) for i in range(3) for j in range(3))] = 4
+
+.. jupyter-execute::
+
+    # Before finalizing the system
+
+    sites = list(builder.sites())  # sites() doe *not* return a list
 
 The ``sites()`` method returns an *iterator* over the system sites, and in the
 above example we create a list from the contents of this iterator, which
@@ -300,9 +442,11 @@ well be returned in a different order.
 After finalization, when we are dealing with a ``System``, the sites themselves
 are stored in a list, which can be accessed via the ``sites`` attribute:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_sites2
-    :end-before: #HIDDEN_END_sites2
+.. jupyter-execute::
+
+    # After finalizing the system
+    syst = builder.finalized()
+    sites = syst.sites  # syst.sites is an actual list
 
 The order of sites in a ``System`` is fixed, and also defines the ordering of
 the system Hamiltonian, system wavefunctions etc. (see `How does Kwant order components of an individual wavefunction?`_ for details).
@@ -310,9 +454,9 @@ the system Hamiltonian, system wavefunctions etc. (see `How does Kwant order com
 ``System`` also contains the inverse mapping, ``id_by_site`` which gives us
 the index of a given site within the system:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_sites3
-    :end-before: #HIDDEN_END_sites3
+.. jupyter-execute::
+
+    i = syst.id_by_site[lat(0, 2)]  # we want the id of the site lat(0, 2)
 
 
 How to use different lattices for the scattering region and a lead?
@@ -322,19 +466,35 @@ which we want to connect to leads that contain sites from a square lattice.
 
 First we construct the central system:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_different_lattice1
-    :end-before: #HIDDEN_END_different_lattice1
+.. jupyter-execute::
+
+    # Define the scattering Region
+    L = 5
+    W = 5
+
+    lat = kwant.lattice.honeycomb(norbs=1)
+    subA, subB = lat.sublattices
 
-.. image:: /code/figure/faq_different_lattice1.*
+    syst = kwant.Builder()
+    syst[(subA(i, j) for i in range(L) for j in range(W))] = 4
+    syst[(subB(i, j) for i in range(L) for j in range(W))] = 4
+    syst[lat.neighbors()] = -1
+
+    kwant.plot(syst);
 
 and the lead:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_different_lattice2
-    :end-before: #HIDDEN_END_different_lattice2
+.. jupyter-execute::
+
+    # Create a lead
+    lat_lead = kwant.lattice.square(norbs=1)
+    sym_lead1 = kwant.TranslationalSymmetry((0, 1))
 
-.. image:: /code/figure/faq_different_lattice2.*
+    lead1 = kwant.Builder(sym_lead1)
+    lead1[(lat_lead(i, 0) for i in range(2, 7))] = 4
+    lead1[lat_lead.neighbors()] = -1
+
+    kwant.plot(syst);
 
 We cannot simply use `~kwant.builder.Builder.attach_lead` to attach this lead to the
 system with the honeycomb lattice because Kwant does not know how sites from
@@ -343,19 +503,22 @@ these two lattices should be connected.
 We must first add a layer of sites from the square lattice to the system and manually
 add the hoppings from these sites to the sites from the honeycomb lattice:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_different_lattice3
-    :end-before: #HIDDEN_END_different_lattice3
+.. jupyter-execute::
+
+    syst[(lat_lead(i, 5) for i in range(2, 7))] = 4
+    syst[lat_lead.neighbors()] = -1
+
+    # Manually attach sites from graphene to square lattice
+    syst[((lat_lead(i+2, 5), subB(i, 4)) for i in range(5))] = -1
 
-.. image:: /code/figure/faq_different_lattice3.*
+    kwant.plot(syst);
 
 ``attach_lead()`` will now be able to attach the lead:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_different_lattice4
-    :end-before: #HIDDEN_END_different_lattice4
+.. jupyter-execute::
 
-.. image:: /code/figure/faq_different_lattice4.*
+    syst.attach_lead(lead1)
+    kwant.plot(syst);
 
 
 How to cut a finite system out of a system with translational symmetries?
@@ -369,30 +532,50 @@ will be repeated in the desired shape to create the final system.
 
 For example, say we want to create a simple model on a cubic lattice:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_fill1
-    :end-before: #HIDDEN_END_fill1
+.. jupyter-execute::
+
+    # Create 3d model.
+    cubic = kwant.lattice.cubic(norbs=1)
+    sym_3d = kwant.TranslationalSymmetry([1, 0, 0], [0, 1, 0], [0, 0, 1])
+    model = kwant.Builder(sym_3d)
+    model[cubic(0, 0, 0)] = 4
+    model[cubic.neighbors()] = -1
 
 We have now created our "template" ``Builder`` which has 3 translational
 symmetries. Next we will fill a system with no translational symmetries with
 sites and hoppings from the template inside a cuboid:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_fill2
-    :end-before: #HIDDEN_END_fill2
+.. jupyter-execute::
+
+    # Build scattering region (white).
+    def cuboid_shape(site):
+        x, y, z = abs(site.pos)
+        return x < 4 and y < 10 and z < 3
 
-.. image:: /code/figure/faq_fill2.*
+    cuboid = kwant.Builder()
+    cuboid.fill(model, cuboid_shape, (0, 0, 0));
+    kwant.plot(cuboid);
 
 We can then use the original template to create a lead, which has 1 translational
 symmetry. We can then use this lead as a template to fill another section of
 the system with a cylinder of sites and hoppings:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_fill3
-    :end-before: #HIDDEN_END_fill3
+.. jupyter-execute::
+
+    # Build electrode (black).
+    def electrode_shape(site):
+        x, y, z = site.pos - (0, 5, 2)
+        return y**2 + z**2 < 2.3**2
 
-.. image:: /code/figure/faq_fill3.*
+    electrode = kwant.Builder(kwant.TranslationalSymmetry([1, 0, 0]))
+    electrode.fill(model, electrode_shape, (0, 5, 2))  # lead
 
+    # Scattering region
+    cuboid.fill(electrode, lambda s: abs(s.pos[0]) < 7, (0, 5, 4))
+
+    cuboid.attach_lead(electrode)
+
+    kwant.plot(cuboid);
 
 How does Kwant order the propagating modes of a lead?
 =====================================================
@@ -402,16 +585,44 @@ achieved with the `~kwant.system.InfiniteSystem.modes` method, which returns a
 pair of objects, the first of which contains the propagating modes of the
 system in a `~kwant.physics.PropagatingModes` object:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_pm
-    :end-before: #HIDDEN_END_pm
+.. jupyter-execute::
+
+    lat = kwant.lattice.square(norbs=1)
+
+    lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
+    lead[(lat(0, i) for i in range(3))] = 4
+    lead[lat.neighbors()] = -1
+
+    flead = lead.finalized()
+
+    E = 2.5
+    prop_modes, _ = flead.modes(energy=E)
 
 ``PropagatingModes`` contains the wavefunctions, velocities and momenta of the
 modes at the requested energy (2.5 in this example).  In order to understand
 the order in which these quantities are returned it is often useful to look at
 the a section of the band structure for the system in question:
 
-.. image:: /code/figure/faq_pm1.*
+.. jupyter-execute::
+    :hide-code:
+
+    def plot_and_label_modes(lead, E):
+        # Plot the different modes
+        pmodes, _ = lead.modes(energy=E)
+        kwant.plotter.bands(lead, show=False)
+        for i, k in enumerate(pmodes.momenta):
+            plt.plot(k, E, 'ko')
+            plt.annotate(str(i), xy=(k, E), xytext=(-5, 8),
+                         textcoords='offset points',
+                         bbox=dict(boxstyle='round,pad=0.1',fc='white', alpha=0.7))
+        plt.plot([-3, 3], [E, E], 'r--')
+        plt.ylim(E-1, E+1)
+        plt.xlim(-2, 2)
+        plt.xlabel("momentum")
+        plt.ylabel("energy")
+        plt.show()
+
+    plot_and_label_modes(flead, E)
 
 On the above band structure we have labelled the 4 modes in the order
 that they appear in the output of ``modes()`` at energy 2.5. Note that
@@ -425,7 +636,22 @@ the modes are sorted in the following way:
 For more complicated systems and band structures this can lead to some
 unintuitive orderings:
 
-.. image:: /code/figure/faq_pm2.*
+.. jupyter-execute::
+    :hide-code:
+
+    lat = kwant.lattice.square(norbs=2)
+
+    s0 = np.eye(2)
+    sz = np.array([[1, 0], [0, -1]])
+
+    lead2 = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
+
+    lead2[(lat(0, i) for i in range(2))] = np.diag([1.8, -1])
+    lead2[lat.neighbors()] = -1 * sz
+
+    flead2 = lead2.finalized()
+
+    plot_and_label_modes(flead2, 1)
 
 
 How does Kwant order scattering states?
@@ -452,10 +678,38 @@ When all the site families present in a system have only 1 degree of freedom
 per site (i.e.  the all the onsites are scalars) then the index into a
 wavefunction defined over the system is exactly the site index:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_ord1
-    :end-before: #HIDDEN_END_ord1
-.. literalinclude:: /code/figure/faq_ord1.txt
+.. jupyter-execute::
+    :hide-code:
+
+    def circle(R):
+        return lambda r: np.linalg.norm(r) < R
+
+
+    def make_system(lat):
+        norbs = lat.norbs
+        syst = kwant.Builder()
+        syst[lat.shape(circle(3), (0, 0))] = 4 * np.eye(norbs)
+        syst[lat.neighbors()] = -1 * np.eye(norbs)
+
+        lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))
+        lead[(lat(0, i) for i in range(-1, 2))] = 4 * np.eye(norbs)
+        lead[lat.neighbors()] = -1 * np.eye(norbs)
+
+        syst.attach_lead(lead)
+        syst.attach_lead(lead.reversed())
+
+        return syst.finalized()
+
+.. jupyter-execute::
+
+    lat = kwant.lattice.square(norbs=1)
+    syst = make_system(lat)
+    scattering_states = kwant.wave_function(syst, energy=1)
+    wf = scattering_states(0)[0]  # scattering state from lead 0 incoming in mode 0
+
+    idx = syst.id_by_site[lat(0, 0)]  # look up index of site
+
+    print('wavefunction on lat(0, 0): ', wf[idx])
 
 We see that the wavefunction on a single site is a single complex number, as
 expected.
@@ -467,10 +721,20 @@ to one another.  In the case where all site families in the system have the
 wavefunction into a matrix, where the row number indexes the site, and the
 column number indexes the degree of freedom on that site:
 
-.. literalinclude:: /code/include/faq.py
-    :start-after: #HIDDEN_BEGIN_ord2
-    :end-before: #HIDDEN_END_ord2
-.. literalinclude:: /code/figure/faq_ord2.txt
+.. jupyter-execute::
+
+    lat = kwant.lattice.square(norbs=2)
+    syst = make_system(lat)
+    scattering_states = kwant.wave_function(syst, energy=1)
+    wf = scattering_states(0)[0]  # scattering state from lead 0 incoming in mode 0
+
+    idx = syst.id_by_site[lat(0, 0)]  # look up index of site
+
+    # Group consecutive degrees of freedom from 'wf' together; these correspond
+    # to degrees of freedom on the same site.
+    wf = wf.reshape(-1, 2)
+
+    print('wavefunction on lat(0, 0): ', wf[idx])
 
 We see that the wavefunction on a single site is a *vector* of 2 complex numbers,
 as we expect.