doc: add examples of current plotting to tutorials

doc/source/images/closed_system.py.diff

@@ -8,7 +8,7 @@
  from cmath import exp
  import numpy as np
  import kwant
-@@ -68,29 +69,39 @@
+@@ -68,40 +69,56 @@
  
          energies.append(ev)
  
@@ -48,6 +48,24 @@
 +            fig_size=size, dpi=_defs.dpi)
  
  
+ 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]
+ 
+     # Calculate and plot the local current of the 10th eigenmode.
+     J = kwant.operator.Current(syst)
+     current = J(evecs[:, 9])
+-    kwant.plotter.current(syst, current, colorbar=False)
++    for extension in ('pdf', 'png'):
++        kwant.plotter.current(
++            syst, current, colorbar=False,
++            file="closed_system_current." + extension,
++            fig_size=size, dpi=_defs.dpi)
+ 
+ 
  def main():
      syst = make_system()
  

doc/source/images/plot_qpc.py.diff

+--- original
++++ modified
+@@ -1,11 +1,19 @@
+ import math
+ import cmath
+ 
++import _defs
+ import matplotlib.pyplot as plt
+ 
+ import kwant
+ 
+ 
++figsize = (_defs.figwidth_in, _defs.figwidth_in)
++def save_plot(fname):
++    for extension in ('.pdf', '.png'):
++        plt.savefig(fname + extension,
++                    size=figsize, dpi=_defs.dpi, bbox_inches='tight')
++
++
+ def hopping(sitei, sitej, phi):
+     xi, yi = sitei.pos
+     xj, yj = sitej.pos
+@@ -42,7 +50,9 @@
+ def main():
+     syst = make_system()
+ 
+-    kwant.plot(syst)
++    fig, ax = plt.subplots(1, 1)
++    kwant.plot(syst, ax=ax)
++    save_plot('plot_qpc_syst')
+ 
+     params = dict(phi=1/40, salt="")
+     psi = kwant.wave_function(syst, energy=0.15, params=params)(0)
+@@ -52,9 +62,13 @@
+     density = sum(D(p) for p in psi)
+     current = sum(J(p) for p in psi)
+ 
+-    kwant.plotter.map(syst, density, method='linear')
+-
+-    kwant.plotter.current(syst, current)
++    fig, ax = plt.subplots(1, 1, figsize=figsize)
++    kwant.plotter.map(syst, density, method='linear', ax=ax)
++    save_plot('plot_qpc_density')
++
++    fig, ax = plt.subplots(1, 1, figsize=figsize)
++    kwant.plotter.current(syst, current, ax=ax)
++    save_plot('plot_qpc_current')
+ 
+ 
+ if __name__ == '__main__':

doc/source/tutorial/closed_system.py

@@ -29,7 +29,7 @@ def make_system(a=1, t=1.0, r=10):
     # `a` is the lattice constant (by default set to 1 for simplicity).
 
 #HIDDEN_BEGIN_qlyd
-    lat = kwant.lattice.square(a)
+    lat = kwant.lattice.square(a, norbs=1)
 
     syst = kwant.Builder()
 
@@ -93,6 +93,19 @@ def plot_wave_function(syst):
 #HIDDEN_END_wave
 
 
+#HIDDEN_BEGIN_current
+def plot_current(syst):
+    # Calculate the wave functions in the system.
+    ham_mat = syst.hamiltonian_submatrix(sparse=True)
+    evecs = sla.eigsh(ham_mat, k=20, which='SM')[1]
+
+    # Calculate and plot the local current of the 10th eigenmode.
+    J = kwant.operator.Current(syst)
+    current = J(evecs[:, 9])
+    kwant.plotter.current(syst, current, colorbar=False)
+#HIDDEN_END_current
+
+
 def main():
     syst = make_system()
 
@@ -112,6 +125,7 @@ def main():
         # 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.")

doc/source/tutorial/plot_qpc.py

