doc: replace the figure creation scripts

...by ones closely based on the actual tutorial examples.

doc/source/images/tutorial2a.py

@@ -11,113 +11,96 @@
 #  - Numpy matrices as values in Builder
 
 import kwant
+
+# For plotting
+import pylab
+
+# For matrix support
 import numpy
 
 import latex, html
 
-# define sigma-matrices for convenience
+# define Pauli-matrices for convenience
 sigma_0 = numpy.eye(2)
 sigma_x = numpy.array([[0, 1], [1, 0]])
 sigma_y = numpy.array([[0, -1j], [1j, 0]])
 sigma_z = numpy.array([[1, 0], [0, -1]])
 
-# First, define the tight-binding system
-
-sys = kwant.Builder()
-
-# Here, we are only working with square lattices
-
-# for simplicity, take lattice constant = 1
-a = 1
-lat = kwant.lattice.Square(a)
-
-t = 1.0
-alpha = 0.5
-e_z = 0.08
-W = 10
-L = 30
-
-# Define the scattering region
-
-def rectangle(pos):
-    (x, y) = pos
-    return ( -0.5 < x < L - 0.5 ) and ( -0.5 < y < W - 0.5 )
-
-sys[lat.shape(rectangle, (0, 0))] = 4 * t * sigma_0 + e_z * sigma_z
-# hoppings in x-direction
-sys[sys.possible_hoppings((1, 0), lat, lat)] = - t * sigma_0 - \
-    1j * alpha * sigma_y
-# hoppings in y-directions
-sys[sys.possible_hoppings((0, 1), lat, lat)] = - t * sigma_0 + \
-    1j * alpha * sigma_x
-
-# Then, define the leads:
-
-# First the lead to the left
-
-# (Note: in the current version, TranslationalSymmetry takes a
-# realspace vector)
-sym_lead0 = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
-lead0 = kwant.Builder(sym_lead0)
-lead0.default_site_group = lat
-
-def lead_shape(pos):
-    (x, y) = pos
-    return (-1 < x < 1) and ( -0.5 < y < W - 0.5 )
-
-lead0[lat.shape(lead_shape, (0, 0))] = 4 * t * sigma_0 + e_z * sigma_z
-# hoppings in x-direction
-lead0[lead0.possible_hoppings((1, 0), lat, lat)] = - t * sigma_0 - \
-    1j * alpha * sigma_y
-# hoppings in y-directions
-lead0[lead0.possible_hoppings((0, 1), lat, lat)] = - t * sigma_0 + \
-    1j * alpha * sigma_x
-
-# Then the lead to the right
-# there we can use a special function that simply reverses the direction
-
-lead1 = lead0.reversed()
-
-# Then attach the leads to the system
-
-sys.attach_lead(lead0)
-sys.attach_lead(lead1)
-
-# finalize the system
-
-fsys = sys.finalized()
-
-# Now that we have the system, we can compute conductance
-
-energies = []
-data = []
-for ie in xrange(100):
-    energy = ie * 0.01 - 0.3
-
-    # compute the scattering matrix at energy energy
-    smatrix = kwant.solvers.sparse.solve(fsys, energy)
-
-    # compute the transmission probability from lead 0 to
-    # lead 1
-    energies.append(energy)
-    data.append(smatrix.transmission(1, 0))
-
-# Use matplotlib to write output
-# We should see conductance steps
-import pylab
 
-pylab.plot(energies, data)
-pylab.xlabel("energy [in units of t]",
-             fontsize=latex.mpl_label_size)
-pylab.ylabel("conductance [in units of e^2/h]",
-             fontsize=latex.mpl_label_size)
-fig = pylab.gcf()
-pylab.setp(fig.get_axes()[0].get_xticklabels(),
-           fontsize=latex.mpl_tick_size)
-pylab.setp(fig.get_axes()[0].get_yticklabels(),
-           fontsize=latex.mpl_tick_size)
-fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
-fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-fig.savefig("tutorial2a_result.pdf")
-fig.savefig("tutorial2a_result.png",
-            dpi=(html.figwidth_px/latex.mpl_width_in))
+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)
+
+    sys = kwant.Builder()
+    sys.default_site_group = lat
+
+    #### Define the scattering region. ####
+    sys[((x, y) for x in range(L) for y in range(W))] = 4 * t * sigma_0 + \
+        e_z * sigma_z
+    # hoppings in x-direction
+    sys[sys.possible_hoppings((1, 0), lat, lat)] = - t * sigma_0 - \
+        1j * alpha * sigma_y
+    # hoppings in y-directions
+    sys[sys.possible_hoppings((0, 1), lat, lat)] = - t * sigma_0 + \
+        1j * alpha * sigma_x
+
+    #### Define the leads. ####
+    # left lead
+    sym_lead0 = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
+    lead0 = kwant.Builder(sym_lead0)
+    lead0.default_site_group = lat
+
+    lead0[((0, j) for j in xrange(W))] = 4 * t * sigma_0 + e_z * sigma_z
+    # hoppings in x-direction
+    lead0[lead0.possible_hoppings((1, 0), lat, lat)] = - t * sigma_0 - \
+        1j * alpha * sigma_y
+    # hoppings in y-directions
+    lead0[lead0.possible_hoppings((0, 1), lat, lat)] = - t * sigma_0 + \
+        1j * alpha * sigma_x
+
+    # Then the lead to the right
+    # (again, obtained using reverse()
+    lead1 = lead0.reversed()
+
+    #### Attach the leads and return the finalized system. ####
+    sys.attach_lead(lead0)
+    sys.attach_lead(lead1)
+
+    return sys.finalized()
+
+def plot_conductance(fsys, energies):
+    # Compute conductance
+    data = []
+    for energy in energies:
+        smatrix = kwant.solve(fsys, energy)
+        data.append(smatrix.transmission(1, 0))
+
+    pylab.plot(energies, data)
+    pylab.xlabel("energy [in units of t]",
+                 fontsize=latex.mpl_label_size)
+    pylab.ylabel("conductance [in units of e^2/h]",
+                 fontsize=latex.mpl_label_size)
+    fig = pylab.gcf()
+    pylab.setp(fig.get_axes()[0].get_xticklabels(),
+               fontsize=latex.mpl_tick_size)
+    pylab.setp(fig.get_axes()[0].get_yticklabels(),
+               fontsize=latex.mpl_tick_size)
+    fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
+    fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
+    fig.savefig("tutorial2a_result.pdf")
+    fig.savefig("tutorial2a_result.png",
+                dpi=(html.figwidth_px/latex.mpl_width_in))
+
+
+def main():
+    fsys = make_system()
+
+    # We should see non-monotonic conductance steps.
+    plot_conductance(fsys, energies=[0.01 * i - 0.3 for i in xrange(100)])
+
+
+# Call the main function if the script gets executed (as opposed to imported).
+# See <http://docs.python.org/library/__main__.html>.
+if __name__ == '__main__':
+    main()

doc/source/images/tutorial2b.py

@@ -8,108 +8,102 @@
 
 import kwant
 
-import latex, html
-
-# First, define the tight-binding system
-
-sys = kwant.Builder()
-
-# Here, we are only working with square lattices
-
-# for simplicity, take lattice constant = 1
-a = 1
-lat = kwant.lattice.Square(a)
-
-t = 1.0
-alpha = 0.5
-e_z = 0.08
-W = 10
-L = 30
-
-# Define the scattering region
-
-def rectangle(pos):
-    (x, y) = pos
-    return ( -0.5 < x < L - 0.5 ) and ( -0.5 < y < W - 0.5 )
-
-def potential(site):
-    (x, y) = site.pos
-    if 10 < x < 20:
-        return pot
-    else:
-        return 0
-
-def onsite(site):
-    return 4 * t + potential(site)
-
-sys[lat.shape(rectangle, (0, 0))] = onsite
-for hopping in lat.nearest:
-    sys[sys.possible_hoppings(*hopping)] = - t
-
-# Then, define the leads:
-
-# First the lead to the left
-
-# (Note: in the current version, TranslationalSymmetry takes a
-# realspace vector)
-sym_lead0 = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
-lead0 = kwant.Builder(sym_lead0)
-lead0.default_site_group = lat
-
-def lead_shape(pos):
-    (x, y) = pos
-    return (-1 < x < 1) and ( -0.5 < y < W - 0.5 )
-
-lead0[lat.shape(lead_shape, (0, 0))] = 4 * t
-for hopping in lat.nearest:
-    lead0[lead0.possible_hoppings(*hopping)] = - t
-
-# Then the lead to the right
-# there we can use a special function that simply reverses the direction
-
-lead1 = lead0.reversed()
-
-# Then attach the leads to the system
-
-sys.attach_lead(lead0)
-sys.attach_lead(lead1)
-
-# finalize the system
-
-fsys = sys.finalized()
-
-# Now that we have the system, we can compute conductance
-
-energy = 0.2
-wellpot = []
-data = []
-for ipot in xrange(100):
-    pot = - ipot * 0.01
-
-    # compute the scattering matrix at energy energy
-    smatrix = kwant.solvers.sparse.solve(fsys, energy)
-
-    # compute the transmission probability from lead 0 to
-    # lead 1
-    wellpot.append(-pot)
-    data.append(smatrix.transmission(1, 0))
-
-# Use matplotlib to write output
-# We should see conductance steps
+# For plotting
 import pylab
 
