Commit b4ab005a authored by Anton Akhmerov's avatar Anton Akhmerov Committed by Christoph Groth
Browse files

add polar decomposition to sm.py

parent 49add269
......@@ -2,7 +2,7 @@
from __future__ import division
import numpy as np, numpy.matlib as ml
import scipy.linalg as la
import scipy.linalg as la, scipy.sparse as sp
class SMatrix(object):
"""A scattering matrix.
......@@ -45,9 +45,92 @@ class SMatrix(object):
s.rp = m.m12 * s.tp
return s
@classmethod
def from_polar(cls, tran, u, v, up, vp):
"""Make a SMatrix from its polar decomposition."""
s = cls()
tt = np.mat(np.diag(np.sqrt(tran)))
rr = np.mat(np.diag(np.sqrt(1 - tran)))
s.tp = u * tt * vp
s.t = v * tt * up
s.r = - u * rr * up
s.rp = v * rr * vp
return s
def to_matrix(self):
return np.bmat([[self.r, self.tp], [self.t, self.rp]])
def to_polar(self):
"""Calculate the polar decomposition of a SMatrix.
Return a 5-tuple consisting of the list of transmission eigenvalues
tran, and 4 unitary matrices U, V, U', V', as defined in
arXiv:cond-mat/9612179 Eq. 1.24. The function tries to use all the
available blocks of the scattering matrix and returns unit matrices
instead of the elements of decomposition which are unrecoverable.
This function does not use LU deconmposition for efficiency yet, so it
may be optimized.
"""
assert (self.t != None and self.tp != None and self.r != None and
self.rp != None)
(u, tran, vp) = la.svd(self.tp)
tran = tran**2
(u, vp) = (np.mat(u), np.mat(vp))
tt = np.mat(np.diag(np.sqrt(tran)))
rr = np.mat(np.diag(np.sqrt(1 - tran)))
v = self.rp * vp.H * rr.I
up = - rr.I * u.H * self.r
(tran, u, up, v, vp) = 5 * (None,)
if self.r != None:
(u, r, up) = la.svd(self.r)
ri = 1/r
tran = 1 - r**2
n = tran.shape[0]
ti = 1/np.sqrt(tran)
(u, up) = (np.mat(-u), np.mat(up))
tti = sp.csc_matrix(sp.lil_diags([ti], [0], (n, n)))
rri = sp.csc_matrix(sp.lil_diags([ri], [0], (n, n)))
if self.tp != None:
if self.r != None:
vp = tti * u.H * self.tp
else:
(u, t, vp) = la.svd(self.tp)
tran = t**2
n = tran.shape[0]
ti = 1/t
ri = 1/np.sqrt(1 - tran)
(u, vp) = (np.mat(u), np.mat(vp))
tti = sp.csc_matrix(sp.lil_diags([ti], [0], (n, n)))
rri = sp.csc_matrix(sp.lil_diags([ri], [0], (n, n)))
if self.rp != None:
if self.tp != None:
v = self.rp * vp.H * rri
else:
(v, r, vp) = la.svd(self.rp)
ri = 1/r
tran = 1 - r**2
n = tran.shape[0]
ti = 1/np.sqrt(tran)
(v, vp) = (np.mat(v), np.mat(vp))
tti = sp.csc_matrix(sp.lil_diags([ti], [0], (n, n)))
rri = sp.csc_matrix(sp.lil_diags([ri], [0], (n, n)))
if self.t != None:
if self.r == None and self.rp == None:
(v, t, up) = la.svd(self.t)
tran = t**2
(u, vp) = (np.mat(u), np.mat(vp))
elif self.r == None and self.rp != None:
up = tti * v.H * self.t
elif self.r != None and self.rp == None:
v = self.t * up.H * tti
elif self.tp == None:
up = tti * v.H * self.t
vp = rri * v.H * self.r
return tran, u, v, up, vp
def cond_noise(self):
"""Calculate conductance and noise.
......@@ -88,6 +171,7 @@ class SMatrix(object):
return noise.real
class MMatrix(object):
"""A transfer matrix.
......
......@@ -74,3 +74,14 @@ def test_get_wave_function():
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
def test_polar():
"""Test polar_decomposition and from_polar methods of SMatrix."""
s = make_some_circular_s(4, 2)
dec = s.to_polar()
s2 = SMatrix.from_polar(*dec)
for i in dec[0]:
assert 0 <= i <= 1
for i in dec[1:]:
assert arrays_almost_equal(i.H * i, ml.identity(4), 1e-13)
assert arrays_almost_equal(s.to_matrix(), s2.to_matrix(), 1e-13)
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