+import math
+import cmath
+
+import matplotlib.pyplot as plt
+
+import kwant
+
+
+#HIDDEN_BEGIN_syst
+def hopping(sitei, sitej, phi):
+    xi, yi = sitei.pos
+    xj, yj = sitej.pos
+    return -cmath.exp(-0.5j * phi * (xi - xj) * (yi + yj))
+
+
+def onsite(site, salt):
+    return 0.05 * kwant.digest.gauss(site.tag, salt) + 4
+
+
+def make_system(L=50):
+    def central_region(pos):
+        x, y = pos
+        return -2*L < x < 2*L and \
+            abs(y) < L - 37.5 * math.exp(-x**2 / 12**2)
+
+    lat = kwant.lattice.square(norbs=1)
+    syst = kwant.Builder()
+
+    syst[lat.shape(central_region, (0, 0))] = onsite
+    syst[lat.neighbors()] = hopping
+
+    sym = kwant.TranslationalSymmetry((-1, 0))
+    lead = kwant.Builder(sym)
+    lead[(lat(0, y) for y in range(-L + 1, L))] = 4
+    lead[lat.neighbors()] = hopping
+
+    syst.attach_lead(lead)
+    syst.attach_lead(lead.reversed())
+
+    return syst.finalized()
+#HIDDEN_END_syst
+
+
+def main():
+    syst = make_system()
+
+    kwant.plot(syst)
+
+#HIDDEN_BEGIN_wf
+    params = dict(phi=1/40, salt="")
+    psi = kwant.wave_function(syst, energy=0.15, params=params)(0)
+    J = kwant.operator.Current(syst).bind(params=params)
+    D = kwant.operator.Density(syst).bind(params=params)
+    # Calculate density and current due to modes from the left lead
+    density = sum(D(p) for p in psi)
+    current = sum(J(p) for p in psi)
+#HIDDEN_END_wf
+
+#HIDDEN_BEGIN_density
+    kwant.plotter.map(syst, density, method='linear')
+#HIDDEN_END_density
+
+#HIDDEN_BEGIN_current
+    kwant.plotter.current(syst, current)
+#HIDDEN_END_current
+
+
+if __name__ == '__main__':
+    main()

doc/source/tutorial/tutorial3.rst

@@ -120,6 +120,19 @@ The last two arguments to `~kwant.plotter.map` are optional.  The first prevents
 a colorbar from appearing.  The second, ``oversampling=1``, makes the image look
 better for the special case of a square lattice.
 
+
+As our model breaks time reversal symmetry (because of the applied magnetic
+field) we can also see an intereting property of the eigenstates, namely
+that they can have *non-zero* local current. We can calculate the local
+current due to a state by using `kwant.operator.Current` and plotting
+it using `kwant.plotter.current`:
+
+.. literalinclude:: closed_system.py
+    :start-after: #HIDDEN_BEGIN_current
+    :end-before: #HIDDEN_END_current
+
+.. image:: ../images/closed_system_current.*
+
 .. specialnote:: Technical details
 
   - `~kwant.system.System.hamiltonian_submatrix` can also return a sparse

doc/source/tutorial/tutorial6.rst

@@ -226,3 +226,61 @@ crystal lattices out there!
       3d module)
     - Plotting hoppings in 3D is inherently much slower than plotting sites.
       Hence, this is not done by default.
+
+
+Interpolated density and current: QPC with disorder
+...................................................
+
+.. seealso::
+    The complete source code of this example can be found in
+    :download:`tutorial/plot_qpc.py <../../../tutorial/plot_qpc.py>`
+
+In the above examples we saw some useful methods for plotting systems where
+single-site resolution is required. Sometimes, however, having single-site
+precision is a hinderance, rather than a help, and looking at *averaged*
+quantities is more useful. This is particularly important in systems with a
+large number of sites, and systems that are discretizations of continuum
+models.
+
+Here we will show how to plot interpolated quantities using `kwant.plotter.map`
+and `kwant.plotter.current` using the example of a quantum point contact (QPC)
+with a perpendicular magnetic field and disorder:
+
+.. literalinclude:: plot_qpc.py
+    :start-after: #HIDDEN_BEGIN_syst
+    :end-before: #HIDDEN_END_syst
+
+.. image:: ../images/plot_qpc_syst.*
+
+Now we will compute the density of particles and current due to states
+originating in the left lead with energy 0.15.
+
+.. literalinclude:: plot_qpc.py
+    :start-after: #HIDDEN_BEGIN_wf
+    :end-before: #HIDDEN_END_wf
+
+We can then plot the density using `~kwant.plotter.map`:
+
+.. literalinclude:: plot_qpc.py
+    :start-after: #HIDDEN_BEGIN_density
+    :end-before: #HIDDEN_END_density
+
+.. image:: ../images/plot_qpc_density.*
+
+We pass ``method='linear'`` to ``map``, which produces a smoother plot than the
+default style. We see that density is concentrated on the edges of the sample,
+as we expect due to Hall effect induced by the perpendicular magnetic field.
+
+Plotting the current in the system will enable us to make even more sense
+of what is going on:
+
+.. literalinclude:: plot_qpc.py
+    :start-after: #HIDDEN_BEGIN_current
+    :end-before: #HIDDEN_END_current
+
+.. image:: ../images/plot_qpc_current.*
+
+We can now clearly see that current enters the system from the lower-left edge
+of the system (this matches our intuition, as we calculated the current for
+scattering states originating in the left-hand lead), and is backscattered by
+the restriction of the QPC in the centre.