-pylab.plot(wellpot, data)
-pylab.xlabel("well depth [in units of t]",
-             fontsize=latex.mpl_label_size)
-pylab.ylabel("conductance [in units of e^2/h]",
-             fontsize=latex.mpl_label_size)
-fig = pylab.gcf()
-pylab.setp(fig.get_axes()[0].get_xticklabels(),
-           fontsize=latex.mpl_tick_size)
-pylab.setp(fig.get_axes()[0].get_yticklabels(),
-           fontsize=latex.mpl_tick_size)
-fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
-fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-fig.savefig("tutorial2b_result.pdf")
-fig.savefig("tutorial2b_result.png",
-            dpi=(html.figwidth_px/latex.mpl_width_in))
+import latex, html
+
+# global variable governing the behavior of potential() in
+# make_system()
+pot = 0
+
+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)
+
+    sys = kwant.Builder()
+    sys.default_site_group = lat
+
+    #### Define the scattering region. ####
+    # Potential profile
+    def potential(site):
+        (x, y) = site.pos
+        if (L - L_well) / 2 < x < (L + L_well) / 2:
+            # The potential value is provided using a global variable
+            return pot
+        else:
+            return 0
+
+    def onsite(site):
+        return 4 * t + potential(site)
+
+    sys[((x, y) for x in range(L) for y in range(W))] = onsite
+    for hopping in lat.nearest:
+        sys[sys.possible_hoppings(*hopping)] = -t
+
+    #### Define the leads. ####
+    # First the lead to the left, ...
+    # (Note: in the current version, TranslationalSymmetry takes a
+    # realspace vector)
+    sym_lead0 = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
+    lead0 = kwant.Builder(sym_lead0)
+    lead0.default_site_group = lat
+
+    lead0[((0, j) for j in xrange(W))] = 4 * t
+    for hopping in lat.nearest:
+        lead0[lead0.possible_hoppings(*hopping)] = - t
+
+    # ... then the lead to the right.  We use a method that returns a copy of
+    # `lead0` with its direction reversed.
+    lead1 = lead0.reversed()
+
+    #### Attach the leads and return the finalized system. ####
+    sys.attach_lead(lead0)
+    sys.attach_lead(lead1)
+
+    return sys.finalized()
+
+def plot_conductance(fsys, energy, welldepths):
+    # We specify that we want to not only read, but also write to a
+    # global variable.
+    global pot
+
+    # Compute conductance
+    data = []
+    for welldepth in welldepths:
+        # Set the global variable that defines the potential well depth
+        pot = -welldepth
+
+        smatrix = kwant.solve(fsys, energy)
+        data.append(smatrix.transmission(1, 0))
+
+    pylab.plot(welldepths, data)
+    pylab.xlabel("well depth [in units of t]",
+                 fontsize=latex.mpl_label_size)
+    pylab.ylabel("conductance [in units of e^2/h]",
+                 fontsize=latex.mpl_label_size)
+    fig = pylab.gcf()
+    pylab.setp(fig.get_axes()[0].get_xticklabels(),
+               fontsize=latex.mpl_tick_size)
+    pylab.setp(fig.get_axes()[0].get_yticklabels(),
+               fontsize=latex.mpl_tick_size)
+    fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
+    fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
+    fig.savefig("tutorial2b_result.pdf")
+    fig.savefig("tutorial2b_result.png",
+                dpi=(html.figwidth_px/latex.mpl_width_in))
+
+
+def main():
+    fsys = make_system()
+
+    # We should see conductance steps.
+    plot_conductance(fsys, energy=0.2,
+                     welldepths=[0.01 * i for i in xrange(100)])
+
+
+# Call the main function if the script gets executed (as opposed to imported).
+# See <http://docs.python.org/library/__main__.html>.
+if __name__ == '__main__':
+    main()

