add 2-band model example and include links in whatsnew

doc/source/images/plot_qahe.py.diff

+--- original
++++ modified
+@@ -11,6 +11,7 @@
+ # + Use of `kwant.operator` to compute local current
+ # + Use of `kwant.plotter.current` to plot local current
+ 
++import _defs
+ import math
+ import matplotlib.pyplot
+ import kwant
+@@ -60,7 +61,10 @@
+     psi = kwant.wave_function(syst, energy=params['m'], params=params)(0)
+     J = kwant.operator.Current(syst).bind(params=params)
+     current = sum(J(p) for p in psi)
+-    kwant.plotter.current(syst, current)
++    for extension in ('pdf', 'png'):
++        kwant.plotter.current(syst, current,
++                              file="plot_qahe_current." + extension,
++                              dpi=_defs.dpi)
+ 
+ 
+ if __name__ == '__main__':

doc/source/pre/whatsnew/1.3.rst

@@ -24,9 +24,18 @@ with continuum models and for discretizing them into tight-binding models. It
 aims at providing a handy interface to convert symbolic Hamiltonians into a
 builder with N-D translational symmetry that can be use to calculate
 tight-binding band structures or construct systems with different/lower
-symmetry.
+symmetry. For example in just a few lines we can construct a two-band model that exhibits
+a quantum anomalous spin Hall phase:
 
-See :doc:`../../tutorial/discretize`
+.. literalinclude:: ../../tutorial/plot_qahe.py
+    :start-after: HIDDEN_BEGIN_model
+    :end-before: HIDDEN_END_model
+
+From: :download:`QAHE example script <../../tutorial/plot_qahe.py>`
+
+See the tutorial: :doc:`../../tutorial/discretize`
+
+See the reference documentation: :doc:`../../reference/kwant.continuum`
 
 Calculating charges and currents using the operator module
 ----------------------------------------------------------
@@ -57,10 +66,17 @@ Quantities defined on system hoppings (e.g. currents calculated using
 `~kwant.operator.current`) can be directly plotted as a streamplot over the
 system using `kwant.plotter.current`. This is similar to how
 `kwant.plotter.map` can be used to plot quantities defined on sites.
+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:: ../../tutorial/plot_qahe.py
+    :start-after: HIDDEN_BEGIN_current
+    :end-before: HIDDEN_END_current
 
-See :doc:`../../tutorial/plotting`
+.. image:: ../../images/plot_qahe_current.*
 
-.. image:: ../../images/plot_qpc_current.*
+From: :download:`QAHE example script <../../tutorial/plot_qahe.py>`
 
 Scattering states with discrete symmetries and conservation laws
 ----------------------------------------------------------------

doc/source/tutorial/plot_qahe.py

+# 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
+
+
+# 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)
+#HIDDEN_END_current
+
+
+if __name__ == '__main__':
+    main()