Commit 344281eb authored by Lucas Laplanche's avatar Lucas Laplanche
Browse files

Merge remote-tracking branch 'origin/master'

# Conflicts:
#	optic.py
parents 7b3fdf6a e52aaf78
import numpy as np
import pandas as pd
from scipy.interpolate import splrep, splev
from tqdm import tqdm
......@@ -12,6 +11,7 @@ import plot as plt
import quantum as qt
import scipy.io
import super_lattice_structure as st
import transfer_matrix_method as tmm
import vertex_70 as vtx
......@@ -33,29 +33,6 @@ def refra(bypass_dbr=True, electric_field=0., wavelength=850e-9):
plt.plot_refra(sl)
def dbr_refra_doping(electric_field=0., wavelength=850e-9):
sl = pd.DataFrame(columns = ['name', 'thickness', 'al', 'na', 'nd', 'refractive_index', 'ga6', 'ga11', 'al5', 'al12'])
sl = sl.append(pd.DataFrame([['air', 1000e-9, 0., 0., 0., False, False, False, False]],
columns=['name', 'thickness', 'al', 'na', 'nd', 'ga6', 'ga11', 'al5', 'al12']),
ignore_index=True)
l_dbr_low_al = 850e-9 / (4 * op.algaas_refractive_index(0.15, 0., 850e-9))
l_dbr_high_al = 850e-9 / (4 * op.algaas_refractive_index(0.9, 0., 850e-9))
sl = sl.append(st.structure_dbr(period=20,
grading_type='linear slope',
l_dbr_low_al=l_dbr_low_al,
l_dbr_high_al=l_dbr_high_al))
sl = sl.append(pd.DataFrame([['substrate', 600e-6, 0., 0., 0., False, False, False, False]],
columns=['name', 'thickness', 'al', 'na', 'nd', 'ga6', 'ga11', 'al5', 'al12']),
ignore_index=True)
sl = op.algaas_super_lattice_refractive_index(sl, electric_field, 850e-9, lengyel=True)
plt.plot_refra_doping(sl)
def al_doping(bypass_dbr=True):
......@@ -66,11 +43,11 @@ def al_doping(bypass_dbr=True):
def reflectivity(bypass_dbr=False,
start_wavelength=850e-9,
stop_wavelength=851e-9,
start_wavelength=849e-9,
stop_wavelength=852e-9,
electric_field=0.,
n_points=100,
l_eam_clad=5e-9,
l_eam_clad=15e-9,
l_vcsel_clad=15e-9,
plot=True):
wavelength = np.linspace(start_wavelength, stop_wavelength, num=n_points)
......@@ -79,7 +56,7 @@ def reflectivity(bypass_dbr=False,
# wavelength in [m]
# wavelength must be a numpy array
for i in tqdm(range(len(wavelength))):
sl = st.structure_eam_vcsel(bypass_dbr=bypass_dbr, vcsel_only=False, eam_only=False, grading_type='linear digital', mqw_alloy_type='digital',
sl = st.structure_eam_vcsel(bypass_dbr=bypass_dbr, vcsel_only=True, eam_only=False, grading_type='linear digital', mqw_alloy_type='digital',
l_eam_clad=l_eam_clad, l_vcsel_clad=l_vcsel_clad)
sl = op.algaas_super_lattice_refractive_index(sl, electric_field, wavelength[i], lengyel=True)
......@@ -87,7 +64,7 @@ def reflectivity(bypass_dbr=False,
n = sl['refractive_index'].to_numpy(dtype=np.complex128)
d = sl['thickness'].to_numpy(dtype=float)
r[i] = op.reflection(n, d, wavelength[i])
r[i] = tmm.reflection(n, d, wavelength[i])
if plot:
plt.plot_reflectivity(wavelength, r)
......@@ -112,7 +89,7 @@ def dbr_reflectivity(start_wavelength=700e-9,
n = sl['refractive_index'].to_numpy(dtype=np.complex128)
d = sl['thickness'].to_numpy(dtype=float)
r[i] = op.reflection(n, d, wavelength[i])
r[i] = tmm.reflection(n, d, wavelength[i])
if plot:
plt.plot_reflectivity(wavelength, r)
......@@ -149,22 +126,16 @@ def reflectivity_heatmap(bypass_dbr=True,
def electromagnetic_amplitude(bypass_dbr=False, electric_field=0., wavelength=850e-9, algo='transmission'):
l_eam_clad = 15e-9
def electromagnetic_amplitude(bypass_dbr=False, electric_field=0., wavelength=850e-9):
l_eam_clad = 8e-9
l_vcsel_clad = 15e-9
sl = st.structure_eam_vcsel(bypass_dbr=bypass_dbr, vcsel_only=False, grading_type='none', mqw_alloy_type='none', l_eam_clad = l_eam_clad, l_vcsel_clad=l_vcsel_clad)
sl = st.structure_eam_vcsel(bypass_dbr=bypass_dbr, vcsel_only=False, grading_type='linear digital', mqw_alloy_type='digital', l_eam_clad = l_eam_clad, l_vcsel_clad=l_vcsel_clad)
sl = op.algaas_super_lattice_refractive_index(sl, electric_field, wavelength)
sl = pt.cut_in_equal_layers_thickness(sl, 1e-8)
# scattering matrix or transmission matrix for transfer method
if 'scattering' in algo:
em = op.em_amplitude_smm(sl, wavelength)
elif 'transmission' in algo:
em = op.em_amplitude_tmm(sl, wavelength)
elif 'potato' in algo:
em = op.em_amplitude_potato(sl, wavelength)
else:
em = op.em_amplitude_tmm(sl, wavelength)
em = tmm.em_amplitude_scattering_matrix(sl, wavelength)
# insert the electromagnetic amplitude within the super lattice
sl.insert(sl.shape[1], 'electromagnetic_amplitude', value=em)
......@@ -179,11 +150,11 @@ def vcsel_electromagnetic_resonnance(electric_field=0.,
max_em_amp = np.zeros(n_points)
for i in tqdm(range(len(wavelength))):
sl = st.structure_eam_vcsel(grading_type='linear slope', mqw_alloy_type='mean', l_eam_clad = 15.5e-9, l_vcsel_clad=8.5e-9)
sl = st.structure_eam_vcsel(grading_type='linear digital', mqw_alloy_type='digital', l_eam_clad = 8e-9, l_vcsel_clad=15e-9)
sl = op.algaas_super_lattice_refractive_index(sl, electric_field, wavelength[i], lengyel=False)
sl['refractive_index'] = sl['refractive_index'].apply(np.real).astype(float)
sl = pt.cut_in_equal_layers_thickness(sl, 1e-8)
em = op.em_amplitude_smm(sl, wavelength[i], normalize=False)
em = tmm.em_amplitude_scattering_matrix(sl, wavelength[i], normalize=False)
max_em_amp[i] = np.max(em)
# normalize
......@@ -210,9 +181,9 @@ def clad_coupling_electromagnetic_amplitude(bypass_dbr=False, electric_field=0.,
#sl2 = sl2[::-1].reset_index(drop=True)
#sl3 = sl3[::-1].reset_index(drop=True)
em1 = op.em_amplitude_smm(sl1, wavelength)
em2 = op.em_amplitude_smm(sl2, wavelength)
em3 = op.em_amplitude_smm(sl3, wavelength)
em1 = tmm.em_amplitude_scattering_matrix(sl1, wavelength)
em2 = tmm.em_amplitude_scattering_matrix(sl2, wavelength)
em3 = tmm.em_amplitude_scattering_matrix(sl3, wavelength)
sl1.insert(sl1.shape[1], 'electromagnetic_amplitude', value=em1)
sl2.insert(sl2.shape[1], 'electromagnetic_amplitude', value=em2)
......@@ -232,7 +203,7 @@ def animation_clad_coupling(bypass_dbr=False, electric_field=0., wavelength=850e
sl = st.structure_eam_vcsel(bypass_dbr=bypass_dbr, vcsel_only=True, grading_type='none', mqw_alloy_type='none', l_vcsel_clad=cl)
sl = op.algaas_super_lattice_refractive_index(sl, electric_field, wavelength)
sl = pt.cut_in_equal_layers_thickness(sl, 5e-9)
em = op.em_amplitude_smm(sl, wavelength)
em = tmm.em_amplitude_scattering_matrix(sl, wavelength)
sl.insert(sl.shape[1], 'electromagnetic_amplitude', value=em)
plt.plot_refra_em(sl, str(cl))
......@@ -241,7 +212,7 @@ 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)
em = op.em_amplitude_smm(sl, wavelength)
em = tmm.em_amplitude_scattering_matrix(sl, wavelength)
sl.insert(sl.shape[1], 'electromagnetic_amplitude', value=em)
plt.plot_refra_em(sl)
......@@ -288,9 +259,9 @@ def ftir_theory_comparison_reflectivity(filename='eam_2', electric_field=0., doc
wavelength_theory = np.linspace(750e-9, 950e-9, num=500)
if 'eam_2' in filename:
r_theory = op.reflectivity(True, electric_field, wavelength_theory)
r_theory = tmm.reflectivity(True, electric_field, wavelength_theory)
else:
r_theory = op.reflectivity(False, electric_field, wavelength_theory)
r_theory = tmm.reflectivity(False, electric_field, wavelength_theory)
# plot
plt.plot_mult_reflectivity(wavelength, r, wavelength_theory, r_theory)
......@@ -298,7 +269,7 @@ def ftir_theory_comparison_reflectivity(filename='eam_2', electric_field=0., doc
def mqw_psi(lz=0.01e-9, electric_field=0., wavelength=850e-9):
def mqw_psi(lz=0.01e-9):
sl = st.structure_eam_vcsel(eam_mqw_only=True, amount_eam_qw=5)
sl = pt.cut_in_equal_layers_thickness(sl, lz)
v_e, v_hh, eig_e, psi_e, eig_hh, psi_hh = qt.solve_schrodinger(sl)
......
import numpy as np
from scipy.linalg import expm, fractional_matrix_power
from tqdm import tqdm
import pandas_tools as pt
import super_lattice_structure as st
from globals import C, EPS0, HEV, N0, NGAAS, T
from globals import C, HEV, N0, T
from sqw_lengyel_absorption import gaas_sqw_absorption_at_wavelength
......@@ -90,8 +86,6 @@ def algaas_super_lattice_refractive_index(super_lattice, electric_field, wavelen
return super_lattice
def algaas_refractive_index(al, electric_field, wavelength, quantum_well=False, quantum_well_thickness=0., confinement_barrier_al=0., na=0., nd=0.):
# layer is a quantum well -> Lengyel model
if quantum_well:
......@@ -124,129 +118,6 @@ def algaas_refractive_index(al, electric_field, wavelength, quantum_well=False,
def element_matrix(n, d, wavelength):
# equation (1.11)
k0 = 2*np.pi/wavelength
beta = k0*n*d
co = np.cos(beta)
si = np.sin(beta)
return np.array([[co, 1j*si/n], [1j*si*n, co]])
def reflection(n, d, wavelength):
# element matrix
m = np.identity(2)
for i in range(np.size(n)):
m = np.matmul(m, element_matrix(n[i], d[i], wavelength))
# electromagnetic field element amplitude
[e_m, h_m] = np.matmul(m, np.array([1, NGAAS])) # eta_s = NGAAS since theta = 0
# reflection
return np.abs( ((N0*e_m -h_m) / (N0*e_m +h_m))**2 )
def reflection_with_intermediary_values(n, d, wavelength):
# reflection array
r = np.zeros(np.size(n))
# element matrix
m = np.identity(2)
for i in range(np.size(n)):
m = np.matmul(m, element_matrix(n[i], d[i], wavelength))
# electromagnetic field element amplitude
[e_m, h_m] = np.matmul(m, np.array([1, NGAAS])) # eta_s = NGAAS since theta = 0
# reflection
r[i] = np.abs(((N0 * e_m - h_m) / (N0 * e_m + h_m)) ** 2)
return r
def absorption(n, d, wavelength):
# element matrix
m = np.identity(2)
for i in range(np.size(n)):
m = np.matmul(m, element_matrix(n[i], d[i], wavelength))
# electromagnetic field element amplitude
[e_m, h_m] = np.matmul(m, np.array([1, NGAAS])) # eta_s = NGAAS since theta = 0
# absorption
return np.imag( ((N0*e_m -h_m) / (N0*e_m +h_m))**2 )
def reflectivity(bypass_dbr,
electric_field,
wavelength):
r = np.zeros(len(wavelength))
# wavelength in [m]
# wavelength must be a numpy array
for i in tqdm(range(len(wavelength))):
sl = st.structure_eam_vcsel(bypass_dbr=bypass_dbr, eam_only=True)
sl = algaas_super_lattice_refractive_index(sl, electric_field, wavelength[i])
n = sl['refractive_index'].to_numpy(dtype=np.complex128)
d = sl['thickness'].to_numpy(dtype=np.complex128)
r[i] = reflection(n, d, wavelength[i])
return r
def reflectivity_meshgrid(bypass_dbr, electric_field, wavelength, temperature):
# get the dimension of the dataframe
df = st.structure_eam_vcsel(eam_only=True, bypass_dbr=bypass_dbr)
dim = df.shape[0]
# create empty reflection map
r_meshgrid = np.zeros((len(wavelength), dim))
# add the depth column
pt.add_depth_column(df)
# extract the depth column
depth = df['depth'].to_numpy()
growth_speed = 1. # [um/h]
growth_speed *= 1e-6 # [m/h]
growth_speed /= 3600. # [m/s]
# calculate the layer epitaxial time
time = growth_speed * depth
# wavelength in [m]
# wavelength must be a numpy array
for i in tqdm(range(len(wavelength))):
sl = st.structure_eam_vcsel(eam_only=True, bypass_dbr=bypass_dbr)
sl = algaas_super_lattice_refractive_index(sl, electric_field, wavelength[i], temperature)
n = sl['refractive_index'].to_numpy(dtype=np.complex128)
d = sl['thickness'].to_numpy(dtype=np.complex128)
# reflection array
r = reflection_with_intermediary_values(n, d, wavelength[i])
# insert the reflection array
# into the corresponding wavelength column of the reflection map
fill_array_row(r_meshgrid, r, i)
return time, r_meshgrid
def afromovitz_real_algaas_refractive_index(al, wavelength):
# Refractive index of AlGaAs for all compositions
# as a function of temperature and for wavelengths
......@@ -347,267 +218,6 @@ def ioffe_algaas_permittivity(al, high_frequency=True):
def transmission_matrix(super_lattice, m, wavelength):
e = np.zeros(super_lattice.shape[0])
for k in range(super_lattice.shape[0] -1):
i, m = interface_intensity(super_lattice, m, k)
i, m = core_intensity(super_lattice, wavelength, m, k)
e[k] = i
return e, m
def interface_intensity(super_lattice, m, idx):
n = np.real(super_lattice['refractive_index'].to_numpy(dtype=np.complex128))
ra = n[idx] / n[idx +1]
di = np.array([
[0.5 * (1. + ra), 0.5 * (1. - ra)],
[0.5 * (1. - ra), 0.5 * (1. + ra)]
])
m = di @ m
# electric field of the forward plane wave
ep = m[0, 0]
# electric field of the backward plane wave
em = m[1, 0]
i = np.abs(ep +em) ** 2
return i, m
def core_intensity(super_lattice, wavelength, m, idx):
n = np.real(super_lattice['refractive_index'].to_numpy(dtype=np.complex128))[idx + 1]
l = super_lattice['thickness'].to_numpy(dtype=float)[idx + 1]
dp = np.array([
[np.exp(2j * np.pi * n * l / wavelength), 0.],
[0., np.exp(-2j * np.pi * n * l / wavelength)]
])
m = dp @ m
# electric field of the forward plane wave
ep = m[0, 0]
# electric field of the backward plane wave
em = m[1, 0]
i = np.abs(ep + em) ** 2
return i, m
def em_amplitude_tmm(super_lattice, wavelength, normalize=True):
# compute the amplitude of the electromagnetic field through the structure
# using transfer matrix method
# very likely to diverge
# rather use scattering matrix method
# initial values
m = np.identity(2)
# refractive index of the input medium
n_ini = np.real(super_lattice.at[0, 'refractive_index'])
# injected power [W]
p_inj = 1.
# pump area [cm2]
p_area = 1.
# input field in [V/cm]
e_ini = np.sqrt(2. * p_inj / (EPS0 * C * n_ini * p_area))
# transmission matrix
i, tm = transmission_matrix(super_lattice, m, wavelength=wavelength)
tm = np.array([1., -tm[1, 0] / tm[1, 1]])
m_ini = e_ini * tm
m = np.identity(2) * m_ini
# electric field intensity
e, m = transmission_matrix(super_lattice, m, wavelength=wavelength)
# normalize
if normalize:
e = e / np.max(e)
return e
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)
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(xbax, b)
e = np.matmul(xbax, a) - b
f = np.matmul(xi, a - bab)
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]])
return s
def redheffer_star_product(sg, si):
ig = np.matmul(si[0, 0], sg[1, 1])
gi = np.matmul(sg[1, 1], si[0, 0])
d = np.matmul(sg[0, 1], np.linalg.inv(np.identity(2) -ig))
f = np.matmul(si[1, 0], np.linalg.inv(np.identity(2) -gi))
s = np.array([[
sg[0, 0] +np.matmul(np.matmul(d, si[0, 0]), sg[1, 0]),
np.matmul(d, si[0, 1])
], [
np.matmul(f, sg[1, 0]),
si[1, 1] +np.matmul(np.matmul(f, sg[1, 1]), si[0, 1])
]])
return s
def em_amplitude_smm(super_lattice, wavelength, normalize=True):
# compute the amplitude of the electromagnetic field through the structure
# using scattering matrix method
# based on 'TMM using scattering matrices' by EMPossible
# 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
# compute gap medium parameters
# here since theta = 0, kx = ky = 0
q_global = np.array([[0., 1.],
[-1., 0.]])
eigen_g = 1j * np.identity(2)
v_global = np.matmul(q_global, np.linalg.inv(eigen_g))
# initialize global scattering matrix
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(n_layers)
# empty matrices
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 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[i] = 1j * k_zt_i * np.identity(2)
q_i = np.array([[0., k_zt_i**2],
[-k_zt_i**2, 0.]])
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
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 = 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 range(n_layers):
# iterator from n to 0
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]]])
# external mode coefficients
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
# 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()))**2
# normalize
if normalize:
e = e / np.max(e)
return e
def confinement_barrier_mean_al_content(super_lattice, i):
name = super_lattice.at[i +1, 'name']
al = 0.
......@@ -624,104 +234,3 @@ def confinement_barrier_mean_al_content(super_lattice, i):
return al
def fill_array_row(m, v, i):
# m is a matrix of shape [len(v), ?]
for x in range(len(v)):
m[i, x] = v[x]