doc/source/images/tutorial2c.py

@@ -11,204 +11,177 @@
 
 from cmath import exp
 from math import pi
-import kwant
-
-import latex, html
-
-# First, define the tight-binding system
-
-sys = kwant.Builder()
-
-# Here, we are only working with square lattices
-
-# for simplicity, take lattice constant = 1
-a = 1
-lat = kwant.lattice.Square(a)
-
-t = 1.0
-W = 10
-r1 = 10
-r2 = 20
-
-# Define the scattering region
-# Now, we aim for a more compelx shape, namely a ring (or annulus)
-
-def ring(pos):
-    (x, y) = pos
-    rsq = x**2 + y**2
-    return ( r1**2 < rsq < r2**2)
-
-sys[lat.shape(ring, (0, 11))] = 4 * t
-for hopping in lat.nearest:
-    sys[sys.possible_hoppings(*hopping)] = - t
-
-# 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 Builder repeatedly,
-# we define the modified hoppings as a function that takes the flux
-# as a global variable.
-
-def fluxphase(site1, site2):
-    return exp(1j * phi)
-
-# Now go through all the hoppings and modify those in the lower
-# arm of the ring that go from x=0 to x=1
-
-for (site1, site2) in sys.hoppings():
-    ix1, iy1 = site1.tag
-    ix2, iy2 = site2.tag
-
-    hopx = tuple(sorted((ix1, ix2)))
-
-    if hopx == (0, 1) and iy1 == iy2 and iy1 < 0:
-        sys[lat(hopx[1], iy1), lat(hopx[0], iy1)] = fluxphase
-
-# Then, define the leads:
-
-# First the lead to the left
-
-# (Note: in the current version, TranslationalSymmetry takes a
-# realspace vector)
-sym_lead0 = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
-lead0 = kwant.Builder(sym_lead0)
-lead0.default_site_group = lat
-
-def lead_shape(pos):
-    (x, y) = pos
-    return (-1 < x < 1) and ( -W/2 < y < W/2  )
-
-lead0[lat.shape(lead_shape, (0, 0))] = 4 * t
-for hopping in lat.nearest:
-    lead0[lead0.possible_hoppings(*hopping)] = - t
-
-# Then the lead to the right
-# there we can use a special function that simply reverses the direction
-
-lead1 = lead0.reversed()
-
-# Then attach the leads to the system
-
-sys.attach_lead(lead0)
-sys.attach_lead(lead1)
-
-# finalize the system
-
-fsys = sys.finalized()
-
-# and plot it, to make sure it's proper
-
-kwant.plot(fsys, "tutorial2c_sys.pdf", width=latex.figwidth_pt)
-kwant.plot(fsys, "tutorial2c_sys.png", width=html.figwidth_px)
 
-# Now that we have the system, we can compute conductance
-
-energy = 0.15
-phases = []
-data = []
-for iphi in xrange(100):
-    phi = iphi * 0.01 * 3 * 2 * pi
-
-    # compute the scattering matrix at energy energy
-    smatrix = kwant.solve(fsys, energy)
-
-    # compute the transmission probability from lead 0 to
-    # lead 1
-    phases.append(phi / (2 * pi))
-    data.append(smatrix.transmission(1, 0))
+import kwant
 
-# Use matplotlib to write output
-# We should see conductance steps
+# For plotting
 import pylab
 
