Thanks a lot for you package which is rather intuitive in using!

I could not understand what kwant.wave_function gives as an output:
-- how to distinguish between the modes which are propagating in the leads and which are decaying in the leads?
kwant.wave_function gives an answer at any energy.

(1)For example, how do get the density in some transport setup? Is this example correct?

def Fermi(energy, T, mu):
f = 1./(1. + math.exp ( (energy - mu) / T) )
return f

wf = kwant.wave_function(sys,energy)

Tmu = numpy.array([[Tleft,muleft],[Tright,muright]])

#assuming that this range of energies encompasses well the chemical potential range
de = 0.05
nenergies = 100
energies=[de * (i - nenergies/2) for i in xrange(nenergies)]

for energy in energies:

d = d * de

(2) Another example which I do not understand properly is  about  transport through a translationally invariant 2d bar.
From physical reasoning, I would expect that the current through each link in the direction of the current is the same.
So I tried to apply the formula (11) from the NJP-introduction to KWANT.
But I get a very different looking answers, comparing to scattering matrix approach, as well as current throught different links in the direction of the current is different. Why are these two answers so different and one is not oscillating, and another is so much depends on energy?
(In 1d the comparission works well)

def Fermi(energy, T, mu):
f = 1./(1. + math.exp ( (energy - mu) / T) )
return f

muleft = -2.
muright = 2.
mu = 0

Tleft = 0.1
Tright = 0.1

#potentials and temepratures of all leads
muTAll = []
muTAll.append([muleft,Tleft])
muTAll.append([muright,Tright])

energies = []
data = []
nenergies = 300
de = 0.032424

CurrentWave = []
CurrentLandauer = []

t = 1.
W = 6
L = 10

numsite = (L-1)* W

sys = BuildSystem(t, W, L,mu,mu,mu)
#beginning and end of stm link in terms of program;
x=[]
for i1 in xrange(numsite):
x.append(sys.sites[i1].pos)
beg=[]
for i1 in xrange(numsite):
if x[i1] == [1,1]:
beg.append(i1)
end=[]
for i1 in xrange(numsite):
if x[i1] == [2,1]:
end.append(i1)
print beg, end

for imu in xrange(1):
sysFinal = BuildSystem(t, W, L, mu, mu, mu)
for ienergy in xrange(nenergies):
energy = de * ienergy - nenergies * de / 2.
energies.append(energy)
wf = kwant.wave_function(sysFinal,energy)
smatrix = kwant.smatrix(sysFinal, energy)

#wave function based calculation which can be used for any correlation function
wfInside = numpy.zeros(numsite)
#wave function inside of the system
wave = wf(p)
lenwave = wave.shape [0]
for ipsi in xrange(lenwave):
wfInside = wfInside +  Fermi(energy, muTAll[p][1],muTAll[p][0] ) * wave[ipsi]

ProbeCurrent =  wfInside[beg] * wfInside[end].conjugate() - wfInside[end] * wfInside[beg].conjugate()
CurrentWave.append( ProbeCurrent[0].imag)
#Getting wave function is done

#transmission based calculation
t10 = smatrix.transmission(1, 0)
data.append(t10)
t2 = t10 * (Fermi(energy, muTAll[1][1],muTAll[1][0] ) - Fermi(energy, muTAll[0][1],muTAll[0][0] ))
CurrentLandauer.append(t2)

plt.figure()
plt.plot(energies, CurrentLandauer)
plt.plot(energies, CurrentWave)
plt.xlabel("energy [t]")
plt.ylabel("conductance [e^2/h]")
plt.show()

def BuildSystem(t, W, L, mu, muleft, muright):

sys = kwant.Builder()

#DEFINING THE SYSTEM
a=1.
lat = kwant.lattice.square(a)

for i in xrange(L-1):
for j in xrange(W):
sys[lat(i,j)] = mu

if j > 0:
sys[lat(i,j), lat(i,j-1)] = -t

if i > 0:
sys[lat(i,j), lat(i-1,j)] = -t

#lead is pointing away from scattering region

for j in xrange(W):
if j > 0:
#modeling large band-width in the leads, that is why 10. t

for j in xrange(W):
numx = L-2