improve the tutorial script and image buiding machinery

The basic idea remains the same, but should work much better now. The difference is that images/example.py.diff is now a patch with complete context and becomes the authoritative source for both the visible flavor of an example as well as for its figure-generating variant. Both tutorial/example.py and images/example.py are extracted from this file by 'make html'. Thanks to the complete context the diffs are quite readable and may be modified directly. Alternatively, one may also modify the generated scripts. When tutorial/example.py has been modified, it will be patched and saved as images/example.py. The patching is done using the tool 'wiggle' that works much better than 'patch'. If a conflict occurs, conflict markers are added to the output file and its dated back to the dawn of time (i.e. 1970) in order to mark the conflicts as not yet resolved. After resolving, 'make html' is simply run again. Upon a successful execution of the figure-generating script the diff gets recreated automatically, 'generate-diffs' no longer exists.

improve the tutorial script and image buiding machinery
001dafbb · Christoph Groth · 503acc43 · 001dafbb · 001dafbb · 001dafbb
Commit 001dafbb authored 7 years ago by Christoph Groth
--- a/.gitignore
+++ b/.gitignore
@@ -8,6 +8,7 @@
 /build
 /dist
 /doc/build
+/doc/source/tutorial/*.py
 /doc/source/reference/generated/
 /doc/source/images/*.png
 /doc/source/images/*.pdf

--- a/INSTALL.rst
+++ b/INSTALL.rst
@@ -171,6 +171,24 @@ Because of some quirks of how Sphinx works, it might be necessary to execute
 done, Sphinx may mistakenly use PNG files for PDF output or other problems may
 appear.

+When ``make html`` is run, modified tutorial example scripts are executed to
+update any figures that might have changed.  The machinery behind this works as
+follows.  The canonical source for a tutorial script, say ``graphene.py`` is
+the file ``doc/source/images/graphene.py.diff``.  This diff file contains the
+information to recreate two versions of ``graphene.py``: a version that is
+presented in the documentation (``doc/source/tutorial/graphene.py``), and a
+version that is used to generate the figures for the documentation
+(``doc/source/images/graphene.py``).  Both versions are related but differ
+e.g. in the details of the plotting.  When ``make html`` is run, both versions
+are extracted form the diff file.
+
+The diff file may be modified directly.  Another possible way of working is to
+directly modify either the tutorial script or the figure generation script.
+Then ``make html`` will use the command line tool `wiggle
+<http://neil.brown.name/wiggle/>`_ to propagate the modifications accordingly.
+This will often just work, but may sometimes result in conflicts, in which case
+a message will be printed.  The conflicts then have to be resolved much like
+with a version control system.

 ****************************
 Hints for specific platforms

--- a/doc/Makefile
+++ b/doc/Makefile
@@ -35,6 +35,7 @@ GENERATEDPDF    = $(patsubst %.svg,%.pdf,$(IMAGESOURCES))
 # deleting the corresponding flag file.
 SCRIPTS = $(patsubst source/images/%.diff,%,$(wildcard source/images/*.py.diff))
 FLAGS = $(patsubst %.py,source/images/.%_flag,$(SCRIPTS))
+INCLUDES = $(patsubst %,source/tutorial/%,$(SCRIPTS))

 .PHONY: help clean realclean html dirhtml pickle json htmlhelp qthelp latex changes linkcheck doctest

@@ -59,37 +60,38 @@ clean:

 realclean: clean
 	-rm -f $(FLAGS)
+	-rm -f $(INCLUDES)
 	-rm -f $(patsubst %,source/images/%,$(SCRIPTS))
 	-rm -f $(patsubst %.py,source/images/%_*.png,$(SCRIPTS))
 	-rm -f $(patsubst %.py,source/images/%_*.pdf,$(SCRIPTS))

-html:	$(FLAGS)
+html:	$(FLAGS) $(INCLUDES)
 	$(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html
 	@echo
 	@echo "Build finished. The HTML pages are in $(BUILDDIR)/html."

-dirhtml: $(FLAGS)
+dirhtml: $(FLAGS) $(INCLUDES)
 	$(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) $(BUILDDIR)/dirhtml
 	@echo
 	@echo "Build finished. The HTML pages are in $(BUILDDIR)/dirhtml."

-pickle: $(FLAGS)
+pickle: $(FLAGS) $(INCLUDES)
 	$(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) $(BUILDDIR)/pickle
 	@echo
 	@echo "Build finished; now you can process the pickle files."

-json:   $(FLAGS)
+json:   $(FLAGS) $(INCLUDES)
 	$(SPHINXBUILD) -b json $(ALLSPHINXOPTS) $(BUILDDIR)/json
 	@echo
 	@echo "Build finished; now you can process the JSON files."

-htmlhelp: $(FLAGS)
+htmlhelp: $(FLAGS) $(INCLUDES)
 	$(SPHINXBUILD) -b htmlhelp $(ALLSPHINXOPTS) $(BUILDDIR)/htmlhelp
 	@echo
 	@echo "Build finished; now you can run HTML Help Workshop with the" \
 	      ".hhp project file in $(BUILDDIR)/htmlhelp."

-qthelp: $(FLAGS)
+qthelp: $(FLAGS) $(INCLUDES)
 	$(SPHINXBUILD) -b qthelp $(ALLSPHINXOPTS) $(BUILDDIR)/qthelp
 	@echo
 	@echo "Build finished; now you can run "qcollectiongenerator" with the" \
@@ -98,7 +100,7 @@ qthelp: $(FLAGS)
 	@echo "To view the help file:"
 	@echo "# assistant -collectionFile $(BUILDDIR)/qthelp/kwant.qhc"

-latex:  $(GENERATEDPDF) $(FLAGS)
+latex:  $(GENERATEDPDF) $(FLAGS) $(INCLUDES)
 	$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
 	@echo
 	@echo "Build finished; the LaTeX files are in $(BUILDDIR)/latex."
@@ -124,19 +126,48 @@ doctest:
 %.pdf: %.svg
 	rsvg-convert -f pdf -o $@ $<

-# Make the image generation scripts by patching tutorial scipts.
+# Make tutorial scripts by extracting the (complete!) context of the "patches".
+# We make sure not to use 'wiggle' here.
 .SECONDARY:
-%.py: %.py.diff
-	@grep -v '^#HIDDEN' source/tutorial/$(notdir $@) >$@
-	@patch $@ $<
-
-# The image generation scripts depend on their unpatched originals
+source/tutorial/%.py: source/images/%.py.diff
+	@sed -n '/^[- ]/ s/^.//p' <$< >$@
+	@touch -r $< $@
+
+# Make the image generation scripts by patching tutorial scripts.  If the
+# tutorial scripts haven't been modified, don't patch but directly extract the
+# image generation scripts.  This means that 'wiggle' is only needed when the
+# tutorial scripts have been modified.
+.SECONDARY:
+source/images/%.py: source/tutorial/%.py
+	@if [ $< -nt $@.diff ]; then \
+	    cp $< $@; \
+	    rm -f $@.porig; \
+	    if ! wiggle --replace $@ $@.diff; then \
+	        command -v wiggle >/dev/null 2>&1 && \
+	        echo "Resolve conflicts by editing the files named below"; \
+	        touch -d@0 $@; \
+	        exit 1; \
+	    fi \
+	else \
+	    sed -n '/^[+ ]/ s/^.//p' <$@.diff >$@; \
+	    touch -r $@.diff $@; \
+	fi
+
+# Make the image generation scripts also depend on the diffs.
 define makedep
-source/images/$(1): source/tutorial/$(1)
+source/images/$(1): source/images/$(1).diff
 endef
 $(foreach name,$(SCRIPTS),$(eval $(call makedep,$(name))))

-# Generation of images
-.%_flag: %.py
+# Run an image generation script.  When successful, and if the script is newer
+# than the corresponding diff, recreate the latter.  Note that the
+# corresponding tutorial script cannot be newer, since if it is, the image
+# generation script is generated from it by patching.
+source/images/.%_flag: source/images/%.py
 	cd $(dir $<) && python3 $(notdir $<)
+	@if [ ! -f $<.diff -o $< -nt $<.diff ]; then \
+	    wiggle --diff --lines source/tutorial/$(notdir $<) $< >$<.diff; \
+	    touch -r $< $<.diff; \
+	fi
+	@rm -f $<.porig
 	@touch $@
--- a/doc/source/images/ab_ring.py.diff
+++ b/doc/source/images/ab_ring.py.diff
--- original
-+++ modified
-@@ -12,6 +12,7 @@
+@@ -1,127 +1,196 @@
+ # Tutorial 2.3.3. Nontrivial shapes
+ # =================================
+ #
+ # Physics background
+ # ------------------
+ #  Flux-dependent transmission through a quantum ring
+ #
+ # Kwant features highlighted
+ # --------------------------
+ #  - More complex shapes with lattices
+ #  - Allows for discussion of subtleties of `attach_lead` (not in the
 #    example, but in the tutorial main text)
 #  - Modifcations of hoppings/sites after they have been added
 
@@ -8,7 +17,79 @@
 from cmath import exp
 from math import pi
 
-@@ -81,6 +82,50 @@
+ import kwant
+ 
+ # For plotting
+ from matplotlib import pyplot
+ 
+ 
+ #HIDDEN_BEGIN_eusz
+ def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
+     # 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. ####
+     # Now, we aim for a more complex shape, namely a ring (or annulus)
+     def ring(pos):
+         (x, y) = pos
+         rsq = x ** 2 + y ** 2
+         return (r1 ** 2 < rsq < r2 ** 2)
+ #HIDDEN_END_eusz
+ 
+     # and add the corresponding lattice points using the `shape`-function
+ #HIDDEN_BEGIN_lcak
+     syst[lat.shape(ring, (0, r1 + 1))] = 4 * t
+     syst[lat.neighbors()] = -t
+ #HIDDEN_END_lcak
+ 
+     # In order to introduce a flux through the ring, we introduce a phase on
+     # the hoppings on the line cut through one of the arms.  Since we want to
+     # change the flux without modifying the Builder instance repeatedly, we
+     # define the modified hoppings as a function that takes the flux as its
+     # parameter phi.
+ #HIDDEN_BEGIN_lvkt
+     def hopping_phase(site1, site2, phi):
+         return -t * exp(1j * phi)
+ 
+     def crosses_branchcut(hop):
+         ix0, iy0 = hop[0].tag
+ 
+         # builder.HoppingKind with the argument (1, 0) below
+         # returns hoppings ordered as ((i+1, j), (i, j))
+         return iy0 < 0 and ix0 == 1  # ix1 == 0 then implied
+ 
+     # Modify only those hopings in x-direction that cross the branch cut
+     def hops_across_cut(syst):
+         for hop in kwant.builder.HoppingKind((1, 0), lat, lat)(syst):
+             if crosses_branchcut(hop):
+                 yield hop
+     syst[hops_across_cut] = hopping_phase
+ #HIDDEN_END_lvkt
+ 
+     #### Define the leads. ####
+     # left lead
+ #HIDDEN_BEGIN_qwgr
+     sym_lead = kwant.TranslationalSymmetry((-a, 0))
+     lead = kwant.Builder(sym_lead)
+ 
+     def lead_shape(pos):
+         (x, y) = pos
+         return (-W / 2 < y < W / 2)
+ 
+     lead[lat.shape(lead_shape, (0, 0))] = 4 * t
+     lead[lat.neighbors()] = -t
+ #HIDDEN_END_qwgr
+ 
+     #### Attach the leads and return the system. ####
+ #HIDDEN_BEGIN_skbz
+     syst.attach_lead(lead)
+     syst.attach_lead(lead.reversed())
+ #HIDDEN_END_skbz
+ 
     return syst
 
 
@@ -59,7 +140,9 @@
 def plot_conductance(syst, energy, fluxes):
     # compute conductance
 
-@@ -90,18 +135,31 @@
+     normalized_fluxes = [flux / (2 * pi) for flux in fluxes]
+     data = []
+     for flux in fluxes:
         smatrix = kwant.smatrix(syst, energy, args=[flux])
         data.append(smatrix.transmission(1, 0))
 
@@ -96,7 +179,9 @@
 
     # Finalize the system.
     syst = syst.finalized()
-@@ -111,6 +169,17 @@
+ 
+     # We should see a conductance that is periodic with the flux quantum
+     plot_conductance(syst, energy=0.15, fluxes=[0.01 * i * 3 * 2 * pi
                                                 for i in range(100)])
 
 
@@ -114,3 +199,4 @@
 # 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()
--- a/doc/source/images/band_structure.py.diff
+++ b/doc/source/images/band_structure.py.diff
--- original
-+++ modified
-@@ -9,6 +9,7 @@
+@@ -1,52 +1,62 @@
+ # Tutorial 2.4.1. Band structure calculations
+ # ===========================================
+ #
+ # Physics background
+ # ------------------
+ #  band structure of a simple quantum wire in tight-binding approximation
+ #
+ # Kwant features highlighted
 # --------------------------
 #  - Computing the band structure of a finalized lead.
 
@@ -8,8 +14,31 @@
 import kwant
 
 # For plotting.
-@@ -36,10 +37,19 @@
+ from matplotlib import pyplot
 
+ #HIDDEN_BEGIN_zxip
+ def make_lead(a=1, t=1.0, W=10):
+     # Start with an empty lead with a single square lattice
+     lat = kwant.lattice.square(a)
+ 
+     sym_lead = kwant.TranslationalSymmetry((-a, 0))
+     lead = kwant.Builder(sym_lead)
+ 
+     # build up one unit cell of the lead, and add the hoppings
+     # to the next unit cell
+     for j in range(W):
+         lead[lat(0, j)] = 4 * t
+ 
+         if j > 0:
+             lead[lat(0, j), lat(0, j - 1)] = -t
+ 
+         lead[lat(1, j), lat(0, j)] = -t
+ 
+     return lead
+ #HIDDEN_END_zxip
+ 
+ 
+ #HIDDEN_BEGIN_pejz
 def main():
     lead = make_lead().finalized()
 -    kwant.plotter.bands(lead, show=False)
@@ -29,6 +58,10 @@
 +    for extension in ('pdf', 'png'):
 +        fig.savefig("band_structure_result." + extension, dpi=_defs.dpi)
 +
+ #HIDDEN_END_pejz
 
 
 # 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()
--- a/doc/source/images/closed_system.py.diff
+++ b/doc/source/images/closed_system.py.diff
--- original
-+++ modified
-@@ -11,6 +11,7 @@
+@@ -1,140 +1,157 @@
+ # Tutorial 2.4.2. Closed systems
+ # ==============================
+ #
+ # Physics background
+ # ------------------
+ #  Fock-darwin spectrum of a quantum dot (energy spectrum in
+ #  as a function of a magnetic field)
+ #
+ # Kwant features highlighted
+ # --------------------------
 #  - Use of `hamiltonian_submatrix` in order to obtain a Hamiltonian
 #    matrix.
 
@@ -8,7 +16,62 @@
 from cmath import exp
 import numpy as np
 import kwant
-@@ -68,24 +69,39 @@
+ 
+ # For eigenvalue computation
+ #HIDDEN_BEGIN_tibv
+ import scipy.sparse.linalg as sla
+ #HIDDEN_END_tibv
+ 
+ # For plotting
+ from matplotlib import pyplot
+ 
+ 
+ def make_system(a=1, t=1.0, r=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).
+ 
+ #HIDDEN_BEGIN_qlyd
+     lat = kwant.lattice.square(a, norbs=1)
+ 
+     syst = kwant.Builder()
+ 
+     # Define the quantum dot
+     def circle(pos):
+         (x, y) = pos
+         rsq = x ** 2 + y ** 2
+         return rsq < r ** 2
+ 
+     def hopx(site1, site2, B=0):
+         # The magnetic field is controlled by the parameter B
+         y = site1.pos[1]
+         return -t * exp(-1j * B * y)
+ 
+     syst[lat.shape(circle, (0, 0))] = 4 * t
+     # hoppings in x-direction
+     syst[kwant.builder.HoppingKind((1, 0), lat, lat)] = hopx
+     # hoppings in y-directions
+     syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = -t
+ 
+     # It's a closed system for a change, so no leads
+     return syst
+ #HIDDEN_END_qlyd
+ 
+ 
+ #HIDDEN_BEGIN_yvri
+ def plot_spectrum(syst, Bfields):
+ 
+     # In the following, we compute the spectrum of the quantum dot
+     # using dense matrix methods. This works in this toy example, as
+     # the system is tiny. In a real example, one would want to use
+     # sparse matrix methods
+ 
+     energies = []
+     for B in Bfields:
+         # Obtain the Hamiltonian as a dense matrix
+         ham_mat = syst.hamiltonian_submatrix(args=[B], sparse=True)
+ 
+         # we only calculate the 15 lowest eigenvalues
+         ev = sla.eigsh(ham_mat, k=15, which='SM', return_eigenvectors=False)
 
         energies.append(ev)
 
@@ -29,8 +92,10 @@
 +    fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
 +    for extension in ('pdf', 'png'):
 +        fig.savefig("closed_system_result." + extension, dpi=_defs.dpi)
+ #HIDDEN_END_yvri
 
 
+ #HIDDEN_BEGIN_wave
 def plot_wave_function(syst):
 +    size = (_defs.figwidth_in, _defs.figwidth_in)
 +
@@ -46,15 +111,17 @@
 +            syst, np.abs(evecs[:, 9])**2, colorbar=False, oversampling=1,
 +            file="closed_system_eigenvector." + extension,
 +            fig_size=size, dpi=_defs.dpi)
+ #HIDDEN_END_wave
 
 
+ #HIDDEN_BEGIN_current
 def plot_current(syst):
 +    size = (_defs.figwidth_in, _defs.figwidth_in)
 +
     # Calculate the wave functions in the system.
     ham_mat = syst.hamiltonian_submatrix(sparse=True)
     evecs = sla.eigsh(ham_mat, k=20, which='SM')[1]
-@@ -93,15 +109,16 @@
+ 
     # Calculate and plot the local current of the 10th eigenmode.
     J = kwant.operator.Current(syst)
     current = J(evecs[:, 9])
@@ -64,6 +131,7 @@
 +            syst, current, colorbar=False,
 +            file="closed_system_current." + extension,
 +            fig_size=size, dpi=_defs.dpi)
+ #HIDDEN_END_current
 
 
 def main():
@@ -75,3 +143,26 @@
     # Finalize the system.
     syst = syst.finalized()
 
+     # The following try-clause can be removed once SciPy 0.9 becomes uncommon.
+     try:
+         # We should observe energy levels that flow towards Landau
+         # level energies with increasing magnetic field.
+         plot_spectrum(syst, [iB * 0.002 for iB in range(100)])
+ 
+         # Plot an eigenmode of a circular dot. Here we create a larger system for
+         # better spatial resolution.
+         syst = make_system(r=30).finalized()
+         plot_wave_function(syst)
+         plot_current(syst)
+     except ValueError as e:
+         if e.message == "Input matrix is not real-valued.":
+             print("The calculation of eigenvalues failed because of a bug in SciPy 0.9.")
+             print("Please upgrade to a newer version of SciPy.")
+         else:
+             raise
+ 
+ 
+ # 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()
--- a/doc/source/images/discretize.py.diff
+++ b/doc/source/images/discretize.py.diff
--- original
-+++ modified
-@@ -9,6 +9,7 @@
+@@ -1,222 +1,236 @@
+ # Tutorial 2.9. Processing continuum Hamiltonians with discretize
+ # ===============================================================
+ #
+ # Physics background
+ # ------------------
+ #  - tight-binding approximation of continuous Hamiltonians
+ #
+ # Kwant features highlighted
 # --------------------------
 #  - kwant.continuum.discretize
 
@@ -8,7 +14,11 @@
 
 import kwant
 import scipy.sparse.linalg
-@@ -20,10 +21,19 @@
+ import scipy.linalg
+ import numpy as np
+ 
+ # For plotting
+ import matplotlib as mpl
 from matplotlib import pyplot as plt
 
 
@@ -21,34 +31,120 @@
 +
 +
 def stadium_system(L=20, W=20):
+ #HIDDEN_BEGIN_template
     hamiltonian = "k_x**2 + k_y**2 + V(x, y)"
     template = kwant.continuum.discretize(hamiltonian)
 -    print(template)
 +    with open('discretizer_verbose.txt', 'w') as f:
 +        print(template, file=f)
+ #HIDDEN_END_template
 
+ #HIDDEN_BEGIN_fill
     def stadium(site):
         (x, y) = site.pos
-@@ -44,7 +54,7 @@
+         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()
+ #HIDDEN_END_fill
+     return syst
+ 
+ 
+ #HIDDEN_BEGIN_plot_eigenstate
+ 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.map(syst, abs(evecs[:, n])**2, show=False)
+ #HIDDEN_END_plot_eigenstate
 -    plt.show()
 +    save_figure('discretizer_gs')
 
 
 def qsh_system(a=20, L=2000, W=1000):
-@@ -91,7 +101,8 @@
+ #HIDDEN_BEGIN_define_qsh
+     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)
+     """
+ 
+     template = kwant.continuum.discretize(hamiltonian, grid_spacing=a)
+ #HIDDEN_END_define_qsh
+ 
+ #HIDDEN_BEGIN_define_qsh_build
+     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()
+ #HIDDEN_END_define_qsh_build
+     return syst
+ 
+ 
+ def analyze_qsh():
+     params = dict(A=3.65, B=-68.6, D=-51.1, M=-0.01, C=0)
+     syst = qsh_system()
+ 
+ #HIDDEN_BEGIN_plot_qsh_band
+     kwant.plotter.bands(syst.leads[0], params=params,
+                         momenta=np.linspace(-0.3, 0.3, 201), show=False)
+ #HIDDEN_END_plot_qsh_band
+ 
+     plt.grid()
+     plt.xlim(-.3, 0.3)
     plt.ylim(-0.05, 0.05)
     plt.xlabel('momentum [1/A]')
     plt.ylabel('energy [eV]')
 -    plt.show()
 +    save_figure('discretizer_qsh_band')
+ #HIDDEN_BEGIN_plot_qsh_wf
 +
     # get scattering wave functions at E=0
     wf = kwant.wave_function(syst, energy=0, params=params)
 
-@@ -119,7 +130,7 @@
+     # 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) = plt.subplots(1, 2, sharey=True, figsize=(16, 4))
+     kwant.plotter.map(syst, wf_sqr, ax=ax1)
+     kwant.plotter.map(syst, rho_sz, ax=ax2)
+ #HIDDEN_END_plot_qsh_wf
+     ax = ax1
+     im = [obj for obj in ax.get_children()
+           if isinstance(obj, mpl.image.AxesImage)][0]
+     fig.colorbar(im, ax=ax)
+ 
+     ax = ax2
+     im = [obj for obj in ax.get_children()
+           if isinstance(obj, mpl.image.AxesImage)][0]
+     fig.colorbar(im, ax=ax)
 
     ax1.set_title('Probability density')
     ax2.set_title('Spin density')
@@ -57,7 +153,46 @@
 
 
 def lattice_spacing():
-@@ -160,7 +171,7 @@
+ #HIDDEN_BEGIN_ls_def
+     hamiltonian = "k_x**2 * identity(2) + alpha * k_x * sigma_y"
+     params = dict(alpha=.5)
+ #HIDDEN_END_ls_def
+ 
+     def plot(ax, a=1):
+ #HIDDEN_BEGIN_ls_hk_cont
+         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]
+ #HIDDEN_END_ls_hk_cont
+ 
+ #HIDDEN_BEGIN_ls_hk_tb
+         template = kwant.continuum.discretize(hamiltonian, grid_spacing=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]
+ #HIDDEN_END_ls_hk_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()
+ 
+     _, (ax1, ax2) = plt.subplots(1, 2, sharey=True, figsize=(12, 4))
 
     plot(ax1, a=1)
     plot(ax2, a=.25)
@@ -66,25 +201,45 @@
 
 
 def substitutions():
-@@ -173,15 +184,18 @@
+ #HIDDEN_BEGIN_subs_1
+     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),
     )
 
 -    print(e[0] == e[1] == e[2])
 +    with open('discretizer_subs_1.txt', 'w') as f:
 +        print(e[0] == e[1] == e[2], file=f)
+ #HIDDEN_END_subs_1
 
+ #HIDDEN_BEGIN_subs_2
     subs = {'A': 'A(x) + B', 'V': 'V(x) + V_0', 'C': 5}
 -    print(sympify('k_x * A * k_x + V + C', locals=subs))
 +    with open('discretizer_subs_2.txt', 'w') as f:
 +        print(sympify('k_x * A * k_x + V + C', locals=subs), file=f)
+ #HIDDEN_END_subs_2
 
 
 def main():
+ #HIDDEN_BEGIN_symbolic_discretization
     template = kwant.continuum.discretize('k_x * A(x) * k_x')
 -    print(template)
 +    with open('discretizer_intro_verbose.txt', 'w') as f:
 +        print(template, file=f)
+ #HIDDEN_END_symbolic_discretization
 
     syst = stadium_system()
     plot_eigenstate(syst)
+ 
+     analyze_qsh()
+     lattice_spacing()
+     substitutions()
+ 
+ # 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()
--- a/doc/source/images/faq.py.diff
+++ b/doc/source/images/faq.py.diff
--- original
-+++ modified
-@@ -3,6 +3,7 @@
+@@ -1,450 +1,479 @@
+ # Frequently asked questions
+ #
 # This script is a disorganized collection of code snippets.  As a whole, it is
 # not meant as an example of good programming practice.
 
@@ -8,7 +8,8 @@
 import kwant
 import numpy as np
 import tinyarray
-@@ -11,6 +12,12 @@
+ import matplotlib
+ from matplotlib import pyplot as plt
 matplotlib.rcParams['figure.figsize'] = (3.5, 3.5)
 
 
@@ -20,32 +21,72 @@
 +
 #### What is a Site?
 
+ #HIDDEN_BEGIN_site
 a = 1
-@@ -21,6 +28,7 @@
+ lat = kwant.lattice.square(a)
+ syst = kwant.Builder()
+ 
+ syst[lat(1, 0)] = 4
 syst[lat(1, 1)] = 4
 
 kwant.plot(syst)
 +save_figure("faq_site")
+ #HIDDEN_END_site
 
 
 #### What is a hopping?
-@@ -34,6 +42,7 @@
+ 
+ a = 1
+ lat = kwant.lattice.square(a)
+ syst = kwant.Builder()
+ 
+ syst[lat(1, 0)] = 4
+ syst[lat(1, 1)] = 4
+ #HIDDEN_BEGIN_hopping
 syst[(lat(1, 0), lat(1, 1))] = 1j
+ #HIDDEN_END_hopping
 
 kwant.plot(syst)
 +save_figure("faq_hopping")
 
 
 #### What is a lattice?
-@@ -50,6 +59,7 @@
+ 
+ #HIDDEN_BEGIN_lattice_monatomic
+ # Two monatomic lattices
+ primitive_vectors = [(1, 0), (0, 1)]
+ lat_a = kwant.lattice.Monatomic(primitive_vectors, offset=(0, 0))
+ lat_b = kwant.lattice.Monatomic(primitive_vectors, offset=(0.5, 0.5))
+ # lat1 is equivalent to kwant.lattice.square()
+ 
+ syst = kwant.Builder()
+ 
+ syst[lat_a(0, 0)] = 4
 syst[lat_b(0, 0)] = 4
 
 kwant.plot(syst)
 +save_figure("faq_lattice")
+ #HIDDEN_END_lattice_monatomic
 
+ #HIDDEN_BEGIN_lattice_polyatomic
 # One polyatomic lattice containing two sublattices
 lat = kwant.lattice.Polyatomic([(1, 0), (0, 1)], [(0, 0), (0.5, 0.5)])
-@@ -70,6 +80,7 @@
+ sub_a, sub_b = lat.sublattices
+ #HIDDEN_END_lattice_polyatomic
+ 
+ #### How to make a hole in a system?
+ 
+ #HIDDEN_BEGIN_hole
+ # Define the lattice and the (empty) system
+ a = 2
+ lat = kwant.lattice.cubic(a)
+ 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)
@@ -53,47 +94,143 @@
 
 # Delete sites to create a hole
 
-@@ -81,6 +92,7 @@
+ 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)
 +save_figure("faq_hole2")
+ #HIDDEN_END_hole
 
 
 #### How can we get access to the sites of our system?
-@@ -127,7 +139,7 @@
+ 
+ builder = kwant.Builder()
+ lat = kwant.lattice.square()
+ builder[(lat(i, j) for i in range(3) for j in range(3))] = 4
+ #HIDDEN_BEGIN_sites1
+ # Before finalizing the system
+ 
+ sites = list(builder.sites())  # sites() doe *not* return a list
+ #HIDDEN_END_sites1
+ #HIDDEN_BEGIN_sites2
+ # After finalizing the system
+ syst = builder.finalized()
+ sites = syst.sites  # syst.sites is an actual list
+ #HIDDEN_END_sites2
+ #HIDDEN_BEGIN_sites3
+ i = syst.id_by_site[lat(0, 2)]  # we want the id of the site lat(0, 2)
+ #HIDDEN_END_sites3
+ 
+ 
+ #### How to plot a polyatomic lattice with different colors?
+ 
+ #HIDDEN_BEGIN_colors1
+ lat = kwant.lattice.kagome()
+ syst = kwant.Builder()
+ 
+ a, b, c = lat.sublattices  # The kagome lattice has 3 sublattices
+ #HIDDEN_END_colors1
+ 
+ #HIDDEN_BEGIN_colors2
+ # 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
 
 ## Plot the system.
 plot_system(syst)
-
 +save_figure("faq_colors")
+ #HIDDEN_END_colors2
 
+-
 #### How to create all hoppings in a given direction using Hoppingkind?
 
-@@ -141,6 +153,7 @@
+ # Monatomic lattice
+ 
+ #HIDDEN_BEGIN_direction1
+ # Create hopping between neighbors with HoppingKind
+ a = 1
+ syst = kwant.Builder()
+ lat = kwant.lattice.square(a)
+ 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)
 +save_figure("faq_direction1")
+ #HIDDEN_END_direction1
 
 # Polyatomic lattice
 
-@@ -171,12 +184,14 @@
+ lat = kwant.lattice.kagome()
+ syst = kwant.Builder()
+ 
+ a, b, c = lat.sublattices  # The kagome lattice has 3 sublattices
+ 
+ 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_size=0.15, site_lw=0.05, site_color=family_color)
+ 
+ 
+ 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
+ 
+ 
+ #HIDDEN_BEGIN_direction2
 # equivalent to syst[kwant.builder.HoppingKind((0, 1), b)] = -1
 syst[kwant.builder.HoppingKind((0, 1), b, b)] = -1
+ #HIDDEN_END_direction2
 plot_system(syst)
 +save_figure("faq_direction2")
 # Delete the hoppings previously created
 del syst[kwant.builder.HoppingKind((0, 1), b, b)]
+ #HIDDEN_BEGIN_direction3
 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
+ #HIDDEN_END_direction3
 plot_system(syst)
 +save_figure("faq_direction3")
 
 
 #### How to create the hoppings between adjacent sites?
-@@ -190,11 +205,13 @@
+ 
+ # Monatomic lattice
+ 
+ #HIDDEN_BEGIN_adjacent1
+ # Create hoppings with lat.neighbors()
+ syst = kwant.Builder()
+ lat = kwant.lattice.square()
+ syst[(lat(i, j) for i in range(3) for j in range(3))] = 4
 
 syst[lat.neighbors()] = -1  # Equivalent to lat.neighbors(1)
 kwant.plot(syst)
@@ -104,66 +241,175 @@
 
 kwant.plot(syst)
 +save_figure("faq_adjacent2")
+ #HIDDEN_END_adjacent1
 
 # Polyatomic lattice
 
-@@ -213,9 +230,11 @@
+ #HIDDEN_BEGIN_FAQ6
+ 
+ # Hoppings using .neighbors()
+ #HIDDEN_BEGIN_adjacent2
+ # Create the system
+ lat = kwant.lattice.kagome()
+ 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  # red
+         syst[b(i, j)] = 4  # green
+         syst[c(i, j)] = 4  # blue
 
 syst[lat.neighbors()] = -1
+ #HIDDEN_END_adjacent2
 plot_system(syst)
 +save_figure("faq_adjacent3")
 del syst[lat.neighbors()]  # Delete the hoppings previously created
+ #HIDDEN_BEGIN_adjacent3
 syst[a.neighbors()] = -1
+ #HIDDEN_END_adjacent3
 plot_system(syst)
 +save_figure("faq_adjacent4")
 del syst[a.neighbors()]  # Delete the hoppings previously created
 
 
-@@ -253,6 +272,7 @@
+ syst[lat.neighbors(2)] = -1
+ plot_system(syst)
+ del syst[lat.neighbors(2)]
+ 
+ 
+ #### How to create a lead with a lattice different from the scattering region?
+ 
+ # Plot sites from different families in different colors
+ 
+ def plot_system(syst):
+ 
+     def family_color(site):
+         if site.family == subA:
+             return 'blue'
+         if site.family == subB:
+             return 'yellow'
+         else:
+             return 'green'
+ 
+     kwant.plot(syst, site_lw=0.1, site_color=family_color)
+ 
+ 
+ #HIDDEN_BEGIN_different_lattice1
+ # Define the scattering Region
+ L = 5
+ W = 5
+ 
+ lat = kwant.lattice.honeycomb()
+ subA, subB = lat.sublattices
+ 
+ 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
+ #HIDDEN_END_different_lattice1
 plot_system(syst)
 +save_figure("faq_different_lattice1")
 
+ #HIDDEN_BEGIN_different_lattice2
 # Create a lead
 lat_lead = kwant.lattice.square()
-@@ -262,6 +282,7 @@
+ sym_lead1 = kwant.TranslationalSymmetry((0, 1))
+ 
+ lead1 = kwant.Builder(sym_lead1)
 lead1[(lat_lead(i, 0) for i in range(2, 7))] = 4
 lead1[lat_lead.neighbors()] = -1
+ #HIDDEN_END_different_lattice2
 plot_system(lead1)
 +save_figure("faq_different_lattice2")
 
+ #HIDDEN_BEGIN_different_lattice3
 syst[(lat_lead(i, 5) for i in range(2, 7))] = 4
 syst[lat_lead.neighbors()] = -1
-@@ -269,9 +290,11 @@
+ 
 # Manually attach sites from graphene to square lattice
 syst[((lat_lead(i+2, 5), subB(i, 4)) for i in range(5))] = -1
+ #HIDDEN_END_different_lattice3
 plot_system(syst)
 +save_figure("faq_different_lattice3")
 
+ #HIDDEN_BEGIN_different_lattice4
 syst.attach_lead(lead1)
+ #HIDDEN_END_different_lattice4
 plot_system(syst)
 +save_figure("faq_different_lattice4")
 
 
 #### How to cut a finite system out of a system with translationnal symmetries?
-@@ -291,6 +314,7 @@
+ 
+ #HIDDEN_BEGIN_fill1
+ # Create 3d model.
+ cubic = kwant.lattice.cubic()
+ 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
+ #HIDDEN_END_fill1
+ 
+ #HIDDEN_BEGIN_fill2
+ # Build scattering region (white).
+ def cuboid_shape(site):
+     x, y, z = abs(site.pos)
+     return x < 4 and y < 10 and z < 3
+ 
 cuboid = kwant.Builder()
 cuboid.fill(model, cuboid_shape, (0, 0, 0));
+ #HIDDEN_END_fill2
 kwant.plot(cuboid);
 +save_figure("faq_fill2")
 
+ #HIDDEN_BEGIN_fill3
 # Build electrode (black).
 def electrode_shape(site):
-@@ -305,6 +329,7 @@
+     x, y, z = site.pos - (0, 5, 2)
+     return y**2 + z**2 < 2.3**2
+ 
+ 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)
+ #HIDDEN_END_fill3
 kwant.plot(cuboid);
 +save_figure("faq_fill3")
 
 
 #### How does Kwant order the propagating modes of a lead?
-@@ -337,6 +362,7 @@
+ 
+ #HIDDEN_BEGIN_pm
+ lat = kwant.lattice.square()
+ 
+ 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)
+ #HIDDEN_END_pm
+ 
+ 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)
@@ -171,7 +417,14 @@
 
 # More involved example
 
-@@ -351,6 +377,7 @@
+ 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)
@@ -179,20 +432,52 @@
 
 
 #### How does Kwant order components of an individual wavefunction?
-@@ -381,7 +408,8 @@
+ 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()
+ 
+ 
+ #HIDDEN_BEGIN_ord1
+ 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])
 +with open('faq_ord1.txt', 'w') as f:
 +    print('wavefunction on lat(0, 0): ', wf[idx], file=f)
+ #HIDDEN_END_ord1
 
+ #HIDDEN_BEGIN_ord2
 lat = kwant.lattice.square(norbs=2)
 syst = make_system(lat)
-@@ -394,4 +422,5 @@
+ 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])
 +with open('faq_ord2.txt', 'w') as f:
 +    print('wavefunction on lat(0, 0): ', wf[idx], file=f)
+ #HIDDEN_END_ord2
--- a/doc/source/images/generate-diffs.sh
+++ b/doc/source/images/generate-diffs.sh
-# !/bin/sh
-
-# This script regenerates the .diff files in this directory.  It's these files
-# that are kept under vesion control instead of the scripts themselves.
-
-for f in [a-zA-Z]*.py; do
-    # We use custom labels to suppress the time stamps which are unnecessary
-    # here and would only lead to noise in version control.
-    grep -v '#HIDDEN' ../tutorial/$f |
-    diff -u --label original --label modified - $f >$f.diff_
-    if cmp $f.diff_ $f.diff >/dev/null 2>&1; then
-        rm $f.diff_
-    else
-        mv $f.diff_ $f.diff
-    fi
-done
--- a/doc/source/images/graphene.py.diff
+++ b/doc/source/images/graphene.py.diff
--- original
-+++ modified
-@@ -10,6 +10,7 @@
+@@ -1,179 +1,203 @@
+ # Tutorial 2.5. Beyond square lattices: graphene
+ # ==============================================
+ #
+ # Physics background
+ # ------------------
+ #  Transport through a graphene quantum dot with a pn-junction
+ #
+ # Kwant features highlighted
+ # --------------------------
 #  - Application of all the aspects of tutorials 1-3 to a more complicated
 #    lattice, namely graphene
 
@@ -8,7 +15,102 @@
 from math import pi, sqrt, tanh
 
 import kwant
-@@ -96,22 +97,40 @@
+ 
+ # For computing eigenvalues
+ import scipy.sparse.linalg as sla
+ 
+ # For plotting
+ from matplotlib import pyplot
+ 
+ 
+ # Define the graphene lattice
+ sin_30, cos_30 = (1 / 2, sqrt(3) / 2)
+ #HIDDEN_BEGIN_hnla
+ graphene = kwant.lattice.general([(1, 0), (sin_30, cos_30)],
+                                  [(0, 0), (0, 1 / sqrt(3))])
+ a, b = graphene.sublattices
+ #HIDDEN_END_hnla
+ 
+ 
+ #HIDDEN_BEGIN_shzy
+ def make_system(r=10, w=2.0, pot=0.1):
+ 
+     #### Define the scattering region. ####
+     # circular scattering region
+     def circle(pos):
+         x, y = pos
+         return x ** 2 + y ** 2 < r ** 2
+ 
+     syst = kwant.Builder()
+ 
+     # w: width and pot: potential maximum of the p-n junction
+     def potential(site):
+         (x, y) = site.pos
+         d = y * cos_30 + x * sin_30
+         return pot * tanh(d / w)
+ 
+     syst[graphene.shape(circle, (0, 0))] = potential
+ #HIDDEN_END_shzy
+ 
+     # specify the hoppings of the graphene lattice in the
+     # format expected by builder.HoppingKind
+ #HIDDEN_BEGIN_hsmc
+     hoppings = (((0, 0), a, b), ((0, 1), a, b), ((-1, 1), a, b))
+ #HIDDEN_END_hsmc
+ #HIDDEN_BEGIN_bfwb
+     syst[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
+ #HIDDEN_END_bfwb
+ 
+     # Modify the scattering region
+ #HIDDEN_BEGIN_efut
+     del syst[a(0, 0)]
+     syst[a(-2, 1), b(2, 2)] = -1
+ #HIDDEN_END_efut
+ 
+     #### Define the leads. ####
+ #HIDDEN_BEGIN_aakh
+     # left lead
+     sym0 = kwant.TranslationalSymmetry(graphene.vec((-1, 0)))
+ 
+     def lead0_shape(pos):
+         x, y = pos
+         return (-0.4 * r < y < 0.4 * r)
+ 
+     lead0 = kwant.Builder(sym0)
+     lead0[graphene.shape(lead0_shape, (0, 0))] = -pot
+     lead0[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
+ 
+     # The second lead, going to the top right
+     sym1 = kwant.TranslationalSymmetry(graphene.vec((0, 1)))
+ 
+     def lead1_shape(pos):
+         v = pos[1] * sin_30 - pos[0] * cos_30
+         return (-0.4 * r < v < 0.4 * r)
+ 
+     lead1 = kwant.Builder(sym1)
+     lead1[graphene.shape(lead1_shape, (0, 0))] = pot
+     lead1[[kwant.builder.HoppingKind(*hopping) for hopping in hoppings]] = -1
+ #HIDDEN_END_aakh
+ 
+ #HIDDEN_BEGIN_kmmw
+     return syst, [lead0, lead1]
+ #HIDDEN_END_kmmw
+ 
+ 
+ #HIDDEN_BEGIN_zydk
+ def compute_evs(syst):
+     # Compute some eigenvalues of the closed system
+     sparse_mat = syst.hamiltonian_submatrix(sparse=True)
+ 
+     evs = sla.eigs(sparse_mat, 2)[0]
+     print(evs.real)
+ #HIDDEN_END_zydk
+ 
+ 
+ def plot_conductance(syst, energies):
+     # Compute transmission as a function of energy
+     data = []
+     for energy in energies:
         smatrix = kwant.smatrix(syst, energy)
         data.append(smatrix.transmission(0, 1))
 
@@ -56,8 +158,17 @@
 +        fig.savefig("graphene_bs." + extension, dpi=_defs.dpi)
 
 
+ #HIDDEN The part of the following code block which begins with family_colors
+ #HIDDEN is included verbatim in the tutorial text because nested code examples
+ #HIDDEN are not supported.  Remember to update the tutorial text when you
+ #HIDDEN modify this block.
+ #HIDDEN_BEGIN_itkk
 def main():
-@@ -123,8 +142,11 @@
+     pot = 0.1
+     syst, leads = make_system(pot=pot)
+ 
+     # To highlight the two sublattices of graphene, we plot one with
+     # a filled, and the other one with an open circle:
     def family_colors(site):
         return 0 if site.family == a else 1
 
@@ -68,10 +179,14 @@
 +        kwant.plot(syst, site_color=family_colors, site_lw=0.1, colorbar=False,
 +                   file="graphene_syst1." + extension,
 +                   fig_size=size, dpi=_defs.dpi)
+ #HIDDEN_END_itkk
 
     # Compute some eigenvalues.
+ #HIDDEN_BEGIN_jmbi
     compute_evs(syst.finalized())
-@@ -133,9 +155,11 @@
+ #HIDDEN_END_jmbi
+ 
+     # Attach the leads to the system.
     for lead in leads:
         syst.attach_lead(lead)
 
@@ -86,3 +201,17 @@
 
     # Finalize the system.
     syst = syst.finalized()
+ 
+     # Compute the band structure of lead 0.
+     momenta = [-pi + 0.02 * pi * i for i in range(101)]
+     plot_bandstructure(syst.leads[0], momenta)
+ 
+     # Plot conductance.
+     energies = [-2 * pot + 4. / 50. * pot * i for i in range(51)]
+     plot_conductance(syst, energies)
+ 
+ 
+ # 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()
--- a/doc/source/images/kernel_polynomial_method.py.diff
+++ b/doc/source/images/kernel_polynomial_method.py.diff
--- original
-+++ modified
-@@ -11,6 +11,9 @@
+@@ -1,219 +1,242 @@
+ # Tutorial 2.8. Calculating spectral density with the Kernel Polynomial Method
+ # ============================================================================
+ #
+ # Physics background
+ # ------------------
+ #  - Chebyshev polynomials, random trace approximation, spectral densities.
+ #
+ # Kwant features highlighted
+ # --------------------------
+ #  - kpm module,kwant operators.
 
 import scipy
 
@@ -10,7 +18,27 @@
 # For plotting
 from matplotlib import pyplot as plt
 
-@@ -36,13 +39,13 @@
+ #HIDDEN_BEGIN_sys1
+ # necessary imports
+ import kwant
+ import numpy as np
+ 
+ 
+ # define the system
+ def make_syst(r=30, t=-1, a=1):
+     syst = kwant.Builder()
+     lat = kwant.lattice.honeycomb(a, norbs=1)
+ 
+     def circle(pos):
+         x, y = pos
+         return x ** 2 + y ** 2 < r ** 2
+ 
+     syst[lat.shape(circle, (0, 0))] = 0.
+     syst[lat.neighbors()] = t
+     syst.eradicate_dangling()
+ 
+     return syst
+ #HIDDEN_END_sys1
 
 
 # Plot several density of states curves on the same axes.
@@ -26,7 +54,14 @@
     plt.clf()
 
 
-@@ -57,10 +60,18 @@
+ def site_size_conversion(densities):
+     return 3 * np.abs(densities) / max(densities)
+ 
+ 
+ # Plot several local density of states maps in different subplots
+ def plot_ldos(fsyst, axes, titles_to_data, file_name=None):
+     for ax, (title, ldos) in zip(axes, titles_to_data):
+         site_size = site_size_conversion(ldos)  # convert LDoS to sizes
         kwant.plot(fsyst, site_size=site_size, site_color=(0, 0, 1, 0.3), ax=ax)
         ax.set_title(title)
         ax.set(adjustable='box-forced', aspect='equal')
@@ -44,9 +79,21 @@
 +
 +
 def simple_dos_example():
+ #HIDDEN_BEGIN_kpm1
     fsyst = make_syst().finalized()
 
-@@ -74,18 +85,24 @@
+     spectrum = kwant.kpm.SpectralDensity(fsyst)
+ #HIDDEN_END_kpm1
+ 
+ #HIDDEN_BEGIN_kpm2
+     energies, densities = spectrum()
+ #HIDDEN_END_kpm2
+ 
+ #HIDDEN_BEGIN_kpm3
+     energy_subset = np.linspace(0, 2)
+     density_subset = spectrum(energy_subset)
+ #HIDDEN_END_kpm3
+ 
     plot_dos([
         ('densities', (energies, densities)),
         ('density subset', (energy_subset, density_subset)),
@@ -59,11 +106,14 @@
 def dos_integrating_example(fsyst):
     spectrum = kwant.kpm.SpectralDensity(fsyst)
 
+ #HIDDEN_BEGIN_int1
 -    print('identity resolution:', spectrum.integrate())
 +    with open('kpm_normalization.txt', 'w') as f:
 +        with redirect_stdout(f):
 +            print('identity resolution:', spectrum.integrate())
+ #HIDDEN_END_int1
 
+ #HIDDEN_BEGIN_int2
     # Fermi energy 0.1 and temperature 0.2
     fermi = lambda E: 1 / (np.exp((E - 0.1) / 0.2) + 1)
 
@@ -71,10 +121,19 @@
 +    with open('kpm_total_states.txt', 'w') as f:
 +        with redirect_stdout(f):
 +            print('number of filled states:', spectrum.integrate(fermi))
+ #HIDDEN_END_int2
 
 
 def increasing_accuracy_example(fsyst):
-@@ -99,7 +116,9 @@
+     spectrum = kwant.kpm.SpectralDensity(fsyst)
+     original_dos = spectrum()  # get unaltered DoS
+ 
+ #HIDDEN_BEGIN_acc1
+     spectrum.add_moments(energy_resolution=0.03)
+ #HIDDEN_END_acc1
+ 
+     increased_resolution_dos = spectrum()
+ 
     plot_dos([
         ('density', original_dos),
         ('higher energy resolution', increased_resolution_dos),
@@ -83,9 +142,13 @@
 +     file_name='kpm_dos_acc'
 +    )
 
+ #HIDDEN_BEGIN_acc2
     spectrum.add_moments(100)
     spectrum.add_vectors(5)
-@@ -109,7 +128,9 @@
+ #HIDDEN_END_acc2
+ 
+     increased_moments_dos = spectrum()
+ 
     plot_dos([
         ('density', original_dos),
         ('higher number of moments', increased_moments_dos),
@@ -96,7 +159,36 @@
 
 
 def operator_example(fsyst):
-@@ -139,7 +160,9 @@
+ #HIDDEN_BEGIN_op1
+     # identity matrix
+     matrix_op = scipy.sparse.eye(len(fsyst.sites))
+     matrix_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=matrix_op)
+ #HIDDEN_END_op1
+ 
+ #HIDDEN_BEGIN_op2
+     # 'sum=True' means we sum over all the sites
+     kwant_op = kwant.operator.Density(fsyst, sum=True)
+     operator_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=kwant_op)
+ #HIDDEN_END_op2
+ 
+     plot_dos([
+         ('identity matrix', matrix_spectrum()),
+         ('kwant.operator.Density', operator_spectrum()),
+     ])
+ 
+ 
+ def ldos_example(fsyst):
+ #HIDDEN_BEGIN_op3
+     # 'sum=False' is the default, but we include it explicitly here for clarity.
+     kwant_op = kwant.operator.Density(fsyst, sum=False)
+     local_dos = kwant.kpm.SpectralDensity(fsyst, operator=kwant_op)
+ #HIDDEN_END_op3
+ 
+ #HIDDEN_BEGIN_op4
+     zero_energy_ldos = local_dos(energy=0)
+     finite_energy_ldos = local_dos(energy=1)
+ 
+     _, axes = plt.subplots(1, 2, figsize=(12, 7))
     plot_ldos(fsyst, axes,[
         ('energy = 0', zero_energy_ldos),
         ('energy = 1', finite_energy_ldos),
@@ -104,6 +196,57 @@
 +     ],
 +     file_name='kpm_ldos'
 +    )
+ #HIDDEN_END_op4
 
 
 def vector_factory_example(fsyst):
+     spectrum = kwant.kpm.SpectralDensity(fsyst)
+ #HIDDEN_BEGIN_fact1
+     # construct vectors with n random elements -1 or +1.
+     def binary_vectors(n):
+         return np.rint(np.random.random_sample(n)) * 2 - 1
+ 
+     custom_factory = kwant.kpm.SpectralDensity(fsyst,
+                                                vector_factory=binary_vectors)
+ #HIDDEN_END_fact1
+     plot_dos([
+         ('default vector factory', spectrum()),
+         ('binary vector factory', custom_factory()),
+     ])
+ 
+ 
+ def bilinear_map_operator_example(fsyst):
+ #HIDDEN_BEGIN_blm
+     rho = kwant.operator.Density(fsyst, sum=True)
+ 
+     # sesquilinear map that does the same thing as `rho`
+     def rho_alt(bra, ket):
+         return np.vdot(bra, ket)
+ 
+     rho_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=rho)
+     rho_alt_spectrum = kwant.kpm.SpectralDensity(fsyst, operator=rho_alt)
+ #HIDDEN_END_blm
+ 
+     plot_dos([
+         ('kwant.operator.Density', rho_spectrum()),
+         ('bilinear operator', rho_alt_spectrum()),
+     ])
+ 
+ 
+ def main():
+     simple_dos_example()
+ 
+     fsyst = make_syst().finalized()
+ 
+     dos_integrating_example(fsyst)
+     increasing_accuracy_example(fsyst)
+     operator_example(fsyst)
+     ldos_example(fsyst)
+     vector_factory_example(fsyst)
+     bilinear_map_operator_example(fsyst)
+ 
+ 
+ # 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()
--- a/doc/source/images/magnetic_texture.py.diff
+++ b/doc/source/images/magnetic_texture.py.diff
--- original
-+++ modified
-@@ -11,6 +11,9 @@
+@@ -1,245 +1,268 @@
+ # Tutorial 2.7. Spin textures
+ # ===========================
+ #
+ # Physics background
+ # ------------------
+ #  - Spin textures
+ #  - Skyrmions
+ #
+ # Kwant features highlighted
+ # --------------------------
 #  - operators
 #  - plotting vector fields
 
@@ -10,8 +18,73 @@
 from math import sin, cos, tanh, pi
 import itertools
 import numpy as np
-@@ -79,7 +82,13 @@
+ import tinyarray as ta
+ import matplotlib.pyplot as plt
+ 
+ import kwant
+ 
+ sigma_0 = ta.array([[1, 0], [0, 1]])
+ sigma_x = ta.array([[0, 1], [1, 0]])
+ sigma_y = ta.array([[0, -1j], [1j, 0]])
+ sigma_z = ta.array([[1, 0], [0, -1]])
+ 
+ # vector of Pauli matrices σ_αiβ where greek
+ # letters denote spinor indices
+ sigma = np.rollaxis(np.array([sigma_x, sigma_y, sigma_z]), 1)
+ 
+ #HIDDEN_BEGIN_model
+ def field_direction(pos, r0, delta):
+     x, y = pos
+     r = np.linalg.norm(pos)
+     r_tilde = (r - r0) / delta
+     theta = (tanh(r_tilde) - 1) * (pi / 2)
+ 
+     if r == 0:
+         m_i = [0, 0, 1]
+     else:
+         m_i = [
+             (x / r) * sin(theta),
+             (y / r) * sin(theta),
+             cos(theta),
+         ]
+ 
+     return np.array(m_i)
+ 
+ 
+ def scattering_onsite(site, r0, delta, J):
+     m_i = field_direction(site.pos, r0, delta)
+     return J * np.dot(m_i, sigma)
+ 
+ 
+ def lead_onsite(site, J):
+     return J * sigma_z
+ #HIDDEN_END_model
+ 
+ 
+ #HIDDEN_BEGIN_syst
+ lat = kwant.lattice.square(norbs=2)
+ 
+ def make_system(L=80):
+ 
+     syst = kwant.Builder()
+ 
+     def square(pos):
+         return all(-L/2 < p < L/2 for p in pos)
+ 
+     syst[lat.shape(square, (0, 0))] = scattering_onsite
+     syst[lat.neighbors()] = -sigma_0
+ 
+     lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)),
+                          conservation_law=-sigma_z)
+ 
+     lead[lat.shape(square, (0, 0))] = lead_onsite
+     lead[lat.neighbors()] = -sigma_0
+ 
+     syst.attach_lead(lead)
+     syst.attach_lead(lead.reversed())
+ 
     return syst
+ #HIDDEN_END_syst
 
 
 -def plot_vector_field(syst, params):
@@ -25,7 +98,8 @@
     xmin, ymin = min(s.tag for s in syst.sites)
     xmax, ymax = max(s.tag for s in syst.sites)
     x, y = np.meshgrid(np.arange(xmin, xmax+1), np.arange(ymin, ymax+1))
-@@ -88,28 +97,38 @@
+ 
+     m_i = [field_direction(p, **params) for p in zip(x.flat, y.flat)]
     m_i = np.reshape(m_i, x.shape + (3,))
     m_i = np.rollaxis(m_i, 2, 0)
 
@@ -52,11 +126,11 @@
 +        save_plot(fname)
 +    else:
 +        plt.show()
-+
 
 
 -def plot_currents(syst, currents):
 -    fig, axes = plt.subplots(1, len(currents))
+
 +def plot_currents(syst, currents, fname=None):
 +    fig, axes = plt.subplots(1, len(currents), figsize=(7*len(currents), 7))
     if not hasattr(axes, '__len__'):
@@ -72,34 +146,80 @@
 
 
 def main():
-@@ -119,7 +138,7 @@
+     syst = make_system().finalized()
+ 
+ #HIDDEN_BEGIN_wavefunction
+     params = dict(r0=20, delta=10, J=1)
     wf = kwant.wave_function(syst, energy=-1, params=params)
     psi = wf(0)[0]
+ #HIDDEN_END_wavefunction
 
 -    plot_vector_field(syst, dict(r0=20, delta=10))
 +    plot_vector_field(syst, dict(r0=20, delta=10), fname='mag_field_direction')
 
+ #HIDDEN_BEGIN_ldos
     # even (odd) indices correspond to spin up (down)
     up, down = psi[::2], psi[1::2]
-@@ -144,7 +163,7 @@
+     density = np.abs(up)**2 + np.abs(down)**2
+ #HIDDEN_END_ldos
+ 
+ #HIDDEN_BEGIN_lsdz
+     # spin down components have a minus sign
+     spin_z = np.abs(up)**2 - np.abs(down)**2
+ #HIDDEN_END_lsdz
+ 
+ #HIDDEN_BEGIN_lsdy
+     # spin down components have a minus sign
+     spin_y = 1j * (down.conjugate() * up - up.conjugate() * down)
+ #HIDDEN_END_lsdy
+ 
+ #HIDDEN_BEGIN_lden
+     rho = kwant.operator.Density(syst)
+     rho_sz = kwant.operator.Density(syst, sigma_z)
+     rho_sy = kwant.operator.Density(syst, sigma_y)
+ 
+     # calculate the expectation values of the operators with 'psi'
+     density = rho(psi)
+     spin_z = rho_sz(psi)
+     spin_y = rho_sy(psi)
+ #HIDDEN_END_lden
+ 
+     plot_densities(syst, [
         ('$σ_0$', density),
         ('$σ_z$', spin_z),
         ('$σ_y$', spin_y),
 -    ])
 +    ], fname='spin_densities')
 
+ #HIDDEN_BEGIN_current
     J_0 = kwant.operator.Current(syst)
     J_z = kwant.operator.Current(syst, sigma_z)
-@@ -159,7 +178,7 @@
+     J_y = kwant.operator.Current(syst, sigma_y)
+ 
+     # calculate the expectation values of the operators with 'psi'
+     current = J_0(psi)
+     spin_z_current = J_z(psi)
+     spin_y_current = J_y(psi)
+ #HIDDEN_END_current
+ 
+     plot_currents(syst, [
         ('$J_{σ_0}$', current),
         ('$J_{σ_z}$', spin_z_current),
         ('$J_{σ_y}$', spin_y_current),
 -    ])
 +    ], fname='spin_currents')
 
+ #HIDDEN_BEGIN_following
     def following_m_i(site, r0, delta):
         m_i = field_direction(site.pos, r0, delta)
-@@ -173,7 +192,7 @@
+         return np.dot(m_i, sigma)
+ 
+     J_m = kwant.operator.Current(syst, following_m_i)
+ 
+     # evaluate the operator
+     m_current = J_m(psi, params=dict(r0=25, delta=10))
+ #HIDDEN_END_following
+ 
     plot_currents(syst, [
         (r'$J_{\mathbf{m}_i}$', m_current),
         ('$J_{σ_z}$', spin_z_current),
@@ -107,8 +227,11 @@
 +    ], fname='spin_current_comparison')
 
 
+ #HIDDEN_BEGIN_density_cut
     def circle(site):
-@@ -183,7 +202,9 @@
+         return np.linalg.norm(site.pos) < 20
+ 
+     rho_circle = kwant.operator.Density(syst, where=circle, sum=True)
 
     all_states = np.vstack((wf(0), wf(1)))
     dos_in_circle = sum(rho_circle(p) for p in all_states) / (2 * pi)
@@ -116,10 +239,18 @@
 +    with open('circle_dos.txt', 'w') as f:
 +        with redirect_stdout(f):
 +            print('density of states in circle:', dos_in_circle)
+ #HIDDEN_END_density_cut
 
+ #HIDDEN_BEGIN_current_cut
     def left_cut(site_to, site_from):
         return site_from.pos[0] <= -39 and site_to.pos[0] > -39
-@@ -197,8 +218,10 @@
+ 
+     def right_cut(site_to, site_from):
+         return site_from.pos[0] < 39 and site_to.pos[0] >= 39
+ 
+     J_left = kwant.operator.Current(syst, where=left_cut, sum=True)
+     J_right = kwant.operator.Current(syst, where=right_cut, sum=True)
+ 
     Jz_left = kwant.operator.Current(syst, sigma_z, where=left_cut, sum=True)
     Jz_right = kwant.operator.Current(syst, sigma_z, where=right_cut, sum=True)
 
@@ -129,10 +260,21 @@
 +        with redirect_stdout(f):
 +            print('J_left:', J_left(psi), ' J_right:', J_right(psi))
 +            print('Jz_left:', Jz_left(psi), ' Jz_right:', Jz_right(psi))
+ #HIDDEN_END_current_cut
 
+ #HIDDEN_BEGIN_bind
     J_m = kwant.operator.Current(syst, following_m_i)
     J_z = kwant.operator.Current(syst, sigma_z)
-@@ -214,7 +237,7 @@
+ 
+     J_m_bound = J_m.bind(params=dict(r0=25, delta=10, J=1))
+     J_z_bound = J_z.bind(params=dict(r0=25, delta=10, J=1))
+ 
+     # Sum current local from all scattering states on the left at energy=-1
+     wf_left = wf(0)
+     J_m_left = sum(J_m_bound(p) for p in wf_left)
+     J_z_left = sum(J_z_bound(p) for p in wf_left)
+ #HIDDEN_END_bind
+ 
     plot_currents(syst, [
         (r'$J_{\mathbf{m}_i}$ (from left)', J_m_left),
         (r'$J_{σ_z}$ (from left)', J_z_left),
@@ -141,3 +283,4 @@
 
 
 if __name__ == '__main__':
+     main()
--- a/doc/source/images/plot_graphene.py.diff
+++ b/doc/source/images/plot_graphene.py.diff
--- original
-+++ modified
-@@ -9,6 +9,7 @@
+@@ -1,112 +1,130 @@
+ # Tutorial 2.8.1. 2D example: graphene quantum dot
+ # ================================================
+ #
+ # Physics background
+ # ------------------
+ #  - graphene edge states
+ #
+ # Kwant features highlighted
 # --------------------------
 #  - demonstrate different ways of plotting
 
@@ -8,7 +14,14 @@
 import kwant
 from matplotlib import pyplot
 
-@@ -22,7 +23,7 @@
+ #HIDDEN_BEGIN_makesyst
+ lat = kwant.lattice.honeycomb()
+ a, b = lat.sublattices
+ 
+ def make_system(r=8, t=-1, tp=-0.1):
+ 
+     def circle(pos):
+         x, y = pos
         return x**2 + y**2 < r**2
 
     syst = kwant.Builder()
@@ -17,11 +30,16 @@
     syst[lat.neighbors()] = t
     syst.eradicate_dangling()
     if tp:
-@@ -32,9 +33,11 @@
+         syst[lat.neighbors(2)] = tp
 
+     return syst
+ #HIDDEN_END_makesyst
 
+ 
+ #HIDDEN_BEGIN_plotsyst1
 def plot_system(syst):
 -    kwant.plot(syst)
+ #HIDDEN_END_plotsyst1
 -    # the standard plot is ok, but not very intelligible. One can do
 -    # better by playing wioth colors and linewidths
 +    # standard plot - not very intelligible for this particular situation
@@ -31,8 +49,10 @@
 +                   fig_size=size, dpi=_defs.dpi)
 
     # use color and linewidths to get a better plot
+ #HIDDEN_BEGIN_plotsyst2
     def family_color(site):
-@@ -43,7 +46,11 @@
+         return 'black' if site.family == a else 'white'
+ 
     def hopping_lw(site1, site2):
         return 0.04 if site1.family == site2.family else 0.1
 
@@ -42,23 +62,38 @@
 +        kwant.plot(syst, site_lw=0.1, site_color=family_color,
 +                   hop_lw=hopping_lw, file="plot_graphene_syst2." + extension,
 +                   fig_size=size, dpi=_defs.dpi)
+ #HIDDEN_END_plotsyst2
 
 
+ #HIDDEN_BEGIN_plotdata1
 def plot_data(syst, n):
-@@ -58,7 +65,11 @@
+     import scipy.linalg as la
+ 
+     syst = syst.finalized()
+     ham = syst.hamiltonian_submatrix()
+     evecs = la.eigh(ham)[1]
+ 
+     wf = abs(evecs[:, n])**2
+ #HIDDEN_END_plotdata1
+ 
     # the usual - works great in general, looks just a bit crufty for
     # small systems
 
+ #HIDDEN_BEGIN_plotdata2
 -    kwant.plotter.map(syst, wf, oversampling=10, cmap='gist_heat_r')
 +    size = (_defs.figwidth_in, _defs.figwidth_in)
 +    for extension in ('pdf', 'png'):
 +        kwant.plotter.map(syst, wf, oversampling=10, cmap='gist_heat_r',
 +                          file="plot_graphene_data1." + extension,
 +                          fig_size=size, dpi=_defs.dpi)
+ #HIDDEN_END_plotdata2
 
     # use two different sort of triangles to cleanly fill the space
+ #HIDDEN_BEGIN_plotdata3
     def family_shape(i):
-@@ -68,15 +79,22 @@
+         site = syst.sites[i]
+         return ('p', 3, 180) if site.family == a else ('p', 3, 0)
+ 
     def family_color(i):
         return 'black' if syst.sites[i].family == a else 'white'
 
@@ -70,8 +105,10 @@
 +                   site_size=0.5, hop_lw=0, cmap='gist_heat_r',
 +                   file="plot_graphene_data2." + extension,
 +                   fig_size=size, dpi=_defs.dpi)
+ #HIDDEN_END_plotdata3
 
     # plot by changing the symbols itself
+ #HIDDEN_BEGIN_plotdata4
     def site_size(i):
         return 3 * wf[i] / wf.max()
 
@@ -82,6 +119,23 @@
 +        kwant.plot(syst, site_size=site_size, site_color=(0,0,1,0.3),
 +                   hop_lw=0.1, file="plot_graphene_data3." + extension,
 +                   fig_size=size, dpi=_defs.dpi)
+ #HIDDEN_END_plotdata4
 
 
 def main():
+     # plot the graphene system in different styles
+     syst = make_system()
+ 
+     plot_system(syst)
+ 
+     # compute a wavefunction (number 225) and plot it in different
+     # styles
+     syst = make_system(tp=0)
+ 
+     plot_data(syst, 225)
+ 
+ 
+ # 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()
--- a/doc/source/images/plot_qahe.py.diff
+++ b/doc/source/images/plot_qahe.py.diff
--- original
-+++ modified
-@@ -11,6 +11,7 @@
+@@ -1,71 +1,75 @@
+ # 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
 
@@ -8,8 +16,53 @@
 import math
 import matplotlib.pyplot
 import kwant
-@@ -60,7 +61,10 @@
+ 
+ 
+ # 2 band model exhibiting quantum anomalous Hall effect
+ #HIDDEN_BEGIN_model
+ 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_spacing=a)
+ #HIDDEN_END_model
+ 
+ 
+ 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()
+ 
+ 
+ def main():
+     # 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)
+     kwant.plot(syst)
+ 
+     # 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)
+ #HIDDEN_BEGIN_current
     J = kwant.operator.Current(syst).bind(params=params)
     current = sum(J(p) for p in psi)
 -    kwant.plotter.current(syst, current)
@@ -17,6 +70,8 @@
 +        kwant.plotter.current(syst, current,
 +                              file="plot_qahe_current." + extension,
 +                              dpi=_defs.dpi)
+ #HIDDEN_END_current
 
 
 if __name__ == '__main__':
+     main()
--- a/doc/source/images/plot_zincblende.py.diff
+++ b/doc/source/images/plot_zincblende.py.diff
--- original
-+++ modified
-@@ -9,6 +9,7 @@
+@@ -1,59 +1,67 @@
+ # Tutorial 2.8.2. 3D example: zincblende structure
+ # ================================================
+ #
+ # Physical background
+ # -------------------
+ #  - 3D Bravais lattices
+ #
+ # Kwant features highlighted
 # --------------------------
 #  - demonstrate different ways of plotting in 3D
 
@@ -8,8 +14,30 @@
 import kwant
 from matplotlib import pyplot
 
-@@ -33,7 +34,10 @@
+ #HIDDEN_BEGIN_zincblende1
+ lat = kwant.lattice.general([(0, 0.5, 0.5), (0.5, 0, 0.5), (0.5, 0.5, 0)],
+                             [(0, 0, 0), (0.25, 0.25, 0.25)])
+ a, b = lat.sublattices
+ #HIDDEN_END_zincblende1
+ 
+ #HIDDEN_BEGIN_zincblende2
+ def make_cuboid(a=15, b=10, c=5):
+     def cuboid_shape(pos):
+         x, y, z = pos
+         return 0 <= x < a and 0 <= y < b and 0 <= z < c
+ 
+     syst = kwant.Builder()
+     syst[lat.shape(cuboid_shape, (0, 0, 0))] = None
+     syst[lat.neighbors()] = None
+ 
+     return syst
+ #HIDDEN_END_zincblende2
+ 
+ 
+ def main():
+     # the standard plotting style for 3D is mainly useful for
     # checking shapes:
+ #HIDDEN_BEGIN_plot1
     syst = make_cuboid()
 
 -    kwant.plot(syst)
@@ -17,10 +45,12 @@
 +    for extension in ('pdf', 'png'):
 +        kwant.plot(syst, file="plot_zincblende_syst1." + extension,
 +                   fig_size=size, dpi=_defs.dpi)
+ #HIDDEN_END_plot1
 
     # visualize the crystal structure better for a very small system
+ #HIDDEN_BEGIN_plot2
     syst = make_cuboid(a=1.5, b=1.5, c=1.5)
-@@ -41,8 +45,12 @@
+ 
     def family_colors(site):
         return 'r' if site.family == a else 'g'
 
@@ -32,6 +62,10 @@
 +                   site_color=family_colors,
 +                   file="plot_zincblende_syst2." + extension,
 +                   fig_size=size, dpi=_defs.dpi)
+ #HIDDEN_END_plot2
 
 
 # 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()
--- a/doc/source/images/quantum_well.py.diff
+++ b/doc/source/images/quantum_well.py.diff
--- original
-+++ modified
-@@ -9,6 +9,7 @@
+@@ -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
 
@@ -8,7 +14,51 @@
 import kwant
 
 # For plotting
-@@ -55,19 +56,25 @@
+ 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=0):
+         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, args=[-welldepth])
         data.append(smatrix.transmission(1, 0))
 
@@ -30,6 +80,7 @@
 +    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():
@@ -41,3 +92,12 @@
     # 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()
--- a/doc/source/images/quantum_wire.py.diff
+++ b/doc/source/images/quantum_wire.py.diff
--- original
-+++ modified
-@@ -11,6 +11,7 @@
+@@ -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
 
-@@ -69,7 +70,10 @@
 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()
-@@ -90,8 +94,13 @@
+ #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)
@@ -38,3 +133,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
--- a/doc/source/images/spin_orbit.py.diff
+++ b/doc/source/images/spin_orbit.py.diff
--- original
-+++ modified
-@@ -13,6 +13,7 @@
+@@ -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
 
@@ -8,7 +18,65 @@
 import kwant
 
 # For plotting
-@@ -70,19 +71,24 @@
+ 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(a=1, 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(a)
+ 
+     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
+     # hoppings in y-directions
+     syst[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
+         -t * sigma_0 + 1j * alpha * sigma_x
+ #HIDDEN_END_uxrm
+ 
+     #### Define the left lead. ####
+     lead = kwant.Builder(kwant.TranslationalSymmetry((-a, 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
+     # hoppings in y-directions
+     lead[kwant.builder.HoppingKind((0, 1), lat, lat)] = \
+         -t * sigma_0 + 1j * alpha * sigma_x
+ #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))
 
@@ -40,3 +108,11 @@
     # 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()
--- a/doc/source/images/superconductor.py.diff
+++ b/doc/source/images/superconductor.py.diff
--- original
-+++ modified
-@@ -11,6 +11,7 @@
+@@ -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.
 
@@ -8,7 +16,7 @@
 import kwant
 
 import tinyarray
-@@ -18,6 +19,7 @@
+ import numpy as np
 
 # For plotting
 from matplotlib import pyplot
@@ -16,10 +24,67 @@
 
 tau_x = tinyarray.array([[0, 1], [1, 0]])
 tau_y = tinyarray.array([[0, -1j], [1j, 0]])
-@@ -74,11 +76,19 @@
+ 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, phs=True):
+     # 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)
@@ -36,9 +101,28 @@
 +        fig.savefig("superconductor_transport_result." + extension,
 +                    dpi=_defs.dpi)
 
+ #HIDDEN_BEGIN_pqmp
 def check_PHS(syst):
     # Scattering matrix
-@@ -103,14 +113,13 @@
+     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)
 
@@ -56,3 +140,9 @@
 
     # 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()
--- a/doc/source/tutorial/README
+++ b/doc/source/tutorial/README
-Note for Kwant developers
-------------------------
-
-Say you have modified SCRIPT.py in this directory.  Make sure that the tutorial
-image generation patches still apply by deleting doc/source/images/SCRIPT.py
-and running ``make doc/source/images/SCRIPT.py`` in doc.  Now examine the newly
-created doc/source/images/SCRIPT.py.  If you do not like the result or the
-patch did not apply, edit doc/source/images/SCRIPT.py until you like it.  You
-can run `make html` during your edits to check things.  Finally, even if you
-did not edit the script, execute generate-diffs.sh in doc/source/images.  If
-the patches applied cleanly the diff files will usually stay unchanged.