-pylab.plot(phases, data)
-pylab.xlabel("flux [in units of the flux quantum]",
-             fontsize=latex.mpl_label_size)
-pylab.ylabel("conductance [in units of e^2/h]",
-             fontsize=latex.mpl_label_size)
-fig = pylab.gcf()
-pylab.setp(fig.get_axes()[0].get_xticklabels(),
-           fontsize=latex.mpl_tick_size)
-pylab.setp(fig.get_axes()[0].get_yticklabels(),
-           fontsize=latex.mpl_tick_size)
-fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
-fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-fig.savefig("tutorial2c_result.pdf")
-fig.savefig("tutorial2c_result.png",
-            dpi=(html.figwidth_px/latex.mpl_width_in))
-
-# Finally, some plots needed for the notes
-
-sys = kwant.Builder()
-
-sys[lat.shape(ring, (0, 11))] = 4 * t
-for hopping in lat.nearest:
-    sys[sys.possible_hoppings(*hopping)] = - t
-
-sym_lead0 = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
-lead0 = kwant.Builder(sym_lead0)
-lead0.default_site_group = lat
-
-def lead_shape(pos):
-    (x, y) = pos
-    return (-1 < x < 1) and ( 0.5 * W < y < 1.5 * W )
-
-lead0[lat.shape(lead_shape, (0, W))] = 4 * t
-for hopping in lat.nearest:
-    lead0[lead0.possible_hoppings(*hopping)] = - t
-
-# Then the lead to the right
-# there we can use a special function that simply reverses the direction
-
-lead1 = lead0.reversed()
-
-# Then attach the leads to the system
-
-sys.attach_lead(lead0)
-sys.attach_lead(lead1)
-
-# finalize the system
-
-fsys = sys.finalized()
-
-# and plot it, to make sure it's proper
-
-kwant.plot(fsys, "tutorial2c_note1.pdf", width=latex.figwidth_small_pt)
-kwant.plot(fsys, "tutorial2c_note1.png", width=html.figwidth_small_px)
-
-sys = kwant.Builder()
-
-sys[lat.shape(ring, (0, 11))] = 4 * t
-for hopping in lat.nearest:
-    sys[sys.possible_hoppings(*hopping)] = - t
-
-sym_lead0 = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
-lead0 = kwant.Builder(sym_lead0)
-lead0.default_site_group = lat
-
-def lead_shape(pos):
-    (x, y) = pos
-    return (-1 < x < 1) and ( -W/2 < y < W/2  )
-
-lead0[lat.shape(lead_shape, (0, 0))] = 4 * t
-for hopping in lat.nearest:
-    lead0[lead0.possible_hoppings(*hopping)] = - t
-
-# Then the lead to the right
-# there we can use a special function that simply reverses the direction
-
-lead1 = lead0.reversed()
-
-# Then attach the leads to the system
-
-sys.attach_lead(lead0)
-sys.attach_lead(lead1, lat(0, 0))
-
-# finalize the system
-
-fsys = sys.finalized()
-
-# and plot it, to make sure it's proper
+import latex, html
 
