Commit 3a7470ad authored by Lucas Laplanche's avatar Lucas Laplanche
Browse files

SCATTERING MATRIX METHOD FONCTIONNE !!!!!!!!! BOOOOOUUUUM !!!!!!!!!!

parent 533aa2d5
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ProjectRootManager" version="2" project-jdk-name="Python 3.8 (eam-vcsel)" project-jdk-type="Python SDK" />
<component name="ProjectRootManager" version="2" project-jdk-name="Python 3.9 (eam-vcsel)" project-jdk-type="Python SDK" />
</project>
\ No newline at end of file
......@@ -2,7 +2,7 @@
<module type="PYTHON_MODULE" version="4">
<component name="NewModuleRootManager">
<content url="file://$MODULE_DIR$" />
<orderEntry type="jdk" jdkName="Python 3.8 (eam-vcsel)" jdkType="Python SDK" />
<orderEntry type="inheritedJdk" />
<orderEntry type="sourceFolder" forTests="false" />
</component>
</module>
\ No newline at end of file
......@@ -45,14 +45,14 @@ def electromagnetic_amplitude(bypass_dbr=True, electric_field=0., wavelength=850
sl = op.algaas_super_lattice_refractive_index(sl, electric_field, wavelength)
sl = pt.cut_in_equal_layers_thickness(sl, sl['thickness'].min())
em = op.em_amplitude_smm(sl, wavelength)
sl['electromagnetic_amplitude'] = em.tolist()
sl.insert(sl.shape[1], 'electromagnetic_amplitude', value=em)
plt.plot_refra_em(sl)
def zandberg_electromagnetic_amplitude(electric_field=0., wavelength=1e-6):
sl = st.structure_zandbergen()
sl = op.algaas_super_lattice_refractive_index(sl, electric_field, wavelength)
sl = pt.cut_in_equal_layers_thickness(sl, 1e-9)
sl = pt.cut_in_equal_layers_thickness(sl, 1e-10)
em = op.em_amplitude_smm(sl, wavelength)
sl.insert(sl.shape[1], 'electromagnetic_amplitude', value=em)
plt.plot_refra_em(sl)
......
import numpy as np
from scipy.linalg import expm
from scipy.linalg import expm, fractional_matrix_power
from tqdm import tqdm
......@@ -496,6 +496,9 @@ def em_amplitude_smm(super_lattice, wavelength):
# https://empossible.net/wp-content/uploads/2020/05/Lecture-TMM-Using-Scattering-Matrices.pdf
# number of layers
n_layers = super_lattice.shape[0]
# vacuum wave number
k_0 = 2 * np.pi / wavelength
......@@ -507,65 +510,89 @@ def em_amplitude_smm(super_lattice, wavelength):
v_global = np.matmul(q_global, np.linalg.inv(eigen_g))
# initialize global scattering matrix
s_global = np.array([scattering_matrix_global()])
s_global_ini = scattering_matrix_global()
s_global = np.zeros((n_layers, 2, 2, 2, 2), dtype=np.complex128)
# empty electromagnetic field
e = np.zeros(super_lattice.shape[0])
e = np.zeros(n_layers)
# empty matrices
eigen_i = np.array([])
v_i = np.array([])
eigen_i = np.zeros((n_layers, 2, 2), dtype=np.complex128)
v_i = np.zeros((n_layers, 2, 2), dtype=np.complex128)
# loop forward through layers
for i in tqdm(range(super_lattice.shape[0])):
for i in tqdm(range(n_layers)):
# calculate parameters for layer i
k_zt_i = super_lattice.at[i, 'refractive_index']
l = super_lattice.at[i, 'thickness']
k_z_i = k_0 * k_zt_i
eigen_i = np.append([eigen_i, 1j * k_zt_i * np.identity(2)])
eigen_i[i] = 1j * k_zt_i * np.identity(2)
q_i = np.array([[0., k_zt_i**2],
[-k_zt_i**2, 0.]])
v_i = np.append([v_i, np.matmul(q_i, np.linalg.inv(eigen_i[i]))])
v_i[i] = np.matmul(q_i, np.linalg.inv(eigen_i[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[i], v_global, x_i)
# update global scattering matrix
s_global = np.append(s_global, redheffer_star_product(s_global[i], s_i))
if i==0:
s_global[i] = redheffer_star_product(s_global_ini, s_i)
else:
s_global[i] = redheffer_star_product(s_global[i -1], s_i)
# w for external mode coefficients
wv_g = np.array([[np.identity(2)], [np.identity(2)],
[v_global], [-v_global]])
wv_g = np.array([[np.identity(2), np.identity(2)],
[v_global, -v_global]])
wv_g = concatenate_2x2x2x2_to_4x4(wv_g)
# incident and reflection coefficient
c_inc = np.array([[1], [0]])
c_ref = s_global[n_layers -1][0][0] @ c_inc
# loop backward through layers
for i in tqdm(range(super_lattice.shape[0])):
for i in tqdm(range(n_layers)):
# iterator from n to 0
di = super_lattice.shape[0] -i -1
di = n_layers -i -1
e_kl = eigen_i[di]
e_mkl = fractional_matrix_power(e_kl, -1.)
el_i = np.array([[e_kl, np.zeros((2, 2))],
[np.zeros((2, 2)), e_mkl]])
wv_i = np.array([[np.identity(2)], [np.identity(2)],
[v_i[di]], [-v_i[di]]])
wv_i = np.array([[np.identity(2), np.identity(2)],
[v_i[di], -v_i[di]]])
el_i = np.array([[], [],
[], []])
# external mode coefficients
c_ip = np.invert(s_global[di][0][1]) @ (c_ref -s_global[i][0][0] * c_inc)
c_im = s_global[i][1][0] * c_inc +ee
mc_e = np.array([[c_ip],
[c_im]])
c_im = np.linalg.inv(s_global[di][0][1]) @ (c_ref -s_global[di][0][0] @ c_inc)
c_ip = s_global[di][1][0] @ c_inc +s_global[di][1][1] @ c_im
c_e_i = np.array([c_ip[0],
c_ip[1],
c_im[0],
c_im[1]])
# internal mode coefficients
wv_i = concatenate_2x2x2x2_to_4x4(wv_i)
c_i_i = np.linalg.inv(wv_i) @ wv_g @ c_e_i
mc_i = np.invert(wv_i) @ wv_g @
# internal fields
el_i = concatenate_2x2x2x2_to_4x4(el_i)
psi = wv_i @ el_i @ c_i_i
# electromagnetic field amplitude on x axis
e[di] = np.abs(np.real(psi[0].item()))
# normalize
e = e / np.max(e)
return e
......@@ -597,4 +624,14 @@ def fill_array_row(m, v, i):
m[i, x] = v[x]
return m
def concatenate_2x2x2x2_to_4x4(arr):
m1 = np.concatenate([arr[0][0], arr[0][1]], axis=1)
m2 = np.concatenate([arr[1][0], arr[1][1]], axis=1)
m = np.concatenate([m1, m2], axis=0)
return m
\ No newline at end of file
Markdown is supported
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