Commit a67d712f authored by Anton Akhmerov's avatar Anton Akhmerov
Browse files

add function calculating the wave function in the middle of the system

parent ec021b38
......@@ -213,3 +213,30 @@ def fast_combine(s2, s1):
def min_combine(s2, s1):
"""Combine two scattering matrices, calculating only the r block"""
return combine_generic(s2, s1, ('r'))
def get_wave_function(s2, s1):
"""Calculate wave function in the middle of the system
s2 (the first argument) is the scattering matrix to the right of the
slice. where the wave function has to be calculated. Accordingly, s1 is
the scattering matrix of the region to the left. Current is sent from the
left lead to the right lead. This function is necessary to calculate
things like current or particle density distribution
The basis in which the wave function is given is the current basis, with
each mode carrying unit current. First half of the modes are right-moving,
second half is left-moving'
"""
assert s1.rp != None
assert s2.r != None
assert s1.t != None
n = s1.rp.shape[0]
assert s1.rp.shape == s2.r.shape == (n, n)
assert s1.t.shape[1] == n
psi = np.mat(s1.t.shape[0] * [1]).T
temp = la.lu_factor(ml.identity(n) - s1.rp * s2.r)
psi1 = la.lu_solve(temp, s1.t * psi)
psi2 = s2.r * psi1
return (psi1, psi2)
from __future__ import division
import numpy as np, numpy.matlib as ml
from guts.sm import SMatrix, MMatrix, combine, fast_combine, min_combine
from guts.sm import *
from guts.misc import arrays_almost_equal
def make_some_scattering_matrix(alpha):
......@@ -74,3 +74,12 @@ def test_from_tb():
[0, 0, 0, t * b**3]])
m = MMatrix.from_tb(h, hop)
misc.check_m_det_current(m.to_matrix())
def test_get_wave_function():
from numpy.linalg import norm
s1 = make_some_scattering_matrix(.7)
s2 = make_some_scattering_matrix(.3)
psi1, psi2 = get_wave_function(s1, s2)
s3 = min_combine(s1, s2)
assert abs(2 - norm(s3.r * np.mat([1, 1]).T)**2 -
norm(psi1)**2 + norm(psi2)**2) < 1e-10
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment