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.
+
+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
@@ -23,6 +22,12 @@ Abstract base classes
 .. autosummary::
    :toctree: generated/
 
-   SiteFamily
    Symmetry
    Lead
+
+Functions
+---------
+.. autosummary::
+   :toctree: generated/
+
+   add_peierls_phase

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,6 +19,17 @@ 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/
 

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.