Commit 930d30bf authored by Lucas Laplanche's avatar Lucas Laplanche
Browse files

smm verified still wrong em field

parent ad221c94
......@@ -381,20 +381,28 @@ def em_amplitude_tmm(super_lattice, wavelength):
def scattering_matrix_global():
iden = np.identity(2)
ze = np.zeros((2, 2))
return np.array([[ze, iden], [iden, ze]])
def scattering_matrix_layer(vi, vg, xi):
a = np.identity(2) + np.matmul(np.linalg.inv(vi), vg)
b = np.identity(2) - np.matmul(np.linalg.inv(vi), vg)
xa = np.matmul(xi, a)
xb = np.matmul(xi, b)
bab = np.matmul(np.matmul(b, np.linalg.inv(a)), b)
xbax = np.matmul(np.matmul(xb, np.linalg.inv(a)), xi)
d = a - np.matmul(np.matmul(xb, np.linalg.inv(a)), xb)
e = np.matmul(np.matmul(xb, np.linalg.inv(a)), xa) - b
d = a - np.matmul(xbax, b)
e = np.matmul(xbax, a) - b
f = np.matmul(xi, a - bab)
s_11_22 = np.matmul(np.linalg.inv(d), e)[0, 0]
s_12_21 = np.matmul(np.linalg.inv(d), f)[0, 0]
s_11_22 = np.matmul(np.linalg.inv(d), e)
s_12_21 = np.matmul(np.linalg.inv(d), f)
s = np.array([[s_11_22, s_12_21],
[s_12_21, s_11_22]])
......@@ -403,20 +411,26 @@ def scattering_matrix_layer(vi, vg, xi):
return s
def redheffer_star_product(sa, sb):
d = sa[0, 1] / (1. - sb[0, 0] * sa[1, 1])
f = sb[1, 0] / (1. - sa[1, 1] * sb[0, 0])
def redheffer_star_product(sa, sb, e):
ba = np.matmul(sb[0, 0], sa[1, 1])
ab = np.matmul(sa[1, 1], sb[0, 0])
d = np.matmul(sa[0, 1], np.linalg.inv(np.identity(2) -ba))
f = np.matmul(sb[1, 0], np.linalg.inv(np.identity(2) -ab))
sab = np.array([[
sa[0, 0] + d * sb[0, 0] * sa[1, 0],
d * sb[0, 1]
sa[0, 0] +np.matmul(np.matmul(d, sb[0, 0]), sa[1, 0]),
np.matmul(d, sb[0, 1])
], [
f * sa[1, 0],
sb[1, 1] + f * sa[1, 1] * sb[0, 1]
np.matmul(f, sa[1, 0]),
sb[1, 1] +np.matmul(np.matmul(f, sa[1, 1]), sb[0, 1])
]])
amp = np.abs(sab[0, 0][0, 0] +sab[1, 1][0, 0])**2
e = np.append(e, amp)
return sab
return sab, e
def em_amplitude_smm(super_lattice, wavelength):
......@@ -433,31 +447,35 @@ def em_amplitude_smm(super_lattice, wavelength):
# here since theta = 0, kx = ky = 0
q_global = np.array([[0., 1.],
[-1., 0.]])
v_global = -1j * q_global
eigen_g = 1j * np.identity(2)
v_global = np.matmul(q_global, np.linalg.inv(eigen_g))
# initialize global scattering matrix
s_global = np.array([[0., 1.],
[1., 0.]])
s_global = scattering_matrix_global()
# empty electromagnetic field
e = np.array([])
# main loop through layers
for i in tqdm(range(super_lattice.shape[0])):
# calculate parameters for layer i
n = super_lattice.at[i, 'refractive_index']
k_zt_i = super_lattice.at[i, 'refractive_index']
l = super_lattice.at[i, 'thickness']
k_z = k_0 * n
eigen = 1j * n * np.identity(2)
k_z_i = k_0 * k_zt_i
eigen_i = 1j * k_zt_i * np.identity(2)
q_i = np.array([[0., n**2],
[-n**2, 0.]])
v_i = np.matmul(q_i, np.linalg.inv(eigen))
x_i = expm(1j * k_z * l * np.identity(2))
q_i = np.array([[0., k_zt_i**2],
[-k_zt_i**2, 0.]])
v_i = np.matmul(q_i, np.linalg.inv(eigen_i))
x_i = expm(1j * k_z_i * l * np.identity(2))
# calculate scattering matrix for layer i
s_i = scattering_matrix_layer(v_i, v_global, x_i)
# update global scattering matrix
s_global = redheffer_star_product(s_global, s_i)
s_global, e = redheffer_star_product(s_global, s_i, e)
# compute reflection side scattering matrix
......
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