-kwant.plot(fsys, "tutorial2c_note2.pdf", width=latex.figwidth_small_pt)
-kwant.plot(fsys, "tutorial2c_note2.png", width=html.figwidth_small_px)
+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)
+
+    sys = 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)
+
+    # and add the corresponding lattice points using the `shape`-function
+    sys[lat.shape(ring, (0, 11))] = 4 * t
+    for hopping in lat.nearest:
+        sys[sys.possible_hoppings(*hopping)] = - t
+
+    # 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 Builder repeatedly,
+    # we define the modified hoppings as a function that takes the flux
+    # through the global variable phi.
+    def fluxphase(site1, site2):
+        return exp(1j * phi)
+
+    def crosses_branchcut(hop):
+        ix0, iy0 = hop[0].tag
+
+        # possible_hoppings 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
+    sys[(hop for hop in sys.possible_hoppings((1,0), lat, lat)
+         if crosses_branchcut(hop))] = fluxphase
+
+    #### Define the leads. ####
+    # left lead
+    sym_lead0 = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
+    lead0 = kwant.Builder(sym_lead0)
+
+    def lead_shape(pos):
+        (x, y) = pos
+        return (-1 < x < 1) and ( -W/2 < y < W/2  )
+
+    lead0[lat.shape(lead_shape, (0, 0))] = 4 * t
+    for hopping in lat.nearest:
+        lead0[lead0.possible_hoppings(*hopping)] = - t
+
+    # Then the lead to the right
+    # (again, obtained using reverse()
+    lead1 = lead0.reversed()
+
+    #### Attach the leads and return the finalized system. ####
+    sys.attach_lead(lead0)
+    sys.attach_lead(lead1)
+
+    return sys.finalized()
+
+
+def make_system_note1(a=1, t=1.0, W=10, r1=10, r2=20):
+    lat = kwant.lattice.Square(a)
+    sys = kwant.Builder()
+    def ring(pos):
+        (x, y) = pos
+        rsq = x**2 + y**2
+        return ( r1**2 < rsq < r2**2)
+    sys[lat.shape(ring, (0, 11))] = 4 * t
+    for hopping in lat.nearest:
+        sys[sys.possible_hoppings(*hopping)] = - t
+    sym_lead0 = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
+    lead0 = kwant.Builder(sym_lead0)
+    def lead_shape(pos):
+        (x, y) = pos
+        return (-1 < x < 1) and ( 0.5 * W < y < 1.5 * W )
+    lead0[lat.shape(lead_shape, (0, W))] = 4 * t
+    for hopping in lat.nearest:
+        lead0[lead0.possible_hoppings(*hopping)] = - t
+    lead1 = lead0.reversed()
+    sys.attach_lead(lead0)
+    sys.attach_lead(lead1)
+    return sys.finalized()
+
+
+def make_system_note2(a=1, t=1.0, W=10, r1=10, r2=20):
+    lat = kwant.lattice.Square(a)
+    sys = kwant.Builder()
+    def ring(pos):
+        (x, y) = pos
+        rsq = x**2 + y**2
+        return ( r1**2 < rsq < r2**2)
+    sys[lat.shape(ring, (0, 11))] = 4 * t
+    for hopping in lat.nearest:
+        sys[sys.possible_hoppings(*hopping)] = - t
+    sym_lead0 = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
+    lead0 = kwant.Builder(sym_lead0)
+    def lead_shape(pos):
+        (x, y) = pos
+        return (-1 < x < 1) and ( -W/2 < y < W/2  )
+    lead0[lat.shape(lead_shape, (0, 0))] = 4 * t
+    for hopping in lat.nearest:
+        lead0[lead0.possible_hoppings(*hopping)] = - t
+    lead1 = lead0.reversed()
+    sys.attach_lead(lead0)
+    sys.attach_lead(lead1, lat(0, 0))
+    return sys.finalized()
+
+
+def plot_conductance(fsys, energy, fluxes):
+    # compute conductance
+    # global variable phi controls the flux
+    global phi
+
+    normalized_fluxes = [flux/(2 * pi) for flux in fluxes]
+    data = []
+    for flux in fluxes:
+        phi = flux
+
+        smatrix = kwant.solve(fsys, energy)
+        data.append(smatrix.transmission(1, 0))
+
+    pylab.plot(normalized_fluxes, data)
+    pylab.xlabel("flux [in units of the flux quantum]",
+                 fontsize=latex.mpl_label_size)
+    pylab.ylabel("conductance [in units of e^2/h]",
+                 fontsize=latex.mpl_label_size)
+    fig = pylab.gcf()
+    pylab.setp(fig.get_axes()[0].get_xticklabels(),
+               fontsize=latex.mpl_tick_size)
+    pylab.setp(fig.get_axes()[0].get_yticklabels(),
+               fontsize=latex.mpl_tick_size)
+    fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
+    fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
+    fig.savefig("tutorial2c_result.pdf")
+    fig.savefig("tutorial2c_result.png",
+                dpi=(html.figwidth_px/latex.mpl_width_in))
+
+
+def main():
+    fsys = make_system()
+
+    # Check that the system looks as intended.
+    kwant.plot(fsys, "tutorial2c_sys.pdf", width=latex.figwidth_pt)
+    kwant.plot(fsys, "tutorial2c_sys.png", width=html.figwidth_px)
+
+    # We should see a conductance that is periodic with the flux quantum
+    plot_conductance(fsys, energy=0.15, fluxes=[0.01 * i * 3 * 2 * pi
+                                                for i in xrange(100)])
+
+    # Finally, some plots needed for the notes
+    fsys = make_system_note1()
+    kwant.plot(fsys, "tutorial2c_note1.pdf", width=latex.figwidth_small_pt)
+    kwant.plot(fsys, "tutorial2c_note1.png", width=html.figwidth_small_px)
+    fsys = make_system_note2()
+    kwant.plot(fsys, "tutorial2c_note2.pdf", width=latex.figwidth_small_pt)
+    kwant.plot(fsys, "tutorial2c_note2.png", width=html.figwidth_small_px)
+
+
+# Call the main function if the script gets executed (as opposed to imported).
+# See <http://docs.python.org/library/__main__.html>.
+if __name__ == '__main__':
+    main()

doc/source/images/tutorial3a.py

@@ -7,52 +7,70 @@
 #  - Computing the band structure of a finalized lead.
 
 import kwant
+
 import numpy as np
 from math import pi
 
+# For plotting
+import pylab
+
 import latex, html
 
-a = 1
-lat = kwant.lattice.Square(a)
 
-t = 1.0
-W = 10
+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)
 
-# Define a lead:
+    sym_lead = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
+    lead = kwant.Builder(sym_lead)
+    lead.default_site_group = lat
 
-sym_lead = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
-lead = kwant.Builder(sym_lead)
-lead.default_site_group = lat
+    # build up one unit cell of the lead, and add the hoppings
+    # to the next unit cell
+    for j in xrange(W):
+        lead[(0, j)] = 4 * t
 
-for j in xrange(W):
-    lead[(0, j)] = 4 * t
+        if j > 0:
+            lead[(0, j), (0, j-1)] = - t
 
-    if j > 0:
-        lead[(0, j), (0, j-1)] = - t
+        lead[(1, j), (0, j)] = - t
 
-    lead[(1, j), (0, j)] = - t
+    # return a finalized lead
+    return lead.finalized()
 
-# Now compute the band structure
 
-# Only a finalized lead has the information about bandstructure
-flead = lead.finalized()
+def plot_bandstructure(flead, momenta):
+    # Use the method ``energies`` of the finalized lead to compute
+    # the bandstructure
+    energy_list = [flead.energies(k) for k in momenta]
 
-momenta = np.arange(-pi, pi + .01, 0.02 * pi)
-energy_list = [flead.energies(k) for k in momenta]
+    pylab.plot(momenta, energy_list)
+    pylab.xlabel("momentum [in untis of (lattice constant)^-1]",
+                 fontsize=latex.mpl_label_size)
+    pylab.ylabel("energy [in units of t]",
+                 fontsize=latex.mpl_label_size)
+    fig = pylab.gcf()
+    pylab.setp(fig.get_axes()[0].get_xticklabels(),
+               fontsize=latex.mpl_tick_size)
+    pylab.setp(fig.get_axes()[0].get_yticklabels(),
+               fontsize=latex.mpl_tick_size)
+    fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
+    fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
+    fig.savefig("tutorial3a_result.pdf")
+    fig.savefig("tutorial3a_result.png",
+                dpi=(html.figwidth_px/latex.mpl_width_in))
 
-import pylab
-pylab.plot(momenta, energy_list)
-pylab.xlabel("momentum [in untis of (lattice constant)^-1]",
-             fontsize=latex.mpl_label_size)
-pylab.ylabel("energy [in units of t]",
-             fontsize=latex.mpl_label_size)
-fig = pylab.gcf()
-pylab.setp(fig.get_axes()[0].get_xticklabels(),
-           fontsize=latex.mpl_tick_size)
-pylab.setp(fig.get_axes()[0].get_yticklabels(),
-           fontsize=latex.mpl_tick_size)
-fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
-fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-fig.savefig("tutorial3a_result.pdf")
-fig.savefig("tutorial3a_result.png",
-            dpi=(html.figwidth_px/latex.mpl_width_in))
+
+def main():
+    flead = make_lead()
+
+    # list of momenta at which the bands should be computed
+    momenta = np.arange(-pi, pi + .01, 0.02 * pi)
+
+    plot_bandstructure(flead, momenta)
+
+
+# Call the main function if the script gets executed (as opposed to imported).
+# See <http://docs.python.org/library/__main__.html>.
+if __name__ == '__main__':
+    main()

doc/source/images/tutorial3b.py

@@ -12,83 +12,91 @@
 from cmath import exp
 import kwant
 
-import latex, html
+# For eigenvalue computation
+import scipy.linalg as la
 
-# First, define the tight-binding system
+# For plotting
+import pylab
 
-sys = kwant.Builder()
+import latex, html
 
-# Here, we are only working with square lattices
 
-# for simplicity, take lattice constant = 1
-a = 1
-lat = kwant.lattice.Square(a)
+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).
 
