From 74bcef6adb5a514f328e010bd5ac5a6933ab8879 Mon Sep 17 00:00:00 2001 From: Phoebe Pearce Date: Fri, 19 Jan 2024 11:48:40 +1100 Subject: [PATCH] S4 ellipsometry parameter calculation --- examples/rcwa_ellipsometry.py | 197 ++++++++++++++++++++++++++++++++++ 1 file changed, 197 insertions(+) create mode 100644 examples/rcwa_ellipsometry.py diff --git a/examples/rcwa_ellipsometry.py b/examples/rcwa_ellipsometry.py new file mode 100644 index 0000000..5c7aacd --- /dev/null +++ b/examples/rcwa_ellipsometry.py @@ -0,0 +1,197 @@ +from rayflare.rigorous_coupled_wave_analysis import rcwa_structure +from rayflare.rigorous_coupled_wave_analysis.rcwa import initialise_S +from rayflare.options import default_options +from solcore import material +from solcore.structure import Layer +import numpy as np +import matplotlib.pyplot as plt + +Si = material('Si')() +Air = material('Air')() +size = ((700, 0), (0, 700)) # 700 nm square unit cell (doesn't matter for planar layers) + +layers = [] + +options = default_options() +options.wavelength = np.linspace(300, 1000, 60) * 1e-9 +orders = 2 + +rcwa_strt = rcwa_structure(layers, size, options, Air, Si) + +fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 8)) + +for angle in [45, 55, 65]: + r_s = np.zeros_like(options.wavelength, dtype=complex) + r_p = np.zeros_like(options.wavelength, dtype=complex) + + for wl_ind, wl in enumerate(options.wavelength): + + S = initialise_S(rcwa_strt.size, orders, rcwa_strt.geom_list, rcwa_strt.layers_oc[wl_ind], + rcwa_strt.shapes_oc[wl_ind], rcwa_strt.shapes_names, rcwa_strt.widths, + options.S4_options) + actual_orders = len(S.GetBasisSet()) # get the actual number of orders S4 chose to use + S.SetExcitationPlanewave((angle, 0), 1, 0, 0) + S.SetFrequency(1 / wl) + + (forw, back) = S.GetAmplitudes('layer_1') + + # this isn't clear in the document, but I think this is indexed as: for N orders, + # the first N elements in "forw" are the forward-travelling E-field amplitudes in each + # order, while the elements N:end are the H-field amplitudes (similarly for "back", but + # for the backwards-travelling E and H field amplitudes). + + r_s[wl_ind] = back[0]/forw[0] + + S = initialise_S(rcwa_strt.size, actual_orders, rcwa_strt.geom_list, rcwa_strt.layers_oc[wl_ind], + rcwa_strt.shapes_oc[wl_ind], rcwa_strt.shapes_names, rcwa_strt.widths, + options.S4_options) + S.SetExcitationPlanewave((angle, 0), 0, 1, 0) + S.SetFrequency(1 / wl) + (forw, back) = S.GetAmplitudes('layer_1') + + # E field amplitudes are zero for p-polarised light, need to use H field amplitudes + H_forw = forw[actual_orders:] + H_back = back[actual_orders:] + + r_p[wl_ind] = H_back[0]/H_forw[0] + + # Fresnel equations: + + n1 = np.sqrt(rcwa_strt.layers_oc[:,0]) + n2 = np.sqrt(rcwa_strt.layers_oc[:,1]) + + theta_t = np.arcsin(n1*np.sin(angle*np.pi/180)/n2) + + r_s_F = (n1*np.cos(angle*np.pi/180) - n2*np.cos(theta_t)) / (n1*np.cos(angle*np.pi/180) + n2*np.cos(theta_t)) + r_p_F = (n2*np.cos(angle*np.pi/180) - n1*np.cos(theta_t)) / (n2*np.cos(angle*np.pi/180) + n1*np.cos(theta_t)) + + rho = r_p/r_s + + rho_F = r_p_F/r_s_F + + rho_mag = np.abs(rho) + delta = np.angle(rho) + psi = np.arctan(rho_mag) + + rho_mag_F = np.abs(rho_F) + delta_F = np.angle(rho_F) + psi_F = np.arctan(rho_mag_F) + + ax1.plot(options.wavelength*1e9, 180*rho_mag/np.pi, label=str(angle) + ' (RCWA)') + ax2.plot(options.wavelength*1e9, 180 - delta*180/np.pi) + # not completely sure why I need to do 180 - angle (or take the negative below for rho) + + ax1.plot(options.wavelength*1e9, 180*rho_mag_F/np.pi, '--', label=str(angle) + ' (Fresnel)') + ax2.plot(options.wavelength*1e9, -delta_F*180/np.pi, '--') + +ax1.legend() +ax1.set_ylabel('Psi (degrees)') +ax2.set_ylabel('Delta (degrees)') +ax2.set_xlabel('Wavelength (nm)') +ax1.set_title('Planar Air/Si') +plt.show() + +r_s_65 = r_s +r_p_65 = r_p + +# try it with some kind of grating + +Si = material('Si')() +Air = material('Air')() + +layers = [Layer(1e-6, Si, + geometry=[{"type": "rectangle", "mat": Air, + "center": (500, 500), "halfwidths": (100, 100), "angle": 45}])] + +orders = 60 + +rcwa_strt = rcwa_structure(layers, size, options, Air, Si) + +fig2, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 8)) + +for angle in [45, 55, 65]: + r_s = np.zeros_like(options.wavelength, dtype=complex) + r_p = np.zeros_like(options.wavelength, dtype=complex) + + for wl_ind, wl in enumerate(options.wavelength): + S = initialise_S(rcwa_strt.size, orders, rcwa_strt.geom_list, rcwa_strt.layers_oc[wl_ind], + rcwa_strt.shapes_oc[wl_ind], rcwa_strt.shapes_names, rcwa_strt.widths, + options.S4_options) + actual_orders = len(S.GetBasisSet()) + S.SetExcitationPlanewave((angle, 0), 1, 0, 0) + S.SetFrequency(1 / wl) + + (forw, back) = S.GetAmplitudes('layer_1') + + # this isn't clear in the document, but I think this is indexed as: for N orders, + # the first N elements in "forw" are the forward-travelling E-field amplitudes in each + # order, while the elements N:end are the H-field amplitudes (similarly for "back", but + # for the backwards-travelling E and H field amplitudes). + + r_s[wl_ind] = back[0] / forw[0] + + S = initialise_S(rcwa_strt.size, actual_orders, rcwa_strt.geom_list, rcwa_strt.layers_oc[wl_ind], + rcwa_strt.shapes_oc[wl_ind], rcwa_strt.shapes_names, rcwa_strt.widths, + options.S4_options) + S.SetExcitationPlanewave((angle, 0), 0, 1, 0) + S.SetFrequency(1 / wl) + (forw, back) = S.GetAmplitudes('layer_1') + + # E field amplitudes are zero for p-polarised light, need to use H field amplitudes + H_forw = forw[actual_orders:] + H_back = back[actual_orders:] + + r_p[wl_ind] = H_back[0] / H_forw[0] + + # Fresnel equations: + + n1 = np.sqrt(rcwa_strt.layers_oc[:, 0]) + n2 = np.sqrt(rcwa_strt.layers_oc[:, 2]) + + theta_t = np.arcsin(n1 * np.sin(angle * np.pi / 180) / n2) + + r_s_F = (n1 * np.cos(angle * np.pi / 180) - n2 * np.cos(theta_t)) / ( + n1 * np.cos(angle * np.pi / 180) + n2 * np.cos(theta_t)) + r_p_F = (n2 * np.cos(angle * np.pi / 180) - n1 * np.cos(theta_t)) / ( + n2 * np.cos(angle * np.pi / 180) + n1 * np.cos(theta_t)) + + rho = r_p / r_s + + rho_F = r_p_F / r_s_F + + rho_mag = np.abs(rho) + delta = np.angle(rho) + psi = np.arctan(rho_mag) + + rho_mag_F = np.abs(rho_F) + delta_F = np.angle(rho_F) + psi_F = np.arctan(rho_mag_F) + + ax1.plot(options.wavelength * 1e9, 180 * rho_mag / np.pi, label=str(angle) + ' (RCWA)') + ax2.plot(options.wavelength * 1e9, 180 - delta * 180 / np.pi) + + ax1.plot(options.wavelength * 1e9, 180 * rho_mag_F / np.pi, '--', + label=str(angle) + ' (Fresnel)') + ax2.plot(options.wavelength * 1e9, -delta_F * 180 / np.pi, '--') + +ax1.legend() +ax1.set_ylabel('Psi (degrees)') +ax2.set_ylabel('Delta (degrees)') +ax2.set_xlabel('Wavelength (nm)') +ax1.set_title('Air/Si with grating') +plt.show() + +r_s_65_grating = r_s +r_p_65_grating = r_p + +R_zeroorder_nograting = 0.5*(np.abs(r_s_65)**2 + np.abs(r_p_65)**2) +R_zeroorder_grating = 0.5*(np.abs(r_s_65_grating)**2 + np.abs(r_p_65_grating)**2) + +plt.figure() +plt.plot(options.wavelength*1e9, R_zeroorder_nograting, label='Planar') +plt.plot(options.wavelength*1e9, R_zeroorder_grating, '--', label='Grating') +plt.legend() +plt.xlabel('Wavelength (nm)') +plt.ylabel('R (zeroth order)') +plt.show() +