Dear Khani,

I don't understand which comparaison are you talking about, but I can see that you are not defining the angle theta correctly:
Theta is the angle made by the vector k and the x-axis so, it should be 
k=sqrt(kx**2+ky**2)
Theta=180/pi *arcsin(kx/k)

You should notice that E= vf * |k|, so the expression you wrote is not correct.

I hope this helps,
Adel



On Sun, Sep 6, 2020 at 12:23 PM Khani Hosein <hoseinkhaniphy@gmail.com> wrote:
Dear all,
I want to calculate the transmission vs incident angle in graphene with a potential barrier. I use sys = kwant.wraparound.wraparound(sys)  and  lead = kwant.wraparound.wraparound(lead, keep=0) to construct the periodical system. I almost get the correct results, but I find that the definition of angle is different. Could you please help me to find out the reason. The results comparison is also attached.
Regards,
Khani Hosein
my code is:
import numpy as np
import matplotlib.pyplot as plt
from numpy import pi,sqrt
import kwant
import sys
from numpy import pi
from math import asin
from matplotlib import pyplot

a = 1.
V0 = 0.0
W = 2.
L = 448
t = 3

def Rectangle(pos):
    x,y=pos
    return 0<=x<=L and -0.5<=y<=W

sin_30, cos_30 = (1 / 2, sqrt(3) / 2)
lat = kwant.lattice.general([(1, 0), (sin_30, cos_30)],
                                 [(0, 0), (1/2, 1 / (2*sqrt(3)))])
a,b= lat.sublattices


sys = kwant.Builder(kwant.TranslationalSymmetry(lat.vec((-1, 2))))
sys[lat.shape(Rectangle,(0,0.6))]= 0.2  
sys[lat.neighbors(1)] = -t
sys = kwant.wraparound.wraparound(sys)

def lead_shape(pos):
    x,y=pos
    return    -0.5<=y<=W

lead = kwant.Builder(kwant.TranslationalSymmetry((1, 0), lat.vec((-1, 2))))
lead[lat.shape(lead_shape,(0,0.6))] = 0.
lead[lat.neighbors(1)] = -t

lead = kwant.wraparound.wraparound(lead, keep=0)
kwant.plot(sys)

sys.attach_lead(lead)
sys.attach_lead(lead.reversed())

kwant.plot(sys)
sys = sys.finalized()

ky=0.
momenta = [-2*pi + 0.002 * 2*pi * i for i in range(1001)]
bands = kwant.physics.Bands(sys.leads[0],[ky])
energies = [bands(k) for k in momenta]

pyplot.figure()
pyplot.plot(momenta, energies)
pyplot.xlabel("momentum [(lattice constant)^-1]")
pyplot.ylabel("energy [t]")
pyplot.show()


kys = np.arange(-0.08, 0.08, 0.0002)
kysp =np.arange(-90, 90, 0.225)
energy=0.08
transmission = []
num_prop = []
thetap=[]

for ky in kys:
    theta=180.*asin(ky/energy)/pi
    thetap.append(theta)
    smatrix = kwant.smatrix(sys, energy, [ky])
    transmission.append(smatrix.transmission(1, 0))
    num_prop.append(smatrix.num_propagating(0))
    print (ky)

plt.plot(thetap, transmission)
plt.show()





--
Abbout Adel