-t = 1.0
-r = 10
+    lat = kwant.lattice.Square(a)
 
-# Define the quantum dot
+    sys = kwant.Builder()
 
-def circle(pos):
-    (x, y) = pos
-    rsq = x**2 + y**2
-    return rsq < r**2
+    # Define the quantum dot
+    def circle(pos):
+        (x, y) = pos
+        rsq = x**2 + y**2
+        return rsq < r**2
 
-def hopx(site1, site2):
-    y = site1.pos[1]
-    return - t * exp(-1j * B * y)
+    def hopx(site1, site2):
+        # The magnetic field is controlled by the global variable B
+        y = site1.pos[1]
+        return - t * exp(-1j * B * y)
 
-sys[lat.shape(circle, (0, 0))] = 4 * t
-# hoppings in x-direction
-sys[sys.possible_hoppings((1, 0), lat, lat)] = hopx
-# hoppings in y-directions
-sys[sys.possible_hoppings((0, 1), lat, lat)] = - t
+    sys[lat.shape(circle, (0, 0))] = 4 * t
+    # hoppings in x-direction
+    sys[sys.possible_hoppings((1, 0), lat, lat)] = hopx
+    # hoppings in y-directions
+    sys[sys.possible_hoppings((0, 1), lat, lat)] = - t
 
-# It's a closed system for a change, so no leads
+    # It's a closed system for a change, so no leads
+    return sys.finalized()
 
-# finalize the system
 
-fsys = sys.finalized()
+def plot_spectrum(fsys, Bfields):
+    # global variable B controls the magnetic field
+    global B
 
-# and plot it, to make sure it's proper
+    # 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
 
-kwant.plot(fsys, "tutorial3b_sys.pdf")
-kwant.plot(fsys, "tutorial3b_sys.png")
+    energies = []
+    for Bfield in Bfields:
+        B = Bfield
 
-# 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
+        # Obtain the Hamiltonian as a dense matrix
+        ham_mat = fsys.hamiltonian_submatrix()[0]
 
-import scipy.linalg as la
+        ev = la.eigh(ham_mat, eigvals_only=True)
 
-Bs = []
-energies = []
-for iB in xrange(100):
-    B = iB * 0.002
+        # we only plot the 15 lowest eigenvalues
+        energies.append(ev[:15])
 
-# Obtain the Hamiltonian as a dense matrix
-    ham_mat = fsys.hamiltonian_submatrix()[0]
+    pylab.plot(Bfields, energies)
+    pylab.xlabel("magnetic field [some arbitrary units]",
+                 fontsize=latex.mpl_label_size)
+    pylab.ylabel("energy [in units of t]",
+                 fontsize=latex.mpl_label_size)
+    fig = pylab.gcf()
+    pylab.setp(fig.get_axes()[0].get_xticklabels(),
+               fontsize=latex.mpl_tick_size)
+    pylab.setp(fig.get_axes()[0].get_yticklabels(),
+               fontsize=latex.mpl_tick_size)
+    fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
+    fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
+    fig.savefig("tutorial3b_result.pdf")
+    fig.savefig("tutorial3b_result.png",
+                dpi=(html.figwidth_px/latex.mpl_width_in))
 
-    ev = la.eigh(ham_mat, eigvals_only=True)
 
-    Bs.append(B)
-    energies.append(ev[:15])
+def main():
+    fsys = make_system()
+
+    # We should observe energy levels that flow towards Landau
+    # level energies with increasing magnetic field
+    plot_spectrum(fsys, [iB * 0.002 for iB in xrange(100)])
 
-import pylab
 
-pylab.plot(Bs, energies)
-pylab.xlabel("magnetic field [some arbitrary units]",
-             fontsize=latex.mpl_label_size)
-pylab.ylabel("energy [in units of t]",
-             fontsize=latex.mpl_label_size)
-fig = pylab.gcf()
-pylab.setp(fig.get_axes()[0].get_xticklabels(),
-           fontsize=latex.mpl_tick_size)
-pylab.setp(fig.get_axes()[0].get_yticklabels(),
-           fontsize=latex.mpl_tick_size)
-fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
-fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-fig.savefig("tutorial3b_result.pdf")
-fig.savefig("tutorial3b_result.png",
-            dpi=(html.figwidth_px/latex.mpl_width_in))
+# 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()