From 625cd64f5320cdd08e974079b3bbc5fe7b1db13c Mon Sep 17 00:00:00 2001 From: Edoardo Vergerio Date: Wed, 14 Oct 2026 16:40:05 +0200 Subject: [PATCH] implemented new filtering system based on a Filter Class. Previously, filtering was made expensive by extreme padding, now this has been fixed by choosing a more appropriate padding size. The gradient chain based on unitcells has been implemented but the results have yet to be tested using FDTD, which will be the next step. Added some utilities to visualize the effect of filter radius on given input masks and added all data concerning unit cell simulations. Still to be done is to fix materials in the MEEP simulations, so as to account for true material behavior instead of relying on constant refractive indices. --- filters.py | 256 ++--- mpid.py | 177 +++- single_converter_symm.py | 2 +- single_converter_symm_ego.py | 20 +- single_converter_symm_source_sweep.py | 156 ++- testing_filter_radius.py | 55 ++ unit_cell/param_sweep_center_1570.log | 29 + unit_cell/results_S21_p700nm_h900nm.txt | 25 + unit_cell/unitCell_grads.py | 51 + ...center_1530_h900nm_res50_wvl_diameters.npy | Bin 0 -> 320 bytes ...enter_1530_h900nm_res50_wvl_mode_phase.npy | Bin 0 -> 320 bytes ...center_1530_h900nm_res50_wvl_mode_tran.npy | Bin 0 -> 320 bytes ..._center_1530_h900nm_res50_wvl_mode_wvl.npy | Bin 0 -> 136 bytes ...center_1550_h900nm_res50_wvl_diameters.npy | Bin 0 -> 320 bytes ...enter_1550_h900nm_res50_wvl_mode_phase.npy | Bin 0 -> 320 bytes ...center_1550_h900nm_res50_wvl_mode_tran.npy | Bin 0 -> 320 bytes ..._center_1550_h900nm_res50_wvl_mode_wvl.npy | Bin 0 -> 136 bytes ...center_1570_h900nm_res50_wvl_diameters.npy | Bin 0 -> 320 bytes ...enter_1570_h900nm_res50_wvl_mode_phase.npy | Bin 0 -> 320 bytes ...center_1570_h900nm_res50_wvl_mode_tran.npy | Bin 0 -> 320 bytes ..._center_1570_h900nm_res50_wvl_mode_wvl.npy | Bin 0 -> 136 bytes unit_cell/unit_cell_center_wvl.py | 169 ++++ unit_cell/unit_cell_fcen.pdf | Bin 0 -> 14577 bytes .../unit_cell_fcen_h900nm_res50_wvl_1530.pdf | Bin 0 -> 14499 bytes .../unit_cell_fcen_h900nm_res50_wvl_1550.pdf | Bin 0 -> 14537 bytes unit_cell/unit_cell_meep.ipynb | 902 ++++++++++++++++++ unit_cell/unit_cell_param_sweep.pdf | Bin 0 -> 38133 bytes unit_cell/unit_cell_param_sweep.py | 203 ++++ unit_cell/unit_cell_param_sweep_diameters.npy | Bin 0 -> 320 bytes unit_cell/unit_cell_param_sweep_lineplot.pdf | Bin 0 -> 17123 bytes .../unit_cell_param_sweep_mode_phase.npy | Bin 0 -> 3968 bytes unit_cell/unit_cell_param_sweep_mode_tran.npy | Bin 0 -> 3968 bytes unit_cell/unit_cell_param_sweep_mode_wvl.npy | Bin 0 -> 288 bytes unit_cell/unit_cell_plots.py | 42 + 34 files changed, 1891 insertions(+), 196 deletions(-) create mode 100644 testing_filter_radius.py create mode 100644 unit_cell/param_sweep_center_1570.log create mode 100644 unit_cell/results_S21_p700nm_h900nm.txt create mode 100644 unit_cell/unitCell_grads.py create mode 100644 unit_cell/unit_cell_center_1530_h900nm_res50_wvl_diameters.npy create mode 100644 unit_cell/unit_cell_center_1530_h900nm_res50_wvl_mode_phase.npy create mode 100644 unit_cell/unit_cell_center_1530_h900nm_res50_wvl_mode_tran.npy create mode 100644 unit_cell/unit_cell_center_1530_h900nm_res50_wvl_mode_wvl.npy create mode 100644 unit_cell/unit_cell_center_1550_h900nm_res50_wvl_diameters.npy create mode 100644 unit_cell/unit_cell_center_1550_h900nm_res50_wvl_mode_phase.npy create mode 100644 unit_cell/unit_cell_center_1550_h900nm_res50_wvl_mode_tran.npy create mode 100644 unit_cell/unit_cell_center_1550_h900nm_res50_wvl_mode_wvl.npy create mode 100644 unit_cell/unit_cell_center_1570_h900nm_res50_wvl_diameters.npy create mode 100644 unit_cell/unit_cell_center_1570_h900nm_res50_wvl_mode_phase.npy create mode 100644 unit_cell/unit_cell_center_1570_h900nm_res50_wvl_mode_tran.npy create mode 100644 unit_cell/unit_cell_center_1570_h900nm_res50_wvl_mode_wvl.npy create mode 100644 unit_cell/unit_cell_center_wvl.py create mode 100644 unit_cell/unit_cell_fcen.pdf create mode 100644 unit_cell/unit_cell_fcen_h900nm_res50_wvl_1530.pdf create mode 100644 unit_cell/unit_cell_fcen_h900nm_res50_wvl_1550.pdf create mode 100644 unit_cell/unit_cell_meep.ipynb create mode 100644 unit_cell/unit_cell_param_sweep.pdf create mode 100644 unit_cell/unit_cell_param_sweep.py create mode 100644 unit_cell/unit_cell_param_sweep_diameters.npy create mode 100644 unit_cell/unit_cell_param_sweep_lineplot.pdf create mode 100644 unit_cell/unit_cell_param_sweep_mode_phase.npy create mode 100644 unit_cell/unit_cell_param_sweep_mode_tran.npy create mode 100644 unit_cell/unit_cell_param_sweep_mode_wvl.npy create mode 100644 unit_cell/unit_cell_plots.py diff --git a/filters.py b/filters.py index 5e13b84..f88a1ec 100644 --- a/filters.py +++ b/filters.py @@ -1,148 +1,162 @@ import autograd.numpy as npa +import numpy as np +from scipy.fft import next_fast_len from numpy.typing import ArrayLike as ArrayLikeType import matplotlib.pyplot as plt -def _centered(arr: npa.ndarray, newshape: ArrayLikeType) -> npa.ndarray: - """Formats the output of an FFT to center the zero-frequency component. +class Filter: + def __init__(self, radius, space_grid): + self.radius = radius + self.space_grid = space_grid + self.p = int(np.ceil(radius / self.space_grid.get_resolution())) + self.Fx = next_fast_len(self.space_grid.get_Nx() + self.p) + self.Fy = next_fast_len(self.space_grid.get_Ny() + self.p) + self.H = self._construct_kernel(radius, space_grid) + + def apply_filter(self, input_parameters): + input_parameters_2d = input_parameters.reshape(self.space_grid.get_Nx(), self.space_grid.get_Ny()) + input_parameters_padded = self._edge_pad(input_parameters_2d, (((self.Fx-1)//2, (self.Fx)//2), ((self.Fy-1)//2, (self.Fy)//2))) + #print("input_parameters_padded shape:",input_parameters_padded.shape) + input_parameters_fft = npa.fft.fft2(input_parameters_padded) + + Y = input_parameters_fft * self.H + y = self._centered( + npa.real(npa.fft.ifft2(Y)), + (self.space_grid.get_Nx(), self.space_grid.get_Ny()) + ) + return y.flatten() - A helper function borrowed from SciPy: - https://github.com/scipy/scipy/blob/v1.4.1/scipy/signal/signaltools.py#L263-L270 + def _centered(self, arr: npa.ndarray, newshape: ArrayLikeType) -> npa.ndarray: + """Formats the output of an FFT to center the zero-frequency component. - Args: - arr: output array from an FFT operation. - newshape: 1d array with two elements (integers) specifying the dimensions - of the array to be returned. + A helper function borrowed from SciPy: + https://github.com/scipy/scipy/blob/v1.4.1/scipy/signal/signaltools.py#L263-L270 - Returns: - The input array with the zero-frequency component as the central element. - """ - newshape = npa.asarray(newshape) - currshape = npa.array(arr.shape) - startind = (currshape - newshape) // 2 - endind = startind + newshape - myslice = [slice(startind[k], endind[k]) for k in range(len(endind))] + Args: + arr: output array from an FFT operation. + newshape: 1d array with two elements (integers) specifying the dimensions + of the array to be returned. - return arr[tuple(myslice)] + Returns: + The input array with the zero-frequency component as the central element. + """ + newshape = npa.asarray(newshape) + currshape = npa.array(arr.shape) + startind = (currshape - newshape) // 2 + endind = startind + newshape + myslice = [slice(startind[k], endind[k]) for k in range(len(endind))] -def _edge_pad(arr: npa.ndarray, pad: npa.ndarray) -> npa.ndarray: - """Border-pads the edges of an array. + return arr[tuple(myslice)] - Used to preprocess the design weights prior to convolution with the filter. - Border padding an image will set the value of each padded pixel equal to - the value of the nearest pixel in the image. Used to implement feature- - preserving convolution filters that prevent unwanted edge effects. + def _edge_pad(self, arr: npa.ndarray, pad: npa.ndarray) -> npa.ndarray: + """Border-pads the edges of an array. - Args: - arr: 2d array whose borders contain the values to use for padding - pad: 2x2 array of integers indicating the size - of the borders to pad the array with. + Used to preprocess the design weights prior to convolution with the filter. + Border padding an image will set the value of each padded pixel equal to + the value of the nearest pixel in the image. Used to implement feature- + preserving convolution filters that prevent unwanted edge effects. - Returns: - A 2d array with border padding. - """ - # fill sides - left = npa.tile(arr[0, :], (pad[0][0], 1)) - right = npa.tile(arr[-1, :], (pad[0][1], 1)) - top = npa.tile(arr[:, 0], (pad[1][0], 1)).transpose() - bottom = npa.tile(arr[:, -1], (pad[1][1], 1)).transpose() + Args: + arr: 2d array whose borders contain the values to use for padding + pad: 2x2 array of integers indicating the size + of the borders to pad the array with. - # fill corners - top_left = npa.tile(arr[0, 0], (pad[0][0], pad[1][0])) - top_right = npa.tile(arr[-1, 0], (pad[0][1], pad[1][0])) - bottom_left = npa.tile(arr[0, -1], (pad[0][0], pad[1][1])) - bottom_right = npa.tile(arr[-1, -1], (pad[0][1], pad[1][1])) + Returns: + A 2d array with border padding. + """ + # fill sides + left = npa.tile(arr[0, :], (pad[0][0], 1)) + right = npa.tile(arr[-1, :], (pad[0][1], 1)) + top = npa.tile(arr[:, 0], (pad[1][0], 1)).transpose() + bottom = npa.tile(arr[:, -1], (pad[1][1], 1)).transpose() + + # fill corners + top_left = npa.tile(arr[0, 0], (pad[0][0], pad[1][0])) + top_right = npa.tile(arr[-1, 0], (pad[0][1], pad[1][0])) + bottom_left = npa.tile(arr[0, -1], (pad[0][0], pad[1][1])) + bottom_right = npa.tile(arr[-1, -1], (pad[0][1], pad[1][1])) + + if pad[0][0] > 0 and pad[0][1] > 0 and pad[1][0] > 0 and pad[1][1] > 0: + return npa.concatenate( + ( + npa.concatenate((top_left, top, top_right)), + npa.concatenate((left, arr, right)), + npa.concatenate((bottom_left, bottom, bottom_right)), + ), + axis=1, + ) + elif pad[0][0] == 0 and pad[0][1] == 0 and pad[1][0] > 0 and pad[1][1] > 0: + return npa.concatenate((top, arr, bottom), axis=1) + elif pad[0][0] > 0 and pad[0][1] > 0 and pad[1][0] == 0 and pad[1][1] == 0: + return npa.concatenate((left, arr, right), axis=0) + elif pad[0][0] == 0 and pad[0][1] == 0 and pad[1][0] == 0 and pad[1][1] == 0: + return arr + else: + raise ValueError("At least one of the padding numbers is invalid.") + + def _quarter_to_full_kernel(self, arr: npa.ndarray, pad_to: npa.ndarray) -> npa.ndarray: + """Constructs the full kernel from its nonnegative quadrant. + + Args: + arr: 2d input array representing the nonnegative quadrant of a + filter kernel with C4v symmetry. + pad_to: 1d array with two elements (integers) specifying the size + of the zero padding. + + Returns: + The complete kernel. + """ + pad_size = pad_to - 2 * npa.array(arr.shape) + 1 + + top = npa.zeros((pad_size[0], arr.shape[1])) + bottom = npa.zeros((pad_size[0], arr.shape[1] - 1)) + middle = npa.zeros((pad_to[0], pad_size[1])) + + top_left = arr[:, :] + top_right = npa.flipud(arr[1:, :]) + bottom_left = npa.fliplr(arr[:, 1:]) + bottom_right = npa.flipud( + npa.fliplr(arr[1:, 1:]) + ) # equivalent to flip, but flip is incompatible with autograd - if pad[0][0] > 0 and pad[0][1] > 0 and pad[1][0] > 0 and pad[1][1] > 0: return npa.concatenate( ( npa.concatenate((top_left, top, top_right)), - npa.concatenate((left, arr, right)), + middle, npa.concatenate((bottom_left, bottom, bottom_right)), ), axis=1, ) - elif pad[0][0] == 0 and pad[0][1] == 0 and pad[1][0] > 0 and pad[1][1] > 0: - return npa.concatenate((top, arr, bottom), axis=1) - elif pad[0][0] > 0 and pad[0][1] > 0 and pad[1][0] == 0 and pad[1][1] == 0: - return npa.concatenate((left, arr, right), axis=0) - elif pad[0][0] == 0 and pad[0][1] == 0 and pad[1][0] == 0 and pad[1][1] == 0: - return arr - else: - raise ValueError("At least one of the padding numbers is invalid.") -def _quarter_to_full_kernel(arr: npa.ndarray, pad_to: npa.ndarray) -> npa.ndarray: - """Constructs the full kernel from its nonnegative quadrant. + def _construct_kernel(self,radius,space_grid): + '''Applies a conic filter to the input parameters. The filter is defined as a + conic region of a specified radius, where the values undergo a weighted average. - Args: - arr: 2d input array representing the nonnegative quadrant of a - filter kernel with C4v symmetry. - pad_to: 1d array with two elements (integers) specifying the size - of the zero padding. + Args: + input_parameters: the input parameters to be filtered (2D) + radius: the radius of the conic filter in micrometers + space_grid: a SpaceGrid object containing the size and resolution of the parameters + + Returns: + The filtered parameters (2D) + ''' + Nx = space_grid.get_Nx() + Ny = space_grid.get_Ny() - Returns: - The complete kernel. - """ - pad_size = pad_to - 2 * npa.array(arr.shape) + 1 + #positive quadrant of the filter + xv = npa.arange(0, Nx//2+1) * (space_grid.get_resolution()) + yv = npa.arange(0, Ny//2+1) * (space_grid.get_resolution()) + #note that the filter's center is at (0,0) + Xv, Yv = npa.meshgrid(xv, yv) - top = npa.zeros((pad_size[0], arr.shape[1])) - bottom = npa.zeros((pad_size[0], arr.shape[1] - 1)) - middle = npa.zeros((pad_to[0], pad_size[1])) + h = npa.where(Xv**2 + Yv**2 < radius**2, (1 - npa.sqrt(abs(Xv**2 + Yv**2)) / radius),0) - top_left = arr[:, :] - top_right = npa.flipud(arr[1:, :]) - bottom_left = npa.fliplr(arr[:, 1:]) - bottom_right = npa.flipud( - npa.fliplr(arr[1:, 1:]) - ) # equivalent to flip, but flip is incompatible with autograd + #impulse response is padded to [Fx,Fy] + h = self._quarter_to_full_kernel(h, npa.array([self.Fx + Nx - 1, self.Fy + Ny - 1])) + + print("Fx, Fy:", self.Fx, self.Fy) - return npa.concatenate( - ( - npa.concatenate((top_left, top, top_right)), - middle, - npa.concatenate((bottom_left, bottom, bottom_right)), - ), - axis=1, - ) - -def conic_filter(phase_mask,radius,space_grid): - '''Applies a conic filter to the input phase mask. The filter is defined as a - conic region of a specified radius, where the values are averaged. - - Args: - phase_mask: the input phase mask to be filtered (2D) - radius: the radius of the conic filter in micrometers - space_grid: an object containing the size and resolution of the phase mask - - Returns: - The filtered phase mask (2D) - ''' - Nx = space_grid.get_Nx() - Ny = space_grid.get_Ny() - - #positive quadrant of the filter - xv = npa.arange(0, Nx//2+1) * (space_grid.get_resolution()) - yv = npa.arange(0, Ny//2+1) * (space_grid.get_resolution()) - #note that the filter's center is at (0,0) - Xv, Yv = npa.meshgrid(xv, yv) - - h = npa.where(Xv**2 + Yv**2 < radius**2, (1 - npa.sqrt(abs(Xv**2 + Yv**2)) / radius),0) - - h = _quarter_to_full_kernel(h, 3 * npa.array([Nx, Ny])) - - h = h / npa.sum(h) - #make a plot of the filter kernel - H = npa.fft.fft2(h) - - phase_mask_2d = phase_mask.reshape(Nx,Ny) - phase_mask_padded = _edge_pad(phase_mask_2d, ((Nx, Nx), (Ny, Ny))) - - phase_mask_fft = npa.fft.fft2(phase_mask_padded) - - Y = phase_mask_fft * H - - y = _centered( - npa.real(npa.fft.ifft2(Y)), - (Nx, Ny) - ) - - return y.flatten() \ No newline at end of file + h = h / npa.sum(h) + H = npa.fft.fft2(h) + + return H \ No newline at end of file diff --git a/mpid.py b/mpid.py index 2bf5083..a648be2 100644 --- a/mpid.py +++ b/mpid.py @@ -7,9 +7,9 @@ import pyfftw import nlopt from autograd import tensor_jacobian_product, grad as gd import autograd.numpy as npa -from filters import conic_filter +from filters import Filter -def make_2Dplot_of(given_field:np.ndarray, grid_X:np.ndarray, grid_Y:np.ndarray, plot_zoom_x:float,plot_zoom_y:float,choose_quantity="amplitude",save_name="plot") -> None: +def make_2Dplot_of(given_field:np.ndarray, grid_X:np.ndarray, grid_Y:np.ndarray, plot_zoom_x:float, plot_zoom_y:float, choose_quantity="amplitude", save_name="plot") -> None: '''Makes a plot of the given 2D field in a zoomed region. Saves the plot with the chosen filename to a PDF file. Args: @@ -25,7 +25,6 @@ def make_2Dplot_of(given_field:np.ndarray, grid_X:np.ndarray, grid_Y:np.ndarray, fig, ax = plt.subplots(figsize=(6, 6)) if choose_quantity == "amplitude": ax.pcolormesh(grid_X*1e6,grid_Y*1e6,np.abs(given_field_2d),rasterized=True) - #add a colorbar for amplitude plt.colorbar(ax.pcolormesh(grid_X*1e6,grid_Y*1e6,np.abs(given_field_2d),rasterized=True), ax=ax) elif choose_quantity == "phase": ax.pcolormesh(grid_X*1e6,grid_Y*1e6,(np.angle(given_field_2d)+2*np.pi)%(2*np.pi), cmap='hsv', vmin=0, vmax=2*np.pi, rasterized=True) @@ -595,9 +594,11 @@ class Simulation: self.filter_radius = None if conic_filter_radius is not None: - self.filter_radius = conic_filter_radius + self.filter_radius = conic_filter_radius + self.filter_obj = Filter(conic_filter_radius, self.space_grid) else: self.filter_radius = None + self.filter_obj = None def get_sources(self): return self.sources @@ -632,7 +633,8 @@ class Simulation: def _update_wavelength(self, new_wavelength): ''' Updates the operating wavelength and related vectors (lda_vector and k_vectors) based on the new wavelength. - The new_wavelength should be provided in meters. + The new_wavelength should be provided in meters. If use_unit_cells is True, updating the wavelength + will also update the current active unit cell. Args: new_wavelength (float): The new operating wavelength in meters. ''' @@ -771,9 +773,10 @@ class Simulation: def _normParam_filter(self, norm_param: np.ndarray) -> np.ndarray: '''Applies a conic filter on normalized parameters. Normalized parameters are the ones the optimization algorithm uses and are always in the range [0,1].''' - return conic_filter(norm_param, - self.filter_radius, - self.space_grid,) + # return conic_filter(norm_param, + # self.filter_radius, + # self.space_grid,) + return self.filter_obj.apply_filter(norm_param) def _phase_filter(self, phase_mask: np.ndarray) -> np.ndarray: '''Applies a conic filter on the phase mask. The filtering action is equivalent to a local weighted averaging of the phase. @@ -794,17 +797,18 @@ class Simulation: phase_mask_cos = npa.cos(phase_mask) phase_mask_sin = npa.sin(phase_mask) - phase_mask_filtered_cos = conic_filter( - phase_mask_cos, - self.filter_radius, - self.space_grid, - ) - - phase_mask_filtered_sin = conic_filter( - phase_mask_sin, - self.filter_radius, - self.space_grid, - ) + # phase_mask_filtered_cos = conic_filter( + # phase_mask_cos, + # self.filter_radius, + # self.space_grid, + # ) + phase_mask_filtered_cos = self.filter_obj.apply_filter(phase_mask_cos) + # phase_mask_filtered_sin = conic_filter( + # phase_mask_sin, + # self.filter_radius, + # self.space_grid, + # ) + phase_mask_filtered_sin = self.filter_obj.apply_filter(phase_mask_sin) return npa.arctan2(phase_mask_filtered_sin, phase_mask_filtered_cos) def _forward_propagate(self, source_object: Source, returns_field: bool = False): @@ -1184,7 +1188,8 @@ class Simulation: #after the phase mask is set, we test the same mask across all targets C_all_wavelengths = [] tested_wavelengths = [] - grad_partial = np.zeros_like(grad) + #per-mask gradient w.r.t. the filter output, accumulated over targets (the filter Jacobian is target-independent) + grads_wrt_filtered = None for target in self.targets: @@ -1197,22 +1202,27 @@ class Simulation: full_grads = [g * self.current_unit_cell.get_complete_gradient_from_norm_param(fw) for g, fw in zip(full_grads, filtered_w_masks)] - if self.use_filtering: - #backpropagate through the filter, which acts on w (unit cells) or on the phase (no unit cells) - if self.use_unit_cells: - filter_fun, filter_inputs = self._normParam_filter, w_masks - else: - filter_fun, filter_inputs = self._phase_filter, phase_masks - jac_grads = [tensor_jacobian_product(filter_fun)(inp, g) for inp, g in zip(filter_inputs, full_grads)] + #the gradient depends on wavelength, because phase of the unit cell is wavelength-dependent. + #Each mask accumulates only its own gradient. + if grads_wrt_filtered is None: + grads_wrt_filtered = [np.array(g, dtype=float) for g in full_grads] else: - jac_grads = list(full_grads) + for k, g in enumerate(full_grads): + grads_wrt_filtered[k] += g - if self.space_grid.xy_symmetry: - jac_grads = [self._reduce_symmetrized_gradient(j) for j in jac_grads] - jac_grad = np.concatenate(jac_grads, axis=0) # = (dC / dPhi) * (dw_f / dw), recall Phi = 2pi * w_f + #a single backward pass per mask through the filter, which acts on w (unit cells) or on the phase (no unit cells) + if self.use_filtering: + if self.use_unit_cells: + filter_fun, filter_inputs = self._normParam_filter, w_masks + else: + filter_fun, filter_inputs = self._phase_filter, phase_masks + jac_grads = [tensor_jacobian_product(filter_fun)(inp, g) for inp, g in zip(filter_inputs, grads_wrt_filtered)] + else: + jac_grads = list(grads_wrt_filtered) - #the gradient depends on wavelength, because phase of the unit cell is wavelength-dependent - grad_partial[:] += jac_grad + if self.space_grid.xy_symmetry: + jac_grads = [self._reduce_symmetrized_gradient(j) for j in jac_grads] + grad_partial = np.concatenate(jac_grads, axis=0) # = sum over targets of (dC / dPhi) * (dw_f / dw), recall Phi = 2pi * w_f if grad.size > 0: @@ -1454,64 +1464,125 @@ class Simulation: self.set_space_grid(new_space_grid) - def inject_phase_masks(self, input_phase_masks: list): + def inject_phase_masks(self, input_phase_masks: list, apply_filtering: bool = True): """ - Inject the specified phase masks into the simulation. + Inject the specified phase masks into the simulation, optionally applying filtering based on the `apply_filtering` flag. Args: input_phase_masks (list): A list of phase masks, flattened, to be injected into the simulation. + apply_filtering (bool): Apply filtering to the injected phase masks. Defaults to True. Returns: None """ - #check if the length of the list of input_phase_masks we are trying to inject - #has the same number of phase masks in this simulation object. if len(input_phase_masks) != self.num_phase_masks: raise ValueError(f"Expected {self.num_phase_masks} phase masks, but got {len(input_phase_masks)}.") + for i, phase_mask in enumerate(input_phase_masks): + #set the phase mask attribute with the input phase masks self.phase_mask[i*self.space_grid.get_S():(i+1)*self.space_grid.get_S()] = phase_mask - self.w_norm[i*self.space_grid.get_S():(i+1)*self.space_grid.get_S()] = self._w_given_phase(phase_mask) - if self.use_filtering: + #self.w_norm[i*self.space_grid.get_S():(i+1)*self.space_grid.get_S()] = self._w_given_phase(phase_mask) + + if apply_filtering: #filtering is applied regardless of the global optimization setting + if self.use_unit_cells: + selected_filter = lambda x: self._normParam_filter(x) + else: + selected_filter = lambda x: self._phase_filter(x) if self.num_phase_masks == 2: - filtered_phase_mask_1 = self._phase_filter(self.phase_mask[:self.space_grid.get_S()]) - filtered_phase_mask_2 = self._phase_filter(self.phase_mask[self.space_grid.get_S():]) + filtered_phase_mask_1 = selected_filter(self.phase_mask[:self.space_grid.get_S()]) + filtered_phase_mask_2 = selected_filter(self.phase_mask[self.space_grid.get_S():]) self.filtered_phase_mask = np.concatenate((filtered_phase_mask_1, filtered_phase_mask_2)) elif self.num_phase_masks == 3: - filtered_phase_mask_1 = self._phase_filter(self.phase_mask[:self.space_grid.get_S()]) - filtered_phase_mask_2 = self._phase_filter(self.phase_mask[self.space_grid.get_S():2*self.space_grid.get_S()]) - filtered_phase_mask_3 = self._phase_filter(self.phase_mask[2*self.space_grid.get_S():]) + filtered_phase_mask_1 = selected_filter(self.phase_mask[:self.space_grid.get_S()]) + filtered_phase_mask_2 = selected_filter(self.phase_mask[self.space_grid.get_S():2*self.space_grid.get_S()]) + filtered_phase_mask_3 = selected_filter(self.phase_mask[2*self.space_grid.get_S():]) self.filtered_phase_mask = np.concatenate((filtered_phase_mask_1, filtered_phase_mask_2, filtered_phase_mask_3)) else: self.filtered_phase_mask = self.phase_mask + - def source_sweep(self, distances: np.ndarray = None, input_phase_masks: list = None): + def _update_phase_masks_with_diameters(self, input_diameters_masks: list, apply_filtering: bool = False): """ - This method runs the propagation algorithm for the specified sources and targets in the simulation object. - For every propagation step, it enforces the same phase masks. It can optionally update propagation distances and phase masks. - If not specified, it is assumed the simulation object has already been optimized and that the required quantities are already set. - It is useful for studying the wavelength-dependent performance of the system. + Constructs the phase masks based on the provided diameter masks and the currently active unit cell. Args: - distances (np.ndarray, optional): Array containing the propagation distances. If None, the current distances are used. + input_diameters_masks (list): A list of diameter masks, flattened, to be used for updating the phase masks. + Returns: + None + """ + if len(input_diameters_masks) != self.num_phase_masks: + raise ValueError(f"Expected {self.num_phase_masks} masks, but got {len(input_diameters_masks)}.") + + #check if unit current unit cell is set and not None + if self.current_unit_cell is None: + raise ValueError("Cannot generate phase masks from diameters if no unit cell is provided.") + + mask1 = input_diameters_masks[0] + mask2 = input_diameters_masks[1] if self.num_phase_masks > 1 else None + mask3 = input_diameters_masks[2] if self.num_phase_masks > 2 else None + + phase_mask_1 = self.current_unit_cell.get_phase_at_param(mask1) + phase_mask_2 = self.current_unit_cell.get_phase_at_param(mask2) + phase_mask_3 = self.current_unit_cell.get_phase_at_param(mask3) if mask3 is not None else None + + phase_mask_list = [phase_mask_1, phase_mask_2, phase_mask_3][:self.num_phase_masks] + self.inject_phase_masks(phase_mask_list, apply_filtering=apply_filtering) + + def source_sweep(self, distances: np.ndarray = None, input_phase_masks: list = None, input_diameter_masks: list = None, apply_filtering: bool = True, save_output: bool = False): + """ + Runs the propagation algorithm for the specified sources and targets in the simulation object. + There are two operating modes of this method: if phase masks are provided via the `input_phase_masks` argument, + If 'input_diameter_masks' is provided, phase masks are + updated based on the source wavelength and unit cell data. If none of the two are provided, the last phase masks + One cannot provide both `input_phase_masks` and `input_diameter_masks` simultaneously. + Example: if the simulation contains two targets T1 and T2, with associated sources S1, S2, at two distinct wavelengths L1, L2, passing input_diameters_masks + will update the phase masks for each source by selecting the appropriate unit cell based on the source wavelength and the provided diameter mask. + Unit cell data must be already be present in the simulation object. + It is typically used for studying the wavelength-dependent performance of the system. + Args: + distances (np.ndarray, optional): Array containing the propagation distances. If None, the current distances in the simulation object are used. input_phase_masks (list, optional): A list of phase masks to be injected into the simulation. If None, the current phase masks are used. + input_diameter_masks (list, optional): A list of diameter masks to be used for updating the phase masks based on the source wavelength and unit cell data. + apply_filtering (bool, optional): Controls filtering with higher priority than Simulation settings. Defaults to True. + save_output (bool, optional): Save the output results to a file. Default is False. + Returns: dict: A dictionary mapping source wavelengths to the overlap between the output field and the associated target field. """ + + #setup + if input_diameter_masks is not None and input_phase_masks is not None: + raise ValueError("Cannot choose: both input_phase_masks and input_diameter_masks were provided.") + results_dict = {} if distances is not None: self.update_propagation_distances(distances) if input_phase_masks is not None: - self.inject_phase_masks(input_phase_masks) - + self.inject_phase_masks(input_phase_masks, apply_filtering=apply_filtering) + + if input_diameter_masks is not None: + use_unit_cells = True + else: + use_unit_cells = False + + #simulate for target in self.targets: source = target.get_associated_source() - self._update_wavelength(source.get_wavelength()) + wavelength = source.get_wavelength() + self._update_wavelength(wavelength) + #retrieve the appropriate phase mask from the simulation's active unit cell + if use_unit_cells: + self._update_phase_masks_with_diameters(input_diameter_masks, apply_filtering=apply_filtering) output_field = self._forward_propagate(source, returns_field = True) #compute the overlap between the output field and associated target field target_field = target.get_field(self.get_space_grid(), flatten = True) overlap = np.abs(np.sum(np.conj(target_field) * output_field))**2 results_dict[source.get_wavelength()] = overlap + if save_output: + save_filename = f"{self.results_folder}/output_field_{wavelength*1e6}" + save_filename = save_filename.replace('.', '_') + ".npy" + np.save(save_filename, output_field) return results_dict diff --git a/single_converter_symm.py b/single_converter_symm.py index a952069..d6f3b12 100644 --- a/single_converter_symm.py +++ b/single_converter_symm.py @@ -25,7 +25,7 @@ d = np.array([d0, d1, d2]) ms_max_available_diameter = 122 #[um] output_fiber_core_diameter = 19 #[um] -opt_max_eval = 1000 +opt_max_eval = 600 space_grid = pid.SpaceGrid( side_pixels=N, diff --git a/single_converter_symm_ego.py b/single_converter_symm_ego.py index e55c84d..f9540ba 100644 --- a/single_converter_symm_ego.py +++ b/single_converter_symm_ego.py @@ -30,7 +30,7 @@ mat_refr_index = np.array([n0, n1, n2]) ms_max_available_diameter = 122 #[um] output_fiber_core_diameter = 19 #[um] -opt_max_eval = 140 +opt_max_eval = 550 def setup_optical_system()->pid.Simulation: ''' @@ -47,7 +47,7 @@ def setup_optical_system()->pid.Simulation: ''' central_wavelength_um = 1.55 - wavelengths_to_test_um = [1.530, 1.550, 1.570] #[um] + wavelengths_to_test_um = [1.53, 1.55, 1.57] #[um] source_files = ["fiber_modes/SMF/SMF_LP01_field_distribution_lda_1530_N_1025.npy", "fiber_modes/SMF/SMF_LP01_field_distribution_lda_1550_N_1025.npy", @@ -65,7 +65,20 @@ def setup_optical_system()->pid.Simulation: associated_source=sources[i], wavelength=wvl, field_norm=True) for i, (file, wvl) in enumerate(zip(target_files, wavelengths_to_test_um))] - + + #import unit cell data + unit_cell_diameters_wvl_1530 = np.load("unit_cell/unit_cell_center_1530_h900nm_res50_wvl_diameters.npy") + unit_cell_diameters_wvl_1550 = np.load("unit_cell/unit_cell_center_1550_h900nm_res50_wvl_diameters.npy") + unit_cell_diameters_wvl_1570 = np.load("unit_cell/unit_cell_center_1570_h900nm_res50_wvl_diameters.npy") + unit_cell_phases_wvl_1530 = np.load("unit_cell/unit_cell_center_1530_h900nm_res50_wvl_mode_phase.npy") + unit_cell_phases_wvl_1550 = np.load("unit_cell/unit_cell_center_1550_h900nm_res50_wvl_mode_phase.npy") + unit_cell_phases_wvl_1570 = np.load("unit_cell/unit_cell_center_1570_h900nm_res50_wvl_mode_phase.npy") + #construct the UnitCell objects + unit_cell_1530 = pid.UnitCell(wavelength=1.53, tune_param=unit_cell_diameters_wvl_1530, sim_phase=unit_cell_phases_wvl_1530) + unit_cell_1550 = pid.UnitCell(wavelength=1.55, tune_param=unit_cell_diameters_wvl_1550, sim_phase=unit_cell_phases_wvl_1550) + unit_cell_1570 = pid.UnitCell(wavelength=1.57, tune_param=unit_cell_diameters_wvl_1570, sim_phase=unit_cell_phases_wvl_1570) + unit_cells_list = [unit_cell_1530, unit_cell_1550, unit_cell_1570] + initial_space_grid = pid.SpaceGrid( side_pixels=N, unit_cell_pitch=unit_cell_pitch, @@ -89,6 +102,7 @@ def setup_optical_system()->pid.Simulation: phase_mask_init_size=output_fiber_core_diameter/3, opt_max_eval=opt_max_eval, conic_filter_radius=1.5, + unit_cells=unit_cells_list, results_while_solving = False, save_results = False, ) diff --git a/single_converter_symm_source_sweep.py b/single_converter_symm_source_sweep.py index 16639ea..0f64b74 100644 --- a/single_converter_symm_source_sweep.py +++ b/single_converter_symm_source_sweep.py @@ -17,6 +17,18 @@ ms_max_available_diameter = 122 #[um] output_fiber_core_diameter = 19 #[um] opt_max_eval = 140 +central_wavelength_um = 1.55 #[um] + +wavelengths_to_test_um = [1.530, 1.550, 1.570] #[um] + +source_files = ["fiber_modes/SMF/SMF_LP01_field_distribution_lda_1530_N_1025.npy", + "fiber_modes/SMF/SMF_LP01_field_distribution_lda_1550_N_1025.npy", + "fiber_modes/SMF/SMF_LP01_field_distribution_lda_1570_N_1025.npy"] + +target_files = ["fiber_modes/FMF/FMF_LP21_field_distribution_lda_1530_N_1025.npy", + "fiber_modes/FMF/FMF_LP21_field_distribution_lda_1550_N_1025.npy", + "fiber_modes/FMF/FMF_LP21_field_distribution_lda_1570_N_1025.npy"] + def setup_optical_system()->pid.Simulation: ''' Builds the mode converter optical system using @@ -30,26 +42,27 @@ def setup_optical_system()->pid.Simulation: Returns: float: The best objective value obtained from the optimization. ''' - central_wavelength_um = 1.55 - - wavelengths_to_test_um = [1.530, 1.550, 1.570] #[um] - - source_files = ["fiber_modes/SMF/SMF_LP01_field_distribution_lda_1530_N_1025.npy", - "fiber_modes/SMF/SMF_LP01_field_distribution_lda_1550_N_1025.npy", - "fiber_modes/SMF/SMF_LP01_field_distribution_lda_1570_N_1025.npy"] - sources = [pid.Source(input_field=np.load(file), wavelength=wvl, field_norm=True) for file, wvl in zip(source_files, wavelengths_to_test_um)] - target_files = ["fiber_modes/FMF/FMF_LP21_field_distribution_lda_1530_N_1025.npy", - "fiber_modes/FMF/FMF_LP21_field_distribution_lda_1550_N_1025.npy", - "fiber_modes/FMF/FMF_LP21_field_distribution_lda_1570_N_1025.npy"] - targets = [pid.Target(input_field=np.load(file), associated_source=sources[i], wavelength=wvl, field_norm=True) for i, (file, wvl) in enumerate(zip(target_files, wavelengths_to_test_um))] + + #import unit cell data + unit_cell_diameters_wvl_1530 = np.load("unit_cell/unit_cell_center_1530_h900nm_res50_wvl_diameters.npy") + unit_cell_diameters_wvl_1550 = np.load("unit_cell/unit_cell_center_1550_h900nm_res50_wvl_diameters.npy") + unit_cell_diameters_wvl_1570 = np.load("unit_cell/unit_cell_center_1570_h900nm_res50_wvl_diameters.npy") + unit_cell_phases_wvl_1530 = np.load("unit_cell/unit_cell_center_1530_h900nm_res50_wvl_mode_phase.npy") + unit_cell_phases_wvl_1550 = np.load("unit_cell/unit_cell_center_1550_h900nm_res50_wvl_mode_phase.npy") + unit_cell_phases_wvl_1570 = np.load("unit_cell/unit_cell_center_1570_h900nm_res50_wvl_mode_phase.npy") + #construct the UnitCell objects + unit_cell_1530 = pid.UnitCell(wavelength=1.530, tune_param=unit_cell_diameters_wvl_1530, sim_phase=unit_cell_phases_wvl_1530) + unit_cell_1550 = pid.UnitCell(wavelength=1.550, tune_param=unit_cell_diameters_wvl_1550, sim_phase=unit_cell_phases_wvl_1550) + unit_cell_1570 = pid.UnitCell(wavelength=1.570, tune_param=unit_cell_diameters_wvl_1570, sim_phase=unit_cell_phases_wvl_1570) + unit_cells_list = [unit_cell_1530, unit_cell_1550, unit_cell_1570] initial_space_grid = pid.SpaceGrid( side_pixels=N, @@ -73,7 +86,8 @@ def setup_optical_system()->pid.Simulation: ms_max_avail_diameter=ms_max_available_diameter, phase_mask_init_size=output_fiber_core_diameter/3, opt_max_eval=opt_max_eval, - #conic_filter_radius=1.0, + conic_filter_radius=1.0, + unit_cells=unit_cells_list, results_while_solving = False, save_results = False, ) @@ -84,17 +98,19 @@ def setup_optical_system()->pid.Simulation: if __name__ == "__main__": simulation = setup_optical_system() - results_folder_name = 'results_nf' + results_folder_name = 'results' - phase_mask_1 = np.load(f"{results_folder_name}/filtered_phase_mask_1.npy").flatten() - phase_mask_2 = np.load(f"{results_folder_name}/filtered_phase_mask_2.npy").flatten() + #phase_mask_1 = np.load(f"{results_folder_name}/filtered_phase_mask_1.npy").flatten() + #phase_mask_2 = np.load(f"{results_folder_name}/filtered_phase_mask_2.npy").flatten() + mask1 = np.load(f"{results_folder_name}/mask_1.npy").flatten() + mask2 = np.load(f"{results_folder_name}/mask_2.npy").flatten() d0 = 195 #propagation distance: source - MS1 [um] FILTERED mf (using LP01 of SMF) d1 = 129 #propagation distance: MS1 - MS2 [um] FILTERED mf (using LP01 of SMF) d2 = d0 #propagation distance: MS2 - target [um] d = np.array([d0, d1, d2]) - results = simulation.source_sweep(distances = d, input_phase_masks=[phase_mask_1, phase_mask_2]) + results = simulation.source_sweep(distances = d, input_diameter_masks=[mask1, mask2], save_output=True, apply_filtering=False) #make a plot of the results, Use the keys as x axis and values as y axis fig, ax = plt.subplots(figsize=(6, 6)) ax.plot(np.array(list(results.keys()))*1e9, np.array(list(results.values())), marker='o', linewidth=2) @@ -105,4 +121,108 @@ if __name__ == "__main__": ax.tick_params(axis='both', which='major', labelsize=15) ax.set_aspect('auto') plt.tight_layout() - plt.savefig(f"{results_folder_name}/source_sweep_results.pdf") \ No newline at end of file + plt.savefig(f"{results_folder_name}/source_sweep_results.pdf") + + #for each of the test wavelength, take the fields from the results folder and make plots of them + for i, wvl in enumerate(wavelengths_to_test_um): + save_filename = f"{results_folder_name}/output_field_{wvl}" + save_filename = save_filename.replace('.', '_') + output_field = np.load(save_filename+ ".npy").reshape(simulation.space_grid.get_Nx(), simulation.space_grid.get_Ny()) + + fig, ax = plt.subplots(figsize=(6, 6)) + ax.pcolormesh(simulation.space_grid.get_X()*1e6, simulation.space_grid.get_Y()*1e6, np.abs(output_field) / np.max(np.abs(output_field)), cmap='viridis', rasterized = True) + colorbar = fig.colorbar(ax.collections[0], ax=ax, label='Magnitude [a.u.]') + ax.set_aspect('equal') + ax.tick_params(axis='x', labelsize=20) + ax.tick_params(axis='y', labelsize=20) + #set the tick fontsize to 20 + ax.tick_params(axis='both', which='major', labelsize=20) + ax.set_xlim(-output_fiber_core_diameter/2-10, output_fiber_core_diameter/2+10) + ax.set_ylim(-output_fiber_core_diameter/2-10, output_fiber_core_diameter/2+10) + #overlay a red circle indicating the output fiber core + circle = plt.Circle((0, 0), output_fiber_core_diameter/2, color='red', fill=False, linewidth=2) + ax.add_artist(circle) + ax.set_xlabel('X [um]', fontsize=20) + ax.set_ylabel('Y [um]', fontsize=20) + plt.title(r'Output Field at $\lambda=$' + str(wvl) + r' $\mu$m') + plt.tight_layout() + plt.savefig(save_filename + ".pdf") + plt.close() + + fig, ax = plt.subplots(figsize=(6, 6)) + ax.pcolormesh(simulation.space_grid.get_X()*1e6, simulation.space_grid.get_Y()*1e6, (np.angle(output_field)+2*np.pi)%(2*np.pi), cmap='hsv', vmin=0, vmax=2*np.pi, rasterized = True) + colorbar = fig.colorbar(ax.collections[0], ax=ax, label='Phase [rad]') + ax.set_aspect('equal') + ax.tick_params(axis='x', labelsize=20) + ax.tick_params(axis='y', labelsize=20) + #set the tick fontsize to 20 + ax.tick_params(axis='both', which='major', labelsize=20) + ax.set_xlim(-output_fiber_core_diameter/2-10, output_fiber_core_diameter/2+10) + ax.set_ylim(-output_fiber_core_diameter/2-10, output_fiber_core_diameter/2+10) + #overlay a red circle indicating the output fiber core + circle = plt.Circle((0, 0), output_fiber_core_diameter/2, color='red', fill=False, linewidth=2) + ax.add_artist(circle) + ax.set_xlabel('X [um]', fontsize=20) + ax.set_ylabel('Y [um]', fontsize=20) + plt.title(r'Output Field Phase at $\lambda=$' + str(wvl) + r' $\mu$m') + plt.tight_layout() + plt.savefig(save_filename + "_phase.pdf") + plt.close() + + + target_field_filename = target_files[i] + target_field = np.load(target_field_filename) + fig, axs = plt.subplots(2, 2, figsize=(12, 12)) + axs = axs.flatten() + # First subplot: output field amplitude + im = axs[0].pcolormesh(simulation.space_grid.get_X()*1e6, simulation.space_grid.get_Y()*1e6, np.abs(output_field)/np.max(np.abs(output_field)), cmap='viridis', shading='auto', rasterized=True) + fig.colorbar(im, ax=axs[0], label='Magnitude [a.u.]') + axs[0].set_aspect('equal') + axs[0].set_xlim(-output_fiber_core_diameter/2-10, output_fiber_core_diameter/2+10) + axs[0].set_ylim(-output_fiber_core_diameter/2-10, output_fiber_core_diameter/2+10) + circle = plt.Circle((0, 0), output_fiber_core_diameter/2, color='red', fill=False, linewidth=2) + axs[0].add_artist(circle) + axs[0].set_xlabel('X [um]', fontsize=20) + axs[0].set_ylabel('Y [um]', fontsize=20) + axs[0].tick_params(axis='x', labelsize=20) + axs[0].tick_params(axis='y', labelsize=20) + axs[0].tick_params(axis='both', which='major', labelsize=20) + axs[0].set_title('Output Field Amplitude') + axs[0].set_aspect('equal') + # Second subplot: y cross section slice + y_slice = output_field[:, output_field.shape[1]//2] + y_target_slice = target_field[:, target_field.shape[1]//2] + axs[1].plot(simulation.space_grid.get_Y()[:, 0]*1e6, np.abs(y_slice)/np.max(np.abs(y_slice)), label='Output', c='blue', linewidth=2) + axs[1].plot(simulation.space_grid.get_Y()[:, 0]*1e6, np.abs(y_target_slice)/np.max(np.abs(y_target_slice)), label='Target', c='red', linewidth=2) + axs[1].set_xlabel('Y [um]', fontsize=20) + axs[1].set_ylabel('Magnitude [a.u.]', fontsize=20) + axs[1].set_title('Y Cross Section') + axs[1].set_ylim(0, 1) + axs[1].set_xlim(-output_fiber_core_diameter/2-10, output_fiber_core_diameter/2+10) + axs[1].legend(fontsize=20) + axs[1].set_aspect('auto') + axs[1].tick_params(axis='x', labelsize=20) + axs[1].tick_params(axis='y', labelsize=20) + axs[1].tick_params(axis='both', which='major', labelsize=20) + # Third subplot: x cross section slice + x_slice = output_field[output_field.shape[0]//2, :] + x_target_slice = target_field[target_field.shape[0]//2, :] + axs[2].plot(simulation.space_grid.get_X()[0, :]*1e6, np.abs(x_slice)/np.max(np.abs(x_slice)), label='Output', c='blue', linewidth=2) + axs[2].plot(simulation.space_grid.get_X()[0, :]*1e6, np.abs(x_target_slice)/np.max(np.abs(x_target_slice)), label='Target', c='red', linewidth=2) + axs[2].set_xlabel('X [um]', fontsize=20) + axs[2].set_ylim(0, 1) + axs[2].set_xlim(-output_fiber_core_diameter/2-10, output_fiber_core_diameter/2+10) + axs[2].set_ylabel('Magnitude [a.u.]', fontsize=20) + axs[2].set_title('X Cross Section') + axs[2].legend(fontsize=20) + axs[2].set_aspect('auto') + axs[2].tick_params(axis='x', labelsize=20) + axs[2].tick_params(axis='y', labelsize=20) + axs[2].tick_params(axis='both', which='major', labelsize=20) + #remove completely the fourth subplot axis since we only have three now + fig.delaxes(axs[3]) + + plt.tight_layout() + plt.savefig(save_filename + "_comparison.pdf") + plt.close() + diff --git a/testing_filter_radius.py b/testing_filter_radius.py new file mode 100644 index 0000000..2b3830f --- /dev/null +++ b/testing_filter_radius.py @@ -0,0 +1,55 @@ +#use this file to import some test phase masks and observe +#the effect of changing the filtering radius +import numpy as np +import matplotlib.pyplot as plt +from autograd import tensor_jacobian_product + +class TempGrid: + def __init__(self, Nx, Ny, res): self.Nx, self.Ny, self.res = Nx, Ny, res + def get_Nx(self): return self.Nx + def get_Ny(self): return self.Ny + def get_resolution(self): return self.res + + +grid = TempGrid(Nx=1025, Ny=1025, res=0.7) +radius_to_test = [2.0, 5.0, 11.0] + +input_mask = np.load('results/mask_2.npy') +#make a plot of the input mask +# fig, ax = plt.subplots() +# ax.imshow(input_mask, cmap='inferno') +# #zoom in the center region +# center_x, center_y = grid.get_Nx() // 2, grid.get_Ny() // 2 +# zoom_size = 100 # adjust as needed +# ax.set_xlim(center_x - zoom_size, center_x + zoom_size) +# ax.set_ylim(center_y - zoom_size, center_y + zoom_size) +# plt.savefig('input_mask.pdf') + +#apply the filter +from filters import Filter + +center_x, center_y = grid.get_Nx() // 2, grid.get_Ny() // 2 +zoom_size = 80 + +colormap = 'gray_r' +fig, ax = plt.subplots(1,4, figsize=(12, 3)) +ax[0].imshow(input_mask, cmap=colormap) +ax[0].set_title('Input', fontsize=20) +ax[0].set_xlim(center_x - zoom_size, center_x + zoom_size) +ax[0].set_ylim(center_y - zoom_size, center_y + zoom_size) +ax[0].set_xticks([]) +ax[0].set_yticks([]) +for i, radius in enumerate(radius_to_test, start=1): + filt = Filter(radius, grid) + output_mask = filt.apply_filter(input_mask) + ax[i].imshow(output_mask.reshape(grid.get_Nx(), grid.get_Ny()), cmap=colormap) + ax[i].set_title(f'r={radius} \u03BCm', fontsize=20) + ax[i].set_xlim(center_x - zoom_size, center_x + zoom_size) + ax[i].set_ylim(center_y - zoom_size, center_y + zoom_size) + ax[i].set_xticks([]) + ax[i].set_yticks([]) +plt.tight_layout() +plt.savefig('filtered_masks.pdf') + + + diff --git a/unit_cell/param_sweep_center_1570.log b/unit_cell/param_sweep_center_1570.log new file mode 100644 index 0000000..a4b62c2 --- /dev/null +++ b/unit_cell/param_sweep_center_1570.log @@ -0,0 +1,29 @@ +nohup: ignoring input +Rank 0 is processing diameter 0.100 (index 0) +Rank 1 is processing diameter 0.361 (index 12) +Rank 0 is processing diameter 0.122 (index 1) +Rank 1 is processing diameter 0.383 (index 13) +Rank 0 is processing diameter 0.143 (index 2) +Rank 1 is processing diameter 0.404 (index 14) +Rank 0 is processing diameter 0.165 (index 3) +Rank 1 is processing diameter 0.426 (index 15) +Rank 0 is processing diameter 0.187 (index 4) +Rank 1 is processing diameter 0.448 (index 16) +Rank 0 is processing diameter 0.209 (index 5) +Rank 1 is processing diameter 0.470 (index 17) +Rank 0 is processing diameter 0.230 (index 6) +Rank 0 is processing diameter 0.252 (index 7) +Rank 1 is processing diameter 0.491 (index 18) +Rank 0 is processing diameter 0.274 (index 8) +Rank 1 is processing diameter 0.513 (index 19) +Rank 0 is processing diameter 0.296 (index 9) +Rank 1 is processing diameter 0.535 (index 20) +Rank 0 is processing diameter 0.317 (index 10) +Rank 1 is processing diameter 0.557 (index 21) +Rank 0 is processing diameter 0.339 (index 11) +Rank 1 is processing diameter 0.578 (index 22) +Rank 1 is processing diameter 0.600 (index 23) + +Elapsed run time = 17749.7805 s + +Elapsed run time = 17750.1165 s diff --git a/unit_cell/results_S21_p700nm_h900nm.txt b/unit_cell/results_S21_p700nm_h900nm.txt new file mode 100644 index 0000000..348e66d --- /dev/null +++ b/unit_cell/results_S21_p700nm_h900nm.txt @@ -0,0 +1,25 @@ +% Model: cyljnderUnjtCell_solved.mph +% Versjon: COMSOL 6.1.0.282 +% Date: Aug 7 2026, 11:04 +% Table: Table 3 - Global Evaluatjon 3 +% tPost (nm) R (nm) lambda0 (µm) S21 (1) +900 50 1.55 -0.845669386929+0.5044362611307165j +900 63.157894736842096 1.55 -0.8251549977113691+0.5387569255954204j +900 76.31578947368419 1.55 -0.7941914472021001+0.5854009666643545j +900 89.47368421052632 1.55 -0.7471355423796626+0.6469515347910716j +900 102.6315789473684 1.55 -0.6739731961216431+0.7260451361240156j +900 115.78947368421052 1.55 -0.5559252634519462+0.8236596878659627j +900 128.94736842105263 1.55 -0.35623973597781955+0.9314716653309569j +900 142.1052631578947 1.55 -0.005446778866258663+0.9995911107511558j +900 155.26315789473682 1.55 0.5624693373844357+0.8192539801083378j +900 168.42105263157896 1.55 0.9716453344468975-0.026879486169522106j +900 181.57894736842104 1.55 0.33229983484123915-0.9109913324745208j +900 194.73684210526312 1.55 -0.5138804677648452-0.8375112454139356j +900 207.89473684210526 1.55 -0.8979812315659875-0.41599951560137755j +900 221.05263157894737 1.55 -0.8666228661809928-0.35839703156981834j +900 234.21052631578942 1.55 -0.9706534677777011+0.1672691323977627j +900 247.36842105263153 1.55 -0.875683679549086+0.4301009457755804j +900 260.52631578947364 1.55 -0.7186563387658285+0.6394472785795947j +900 273.6842105263158 1.55 -0.6152456609406248+0.5634675095025711j +900 286.8421052631579 1.55 -0.3669553952916399+0.7919528013254294j +900 300 1.55 -0.07785011501821787+0.7863645334693915j diff --git a/unit_cell/unitCell_grads.py b/unit_cell/unitCell_grads.py new file mode 100644 index 0000000..3b2cf97 --- /dev/null +++ b/unit_cell/unitCell_grads.py @@ -0,0 +1,51 @@ +''' +Utility script to test the gradients functionality of the unit cell class of mpid +Given the simulated transmittance phase values and associated tune parameter, +the unit cell class should return the gradient of the simulated phase with respect to the tuning parameter. +the class also provides a method to compute the gradient at a specific tuning parameter value, through interpolation. +''' + +import numpy as np +from matplotlib import pyplot as plt +from mpid import UnitCell + +wavelength_um_1_53 = 1.53 +wavelength_um_1_55 = 1.55 +wavelength_um_1_57 = 1.57 +unit_cell_diameters = np.load("unit_cell/unit_cell_center_1550_h900nm_res50_wvl_diameters.npy") +unit_cell_phases_wvl_1530 = np.load("unit_cell/unit_cell_center_1530_h900nm_res50_wvl_mode_phase.npy") +unit_cell_phases_wvl_1550 = np.load("unit_cell/unit_cell_center_1550_h900nm_res50_wvl_mode_phase.npy") +unit_cell_phases_wvl_1570 = np.load("unit_cell/unit_cell_center_1570_h900nm_res50_wvl_mode_phase.npy") + +unit_cell_1530 = UnitCell(wavelength=wavelength_um_1_53, tune_param=unit_cell_diameters, sim_phase=unit_cell_phases_wvl_1530) +unit_cell_1550 = UnitCell(wavelength=wavelength_um_1_55, tune_param=unit_cell_diameters, sim_phase=unit_cell_phases_wvl_1550) +unit_cell_1570 = UnitCell(wavelength=wavelength_um_1_57, tune_param=unit_cell_diameters, sim_phase=unit_cell_phases_wvl_1570) + +fig, ax = plt.subplots(figsize=(6, 6)) +plt.rcParams['mathtext.default'] = 'regular' +ax.plot(unit_cell_1530.get_tune_param(), unit_cell_1530.get_gradient(), color='r', label=r"$\lambda=1530\:[nm]$", linewidth=3) +ax.plot(unit_cell_1550.get_tune_param(), unit_cell_1550.get_gradient(), color='g', label=r"$\lambda=1550\:[nm]$", linewidth=3) +ax.plot(unit_cell_1570.get_tune_param(), unit_cell_1570.get_gradient(), color='b', label=r"$\lambda=1570\:[nm]$", linewidth=3) +ax.set_xlabel("diameter [um]", fontsize=20) +ax.set_ylabel(r"$\partial \phi/\partial d$", fontsize=20, rotation='vertical') +ax.yaxis.set_label_coords(-0.11, 0.5) +ax.set_aspect('auto') +ax.set_xticks(np.arange(min(unit_cell_diameters), max(unit_cell_diameters)+0.1, 0.1)) +ax.tick_params(axis='x', labelsize=20) +ax.tick_params(axis='y', labelsize=20) +plt.tight_layout() +plt.savefig("unit_cell_gradient.pdf") + +fig, ax = plt.subplots(figsize=(6, 6)) +ax.plot(unit_cell_1530.get_tune_param(), unit_cell_1530.get_sim_phase(), color='r', label=r"$\lambda=1530\:[nm]$", linewidth=3) +ax.plot(unit_cell_1550.get_tune_param(), unit_cell_1550.get_sim_phase(), color='g', label=r"$\lambda=1550\:[nm]$", linewidth=3) +ax.plot(unit_cell_1570.get_tune_param(), unit_cell_1570.get_sim_phase(), color='b', label=r"$\lambda=1570\:[nm]$", linewidth=3) +ax.set_xlabel("diameter [um]", fontsize=20) +ax.set_ylabel(r"$\phi\:$[rad]", fontsize=20) +ax.tick_params(axis='x', labelsize=20) +ax.tick_params(axis='y', labelsize=20) +ax.legend(fontsize=20, loc='best') +ax.set_aspect('auto') +ax.set_xticks(np.arange(min(unit_cell_diameters), max(unit_cell_diameters)+0.1, 0.1)) +plt.tight_layout() +plt.savefig("unit_cell_phase.pdf") \ No newline at end of file diff --git a/unit_cell/unit_cell_center_1530_h900nm_res50_wvl_diameters.npy b/unit_cell/unit_cell_center_1530_h900nm_res50_wvl_diameters.npy new file mode 100644 index 0000000000000000000000000000000000000000..b46d9806ba6a992787b3c72565a5e36e1f70c729 GIT binary patch literal 320 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh= zXCxM+0{I$7COVor3bhL41Fl&!L13r-=fpefy|wn+t4-y9To-%DKIo^xyX&e)?dxMf ztmF1;4i>X?Z98RuLlvaptUZ%!H19m+3-%JyZl zs(sOZkj!=ajzSRQhCNvH!A<)uS07k1x!$%v)ee%lV-Mm*^WL-n1~RVjzCHVSL5*8k f5A212>aTrvjF(z<^gnx#31wjdF8}Sd z)k_$sSu;5FM6KESm7CFF`qr8i4`rDgw!eLy>^X_qfvIEZjhCxg9j^JDIn4Z#!{Ons zYbuL^Sq>at+4D9|{nvi}c}=Wy4?fwy`&N9c(eXp}?(;XT-^}>TKJ8&^ec$)b_5z1% zY?kGHw!gYyx8O|uANI*Jk3DMr_S;@vuj#u$|9^XrfEiCK4VWFo-7L72e3=|h>=!$E QRf@%->Ws-Dg#soA08C4E2LJ#7 literal 0 HcmV?d00001 diff --git a/unit_cell/unit_cell_center_1530_h900nm_res50_wvl_mode_tran.npy b/unit_cell/unit_cell_center_1530_h900nm_res50_wvl_mode_tran.npy new file mode 100644 index 0000000000000000000000000000000000000000..9cb9dbd5965f4c67d875784698a7f681d06b73c2 GIT binary patch literal 320 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh= zXCxM+0{I$7COVor3bhL41FjP{U28Y6ytn7}3H`z<@ZR3?@ST;x^6%~EMmCye7{9j{ zQ~2gN)$hIi*SzByQbq6WXUyx1Z-42Ij!OH*_xAA> zk2xKr-rLvdc5HC{`qqBdoyjh8+VAY;7i=`R{q>!F3&*r&eqrzJo0$&T{}+5`AI!5( z;`;x0_8L!=-0W_>v(H_;SN(&-JNpX~yDH6F-r7H_UT?p+?Y(`-W)bZr?_St1@O`ss R)86OyC8qr>7n|PM0|1)id@%q3 literal 0 HcmV?d00001 diff --git a/unit_cell/unit_cell_center_1530_h900nm_res50_wvl_mode_wvl.npy b/unit_cell/unit_cell_center_1530_h900nm_res50_wvl_mode_wvl.npy new file mode 100644 index 0000000000000000000000000000000000000000..0d92f0327349dd1599b2001288d4063d8b8ce518 GIT binary patch literal 136 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh= eXCxM+0{I$-I+{8PwF(pfu4<8W?hmVe*aHAR;U97U literal 0 HcmV?d00001 diff --git a/unit_cell/unit_cell_center_1550_h900nm_res50_wvl_diameters.npy b/unit_cell/unit_cell_center_1550_h900nm_res50_wvl_diameters.npy new file mode 100644 index 0000000000000000000000000000000000000000..b46d9806ba6a992787b3c72565a5e36e1f70c729 GIT binary patch literal 320 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh= zXCxM+0{I$7COVor3bhL41Fl&!L13r-=fpefy|wn+t4-y9To-%DKIo^xyX&e)?dxMf ztmF1;4i>X?Z98RuLlvaptUZ%!H19m+3-%JyZl zs(sOZkj!=ajzSRQhCNvH!A<)uS07k1x!$%v)ee%lV-Mm*^WL-n1~RVjzCHVSL5*8k f5A212{h3$nJ38Xsc=4basbpjGm^)I5-@p zE!w9zZ!U*})W+MVXWZpDaN7Kt%<1dw2ef%>-?&w<9VmZ$vM9TS<-imEw4|_3#sm9I z4YcP={MsLt#JE2D_{;q!^BuSvUmUh)k4VZ$J@DSX)Z}BvCe5GrQu7~g+Gfq*VDo#$ zE8PwM?90yVij-_&cHmK))LYoa>TpTI^Q&AqyMwWzp3Z{L90vrp`dj!naU2j=t`9pm Qk?p{vU#h&ln^+D206wB{dH?_b literal 0 HcmV?d00001 diff --git a/unit_cell/unit_cell_center_1550_h900nm_res50_wvl_mode_tran.npy b/unit_cell/unit_cell_center_1550_h900nm_res50_wvl_mode_tran.npy new file mode 100644 index 0000000000000000000000000000000000000000..9e7bb1f001b6fa6ddb4366598d0588618280696e GIT binary patch literal 320 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh= zXCxM+0{I$7COVor3bhL411|Zw7wez?d1rt00K1|m`+NHYlZG!nqVMhJfB03=rSabW z!`7zIX-@C$#s9PhTcy6Y-^U){lhXIz-txhIkq3L;+Y8*O=J$W`-ah(?*DA}G@9i~X zr>*i!e{UbCWx9E-)jRvO=NX?Z98RuLlvaptUZ%!H19m+3-%JyZl zs(sOZkj!=ajzSRQhCNvH!A<)uS07k1x!$%v)ee%lV-Mm*^WL-n1~RVjzCHVSL5*8k f5A212gXSP!Mx6l`bB^3 zYdrRAh93TBKl$^mFgIlehp2T|YoEz7I-FPeZdsDWm^~-+4+7<6NtNpXzFLiY3j!%E=O{9ZL RlQuFs>CQN=QtMT4` zn@P!IafkQz#Z#vFW+%M2*T~WE&1!#d-{7Pv^>xd8`z)J?>i);??PrBAb65TG-hO(x z5Z}HP@9l56G<5u%`_6vdA=zb*@4dA*yuA6{iGX+Z8w^WJ6u!K(XMJ*sQ!?(oea6?i z6KTos?9Y7;w|Q{&o&9WvV+YONytBW(e)IpGLGSFHfAbZrnDf^D-67eOB&nzN+&hoD RmomPzZ>~-;y5RoY9sr^Ke)#|Z literal 0 HcmV?d00001 diff --git a/unit_cell/unit_cell_center_1570_h900nm_res50_wvl_mode_wvl.npy b/unit_cell/unit_cell_center_1570_h900nm_res50_wvl_mode_wvl.npy new file mode 100644 index 0000000000000000000000000000000000000000..8c4e34c7442c6ed58fcc18f2fdeed5bc615705b7 GIT binary patch literal 136 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh= eXCxM+0{I$-I+{8PwF(pfF8S8ifji`W+5-SU0v}2M literal 0 HcmV?d00001 diff --git a/unit_cell/unit_cell_center_wvl.py b/unit_cell/unit_cell_center_wvl.py new file mode 100644 index 0000000..eb971eb --- /dev/null +++ b/unit_cell/unit_cell_center_wvl.py @@ -0,0 +1,169 @@ +import meep as mp +mp.verbosity(0) + +import numpy as np +import matplotlib.pyplot as plt + +from mpi4py import MPI + +comm = MPI.COMM_WORLD +rank = comm.Get_rank() +size = comm.Get_size() + +resolution = 50 # pixels/unit length + +wvl_min = 1.53 +wvl_max = 1.57 +wvl_cen = 1.57 +fmin = 1/wvl_max +fmax = 1/wvl_min +fcen = 1/wvl_cen +df = fmax-fmin +nfreq = 1 + +dpml = 1.0 #PML thickness +dsub = 1.5 + wvl_cen #substrate thickness +dpad = dsub #padding between cell and PML + +k_point = mp.Vector3(0,0,0) #used for PBC + +n_si = 3.45 #silicon refractive index +si = mp.Medium(index=n_si) +n_sio2 = 1.444 +sio2 = mp.Medium(index=n_sio2) +n_air = 1.0 #air refractive index +air = mp.Medium(index=n_air) + +P = 0.7 #unit cell pitch +h = 0.9 #pillar height + +def unitCell_paramsweep(D): + sim_x_size = dpml + dsub + dpad + dpml + sim_z_size = P + sim_y_size = P + + cell_size = mp.Vector3(sim_x_size, sim_y_size, sim_z_size) + + source_position = mp.Vector3(-wvl_cen,0,0) #source position in the middle of the substrate + source = mp.Source( + mp.GaussianSource(fcen, fwidth=df), + component = mp.Ez, + center = source_position, + size = mp.Vector3(0, sim_y_size, sim_z_size), + ) + + pml_layers = [mp.PML(dpml,direction=mp.X)] + + monitor_position = mp.Vector3(h + wvl_cen,0,0) + obs_point = monitor_position + stop_cond = mp.stop_when_fields_decayed(50, mp.Ez, obs_point, 1e-9) + #stop_cond = 10 + sim = mp.Simulation( + resolution = resolution, + cell_size = cell_size, + boundary_layers = pml_layers, + k_point = k_point, + default_material = sio2, + sources = [source], + ) + + flux_mon = sim.add_flux( + [fcen], + mp.FluxRegion(center=monitor_position, size=mp.Vector3(0, sim_y_size, sim_z_size)) + ) + + sim.run(until_after_sources=stop_cond) + input_flux = mp.get_fluxes(flux_mon) + + sim.reset_meep() + + pillar_geo = mp.Cylinder( + material = si, + radius = D/2, + height = h, + center = mp.Vector3(h/2,0,0), + axis=mp.Vector3(1,0,0), + ) + + substrate_geo = mp.Block( + material = sio2, + size = mp.Vector3(dsub+dpml, mp.inf, mp.inf), + center = mp.Vector3(-(dsub+dpml)/2,0,0), + ) + + geometries = [substrate_geo, pillar_geo] + + sim = mp.Simulation( + resolution = resolution, + cell_size = cell_size, + boundary_layers = pml_layers, + geometry = geometries, + k_point = k_point, + sources = [source], + ) + + mode_mon = sim.add_flux( + [fcen], + mp.FluxRegion(center=monitor_position, size=mp.Vector3(0, sim_y_size, sim_z_size)) + ) + + sim.run(until_after_sources=stop_cond) + + freqs = mp.get_eigenmode_freqs(mode_mon) + res = sim.get_eigenmode_coefficients( + mode_mon, [1], eig_parity=mp.ODD_Z + mp.EVEN_Y + ) + coeffs = res.alpha + + mode_wvl = [1/freqs[nf] for nf in range(nfreq)] + mode_tran = [abs(coeffs[0,nf,0])**2 / input_flux[nf] for nf in range(nfreq)] + mode_phase = [np.angle(coeffs[0,nf,0]) for nf in range(nfreq)] + + return mode_wvl, mode_tran, mode_phase + +if rank == 0: + diameters = np.linspace(0.1, 0.6, 24) +else: + diameters = None +diameters = comm.bcast(diameters, root=0) + +local_diameters = np.array_split(diameters, size)[rank] +local_indices = np.array_split(np.arange(len(diameters)), size)[rank] + +local_mode_tran_arr = np.zeros(len(diameters)) +local_mode_phase_arr = np.zeros(len(diameters)) + +for local_i, (global_i, D) in enumerate(zip(local_indices, local_diameters)): + print(f"Rank {rank} is processing diameter {D:.3f} (index {global_i})") + mode_wvl, mode_tran, mode_phase = unitCell_paramsweep(D) + local_mode_tran_arr[global_i] = mode_tran[0] + local_mode_phase_arr[global_i] = mode_phase[0] + +comm.Barrier() +mode_tran_arr = comm.gather(local_mode_tran_arr, root=0) +mode_phase_arr = comm.gather(local_mode_phase_arr, root=0) +final_mode_tran_arr = np.zeros(len(diameters)) +final_mode_phase_arr = np.zeros(len(diameters)) +if rank == 0: + for i in range(size): + final_mode_tran_arr += mode_tran_arr[i] + final_mode_phase_arr += mode_phase_arr[i] + + + np.save("unit_cell_center_wvl_diameters.npy", diameters) + np.save("unit_cell_center_wvl_mode_wvl.npy", mode_wvl) + np.save("unit_cell_center_wvl_mode_tran.npy", final_mode_tran_arr) + np.save("unit_cell_center_wvl_mode_phase.npy", final_mode_phase_arr) + + plt.figure() + plt.subplot(1,2,1) + plt.plot(diameters,final_mode_tran_arr) + plt.xlabel("Pillar Diameter (um)") + plt.ylabel("Transmission") + plt.subplot(1,2,2) + plt.plot(diameters,final_mode_phase_arr) + plt.xlabel("Pillar Diameter (um)") + plt.ylabel("Phase (rad)") + plt.subplots_adjust(wspace=0.5) + plt.savefig("unit_cell_fcen.pdf") + diff --git a/unit_cell/unit_cell_fcen.pdf b/unit_cell/unit_cell_fcen.pdf new file mode 100644 index 0000000000000000000000000000000000000000..574179c0799318948282669124a56aa868113c91 GIT binary patch literal 14577 zcmc(GcRZEv8$TJzUXlISBb@En*?aGqd5mlh5e=&(TgWWCfs9H%7kMT+cnOabMT?5rC`M~;94Uk41fZD&POQa z?>!CJ2Co4NL$}Vz7WMIIWrn z$Ohk zj`aXw-)GeD^akYxAR6AF64Wp*ey$iC+dYk?FpMiHh`qSFFf+sz|S5ji*AS>k#{%_^-wdX$t~oNTTaPCmT7{jTd; zc{$xZo)c^aGmp+)i@dkK(ah!R^2JMYC3Ne!e7O8psFXY#dk^sg5~hH&h{hMUXp>hG zYBTM+++XZfyhGxUVU%uhm8)goR^;{9N@k}0y&%2i+4RD#T1Ej5x8dmU4WUCk$*IxV zOK#kEPa~^=0QrURGm1f#$#5S=?1GUZ_Pbi?+S2arg<}yk-Js)bN6UYm8+_fmI6C z{UD!D?9NH|$C&My94aKHV_QtxX>Pq2^z4qS6I8EoWIz^OFr9K6Bwk|`WUW$n?J6qr zWj)|6xF!`UgqD{fJEtQZ)@L;YsgP-h?0twt$w$pKEP2fY z6ek?dXBUYWnYmlX*3wy7I2tO@PCTA*ZsyB38C0U58 zn)+(#i!mmq(36IdcFzixpbMvtWV_k5*45q@h~`d~)4s_hAz3b8doH+0i9vRCBZA#< zLEkST68WKl@76<&J3jf(iw#H()BUa`k*=J1YkRr)#jPYm`e#v#+qdSmu3S3sI#zbx zXJLbRg-GnH!)5;KNbC8X4{vkt7w4|2qTL~pm6De@?{|c_i-p}}3OfI`q^J1;%}t+~ zjTH+Z;^^@wtJ|M7=4L)e9BJUHLn|3uYkMWtkjPHZ5=pQ_RbV43>wDIcUbhn8uPOSB zu~%Li?;b4CP|}X6?H4V6>lAt=B!r8bNr3<-%d|tIMVF- zqIiME#(xgLc9YQaUDjSVlz5BpoQG@AI|+(fAu^>#>KnH*nB05P8q^0vrtF!k>mTsl5GJqsX; zI8oLisq}PS#Xh2&ZI==yrh($<9ib^orEAyYyieIbjJh@=M9e4|M`eYnnPTxW0dkg8^kz*99Mk?Xvb829T&lsA}(<8Bf~ zE%G_?wHCc@-G7_tRNSMqE89na7Y0WCZ9Dd|^g7j_KX*MQKckAh`@DUhyk{N%e4%{k z`t;rdMCVf{>PU4)HV#c)trdOZw)d(J|nSEvpSBr!jbMBwx=Q!lBu0W@m{`L zJVT_gQX1MDdLv_RTj8(i883|=BripNNiJqTa`@n@JsNYZ3^_L!RNJ6U_M&`ydjyHr z$PJu_T-EN?*^MdZht_O})A*mWdpR>;h3UNw@kRJrQl0TQcF3KP&D)GO$r(u zRQR6fzJmsQQp6=n7#fBCfh6$R6i2iL?EjvvIuCumnHj#6LsMA}yc{5u9<|lle?-B+Nx%v_`8Gc3PTh-0d9T; zkKb1t0xJD?b`*?teO+MG6PwJA6b3z4J4 z@U1!qw)gfMY%T7HqpzM;FTBqQrM2Ccd}`%^JL;*}^SHUw9a40Q2L|568lG9}Xj4q% zm)EpPznplS;f%ewwc&^a1fF-?vN%f?S$u8uWcl?K1!kac&c?g8Y)3IvrIQ-V>TZ0N zi|WJv*K`&GQ^1Azoa7M4H9oiA6MZU2+dmw3+%TwT*AWxQ8443d*t`!vW)%7Ufo_N% zw<2SUOWYkAjW=hWt2|ZyB7}COlOIFW-;j=1%?(EfR(-adVp)AiH$7)q8I_<56BGRC zV~}L*+l4F}eRx&uYL;cT!u)*FfMr*uw4o?>aO?+@{-JcgxV{UPja^*s`6CudgC+%b z4@V&+J6Li8%EtlaT}uIuK>uC&D5Ks~5?=aErdUS*{rWu6l>{#@o`tQ}gORg>kIs~o z$v;Xxn&^1x7mV)g9i=a{v!g&H#MFgM`4%Kuwl5XwJN88% z<(SE}K3q*)cUEV_Jr&Ybiz0K8eIPKP9pX}x!ifnF>QIHdV2K74x;uqW!VstG;{K1C^J}=gCR(SK?a8T>0A* zH$`Aazp@E$)!r#yy{w)p@|6SS@Vgy3;8%x8f8HAz^R^2)G(#>LrNp}cQm8b&-Cq>z1n34CdUxQ zYk3@Nl&mezdFaPQ#(tes+6CVA35Ek8!A}I(S{vxKs0Q76=aQ@>ab(Y zH%Y}Asb_C6=}V7}(tnASIWObP z8nf?bccXhTzLh!U%p6rT)3aSB^Lwqu{K`hME;NU{d^MomeYk@$GqHZ~W1 zsEr)~CF1s^yOfBK`4=Uc=nrtgC?*C)Z!49)&wVX_{mZcaa%&}xU*myxf$jr&AJPPG zzR4+~b98qv6S=4CO zQCyv;)sbd+PvI>AXOiGYzFC}$s+Fp;NvX_VgEP$b(0SW0p&Y$3jnc;of)8ov#exQp zjtgo<&0VvG-?U~c*LRm+nSe94BvjoV(c$+yNcUWU>S&}V%uv!Q_IQbtGl$e$6-`+$ zPp0sbeayBx-sF;R0zLbUH;8xI&mD+9;Q8#G9-1XeheKLOuq0jXbk9?~ZGKvc% zjOuSrj?&r`-{TeXpg~kvLp1z`J#%@Z9*~|t5Um zyB^x5!hKIsZ(Q>3aKWeQ68Qo9vS@*X2Q&Mc)2yx#O9K!l18SwiGi;8Vqz9)2i+}+ zh#tXf)gLD}KQP^Wad#+P&tTi7AXl$z(zdYi?7m=LmUp`RQ-Tlg#4U>OG(#R7m|^~C z_jnl@=uk-Pa?%zaX4JnKqZj2i9uQHYVs&($MLgb^Gsit;{*$dWmkA*EUWYqQj*5Py z77gTfX=FPdZcaG^?zJ;Ek)hyR$dkqTRWo4~w&#e~O^k9J;`YvKd_q?%8n4U?R@ zwzd#%b8;+-woXa&Vy&xjNI2C|+5PW&XlaJB^Pk)@nN6}lzo#2qs2;u9wI59U@kYGAkvrCvahX-r;o=bjA|&D)+h^+85bK#6pp$s#ws z_*6y>v~jYu@^aH;Q}(2HDqUk@ZepyX##=S3fP0TrHs(ti%OA~k6p!2GT8kmb`#7ia z_x6I*UmS8HBv2*`&#YTKdwiuEL?YXxJ!|wVcRFfHIQXm)B24JI?(2E ze-VPwtTf3)dwt1P?(+D%O3tOnFAnj~%We}8mH;G1|3j=o!0U~mM%74S#?vfilrwJM zA=8z!4JwjZg=8@qByMa(p1RUx41+{4<)=WrvUr>`f(a2PD^0?e)t*64B&5o6p8!7u zKr8YeFe_nx6yq`;U_6MO-kh=qo-!IJcOq^N5j7Ohi2OM*f&g}eOZ|f;)Qx4*NO`Fr zgvDN*2;`_n|r#~xPuOP$R6@8JoKu;JI6Go z|A8>`b0g#BgofIh+9lJ-GqY;PL^6#!9>}~u1L!LBG%;5PhMasa;BlJtNg?CyKDo%H zC;O-$iype`YOu%tk=_GHx>@Re+R+r6*nC%QVUI=9uKQspd7HPceTCl3T{qG=gf-E> z>fQVQ-!pE$iIXw#0%Kmq6Yz=(PKrkmVir>noddYmRQO38<3* z_=Ww;R!be+YTX5>)xr&GG#_RShUN8^_`iGU7JDy%c=I8p_Cp20szbtsx9L@B?}0HITQg4z1WgMC&8*D z4y!fqNS(J#wCX7W{+A&9TUw}v_#KOs+b(2M8V7i3H~aRGv)6Si(TB4qIx=fh3)?P4CwU0+3Jrk zi}UXmDZOq!eMdH3Ih#V(X>dAG9IHDLumk4nmajy{vEhdZSc4FX1^0%3-{u=>`+_L8 zY4G4iBvBd5#A?C5_`_dzywurxwQJP^R!W=zCWu=`VqRx%R`S`g$kqo(E)RGEZc+8T z-ct?IBq!*gd!HS(w5>{PG)`B6x<99XU6r`fEhNrbrYvS&@Bs0}{@KmMG-u+A2c6N% z!zMBIwkr8*vG0$s_tywAvh6QVI9Vh=>UZzHj9BKVn!VbaPikf<(eqj1+uc!vrpaZp z^HuaE>-rO{`}PI2gfg2S^p!EWCqs5=T*m6!zWJBAkk&nwR<>;PU#x~2)YOF^Upx@M zwo13~!jzDM5F&JDq>jsSL?*Au4P`6|5xMcE&gWIEd(zn52-uGR-h%y0 z5>PN~ry+q+Y`RlI{7tB1Qo>BkkK3*9ut=IXSdEQYS!!xh+nGhCU2vJ|DqE!6l3LkE z7k#-s=%5*O{#c@}qpu@9`+}Yg?5tW)I%QdrFWj>pKBX3R*rmv%;I;GXHo5Dlabg=W zIgu#ys)Iw~^t7-MX(uvgKI_$D?yDbTWSR}n=byR6cy;J&R?|fJos&*(Py*Wn1OXTX zL&5&8d6Z!TC>;)7uu@zayoF4T2D&giQfZJYS{IT%bjkS-(7|F?+St#Q8!4bRQxc>X z`QF5`Q3l@?uR=ikshrXB(;6dPT+`O z6{B#8qr~tt{WN2YdEdU9@}U>y#2J^~2Vz~v=JkU0oUpUh_m)rQzF7OR=cPDwlz{aJ zz^8u!%Lu6(AfkXGp(JSpT?va?~7jES0@`v^L!KfQC9RSWxYj7lr~ju?D|Z5 zi0SzAY3`5F@70G^53QQP^A|e=%r!;U)d2B9w_2kMpG?O@r-6y8<(NdecO`W;qWmt? z+M42P6VYc!ms#Jfun*E{W}qhsNJ;=HA;87;*C>ZdYw#h3QL6`;x06%HWhEW!hH{Q= zm!EiD6-|sfzX1`|$LO4n@X+zjvJ<-~%V*Vi-|sqK=2pT({#47YtDUm14j#}BJH44H zEu<(DVB-CCQ}>hQ#FKX_Is_kEst$aJ`RU_>(_k?*^WVqy1Jg|3(nE~pQ}V1%PrFLkxueAD$ZY@k8V0E+^JP4 z=DnO$Fh5ziN1PHRJP`E8qO3CEs>Gsqq?sqTdxz_eZHGEsXb%h0tX}cR z6OnweQu?LJ7<%Vwo(@c}AR{c?-;Tencf?z{MtkbUxi7jUo|w$0E;jVBI{t|60tw#K zwqWzQkP+1jwY~2LZe<<&GIM`@eqx|QE!^C!S?*n)L>^}NR?OqhkULe>9Vc6x-6N8( zAJW4_#QkX)BmpFe{)g=YXgpI8QbYGJmX)3UAE%JweUak)R4JL9lB!fLP#@ za#w2r1~lTFB9n=6lK@8qz!(DlFNa|!GQQwaS(x3w1)fCS*e;-95l!2v=k@UxOcXaR zk>*N_e6dlh(%;$FRO`b4?tOieHD{Pq?kl z)qIl1M2^*KZ@K2fh4B6AzJrR{0m+Z_?^~Ow@Mcn_e<~I2tT}Q&F;U;F3O)2zn(wxI zlDaARNxsxorZDkOV!D>G)4j`=UHA^xo{<+eT&>ZxxU)m-A{3E8z#atPBJ3XyhVdct z)LmxPDkKrt42cPWkqJN_7~&rmpMLTKTskrj+g~?0^ z)oq{0fSj$%1zwc*kTFAadJKn6tfn?r*E|39An#MUXS0o5J>?_Ichp($w3Y|M=m>~I z0Hh$0{~(UU_l3Jnc^o{CoLC+Qv}!(||Xhzgh3p1f1GKKtf)1v|eIyX}ywf8QKXhpU}^H^SmS1XONIRy&6&%Xsz#3x zwq`y3HIa^sX+F8j?gNQzkGe)=u;ZKB@yn`|b)1$J5mv+4x~Iu|Sw&u8RxetbY8l-= zeN3m2cCYyKfm5d+rtD8WJ@oV;XAm{KelS1frGfTmBkIxXHLgX`w{1C?`Rg9Il%$zX zOE^k?dF&QA35ev7OR$GdUlH(%S|`3j;Y+nVi)tBJ1`muyg-&%6Fc1N#_D@G+g2qh& z;L+HE0+|N=&anKz8xw)edg9NNY}F#fKt`st$h2^Uk?%&DW>pc*NyBQr4^6M>*I#$| zp5u%Y6Cg98eG>DuYxC{?v49H<3JvU;?eQNYFIrQgi=^etUfI0AlEWK%)Hq!1Zkl;R zoUS0b-c93Tw_}fHfoT;v*6dR`O?QLXXy;~n8~C?$U5<#B2wWAN7f%uBU3+qQE|DFz z{qar8z~!syD2+RIbF=64RIWQ!af&D$PlNSul4hiB?^T!{+3Q1{Rr7fASu-Oq@w#^F z_T%v@etD&cx~^v{`EC`-H6+stmHKU>qVAWwMzun=!n)T#hqdewAI`UpCZI+_m=OIB z`$zp0U0-$ZhsoKB4vyQ519Vn{l~-vGuOF(MIzF6AM@uUFCNhcb_JtmQR|W@fgp>G| zZ_j~HWYUhS9mHARsNREWr1W_Jx~}s;wVqjAl)WO8`lG$dGBr|w?Re)m1 zbw{cs#`0`RSg)t3F8p;dud(iZ{g#g~g|xX3>x&f{;SSnOwOQJE2OW|}Qrk3mk3yws zm>^2I9MH7{jYRfFe`V7wgTl5Cj2p5M!Qpe#y>ytW5$2g~I>|q!a>xp=FR@VHyH~ zt8K$gq3?tZ^u*V`NkXAG=mGx!cV$lr!Et8xcWLcYZKpMej(rtpiw~f zg26HbutE=F>FJ8~042vIY7=lkKXYGCeEA;iTW#adT1dQkew;)9wEr&##B2M1SPL|T zx08DisIPyrn{to~uAor{4xY8*X29J@z;W9NRVRO~Zz_ROLacsVg5{E6A*ZK{vah>0 z1{4KiipBU?;4bu?Lh&+!nnM0sO8EnsfF=a7jufcC|Mfx0Kmpum_w)byfow@3q$t4; z?)!g#{x5U&rTq2&@7IXmc>~Ktf9dmYV}T3czdi^QAO(8+&wC^QN8@#b9}9Pl zQ>GLMb_MzzC<*@hr*DJ$!4HI|d~DxEM7Om$9x$SPPOV_;MbP48N%>4o&D}k7*DhkuQLROzxM%0)Q5|fMX#G!Y7nC6| z$|i2PtdjH3s#c?1TsB3oh-tml%42Qqm}%E@Z`LIOj<1k!J3^{IOW#zvX*t^xpw`qT zdBqDl_Bg;WKm46W3{t?4L-k8_g{{tJQ1f+KUaFPDHbd#{m^33TN5i{Rrmjvq9y(cYh9wrKLGH0s)9Z8J>0hq;Mzzv$d46^q-h44aZACm7~g zb21dgY?%2^?NQj&X*tX`DEX$U^x4hxI)au(Y+GNwx%al_PgHz-63!kJ>gd!WKi>Xo ziKekPfaT<&q?(3Te`Z(ED1TTA^$%O(I_Rx$$5iG}INA(~mNJ1#q0F^lui@ay6(Xp? zF-pcUGqr&sg)B|I8mMKn^CVK8rXx_bhWhoG;AeOTeo60|R%+(j+W#ZHgQ51n><&ie zUAY}RF2yfi!nysAf&VE-{J-UP@JJKS!H+~vAh(05=C|Aqisa+z8yxhl2LDH@1XDPE zpr5G{tVqF8ex^#iv_IVLr!EPn%0Gnv|DP(Q!Kn3%G&q5A(%1sH2mgh2yheXW`At#y zRs5?a@T&RI=6e8wgJJrwh}u2{<_8x2#U4G;mDo$qfZ#6?vzITad=+w4)=^r#gq^-S zdgnM@TBC{KtCHzEbP7ohUuG3I_H0RLA{>clqx)}nqBvq0gn6ul+m4xXn=ZJGoDGON z<0C4)-mWer_@3jeZZ8LCfj8$$57DVsw{zoX_R;cWO`Skhh#sF#C~V6XC1c@IIQnEx zeUPHLz5GSmy9XhBi-RjIf~1dT-f{P;ozop6ucKvxhn1NCTp7*>77uMJC;HZ>H{3E5 zRTCbz#GH5%Kic|kye|jQb%e?@>t5$KD;D=I`VZW~A|6W9as49bvATA&_#F-MJ4F8|az|-jO&3-ZZFJS-YO86y@ z|0MpBVtx{F9`S?2Zw~P71sjO2rz;LNo= zP?XUB{s1Wy*cfmF{-fcx61YFW_d5-SfPuJa7Y#0h2F<>!9h6XC+;Z5}4u%6AyJ=v~ z_?-scfxmO_rhz4Kk^q+m<;3i^MUpP$ak9nbJG#KD6YO;6HP&gpEizfZ&97x2U zeZjrwAN(Pup?~m)f+9ih_O=v{qfaQOCjnl$uJ8%ZMx&`{_7pkv4KvW-k3oAtBIzb$1q@g0U#!Rez?$xha1|yZhinhG~zQs5Y{OWiw{LG1VRQ) MDJZCJs6qLE01Q(A*Z=?k literal 0 HcmV?d00001 diff --git a/unit_cell/unit_cell_fcen_h900nm_res50_wvl_1530.pdf b/unit_cell/unit_cell_fcen_h900nm_res50_wvl_1530.pdf new file mode 100644 index 0000000000000000000000000000000000000000..e440b2cb9ffb9a6e5f85337e73d48f8b75684d1b GIT binary patch literal 14499 zcmc(GcRZEv8$TK8*b-6c5EiBYKfaa0X1}om+bsZ|Z#?(z7=XL)>_c=g9xk5l)_^>C{(h%6u&(Axuy{?48>7V84s!SGFgR|fv)pW0Ky zVC``BVAzit6&#!ZdBG3`CqM~u7;6_B44!QdoIA$S8RVV0kgD#%EXok(k!=oB+`RN+ zmcGc7hBGHpm&9O-!ToxZf0FtWv3(a$xj1^+*$uriyRN2@^F;WTqE2JNJ(E~5KiqQf zrU7d&$Em{?eWfy;P0F!hCoeShj>UJ-Z+t#Ey$rkXiY{d6+4azb$@0!M^Bc(^W90T4 z_TVcTyHi$=Xun8ul(YzD^RY(R-N>6g)H*eEIvp+{c#BlmL-J_(v>&xlA4Tsq>rWLl z(M$()t4lM|loLG_TV)03R$QZBBc;GylQgVXllaA}xVwb~o^==)w6?BH%%_~RW$%5x zdgodr$AzPs@0fe)+)l(?;N?$mZ#Slqm>5`&UEX^9Vde7&+&iwlYV_(?vub9i?{0sD zH)v}^)NjlLyc>{MMJadR*3NjJF-rGtB_TZhj@8o@snw*qhO}F=$|~})=LZ$VKCj{4 zW}Yj(mCd`P3G42fn~ML`Ibiu}i9Zj!n-kQ8%_ezM?z7S`ihJ?X{tQ;-Sx3+A!u!w#7)& zj^+`2s+8ilgT3!e=o^I_pD!3LfDlGccp0=wSS$xVj{-X3y1t${R~vh#g#r z49jlGxfi`HuDt&)8?Den$DD6$60;C3iyK;G;nn5LYZ2O$!+yHDv!s?<049=A*;IhN z@$&4kb4bO=H*eOA&umwg^DIV-i__~p*3zqk0= z^`_i?7Ac%}$fu5gIyR2r8t##0J)NV%YDx}Sf9irSZzj5}u^2u*Ur)yH?(tld?SaZ_ zXLF4OhX#i>18BhNfgR%=Vq$@?3mANIcoJJ#yYWaMQIAD0zQB zskBuQl!~HzGNn_`P71sLWwFKb$sf_U*(Jp4r8n)@mriew zdZAY<=$PgBl6`4tR!M5n*ewSoU+TWDLpL)lwJ+*&v{k}an=z$Swpe8nDQ5C(`;oU> zneQ0{o!gUgYq~C^bU8faG{uY46kf8Om+iMt_mRJ0m>?C8T})d}C0FFq`b0v0z3IZf zD*AG1_<@@viO8F}DvM7&YBPnNwA{y>TO#??jSXe9MZGnACKkR}o=f6vYkVfIN3fIn z+pRd?)MIPjNH8;K$T92Iv$&&1;-5$%iKsb?(o{vp{jgIa)QnwOCEZQUjzN@x$GKH`9s5cy zDT_T}8b3+OEnyUF8@O@(Bi7}n(%^(I7q(WT=?i#n9|sRpj=Ke_ufJ5wtFZE^&~@<* z_K|2_>k(Qn&ZdJ}N3LIS<}7_lH$vG*79FhRb~fovF{?e#NmOugCH$Pgb^wT_xlXBY zHi}OmOosh@iS%vH;*KX3?=g3_o-<`X%POY8X;D!K&Ra|diPn_8dL2T=7cgeh!RB2! z8+6Jf#vbIStaSJKIZAYz@!Vipo6;Uo$s#PwtG;x8e*>7ANj7FPC1*{@_Ca)BDrR4m zL5;Kscbz_6O;<%^^fQiV+Tl(d&(;KasZ*-$nAhi1oHt=fDyArp2qs5@!3v|mRUNprHD3kr?Ld&-c zQZjN4Z9wmu?YDV{doy5X&AscymWZf@SjHDY|AVwzXr@*$!<TFKWEc9h4hyT(OmLgYkl`@LJ;H|@+a0bY(-zJPw!udEInzQdZR^)htEuM{h9~>S zL+UEe#idQhNO!8{9#CNImcS*I>IZudgC*L1t+vZna8?i2>3bWb40cNl$77xy7)KQP zcYzTIq#jhO4@Jpj8+IHt-9ZWB4>LX8P#2#sBg|?@^XLg_b(dP#B|&RlUG8MkI1xB1 zcz}WJI;ktA+(k+Qb`mXT=6WLB{0{eitTqHx?C9=4N_*;C%O# zmKueMkk!yfyr*>YZVN&yk8_;i^BRP2RWq`EFxy~jutRvgd0sL1A>EtKWMAy1(J=95bO7UJGeINGAt4AnGr%A>MI4!XYxrVe z`LZN4xa+;KQ&qt?DQ}rJ8XT*gcEl0tj|bk;pXr$Zr$(m7p0`-#vVDA^OUA4D<5`Oh zjT&|p0q*pEA6|s<2R~n}zz?0O=hZl*m>R6ZAJ8gHU4AL^T>8s#Q5*V`BZ!*YVv({L zexmMWpA9BhR{H2C-)oizMXACBc-EXXVzjaC$b#X%Vu9imgEYz6*_a-K_EIrTeh$yj zk2>A`i7w$?sRsAjkJx1mor&qy$u{j9hLC;3Q4&!;9@_3$3UGwz-<6NjdK^#2$*{RU zl*#pg`cdd|l%wO(xvj<%fipZ$FXt7Ud>S8-!UNw^E&ghVi%dJP*H+-`ox*@rtp|ID z8FEZbsSr^i)#oR$*)ayq3pdm)x`Mq9O=lRLtst$AP#Lm|hqPCrD6GXhxqFn)TbD>4 z=`T)s%d9{- z``X1NKmrFSQf56Od%k?$nFfb>rx{m6e{fypmZSV_0N7*Ww@3DTxvP2Fuyf$Z{9NnJ)atrRc^r#Ir z2r}veY1u8(4p+QTH)VEgM)!=5AhcIA53Pb&8>})#2X1J)v`r{8Vj>mG16Wox48@;X z-lgx)UQSnr9W+POZ98VCjOwbe0+QSLa~!jv3c)4R5`abt6Z)`{QPl!y}XyT~0%M4Sbodwt&&8u-&(nOuC#&zY4=45J z+C^kB$k(-J)h38Ox~&(kHD;FE=_;EXO%}(p(MK_{;?0vtLRn+EP90lvJe=!+XC_Nl z4nHGvb9%fy6Nra5v@zE|`&2By){KRq{Lw}W+QB_Jce$;|JfC7y4$sS$%8JLtGk^6= z(%VDtWVV2^a7@-p9Le@PrJxoH7(6)KqY?G~mJ$4p5nG|U-O1%KI8#GZS?!R@F_#nc zFD0qH0v%wQ!iJ#%d6rg(MBd3LiaR>&_q*7|Y@*^sDLm!w(5<~e`mOoO!Qg`qFBa8A zS%Op!i5-_&OPM?Mxn}GswrnmbH(N-n=Fa#qopJ6%&g1s9h!P`+g3F*o24_$w`DJOj zkw?~V$&5YRA26Es49eM8%h5SL48LIdh;AWg-*eQIb!ICEYt_)E#%4|x0wKT|^t>SK z+~fY)<}W?1qZhuCMwkM7&t3Qejz0(3VNf^<`FF2S*OtJxLutl5KxfBEhlW4;=rA8S zyVRJhV|-oqv_i7Hj)tn>d~^zV3;kUo|3^HxD%Qp~KkjdR)!Lt^rm=0EouSq~Zjy5^ zVxK1`%cAPB37)mi*7r$1W?~GNg!-23M(Mtu{TL8ncI{M)ylFNOHR8`tcBm08@()^!mZob0Q;{##TU`TR zKmBb0Df#YOxxS+{1ZpW#QpCxrc$DX0**@Mx-=?Y?jGd(XY@V+~XieT97uvIwrZDQP zleKzWHScI;XH`c3m-Piq46T*ai&?g;VY0JB$B;z$#2*Xo;1eeD59XHFs$*=0(#)YK z6vl7aam={$XKs&O`iK_hK25$)XP$#Wa3ZM^dT%_x^ydBX`)T7&@$~niGonK+6yC`h zx;=aE*}hD6cP=NZcG{hkaN;>`h;E^Kjf|H;;atc{D#8}AhN(o6v-9qIIHJT z))P_8U*Z{`ye(T5ZHe~hZx7*9c`RcoGpr@ zm{^TUzq6-fmBZtdz4$4ooKtV=oYHmUx;uH9UutPDMb%YRRxRiTUY?QjAsahu_Q7`@g-EbUjXqJsK zul+oE`xBpwob_9`zCs^ntZONp!s)0NJ39vVKl$tYu`DLT8s z+MRO*I%zj`e(4)~`SK?Hs)g-aBI+apeqsNz)shFcTCHFjIX{g`#lDnYpUlU3u8Xg2 zLm#@4ZuWtc`y|1JEkf3u_uu4~_L+Tn$t`B`HWxctP|46ml|^;p{9&Px!n4Aj>?h28 zkr|GO`Ecu|F2+~y->HOf;{@u*@~fV3>d&mDM9rEnUk>u>v9+Q)1J6)0n@7amZ|IIu z>#_2;XZ7wW%^bq%;uM0)26e>WH_Qzk9dueRCwLK0fy5pn(h@6&M|I)RH!rLe#3Ing~XgEP=oCL zHKsZFvU8Pk6znXtgXvUmbB}R(IhH(SF^3FSodQ+f?WdEE<959cIku&Qnv48)Ca%_+ zLPX&pC*5Y(9!mD=4r%(R*9Hnoizvq{FZVxq)X@di5J!%4joW${kxHAmo-y=pa#w2u zpVp8!f|an^Y1|AMpC3DNwh=#P7$#NF%%2z4gYElNiAh8fMW0Wp4bQ^DPS2}$=(-v_RBX%dA7&ORdR02 zd6wY)5Zk2a%(je_*ocw9#?Es$dz`?wK{cFC6Ln%_7wDmTUw9dqlttguPLzS#y<~V> z7QNhYT#&UuT0lR$6Y<6D#hvH`E7EHxtVE>;bwbQcWU}N!KLo6ISMo5i9Vm>tc;n=- z%flyVf#hL1Gr6fxa(Z#WvnhVt9YMXiu?6C@Wej=i>SL_?_PI59GwYwgqIDjkDXx#A z4R7t6eVqYm+*4|3!p88$u)j`Do_FBd!N}DW`ngxS#2l0eN)U=MANXIGw@XK!=WbTd}LI?!X4NeZXclM1Nd7ICY=Kr>D28;#G{J z=hTy?UmxX@7`@;2vW#_40=q2{`w_uguzyJclA6u5WH72tI}pTGhbAP>M@K)vbp0EP zu#UOm$cUkVq9To{USLA1^+bEYJpGo)@;>_Do6R04^k}k1qE#)h77XljYR0e#Igdn8 z!3`|jp$0x7=X2KjhEDcdtG7)O<)~3oV*v@iApNov{eld1upu!^3M($7m0XVEwGec@ z=GCmr*O`j@zoy(DD|~R#(iTc&dw?JUgJ3Aw-!+fYtOKOO!wXib>%Di8vBBVW%r}t) zNRd)%ti3Abv4ixm(B&rfh(aw%)Mi|i*gV%%C>zMLRj>>J?W$1Lok%>fRADF@nbNOj z)o+y)Q@6bU6uPzXDn9A6|8$(fQ9 zh>ml(5!XqWT1|5asYItBYXb+;_(U;nDv{ZV;2D_6KWJB8OM*_8lV%JZ8hg4S zZl%6>_Ji-YJ9Nd*DNHKeXoQ0z(odQ)#|K@enHe1Tz#)}!Ad1DD(j+6bHVl)Xy)Sra zU$wY5?aNK*nmB(ksOC&wkTP{u==yZ?dEL>MlN@WoALRR2POa#{v*ug4^%eQn<-mfy zwpCiGpL9p~C&6Q7OCizpi+R(l@$e7$ATkZdsB7hZT<7EFq|x&h%=$Ebw)+p9RGm?I! zDo@%*Ve>n`(L1TfX zGBDK*e53!%KzplfBDP|zl$366P+QBjev6p(+AM)5s17IR7hTSg1>&q z2ya69%g7pAUJ8xD7!PJe9!&4)vP)|+^M1@#QwlRq?y;)F(jKqOeS`em(QNAx23M+_ z_>b@f3=|FPwl_Ge3W7x<*Wpj1YUPjAU!>KMV0GMEsMt5>cR(K7E0yLJ`&9jjk&X;! zGIipoeEzn|b5Ek9)%D6m``?Li)!N0#>r!6iieK69Blt-`)gW~8@zPCet`k+4Px5N6 zR4Sf%@Qu{^xPKH8dk}$(uzxrhCWOdKt<17zNRlJdWI9AfCIWq6h<{jo>am^p)G4>5 z2jls$9P4{n-p}<*BW3+_G1h5PD-sV1Src-{n4YkdI2vj_f01YXm5-ipzJTQ5mP|l=Zfk z^Ev&CnR`bb6%H*ukY|0+Sm+6(Cn63Jkb*$|gE&Gz7Vb7E0632vTN(u`RervBNT<5} z^x-V#sgDUlW*_tU(}a!)zrJUcvf=sJBwqM}=^N}aTv@)Yrv41qvqfp6h^c@Q_G72n zP5SjUlny5eZ+lU1Ms)rI%$$sLOo#X?K@V-t!r&b`2kx%Vf&urjHP|0^zwDZn1D|QuYPbH@Ze@PRtt%&t?}lQ zqCrct+;=*(%EA1O8n+gFID;u(jsd}I#2WTBn*PZC+1Pf|k7VMHCpU3WH<$Hx$xn`iNq*WW{tCB;vS8(#+rLexZ1My+~&-)I0(7@;ohtv{S0EO;3EYmC5!P* z@|I{}Hxd-fZqQ!TtWf)S|1HD%+ZODV!{Gwl6gqUzLY}v8zB@4DmdYqu$DZ6AxfXcM z2qbz#>}0_k;}1pYoZepAegdrt`gP%|Jd|p8v~z8JpU!|MWh7YBE~VdZ^-(SO)Ks>i8f<$(v2GbJsDYMFZI4SP$A>CUd7DxC-z zjHjm~=bZ|SVXIAjLO=#7l|W@8GmQrD`nryk0G=_Xp;cCbQCuj!0) zUB*N7CN3%zEPej-1T_s;|1-gu^jI@S$wtZdQiJ*TKb?_e5{sYhr#|CKA@)-8xNsKg zjeLH5qdq^|WkWKRGPAb3qx*PX$np<39`M?9C~Gn@bMzUIU0gp}fB57q5xI$AML6&{ z_wN`+tIk;-NKCJpjmOE7Jx5o#srqfciR6VCM8x?#cHmcqzs=>;R(+!0um+P%c>i&I zzC^*#T={-gigM-&^Vp&ICIwC}s2J^j$mxtj(A6k~X!d)q(z+=cIZYp#HpKls{oZp} zkLNo-LNm(I@C5VROfyDWKKb?Y6s`x{e#NPKeq^}7lrChRtMk9LrgK;v`xPJhz7k;&umuz!i| z!_paFTeoxU z@BL-I7orKlpl}#i2q}U9!$i?Ypp*=Oz*nQ;r%<=VxjPW5)r6r?JoEtm|GP2+LGavZ z0Y!8`ks6_-4g`_M0ySs&61!ieeLqPch@6A1Ee0r%!&k(a5Tx>OwX_C`^?)6{E535h z31f@AOmxAc*15OQ66JD7|#Bmd4sSVE|Dex;TvU8T^I1r8hxFKvT%yiywa?6TpN3 z))4^|_`f~~G!%^g?EL+|egIn{2oVtQ!GHhn&;MnP-^Ugeg@S&+|NAxKPu_sS&fog{ z+gQK?_^%HF1r`Ck{ntGb3>PKnh%grZ8m~+d0LX~|(eYRf(L zCvspw-GsHK1%~1GjQCXtCX`A5SC(L400PYh49CM#fQA9X39A$I!xi8^;sn7}2oJ-6 zMg;GH6Z{tptbbrPLFnrMhX0t&2@I?yfCd5_f^fJFbix9*|K5(^T|k5)baDX$!Q=ZC z;B$Z^gzI0v4d{n35P|aD#R0Q_Z}i>GfhYTfTOjKFm>yX9KZO16)`VRdzzp8gkw2LD z5l4P-0mO|TTo3{T7{Hl?G4U}2N#FqQ#XxK!@InYG$Xz~*C;ar;pEVz-e>he%(@}Q` zhSEG(IpIG~cR&bPp+e;_D}4A9;mV-G-E}$YpQ6@|YjbE1a&Z zY~Azz)-~J}>eXNr$ISa#R0|s1H*%dQ}o|z%u2w1*GD$*d3WW1~UzAa`&5riyWaN18$mGev1krNN$%yvR^7nOjI^K z>dWaksh7_h_a`=wda?)DK~nW3j8S@$rQ*_#U*F@h;?6JTp0XBoxpq}$Lw-E zLAS=M>S6ManD2uQ@V~z#bEejYaZ@<&0^=-edXkiYG4ruYdn7ki8qTuy3Qv{gzqoT% zg~#9q+tybnj=hapV*bs&+Qq^Bjq z!0*Qdzrs8CxAd-PD5tNi{6Eq=5NiL+?m%STk=udeQo`aTp4)c~{7*UJ|1GxzN16l< zekO7vxgAI~f8=&RBxeV#r^olI`Jbs0Na2Kmex*vFA_Y(Rl`09+{&c%vx+I<|cM1Rh zKUInWQR_Ep@B-tdF#+R$_%EyzH2PD@?}{R<;@>qvP|eRaKLQXO2-AN>)Mjxo7og}b z)aycf^kW7_#IZa9vy*wHZ_XEsTZjqfu`_f8e+!^bxTmA}CU5cqy=08}ml>&zJzGMG z2n*7g;O^Qs)S(bY-lK-RO}@Gux^uQe5pF@3o%zMqo8^!5d^i-L`uNb{Y^TG^k4P>x z+FlvGypQf^%ESd!34g$3R8CVGKLyJXNv~(`<$I~>n+snhEOwsfn(tk1;30oHy~y!c z?uu$ZWi{P?xKDu&_(+n~!TD3$($UzO#Janh{Bpd52AB)aB8M9nN4wGy?dPZ+e9kc) zb$x9VS){YJMbfq4er2}<{vH4nU|<5=26B`F9E`%d=I>g?XM&$fm4^MklK*>k1R^U9 zG=f+I*br|@h&23vHUmzhe>D5e=)ZyeUn}9aJpPOLTZ;Kb#Cya~4!=9V_g7*dstz`I z*i3l$MCbeaBAx)YBQW`+eik^L|J7L*c;^GZIp70OysMK7&dI?Fd>$bTg9*a~!S*IE?0ZaLth_+a&CUjz(~_jb^b$lWy2-E$z}yZa(hs9pRaMe#uN&pD82pm=8o4TayV z?WDnX%Zox0=_>}^#j^-}*II*$Akn*YDk6fbB}A>_jSFm`?($=166eixFiZ9RQ4RWUJro-Fu=?H zIOOnQ0BRiMjR&AAI0DYy%Nc+g;GFS(00L})1r!t@cuz-SL+L*{sC#)503@ja)cmOZ zal8Wo!2GyX3nXY55pV=Nfcil&z!3;|Urzu{yoErGTyT!Ap3Z>O&tHAL9E|V;zzUpJ zRUKpnA4mY8S{@(^%0K_B{QR#4*!;l`e7gWdzKPuX;z90yun)!J{k;5r9YA@A{DbrB z;T>IZ%3gusj4<#IB?Xg|29PMUq!c(jR0R|a?Cl50{9sSX)6)yQgA<$nsSM)JzqF@| z_jD$>0PvqPs=K;_@&ZtGcTfqccn2>>Jdtfbf-fHD0SV0gVE)MWg%raFOWU^MCOSq+ zYAx}clFz!PdNUU!vn3Hh&T|Cbn(l>_Gb{YlqkVeVE7zUwR?>xaveNPKtP>OR6{!&qnF9$5ceL2J+)R zJ((nm`6##lR=c5bhI!N8%Kpm*30z{Bj4qH~%ML+ZPfMhD8My z{tQU?+pF?VKC)c=>{@$wlbw%c&^4arJR(i5sxvp+@fy!DYer^hboE87CcEUZ6|a^+ znRO2R_v~#W!`czAw3czM1@HRxS5fPHl6=gH#~6-AGx@fDYE`|q7~^@|DmwFHqPAyO zr~CqyZ8OzHEmXI948bXueFJ5HopU~KOL01T_j$``zwC#C^`A#>RZs(Snm-fb?(~4^`WgH(Szm6p} z_p!x03%|hB+`ZAVc|Ga2$E)7gOv1hTDbu)PyEdMuvv-R@cHl2bhwtvNLAMTbF^- zUH1syaQqbf(VO`&jBD&u##grT15cWEk3hO@bED5cikzBK-3aAd>BwdA%GJ|I-q{#> z|21Z52y^yj$K$mmg^6u#-o~l|s}` zh$VF7`}H%nT0Ri9k@o&2Sf-CkSTKwq#^kcIoZBk)!-_L3Ww7-bytL8}YvO*}i6e_& zUo~#`DV;#wQSSz`*o5n^+F+UH>#Rapoi~HVFUF)H=WgD6RuC?@PG_EA8p}f|rAgg2 zI-jd)uQlwh>b^&&wfU`Q&HX62>PC>3%ltk0V%`LvFU$cN@7BE+-fTETs{<1*%LQKmh{^5wS6v{yZE?y!kP zz79Sqf*7c83V!+~vbOe+ImN9hHtjn=`W_$t4$b0#FBBRYO7AA9OJY}Z2QbGA7S0Ec zgi)X3&X1fu7I@2Scquje&{wJpiC=e)exG6TV7@7^(fZOs&KZ|}GGbBMJv&`bZae`Z z5OPERO@R7z{N$sJz?PS%?sB5`%^$)#+uR$tQ*E zJ)4+$ImGSzb|tPP(xa<8W1epoEO?iQ+HYd`Ufo(8-ltgM{0!T6M|zn?O}hF;(D~eE zV49WlXvvf4m!B%GS1 z1Gsh`nLh)GOplw0mT(^8GomY8VPCE>({TjPFh3r{0H+Y zhsI>6_AOIj&ZI?ru)0&exq`L6 zs(SmS-kBrGz;#Kb6p{CK0bi4v%`IvV8N@i=pMFP0=Na&pQJ~y6eQtE`wbD4~e6t-t zoxj<4m&1{@hRk%aPcF?qJ49I@mhHbN&RthJ_ge9*z;f1k&y~YD=S~kjHLVD>-MNoa z?8>>*yk96O%o@1mDYBV{zEeshEeY#qxx@n%ZbRR*VwpPM){Z^hom~56#m{{2v%Za| z(wcY5H|V)2E#~(}1$31WrWZ>2EN#~i!&#|yI;RiY@0s2oWg648saa^x`Cc<^uGBU~ zW!&%i1+L%u7#?TN6~3Sw*S#@j{!)uaydsFr5kl=elv;XUC~)RU8r|6Htli~78nfRv8ZlS&#?31@5_9OE#|{}Pcbc)onR<+MO*DY9(8l0Br~Vb+!6A=&1FM) zs5b|G%r>x2dV!2u$l(GhnZJBC{q#e2mNA(TDn{>W-wlv6YqqNi1J^E!Z|gr9MRRgd?h9$GzD5$c|7^TyF{f& zI)!`GIlJWu={HOw0*3&(u3-C(@+E?O(+YEMvx3=ni4jC(v;_%aVE-*J5{WW~89YJL zvN=T`guJ^;3lWXBKGIy5Qmib=Y0hx}Aw~6L-N$Fd9gK_wGboa!5a@`07Vb+F-n6Qx zY0Y?c72BtK7}b$P=$eD`r1t2QpdhN;W6bo$E8MR|7M@wBCT zRmrA8pmHYzfz#PAb|3wxJ#UyzdnSOaxa_0=yX6B;4^BN+KGF8+nBA&=4Ue{%V0K@S zFw*j4NU%ZJ$45E=x_k<3%?>e5jOuUCzEB=j`gRcG$b5JhS#wJ|P9-M<<6HjKY=UFy z3G?Li(Hr6MI&d+e&mQ`ThMrxhlA$N1Vx^g87vyJV5_`s9oA(~Q{jkM& zH@E$>>n(@Eri7lJEh;(uG$l4u2tipb`*QFAm0?(~li2s$S59UbH0>N>xomAshm4P` z4w&%FPc&=$P@re`IO4?KcR3cvDk!RBwFjM3pj}mH8VA`&f<0OR4rTKD9+?jX%SMT2 zGKf&|=Y4<)-RG4v46ZDUkzo*qr#J6>^`XY1?8B-|*->+D6|tU0ljY8WJy9b|o|HVk zN(&}4v{ZSEG0ilN0&NLvBJdO6xrNuOnhKY$s-=j0--~v5E5+KETW$}zysJxA4kEUP zCUeFwuIJm=oq0;ubY!$fOry4H4UPD3t;rk6BvQc{W;YaNO{Qmj7-_hiyLTDF*=(PS=`S$!>YUJG#m8yfJjt=7Z!Y^3cZazze=%DN z&Si_N+i=Uz95K@7A(*UDi?cCI-P)}uJv7AfElTFRjQvy155T#(!;eK79ZGEm=8nG8 zA6}`_4TT0K*nWAVn?AJ6?)lZ(;9jKXjqEGzd&?Oor@F8+J1xY#N(M8t8UtRv?$Nq? zteq_*p{Doq!fL_V5~I7Qr8OBP5_8QhN<_;1jS{hXJ^SEvW4)raM~Xk@yg7XH+koCe z%MC`adagFXyIi@SQiX26%`9=;H!|_bYX6V`pU~<3W1|;>VkWC^cMYZY7COh~v8dH` z<<+KQ?%y(wF&MKceB`Z?kwBHqvHFB&VkuB4Z3k@)ZiPOn>|o4+CjVTHyvw4j>}}a8 zYHSc8;i#6v(0s)tQG+Qvakc#|c8mj*m+uJLQ~5vj%;cR@xuGJPn8Nx;y0bV|FWTcUsA}zOLPds}!4U=eQ!cT%XVDVmQLJ_evjB z{+v0h__b#2sb~3YdSSkVLCx*)AtuYhd;AAo7?EWbP<5{X*Bt)vN7QGPREO_7yrHu6 z^KD|a?io;Zsa9ojyBBlH`aaW#%ezL=Zyjcnog&t{AcKrSiKhiy42nRb{(D#G8OnKf!5GH;Ajifj28TWcVcGW`TWHC^T3%8) zqMo6K)z=Z9OUR^dXTBp5dSB>z#pm(0PrEx_cJ!s`>TfvY=je8gTV1Y?-Q~~EF|Q*q zA@rmvW=>>(BlHp1JND1k{R=2x+shOVr>zh{2Hk6sy5UYEKA}a*<|k%2#N!Nkvz?P? zzF1l8!vb<2wfSP?=vfAl)3X-1!wfi#$1kn zwstdH7Vqh%)_4|iG!fUezK(R7mg7*c!X?g#wJgP;5cDTx_srjCV(hz+_v{XKD$#sF z;>m*Z2-EjtpH7~%xv1E#W}Qz)jl@%yEozjJ`U@>4C^5AHbkqxt_7{OmM>hLW^6xiq znz%VYVK}L>YyA8g`-Qm5cL}Ejw^kLfKB5ri_J1kGXf=IMf^y-4`iKWMZ~353(f-Ou zRXKg%RzBbp8SNFG&v54rQ5_o;K#}2-czCphPq@@ym|M-DjR`5$avoIX)KgiS$ZX{dn<>s}17~7slOFnClaA5~A$X->aJY+_>L;fQsPg&30dG4DiPT)C>I@|E@m84OH^wzE{|FkLWTe30KDsbSY5P zZFmGY`f^NaSiF_N_*(XN7a-Y;ixaV0D<7=nu8z##;Qi45Qc++=c7u$tWFRr-FMJge zQDXozs6bJ$MRSxu-Z}k*OxI%TlqEAS(?q6IIdSiB)sd!Q?cD)a8U^B%#AED{yOA-n z(o_dB+tO)w98}?WNQNIWpcVBOn3X&~3Nh)A@Gca%w`VMXQ8s<0PGoK0j=FqCk?j*B z$zVr>)L&>q%}^$lnxElOP}Ifoy@5x)JoI2|n*@A6G7m_n$*5Q%*pJ3LOvJN)OJRNZ zrhFOGo)9Y96)B?qKpCezq^W|=LU888^-}J-J#MJ-R^*N0YbzG74H9)rt43$vro6k% z=cni*tLT1N@pYYhwo!8TBVqOz28IjqbybyBAB@7zPN@cqWEiwRlKFTR(2>92z+T}S zaQdU5OEmSf%WSoe<-$HZ+r`i?s@UPEPiga1_YpMBIAsshP%>jwo}-qq%N%vr!=Tgr zjqBIH!|vs*7^o`}uzIB)t^*gV^7e%&amZLzWUU-?-~Ue2B4?>JA-~zdmwy;G>HIce zVUy?P;u`a^ozojK>Ldex;eV5Ase!3h2f&~jqF<@;B(pat_d${O{41xZdp;CvPas-P zB0+?}Y~V;}0m8Q60p zcaUI2P!BI3z{*ZH&kpV%a9_De+C@Soa+GAGC4)ARh`&fI(O4bdE*Ra|6ub2$0fA2w>zHR}VJg>17i!I${P%~0xf=7O&|PJ8XJ6$;a*YPiTYt5$xzQk8~> zgK;36&S&-kq3AMBs6276`LcVM_WRvTYRQ7$m!Jacn&{cMP1EFB2O24LE`Fx9$CR`@ z)pwPc<6oJnYtEzH&b`>(bpP&Sn7%A({J^-Ap9O`ImA9#RV5_fgCvZex%>rJ=>8yV> za(r$)bH!wScCA3Erjh{mLKMkZ zgB*$lbHo2m^9{5-K@{82duTOmM+wK+QvR;EW8XI2)VLpLRjK+c7U6*1P^a{S+|HcL zq}buGmPg00_P7I1;Whm36Lr#5r`6{+Fk>s#3Hq6K3_&-@V!|A~({| z`J$YYGL^@PjQz;qE%@IA0r{hCj8t&CHD?Ia8_N)x9E3GFX}z+^A&Ip$A09S0)6ihB zHV#Y8a+v5UnPXm;THM7Pakb6wkTFBvaDtAVryUE=tga3%@Cq3caS8mE@vBE5OAJ zk6LWyiM?VVk6ufTm!3QDHi{eK-yvR(ggvg%GMY#`v`}G=iOcNMwRyrMdfK@vO43mM zgu>k;wlce=Httyq9v3QS6E3oo79L?Y#%=Q5*&wTbp3U4geiu5 z3#W5lE`OtZB@P=RV?8qP>2JU?QtB3nC}1cURjPXQk&KjU%5nUem&0U|T?>*Q(XgTn z?^9GWoIz z=_rj=r4NH-cx?7zYxb1zvkX%M{U7-ha`wb?*wR|%WYtFFQw?`TEbOY54P<<=2Ky{4 zS_-K#EehA7uZmiE*A`$j@?w(jbHqorz9q#aV?^FuyP%1N$cic;-s@ClkoCoASacE? zD_@9AV4g3kwiFd`nAFk`UmlBy9a`X=U*zd!)=0;Uk&%=PQbK}@>p#65I6($ym95!oAL+>^!Zh&s2*NBI@CqmJ=0q3qU-_l`iEXO4`i&zKM_bV zTQBXDeSL^aD=2y`L;9eCj1SiRRe3Kg$;=s1og`ulTVd#?$&RD#RIv3A{!roj1hp?w zsDAu$KE;EMPQtmu9*#|^b5oK&{>I<3?K!{in!RA>r=X+bipjE%=>A$RVt2(fEf0k| ze5>sI&6$XXSm931%VO>eiTN|*mnp>|XyG2ex27dG;!7pw+{27r`JCGwH?7*$_){CZ zu3=Bz;$S)^#JF_LC09i9FLK;s#KkQzqGR#FoEKTaXVABTwx&?jf` zN~+L1z&z~Ti|m~2X&o&9I8aaA6dAkO*2r)~28HtSyHui>!!+wcCX%Y9;47~rZ-qDm*&naRS`t_ zx$aU`jS_vmU-ahGl&F_3{X^^^-|ETZT-p6HU9rfbVRU3serl*mpMm{%>t z`&$cbrB#J!(keQl5of*6R<+U} zyvEKDW?6Ra$xH0J0Qy>wen8I3A&VdF-e<@fkruhvGD<^>pyQtRrk8(|`T11+zWY}O z7n;;Kn_8~;!| zvFHVfeUh*0?K4;Xzgne8p0a-JxroqG>#S)sJ@9N^$s+dc$ub^+BRp1pCis@kO`UOxszO>bWED}Al61`ID{P01u2s51#2H zV<0k6?XQl;gbW*ez@xEQc^Y+=&4I%`Z?S@FH56YV+!Z1eKzfF>$fR(Yf#+(fMtK3_ z>7y08pBmnv0lTioZcixv<}f%M^$*Oniz zW%CD~FbomvNHwX8(GjB6y=_?N6#R4wm{gYIyl^JFp~H`xY5LuRI)QZ^hvT9}f~BG} z;>m&!mY-dnPT)are14nUbG0-Lt=?okJ#|i3`6jNMS493~D!hA*Iz4q`r~KsLP7j96 z%Kq`^jcoiBD_SiZ{Ug`Ba*L7GUC$TuoXV0asV3!b=(UQ9I$!M?(hOJ+y1Vi)dAHAtJjr0rogOYjZ*J{NnEheJeXU(a! z%WXREjO-G6t|Hp+w&%o}Yk8}MjayK^%KXayM&83SWaK7;6%pX$-2cWf26Y~4U|@R1 zZaGes;Xkq@NZ04MDOD6{7MmROz*STS@ura9Q0Jju^Jlnx>h!0Txia+-TdjtwOs(8Q zwn>92t?K+IVA71cp+|D|!j|LJ6L{*qm5eg=FSmYTTa^v<51Hn37%%p?FT<+JAQT~V z^@1e|_w-AsisjOnbGlb9ensljb8Nk!=_REzdzkK3f5WUB_dXCDEmw|CrRj#Nw^R2l zufG@ceQg~IP02{PbVq%M6_4Yc{dcK%z)~)Y9Qk5gkSpYRuWy0-@buLFVyJQYcICkT z+*?1Y;z*@#P-UDS{>PPysimg2=^<78ah$2Y5zf<(s5xIhf{F{y7eIn~)W!Y0frC%R z!IN@Nf-2t6!PnKB;N=Uze-_0V`P&mnMR&wWz!Xp)=RvH%+dB4-{>nef)g<6B1RRh+ zNg)9^Mg|2IoFS3Knl<7SdN_iwE2*YU5(Xne5Ag4wD;o%u$ejULY6q6GkqYl1P&H4m zGL2Y>_qzb_7YPDYb#-#WgJpEYx;ZP7RDRw#2e70M%;>#|^>yxeC&JH5Usq=r0$4`} z-r{{-y&M4)7*>PvzFq*-%M%oea00v{xcK7n5U3Maiw8iRTmwM5AILA%4z*ncm+0g2cRAR6dc71 zfP?ylf-`_d0on5h%MifIJE)nfBf$lfoEWIF;DBBxp01>_Jot~A#oslKB=h__hyHE< zKMY9H_W!UJXbN|nvmdChzp|T>p98T_Q3e5?wGw9_-bf&b>4XZ-Tl0rXU=XPJuS>8D z5-iYkbx`tjcE^LFK#d4^4^!fW9xjk1Bd96VpCyvNkO^o)5bH>R3jAMRqznuoez*Sq zPd|_?DWntx{1X5F@9+O*jz7o7U|^6x@BjN6`ETC9^3Ol|{M}gK0{Ev75)DX!yZ!e) z3P4~;IwFllye2AB3IuXepm+S$4h7DO26rLx;Ue+3-*qm3h7(l)R5xj@831r%&PZH! z0I5_0yutzC03@0tfFQzBkcI~kq}2)exfKvU{RC;NfTxoG;K_a^X zh@Z2$1K?T$X%KKjkPg?ukv&1%|7b_rU7&{|b@Bp0=kenT+;gBLr0d^%8`KYJAQI)r z76;A#qtTCT4nEl@-GW~4=k(yp|0(Q`ZB5F`KxT+L9rcrmpMK;g7ogwx$py(l03gmJ zjY;$wsGm~%6MYMb7m`y!ZQHX%!mmC1tL6j!7sqNgSUs-@7(>(2q0l*|Jrbx2Z93N( zN!~A)EAvYBx!8*pu@)R*USAvW8Pqzbn!oq#%~zVaoGtC|+H{>8b#?$J7il-_pcP-G zZ!6z6n`-t^ZD^Ie<^~(?_c@vuGOr$q67<@u@~xuGN_)+(@g@^L{o*mpzO=Rxf1Z=h z&@5v)OSG~44Ou0(@6&|5JnI(kMahx`LP`yTBPr`8Dn4W8A%xZ_A6H-#)J`WLChv{@tB#XG`8#+2>~= zJbr<8xaPwnZLdEt)<5v!IIWmiS@(K-bOnvF%~I&UNQtXqcYYjGnZOVjV;DvX3zI^d zXu{thz>_PaU!7gJj9o@b9cwa8s#*n5#cl1nLuHbgkqiTW9uxd_@4z47UBg_}L`&;` zgm=)@{uA9n&%7nJgU6+$#Y;Z6ZyWfpV#NPjYzL1tNgVtN|g3YB08Ck^yFRDu;LM9S|_Ns{*0w)?G1@}Y8@@c;iqr8MZZ{*Z16loRzsR@#5ezp1OfDoXY{!c`0lMMF)i~gccoa#z=z`}|YC=#-LPy}bGEyo>jDoM~}7H*$6t)BenfQ|L0$lauk6 zTQ7*xaO{&m@oZYHm#(qx%FERGM*#=sdKa68sGq)@=X;=fPN$EynrSy8s00h_OSk8m zQ`}HW@T^IzyK_`jRd~P*f9hG>P|N(t<7{NtaeCLF<81rAUs=Rm!+u`h@%V%9x$PV9 zPY0k52S{)mj8W~Rim;55I60GpB~aP{(3CBBv;qAH7kA*5i^;2#Vr4M!e=A%tNt zVbMP{6X=U~g1`VY8Up+014yCKAS!SI{-F_53E~Iv{7Hi&;Y2L(ZyG{|h%~m+Fx&gW z!0dZVJ0uLOR@q7e^QtX0um$*>``_(A*@&oy`1yxFFn=P#hkw&x2r!A*LIV@xZS%q~ zM0E47zHk_NyKFFGmbRrG91h#212{^0Ydf%~)HeRWAh%66gw%F9K*+k4KNu1Y*5dw~ z7Z?%+=2ly2(%a`iA-4BLp}|~iOJ5XbyYA7T61KKO!-+VK`1x%QG-A8FXely%rHNSb z-+iSJ+twPm3W#8LOFJp)ZSzXWY_~ZEPQ=Fl?u$X7x2+8f3dH_f+sP1QxjO+?aL`eKMk_}?@c;xo0aG?er<{$wyjp#HDEa0D96Wwz3gV5Ys51|!#3W_vqg z?z5#YTn4sP4$vR`*_Yspb9Kl2l3qD9at+3V{sn+yy}XF7kpwriJe|A%l4~SI0zU%I Umq2nwa3m6qfd~ny9aV??A7aYz" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.figure()\n", + "plt.subplot(1,2,1)\n", + "plt.pcolormesh(\n", + " mode_wvl,\n", + " diameters,\n", + " mode_tran_arr,\n", + " cmap='hot_r',\n", + " shading='gouraud',\n", + " vmin=0,\n", + " vmax=mode_tran_arr.max(),\n", + ")\n", + "\n", + "plt.axis([wvl_min, wvl_max, diameters[0], diameters[-1]])\n", + "plt.xlabel(\"Wavelength (um)\")\n", + "plt.xticks([t for t in np.linspace(wvl_min, wvl_max, 3)])\n", + "plt.ylabel(\"Pillar Diameter (um)\")\n", + "plt.title(\"Transmission\")\n", + "cbar=plt.colorbar()\n", + "cbar.set_ticks([t for t in np.arange(0, 1.2, 0.2)])\n", + "cbar.set_ticklabels([\"{:.1f}\".format(t) for t in np.linspace(0, 1, 6)])\n", + "\n", + "plt.subplot(1,2,2)\n", + "plt.pcolormesh(\n", + " mode_wvl,\n", + " diameters,\n", + " mode_phase_arr,\n", + " cmap='RdBu',\n", + " shading='gouraud',\n", + " vmin=mode_phase_arr.min(),\n", + " vmax=mode_phase_arr.max(),\n", + ")\n", + "\n", + "plt.axis([wvl_min, wvl_max, diameters[0], diameters[-1]])\n", + "plt.xlabel(\"Wavelength (um)\")\n", + "plt.xticks([t for t in np.linspace(wvl_min, wvl_max, 3)])\n", + "plt.ylabel(\"Pillar Diameter (um)\")\n", + "plt.title(\"phase\")\n", + "cbar=plt.colorbar()\n", + "cbar.set_ticks([t for t in range(-3, 4)])\n", + "cbar.set_ticklabels([\"{:.1f}\".format(t) for t in range(-3,4)])\n", + "\n", + "plt.subplots_adjust(wspace=0.5)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "mp", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.14" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/unit_cell/unit_cell_param_sweep.pdf b/unit_cell/unit_cell_param_sweep.pdf new file mode 100644 index 0000000000000000000000000000000000000000..426c89959089b028808af9f64c00c53a5dbdb0d5 GIT binary patch literal 38133 zcmd442|SeF+dr`O_uBqY1Bua(ed|46b`iZ*14>}{yT6e8=8 zn5=`r*ynfd8#On+-*3uKjvn*E#MFr;oL=9DJc}co~NS)ieKK$#mv^h3dOHwW@Ujz34sYZC>a@A3;XlXL~*JIlpO3` zP-u7p|7p#07cIsq*=48gVvN`Mp#4;K`_iX9-~*uuZ^3;(L1jH$>GT+RSQH$=U&1)y%8e107Z ztb?nwIp7DPAMCGgao*bOn1ct{Q2_iWEGob!jzSBI@QH%u`Q-s);At!hGf$qJy}bkY zCJ0TX76$a=Md_(o*ju?=KnX7FsAO#m_(kz6*#aRbSeQGUw}8mTx;R^y+0lB$55!7h z_lRx2q}!=~?*QH9?`4HK$qL&}YJX;7cB2j3r#f=n`2c4<9=G~uhm+Km6AFhz=)4%c zJ$&-MrxuBqOt|^_4L5%@oa>a#J(NsRUOz)O%a{Mt^cVU4<)kmfNR{2KLUxaaMqRS~ zbecA*?vG7M-!=Vtr^AJh`>Q0)$|FvDa=fYWazqu8yuw@bwmw;RkL&A26}A*sqG$1z zOFMTl-+#H9X7b%Se;je*2SYiwXU{M4b{t`#fmh>pq^wCM+D}&P_GqS)+)Vtpt&z80 z=I2S#0yyl`^pj%M#*4;%lHxh9Ujm?~fGxU28nEGi<8P2_t8Dja| z!$ZI8d8#LROK|;>Z9lSNt?{2QlmvdwMom?>@~!)BvJ>@QpLk}F#k%F;)}6v^t83Jb zhVm7ZOfeo{iAJ9wuB&kqHb_u8tz@90`?eb0d7nsVIox3@KpQ?jTgOJS(`csM@i2s+ z|L)hSc-C*&<~(=n^E10CLW7b;(?^P<-WWuc>szNv#pIdJIo(dNA%5j5FKte9=#ViT z66tbcr+v8%mFCS-r*P?Ftkkum4Ns?gz7%(GyoE!%?HMQM2|l3#XaD|xr(K5KRqW3(P<4t`ZD zdW*SC>}Tq;9|yapbEWQlnkqZk=9+y*+Cl$oRO;( zHtUm?*!D)Tg~Prb|!0 zXPcSXF4#ISzRt8ZcI;~5)ar&Z3*(2MTq90gq(8RR`q9j#s}6lP(5hKaE3SXStF&*75KvBWtfZP zNpr^2oL!dbCQVgmgFh9pJpWV_`#jDU-Q*WZR7gBx%ejg<*JLZ7t@p<_x$jPOb$qU` zDvPy~ee&J?`%}j|3HPERZ!$}v_migs%Ne{3&K>^bG2)xrulG++Swhis<{RY|kJ?Am zg&R}*PM`jPO0m^(z4l2r$@aXW?t#hpqNfa3GWFT)CX%Kca(hcV3}(sY%;(>X4{k9h zKW;tpvg5qt)%7nUv(HN=uU2NXH+{yyj`z0P#CrTN&w;!8cJCsL-l$TFnH=~XYngaN zHYy7AY-ma?7`iktap~lzfXmPHS@q5aKV@3~;;GcG4+Y=k0^MG8OJzropt8rF2HYLc zlAd~fDtf~j!#s({-x=$lhLLo-)+Y4c8^7|-t$%l*Q-^ouTxm{!GhJXo#9-GOuQoJo zs@Eye4?V4JImApWvRvikFso6b|G^N0HJRMGnY<-UjU^7FZOqv5B8@ZJ4>__QUo_%P z6SzNzR*LwU!B;uw%awL++Q_xWDdVKCh)d1)o`>jLr}kZ-+qn;|;B9v`isCPCl~Es= zo1?p(=OS}?$noKK5<9n(LsxyIC*IUmdw7}tglO`m2l;7oH(${8a(CbR=jgg^zkDxg z{H&orK|W+lGtz%IC(=!Ol}Ie6ywut&VESyrq2shiCR1{D-~W1%b0GT+*N3U?7_E(8 zGC{4qyrrIE+q-Q>n{Qw1QUX4^~&Gx;0FF%?+#> zs<=$EhQ~ISRtlp&p>xJpbV&7#@RN|UCk^Wux>_C`D!V8l^{FHEv#-C7OMJg0V{>zT z*PRrn5r6WR*;RDc+19Q=BLwP}(1;cg{|h3uRqf3M>6@pw7agVPBj7xR4iy|HmhwIk z&KCH2T|{ zo}jm)O#Qp6y8<ED zPaHmQx5kqTeP+npOY8DbrK+16o6Ob{^J~QnN?pON$C~7(4vL*;lKPI$%MrgJALlLR zoIa`Fv6J+Lsk28j?aD1x!NZ(mb{e7D_GQ8;4PTxfej1}6E8W{0TCHD}Ca%fN<{CJn zQ&kh?aP9LQ{em)9tAzSfp*WpHqc08oYbY-3R}ek~z?TSLNaSxUL<lkGWIihW}BmQf~5g>+O&{+!%4KD&O3TR=SeopUC>zGDML5f&$XFm4+ym#VM|4No*4+eG+;?pw~`ue z7TQ{JEB#gdiTw^om|CUj+%H=TYVw^9^nGY{Zjb1%W6BaX*4w=OKW{fYVY{BM%h|e0 zdtw#k&COi_yR2IV)x>sQIk89l;IXloKH15<=CAhYebLE@yjpoL+Z!0qBv)`8yug5d zZosePP-h*-aiwx?u$;p859UADoWVL5Z!xN_Q@D_=z-Uu(&CloqqZZ5ub|KsDm5Vy;r|?x)U_q;%LM3~ zv9u=btLhs@Jax9Snh;BG>zuhKe^lwVqK<|tPk%_v+P6%(ygnZ|pS>7spB~xvuI*h- zl$yq@d19PeS-au=g6o@HId%@J?(g9IQhcrd5PK1S<*x4SV@9<^VQ16(tIYikg*>&Y zrmw19v21kmd3fx!d+$!38`^vCT7~zH8yc|cpd^P*uw9d++gzU|hKegwiamRT!$$6w zIbCGnMHWo)elKmUV=M`;5A>=(Z@zGG%*aN6l&89Ct@5`3%}}$lnHjWYkK|4nL$lD{ z!=UZw4)U%h#ws=1=_HIERDH;vS(z1AGesV-2xT}Y)6&b5(6Gj&em{)Pkl+nX z-4$5$Utz)jOsKSus&kpZ#%98HBNV#Zs-IM(k$Q7X&j{?QabzeeD+^F{^}E9}P zool}d<%^Uz`yZLrbK(zwIY4vaasjrRA?3T{710ftyurrwqnY=w87Uet$c-9t-8{v6 zO{>greEo}dn*AF-@HX~-G3lc@JudNuIqF7QZ}^@)muxpEpK4ZoYQB3L&6m-8`z(Hp z)U8$WHG8CfLyAvc%;d_BQlos1i*LEkA2V}OSk?RktNon5tYQRorS&QhK8-TkC(d&&M$v?{Bd8YTF7{ zv?5~u3ujc+%HQ%%fWA*;ol^TVD>i~Ncl>Pgr4bC@fun0T>GZR0=IMye6ewtal=ir= zy)d@jHiD@jBrYWItkN%q(@qs%k4^MG+M4>c=WTMMQJleH^!m?xI}&cBD;z(T`nJ&c z(6jAXy+w0c@`;jjFV0pqs?p1FWnZr=lUb8LYv*>}d1t4xK^I2M`Aluq1N2tCN1-|@ zPcFVZnuC{Ig3>xJ;DtwB3eXuB;J;Nz3)`N#NT z?ytKVy~dJ-##&W;-4+~;pj;E`h8)kebLeg8YZBsX_Qkx3UQcsSe&?qZ=ve{I3jYOu zUzr}6YtfY!7gm{dT{1v5ZPk!_kA77{lb^_NXoZbeK(0cfe<2A)ZA|1^4*E*Zz@+v) z9-0nz>H^ag7mHeS0xEjl?H6R1b4|C*J8o^CirDh$XZonv+Yle_va5$qd^~1$tl_x4 z$Q>c(L_zh4a+}YES&p)Mud%&(#PiCN+a@Yo)? z+TN*4s+VP%uhW`a|LN7tU_}LkBJdt^Ia?>#NU_mj}LW zqOawaeRp1CweeTAO8zLlh#ia#;S7NZ=T*2a^sg=ZOk&yZ8&s*Ahex0~_z~yjJ^P{X! z94855F}He6hl8)US6iOjcuFYlxN$!^ys)GyRIU1)&jn_W>a_TJ7hM;nE9rGQ57hmJYB(R~yUEeT0q(D3Bn1EI5 z=AY>y!{rBgm{a5qpG>SoPZ_to42e9qD(T2MF}XUOtHy@M5)=Z5e92XroLgCTq~7v> zAl2Yd@d(etW^{ioQZAG@I=uJOIBh%bZPHQOkh;eJyG%M-Wu zb%%FbX!G55kID#lHp;%jTDI@&Hck7YIQ+n~b>F(i4hyf38WR34bMDmPyDScb_ga*M z^NQ*(jRU3Ff(lFV-#_X^MC*wpX#`4iMw!KqmD{R>9Db_t)w)Y)Ogf8Gi`ITe-~A(x z*FHUFI+6C^Op@%;8so9$q7z;La|!-U#Fo~4!>~qh{N9C^pQqb#MLX-0H-6H~zBH`pV3^a; z$^&~L>pYdaxuLsw3wAF1Kl%qUndZno~+J_AxR+tabD zm?F)85`yn-eOfaYQ`nqZ>~Cf%utKE;K-n@xDJb$6NKHgDA8-d@0p^YOaJj-^0jM$y zifAPNljC>7E~u{Gzl%vQaQOA^>#17OBGciw#QXPk1+vh(zT-(p3w(Z|qT3O5g!tmL z*o~MPHRCUg-2PTsfqdFZ?lR>^O)=X=MiXFL3=A;oxiC&(4sjY{R;>Stcg?#I(j?e z$*~(8G50TH!mS^KSFYP4to>nCHvJ{xGd?Rd&?XNBEY03U*)e<58(ElEas z4*F(HVA#=;a8l9J-XSl_OU9>7>kKoU==rW(zbj?5J>p6*J)y34h)pJL$E}^F>kZ@X zyb7|2)ZP?8+>|Zh!O%J_FebtMls50w!z(IuS%Kv4H*UI(t(|OR0Yi#4Bw3Q4P(uIP z11FUakrhxpILj=pJL9_FxjRwK>BOrcOoI=z&v5TI?^IFd8Ep=@-aup?9Nvv%QjQjD zUV+mUU=14m7o~|vWWF5;FpAidZ`RXaiwX5A7ufrKHr4NEdcZ1?e<^`)<$N`m#Slf9C0Xi62LHsdxrW-xfb8gK^Ta#i!#0!t||#vcnFU z3Xtim*4@3a{=JN;tJlB_=P$ynp_~QnpA%PoeD~hvW^(mTis(%OpOdTJ)ZKH;bDR2N z&thd%RjtMDoC~O$Oh0rs^;A@Xu$ffO^M5jW0}8Kmy+3~cuq`n(vA6yHY93ka&L3>TLtfvdPS0h@x1XCb-7n--@6s3)4%mEEwb_n%aB#G@VT~vVpf2CXd$ux zoRvg$<|990#|&r&=FBEC?BtH5=^srw$W6;ZQ`z&Qg2`^Z-4|(MQmo~efNeANy5jAO zVIg1FznNvr&3eeU_DO?oS&8*14@&d~S?JTPSBk7f{tP;j%r+ZRmA~|P?@+YI$;3K^ zeO3QtpmU7lHeJ;CBkuQ^7e9rBsOzPR)%+6Q_sS|%QFp!nz6jDbPoD9^s``PQABm66 z_Z`U!mg3SRWh$R4rmQkQ=yPiYeXN8V|Drbw2c%2yw#%mr)3A20(OF^T6)>70`Y*~l z^{`4PQ7I7h1UZL>!@Tw?dOJ>)t}kj5XYO8gjt09<-f*o=XmwCbI6a-7kLb6W%nwZy z$u&_~-pAL=o)-}cb_~vXO?U9gcKXX_okAsG&BR93 z@ytp8Jv!NCNB1Ud?;45ZH6D4y9m~tghc7r6GvPXE7{TXf^uvBwNJa5|UeT$2-v;Fj zu6Ox9-@X6nZo`_B8proW9iMn-e6hgC`r6A-Utak^kMY2h+FQ?ER@&&-v-@`W;ESic z#7E5?%Jg=)9chut&kOW-w}=ni_;%luH*5Wneabq^9=20v7c%9B?+s$gviC zVEY-ht2bWdc^s0V$4KuiwB0kxw*42K%UHnPxUCTio6%eb%t3P@Osv~_!^(_C)<}F% zRJ!bXwkxagu{>?|Uj64jr|Vp@o5D6QA8NB8CF$!P*LoG?b>cqb2AaL@Mr*_GxGBXv7@&+<*xS~|C*4&7bPyW6Al z$pM=y_2smTr82QR%}K|w2F%D>62oZ5v>pr_X-y-sN}rB154OXt9J z-l|DjmKTRsp`ve#AL`_Ku4O+Fshs|R!C&)*+DPHg&E%hN?QibAcKE4Jk_s^)BU*vvx;3;=6Pxyh4qu(C)gzOfX9qS6Oe*82_M5)-Q zhj3HvScX~p-b2#9k%CpzYojA)H%NEZZ?L0}$*gT}DcZ`himXyPTif`=A^s6MyR2n6 z!SZ=n=9*6FH1*fq+*Xgv8jicoc$SkVJxeI7OcD$OmWxqjC44FN7w4Di;i~qE;QW&C z{Oz7sTdSE)Ab`bzKf+d3EQ5>Lix^e%A{%z1OGg^5ZHn(=N9HQF3mR9dz{jYQFJ=-VHO}B3&xoBUAUN zaQu|?R`H158oE=C>%?275Ar35{7`%pQF@Y_CHVB36Y0k9a~n5tw#akW+U#(jwoZR- zU~J=AD?dnPFWM{By8_=UAWyNssFrwdy@9HwH59cs8#le*JpSaVOyq2-;f^*{ZD3ml8VF)K6SaT9Gg0ex83$s7O>XwU^}m(-8L=HUMzd3XzRoZ1g!uSg}^cN ze>_MQ(aN_|1W*{i{Y?9s+pdkI0~>43Q$!zL)xRF@`O%tNRp@6jhqme`^^!3`>BycD za{qHBZ&Q`RtQeK}Bc@^X5wDdv+y%rLw(%d0+aoY~ODSY`fuo#mjK=-fBU>jVd|bVI z*v#7>*?qulQK07x;Cvi=M%e7@xR3l9QnOXGL%ao9bfiV6;J~SSa+h{6R%B0!&FJAi z9%xEG78JRzN>J(T+Un7nUk99j82Ru=+>W@Ht3+eC`+P2Y`C1x*i2H|*j_W;$=d`Y< zA+kvI5ZE8_>qRdYjv%$O&tFu7Z=dlWGs9ZUf093@a7O3M5e18jW~W?r&Fry|+xxn|}piO`|XqeeS zm*SSLJ^%c%`5R=sf&xNa@V^H*nJq`Wd*Irf3f?i$Qne{@lr19HQk{u(q_c05{v9t89iaS6D|gb8yQJ z^g7S0YReg}KNl#4s~sZ_9#I?@Fds(=mb#bisE+z7bw7WIDUQW^DY;_&ANqa zgL)1p?XBUvazLQSyLPaG7SA63SOAp@VR`=xeF0_Inpt6i#QsWia#(Zd_8>+GbVi{a zpf9{aP(2}U=6HNwFaosvrx!nhyNTe|q_w%6y_Ky6Ae>*<#lr3s^g-Rs1LhHh5yw#9 zTU-Q7K#71eQ2~K}`3H>=KtX?&&iVZuAWIZ2N(=r#|04fD>q2AyX8ol<6B84lT{`yf zzoU`sgFBgv^8wR;x0c}2HALpuK!bk)h7=NmWdyGUea46h!{eeLK8iwf0Dr<@UlCw7 z&_ONqj0pJm&n+G5ET@2iFYV|+CJkR`0iX0xLJK!l^-)6bB~{40g&?#Herf?IfIiVK zIEBzcw16E37zsQH_96I)3_M{4^g2I$0diA#f;9@n1$Y%(aMyHx%nk*}hkt610*Zn^ zIiS!B{s8$YdSNd|z}fuMe;gfd4uRd^kp+KW*mlw97xs2W3B!9sJ_t}OOa%`$3y!}HbLro5Vq`$n^?$yr+? zSYQ*TPRZoBGexqUaO7$O|9@Q6{*Sl-Sh3OvA^B$YqqvUwAiPJTWXVjxbmylGr zMeoxOJ99x9p1z$H{w&028=|w5R8%ap1g5Kx{%Tnlc`?{D?I~NS#40OFS$4pOYTwT| z?zqofC$^M^7jX6N;&$E~x_8x$qqk_Zy=#c7_rJ@1SnGVb__f@<&4&!sf}_larL{Ml z*e}0*pU%33Gw-+3@nzn7zRz+;rd?ypb3ItI14*rCE6jjn!KCN(4Mde>4eGqie_@L%PP|60ZeNfxGHF}<&l@j?1QWPHGqowdCycK#ClV)_(W zaP2?o6X*%)gD|&?%KRhCmDA@k=Kt^0r#J|MeG_qknN+wP1Tz9*^JwSJX3T~tkLIo^~=u6h1lms<>zX0 zD{$56Tljsu(NH?r`1vln>>{1D=f$n;j(CF`Pjtp+Xg&`(-(2p1kXbECN=Xm}gK&^- zl!QnEsz+cu8( z@6994hBhaS=8FF$C5R1NC5cM>Y>~dumzFX2%Xf0Fv2^-ZY8giovB;V$KJP0BVNlM{ zzcVbI5;o~qscYiHv1GGa=`gQupD;={Kbv2RxaL;GgANn*ss5%Vc#n<=@-1@&DJJC0 zxEYfYQ;MdktHa)r$?7<_x1M`gVUg4}P}(V<$kP0;CQN^XacIh?b;dAnO?ad6EVAH=bPnD z4PEoLksp_9m#^6Qs&8c89$M6!b6XxbwYV;5dN9S#(Eh9VV8Mb4-D`7s^DgXlH(roq zp>^|ZL*dXuK8L1F1Ur}it?gTjJRgcfwmnaOMZVX7af~r+?f-(K4+P2IBszT$xBtrS z`D;$rMh4G6HQJyl%H`8$yRaA#Mt9ZzvOnxdRDZJ}7WDtGGWfqLU#HJh!dc0D{&@ai zO31q)mzyOMtz0Ezy)WHc!`fyaWUzD3$DkFBt3vR%AOE+khn4UY>{I0Vib~-?M_D#- zO!;pu`M>QCf8uS6`kr8ho}I;g!Ls8YP6D$k9q;!i-3-3DO~yeq_EDuLmpVUxB)<7`Za2|@!XCZ7O6m_Na{lFnu%TVO#F~1wj}JsyuE_GK%63VXD&Vhpg1ex zI?ymO6?eddQ|F?J{EGv|Vh)KsaJe-P!kyAhyc??GQW0uMhRwx=I@M#VmWkaKUW^;D zm@%}N6(Xwk!PMp>sGWN0ZtON++09!cAqb|ngj!blPD|i-MZWd2SknzZXGl#d>F2W5 zZ-Ws;ZLvJ>m)~UOZ+ihmPGXnpvmx97(5aeifSEULRzmr*HB**Tj458K-FSm#*?we6 zwwI{pT@o*UiCZm&QBL+ivWV)c!iiY5Sa<$LC#CZJe_R~Gqf0z;(mfoR#7W(;REnv> z6(M~Klyty0sMHBkB5<{nRIWCe(1YA&(w%nF2N%b+_!JyDWJ!h@(;-BFeTv8EOs6I$ zkpZ~1sF7;P&O11Ak!LQIEi6ao0J^78O{mCsg{_f~+8WO^STs<7w+G5tYU*4t4v0b& zX0gk3J#S@5Ujnu$=D>CYu>pXMOSS|YSH&qW$^s5@QOjJ_752B|)JBXQSmtj&%Z&fd z1Q2xa;^MHN9j1S=mv#VqZ}8(=^nC`6nH1yZLjwSI|L0J;#}0jG3WTfe!3MF#~YWR!!7 z=zoVbT=tUN7F$+}Wj%<3|21$4x;TEp-xq3kxXyQg}-D zeEbc{K;!C9t=f(Ivdjw4EHh|GZJs0_Y9)~p;4xza#g+Usv1tIk++?QUy8lF6R)Pox z2ZhWegtaGTLwF&H8ijYp!|s&2=)16W+pK^g-UD~?w?d@KdN}bQ5Xt0`M@#(4q0|oP z7{A12i-1p4H|#>OB(S>d8iWxf`jliZHbG#y%T7FwnVVAO^O>%9WuJPtK}Le&`cDFc zNKBu(&$s3zswGP)&3Wn5-GliE6G({Z2U0J2K^>_vd5|=4%+%T{apPLznZrO z1u67{UoS6mVJ>+1FPmT1Y&I*Ba5X91r_;FSxVjK ztibTe;3p-|{ZSDpF+j$phO$dj@?2~Ia0sLd=kpp8M>;RrizN95s4x-XF!msiE`^b` zusIk}$$H*=5qid2!YF_rk(A*WnukE^c{gzGzwsvwxt0BTs;%}iClSL+rnf0;YzM8| zc!yH?-<|5eQ1en~F%Nuu5aAC1-=)X}^Ah%loqvQqZ=9rVBxe2$XK7PJmIl>CoR|X= zgTN3x58+;fWho;SF21Y~n`**OG9!{r^YSN`rk2i@lBID4j+sbz!)>PO*UrZkKXU0T zJ{4?`b%>bL%b%`fQ58-m1%WCaTDLZJpHUxewntBh_sAYLX zYFCu;I0vohju3(v9Q!vQGGPO3o~DPXANGTBq#pu=drm^+%pj5ksRr_~$q0nsWq<}6 z5^(Z^1l;AHOr^9=DbuCth5gkFcVoYUi1eHEY*}{E|%cn zr~?#y&^hixwI2g7gbX#fYak0wmLcY4{~$>(F=QV<*tm_T8J%f@&7cbLtSZca=^}=N z;R4f5XpfM^d2r;%&owbg&AtP(O(#dD(|zs%T-oj8e`4_zaMBQ72usDaL^VHbph%FyqbN&R%>ic>r70WNM2vXl5`#-xK0xo}u3J8h` zqeU^`0*a`xFdw*q{$IQsA}A;z0IszBrC*a?z?n+ z*2Hg>!D#uqaIVKWZ3g$-TG?OOAE>BZlO@wN6-u5LP@$$RS zSDR4o|DyR$l#GOVPyNdycV8tImUK2Y`;|_fb-rSgl1O=Fs}v<)U98+R*5~i{OHZ2J zw|5|+oUydM!g&gx^Za_B0b161B<&etnv2)d%(?q1R=-Rr6mycc>Be>jANWn% z*faU0mQ;c-uC%2m(w|6T#cW`1!9%NgCst~}7hk-?a{AtODEsQ-@UQysYcNX>k%>{* zs+KBsBi!A|u?SM+Hu%m;+x~&w!?_u=$3aGD&raZ*s8;qzrFOl#4M`+^lANUaXT2g* z{G>uBd8VTQN7833(5Mc&Ob~{M+w*Kt^Y=8H?_8QQ$%$6delNgp^7{K}RCo5B5`m`q z?qT}FgUIhkUQ=O{ zcZGTIra^M&E=pG9m-P&g1GOCD;4%hj&l6T>hi2^=jfB177yU z?((UZD|faYwi<65fHuWkuvGq|EP{{wSutrc|RzH>Q#>>&4{~#VVioy@EgT!od)X^6RQ^3x4dgl6&5LW8F{{<6!!^} zo-yOpVfD1cGPJkxP0++I~dRJ`6p@K#uoXfxQnHL^&2QdE{J`} zX^jyD8faW%F5kQc9GWi1{48`FZ&8F@VK zK_MnjHepdVAeTFQHDPKUH;t=L(=cxE{?ukZD$=LnSXD>bV`C|X`MJQkl~E~bxmZI@ zwiQ9j%|#1D%G?%al+e(yR+qP-59itpfUA5Rdw_zQ-ia2bu)|SatU%KpT%fv6+u& zhd345XByF z4Z|=87tQY995mLa(~wYKOdBI|e4|*)u!7q_Fh0FA%oUy6R_3o1l z19vztEZBj2MjK>_(=E9YR5U$~zib=?S?(6%ga^WMvn45zdfDzu2vN<4c(B0gHi-B! z?e!6)`oZJ)2(t8TutBMSj~ikkJFXmPgAKA|BErXPSp}uQ26Lr$l?bUCl$yx{smAAc zX(HAsPWOWhxrE}c#7SjFc11zty^z)BCDd{v2ns+SgnYii1oICn!MB8eHX05KiQ{=aBg&-4FL)UpK`t;80a^Y$$-34tlytwabw5t;cX96-Ak>9 z3*JF$@GZ@=8*$Jd)PNl`DVN{D`NM)!dznZ}Lj1R8#vry!Zdgl{sA>^>VTfQ~vsD@c zu|MAX3_(RpzOz0MSHD-~4QeZQUY}T(2)I!$CLyffx!aW(3C2xhF$Wikh9XkDZ&#qDeZ(@IZV1ejv#oHqUlUm$>0pnO4zb0aw@TAcl@(*Sc zKcxApRKnxt6|#|wR6u!hK5}Dc@P3D>^?9N3N;l~k>Sd~j zp4fn(-YsZVh0wRvaRLXB3fg}t0U=#CWAtIj{4&YIh;g=_@9~hCRpF7bT9faUG?>pe z7VSr^jXfD1Tx7`3Ht_ZR7A{~bXLz@al^0Nx&s3Q_^#ffkuROfL_{^XP9km^NVwuQ> z4A8R=?~62KW!!`!c%W{&9+9s4jk0qfZyhUfU5twzFF1-7r%T;h$;_^{D#q@E)eP*V zj)UDl*EXn0mpVC8DYbu67fKnXzQZ{P_M5zu5&-*2B^N=V0HN@OXuE;E2y*UT-!F&+ zTh9L-j**^aGl*r@CwzfydQ9fRB5g2gc{Vi$fh{lg2a+&+n@YE}9S z!eMRtKd6VTY-19a&u;QyFm*HqN@$t?B`^)O4=~*hGO>t=5TBj^MTBSM zFeLgRW_dy?_Sf%5kollc5^{-DTr6^#@T5E_?LNZKB5^>ovdJAFN4$@-F2YrgNQjC; zb=wZQ+(kp2$3(6z(R5t=&4~r}A3TrP&$hk18M(~}JvP&#nj!of-UdkILAspHe~qw@ z2fAKx?q_a8BJ6EPb^_U|%{EpbroU*`2djLhU_bqX5LJ*%HEv5f}DTOVFxCg$vr1*+E7Y5y}HdmN1d- zjZ~}aRZj%xLBqvipQ(p6Dm&H>fmQJABnN7{=z7lbMHW)OH7OFolW2x91O&%?^AYwu zrG6Sx&zP?hVz{X=?_s*lOrq3`YX^dO{Lf90RhS~3Wze`d4T`d1g-!;GB7?(haH3mX z;FRPMVnK)OR0nqws|!GPpynJaHgj1-L}xVjLePZSzjF_Q6I+Q=D3_Gv6EYE;)IT?d z0REWwM}$Daj#vTRtIZZhBiwu3WDR8AZ@v-`KqlK^>LH}FQZYAbouMZq;4qkEPNZf# zwQDGe1#;RXx?bvWtKlc9Ku$g(&*SF13X3IOX!K2JS7X*`1njYnD}vCqU;C?P)Z?Nh z>kD%p0&S*QxdO=erv|ekkC{UzRKkR(=C!&-#1I;mXhDqjtjT~vey<}U;)^9XTY?x% zDaVu|pc&!gx@?HX$OIgWuF)ak(#u8gOuxpX;0J0biQZs-u?^uaBl}x^tH=K=m1DSEg zGlbQq_D7-Xhf(@4`zo7 z#Hsu>h;8|bz15**2#1hKrLPU5JuiNegs^{6rQ_`1tJ+6 zkMOEARc0tz#hvO$*lTK6eHkQI;vQu0C#M+BtfWa2w_YEH1kEVyVe zOlC){PI$hm0djHIe-N}ZPcxi=%>EFSinTG#%rbCbE2WOrdexJvB24neHJ(Rv0l ztMptG>JYRiBi5V`f~{UjDTwJCZo30_1+gX(F+7(h21$jRNKr!=x$YFB8m#wY2-Ub% z`@=mNY!kwKjHG-|!scK?Z4MwtFp0_L1Mrh%zL=5c!B3e@b)0Mt)S-&0zKbQ5+{CFt zWoLFAc96yfdgLNa1EDgsU+7_k!uW|s z5L%UO84Fr;^7+24n>D0_up}=dDgpk6GTk=(rQwJWRrD~Du~ol$lPK}dPiQ$5t(>vGqgD6s9AZzaMo7-#&>01Idt zue(X@^<74?{;H4zPG3YUp5tl{#R5|4P&9&=`dmiaA71rpz8YDJLV){0$vFGT`V_iq zRb%&WxWUN}d!};_aJ{nZH6V&4FJ{Z^a|2M&X8B_a^KBN-pSQs9f!inYLX;#$Awe&(E)D`4XnwZ&bCty=bM7OU!nu~k zf)48EH%*B2i}JV{`-C7cgn`sqEua^lu{5IOBLujxoCY1-VVgl}>_ixuW`klLMZ)D2 zl~H{y_s16IKslfoSBBuf`s;pBT<%eS`KQ-0pD0OV1WJy3&DA>B2EYEK+&=Op$u}ap zxT$Hg--SO~-ww_b4tEj47)&*VYjey1!tm@T74eXIr{wdfp;t8NeN8a%KFgwvBx+8q zCI*f`nppz(UYN+k`r^BV*z=@$I@I~Oco>lmnbPqPC7oDZglBk}2%;fmGvFUW-4AM- zpyc3QLs}=QBAWNfrP(7O4y4jihfK8n-0cL*=!yHW@z>)cWA435vm~H`tv6CuI2vGdH89JpvWlP!^ zI?X;fvXsaoq}_aq4~l0u-K3CZ0cKkaloclS9+sf)>~y3u!o8Y!H^eHeQ3oM7(U2s2 z2N9!3Uur_DWT&PfTZ;X+;nTRKu6Hb>7Q#07pX^71uhRB8z-oTaYh?8Mh29Dsu_ZMP zB34&U+C&5_P37kyOgZ<~F37I%KIe=Pg0ZQG5mGUKGqsOtG*l8H!4Rjq+YvD|C{YTA zLeskty4+jje{ov_=&bbp4pnd^G8{La)F9|^s<&WIH+!u)Ih&sgima23nQ7|3PZoc9 z!lA`a1DTuXyy#(9#otD_P4t>6I4C$E@nXtL*!5&K;Hdei<}a822g)c*96Mgl-91(81`G;buy)WCl>@DPd z8FYv}-`ShOmf0wY&{vR)HN^e1uxT7Z0g1cZEg)J52M|8y>Cp~t+ZqS-TAj!&2h@m1}9m>%cR#Y=fHmxfb+1h;!W z-Fo(?b*E}+59RwL76CM(Wx~Og3%R9o$!ryZo7}J|7((N}A-FMU%Y^EKJ1qKEL~VY*z%VKc(w5lwfE%1w4~RWyxM z+;}`xfLr?5Bd}bkUr+=%Cod1jnE=3OKkS{<0Y&idk3w0xU+93|Ee)5$N zoDQ12gnEO$zDE(m>s?u4`}Cba4Ey27;2z16@+{kZvIx3|P>0*H$PXM~N4sycI`kXm zi@Ub+O3#t&^9>V0t_OvP-_EUz=b7>KA0dD0R@NvHL1=P3T^1{#b>bdSh%|s( zF>}a9JGR~)E-t(eAT%RAZ5#{SGPq>c9Es@g+3t{aY8GJL0Q8)MC`a8(*oBzOQao?)HZ)YCyG3& zU8kD@8=(xMxrqKz0xgdH@;Q{+@sBCV#jS$^)J}$_?{oRr- z62x@SEI;CS0oyS;25TS~lu}f8PQvqK$k&0!yi5U-gD%LxwcLeHo1v%L@O*bf)~bFm z76yp?{E4a{?3%%wffgy~0qm7PnBcoVr z0sbYb3BqxVh;j_0ZAfREg0dEpYeM!JAla_Dll6VE5@hBGK#Jh%tOF zGFKbA<9D9aAVkzcE{XSXzko>LUbyGbeI4KD2u2d`h(M@IhB9&l;bLkg!wPT$IgSXq zn2d@ch;N)L!hEw|;64JrNntmMNS|r=1UHbD1mgfe?cWRQ|pr6!!&ahq~JKF5(Fv>ZOiV4 zObTB>Wm1p#&h7~YLX?>)%RnIgWb$)Qpx+Yy283bujDx+j ze-0&?755FE@u~_Q#i$O>pH`kV{Bs1|^g0vpk?oM0jNmUn7mI^cVBk!R8?cUTpP^x5 zB?5hr`hs&IbeUlyC5DL5_1n%sX7W6lf;_Nh2-pD47h{j3KH1f)E{CEAC?^M(_?Tmh zJfL!E2V+5e!&p5;eKL-Bg-Y>TNy#UuX(tsvYK9`Qj1@9nXlm?%LC=x{)os(`P=Mxd z721bj$v5mX1lRU36(O4-if9QaoDd8L5x2~gDdint`puGYc7*Bo@Ab8YPBHrvE2z8r zpFAIS4nv8u=;GJ$#aM`Uo|C0Z8_ijM`4@V-;^Hgml=Z;t))wB9sH_KGWw-SGb=F7Yc=s zLNTHIq@nDip&X{6D1bpT3Z;%hiKb92D3tG%n=}+7XE16_>7k+2fe%CqMS?O+Lm8)` zoTs74f)Nsh@{qzqL%Bpl`9wo0qoHVmK{AD6M==2(S||@_D0wuLQX0xj8cG@%aiZAJ zP_$_%S>Uq}1xG^(0w1O+6juu69wmT=@|K2z2O|`U8x6%CeDI`Dlqg2frvyq34P`40 zg-ApBOhZuyBfb<0p0W{qSVQrop@`8?SZOFx;A;fsG5C7}rGSQVhlUbPLs<>J(oxiC zD04IvYZ}Ta@RdNJT&5V%P|U#R#}pcxng3VYl}9yo<>A}Tr~#}-)+{5sX3kd;n zdujw+K%_1hw2@7$hAk*ctF_4N zK5A`~cxKL=GxPrV@x3pb74u~b)0CYlJ2-;!nJJyCpq_1yiGXzfH zEdY85szC#c0pKM7FR%o=Ismj7*whyXYB0d5hH($50^lJ46%Yr4yD&C@TL9d^5>8&j zqLUhc%Mdt#G8o-~3V?GMf2NPm2PGIHp9bJ02FCYI0>=S33UN13!nhd}z!(7z0I(kd zU$75=JrMMr+XX-l0NViA3czLnvN7`FCM-H{gs}|t6@U!@tcN%itOXzgf*P8J0DvC^`vDsOCYE{C2S!)m1>;-bi7_&FEY5O;@h9Md;Tsef z@enkQfZ#nK0YJd`@V9!j^PubhzSFNWn~H(C!;SiWM=5 z4#a|lgG--eDT6WsFMqgKab<|$0$#taegd_!T;0&eKI%|JnWiSy z5299fsza!itJI~4GEsa&61Ldtwxj@2UQ_QxloGWQqO6dFO2YDS@5JX3rN258QFcl? z#bM3X&A3DxMXhZu9w~FX=qX+^HGKh&J$)B-h8;4~;Sf%-2*1w>8RtI378y4>kQV4A zUK5^<;EcQIi*Rfn0iH$>sDWd|B=5lv$T*K@K4Mkqxp_Et5P{O}gc$7%!w{q68f(QG z+$w8@0WY*uSt(T3r8q@^jh~f5XfMPgC3ZzlkyeTyv&CU)lAjUf1W6B~6sV7@In>H) z>haXdN$N;y3arst?5SuMFihnv5;dkQ{ z-VTvA3ZbI`_9;bbi`=!ZwSR4+$iYb@VP8weBg$~~d_^|CaP!HDIOEZi3veZKJSbvM z7kj%QcE|{7ggvs?JI_m8=;iB(M0ttu21n$>7`}_x)BC$I;zIXIJhEe0kzIh5Vk2$~ zPT^+rND`JUu0oVX^&-VqygsZxi5lTrj|7)YKNk}I!F9OTaGxWdxcK10-KU)QIldlX zxh?jOXl$$FaQ(1PG1$MyVc<{Cf&a9nL1p~^eoOnHz5TwfqW9I49<;8z8SQ%4(rIk( z=}WG~bd?`$hUmoRykYKdw}~3%Df(6;nk{^@Y5D8s_E#I0>YCoY$XHU}^sctC zP21S^^xLm%(6;xrcb}lGz4uD18d|#z4G&(m-2Xc0Ug*?c?D!tHgOOYHhMUtSbzeQA znKb@d-P@~U16L1pQET6pv$$nX_2*m}XRF_IH@zyZYAAa2y2#1%$b(l(2mUCda^o*| zb!A@cC$Y2l-mdRFpSS1c^WHmWa<13)+&Hn7pWJsl&C1cXgj7~l6nS6?taNH&S7$x?obF}gwf`So?f2mnIr{_(rZLoOQwDtI@YllwG z7&2)!KV`zo!Q)c-!m&x+>cba2eY?sOAD-QO|Kefq>zW>YeRoS^S9@D$M~AVav;EM9 zcL&$E?)$z4-kH+erZ(f-Hlu!X657(<(9-^TT_ton9^lzTOeOI?~kv3_; zGvU0a`nH~iXVq_7AAd8m@^x#s=Gep6Efp`D@4RTfF=gWQ=S|nfkE;mhsks5aa{Mnp zYuq-(?Gl}Hp+>LbdobTum7jZpmb%c-4EH+y_)RfkbS7Uova=-T*uz&x9X(eLpdP;e zvf!R>|D6~4wtSB^#kyR)`>q?$dtW}@rm5X}RhwnwDpNnpy8JZrQjPfHZ|hZ0*5PPt z&sArfef<8>#q={(y}HV@QpH$>}O@!&B%QDL!1?N1&*_9kZ=4z zX3PM(;Qzd}!3p~Bw>D0{;Rq(M{J+@R#z&thz8+Ni^hxYN@;B{%4~?uI5^7-d5Xxkt zGU*7Os!7A4Y5O#DY1(y~a*ZQR8+Ra6Zt&L4qG?BHns}k zzBMGYh+acLBgx05qMhVAVjT&+#K<9|IizfvC|fGwsRn4AGz~meg|wTe3YFSOMen34 zo~lNZ$(u>j4$}1URH@QcQc*lPo`}}c3yCL4=ps%b0gWTZ6VO@gLnPFb=tHO=qqoQv zMAVJwPC(sAu~N}>Vko&lCb~`tC7>b1X+%^*oJK%3gafiGGSLb@mrjdHIq~3XQ1yfKJ3)R==Px<9FnT-tVHK;0sVAdue=`M0rG4>CZe>u}02Q?bW2y zv=1~op2|w&quIt&9hOGRMB`+UJk-+D#e}Eo#7lT)APSE{3MHKX8h1 zH)z#v#+4wyN8Cq2oj)gEk_)%CX3ee&>fD`NEH|V@Z<8Bl9@r>1MCtgt(X^4KQYO;B%fxNN0IW+A;UNlYXshgg?N}jB`mQm4FvsBPJm#}le#T{pzJI-$3Zn?Gf zk7%>4?XRv%zdx*n!6!3lL^C}h2J_H?t|6FGKD6tZoR|jH!27kyKW>ZAD)Y~$!pSMI z0_-nI4`}I=v^-I?JYflCB@4dF>CG5JSt=4Gr3U!;toZbqcVbHXC~F_+?Guw09cx-F zDJG7xL=gQgZ03(DU}}E@35mifm{PgNlmI53$)R)Ue$Wq<>FYn5&UB~K-ADaJ5IsF5 zHqIIvyfCflpAHI#&4!L5agPaoKQ^j;LEWV9GwYWg)H(gKNXX*40tW(?ZWGWRxqY{ zP>>z-Loq4K^5tRz}WW};xmu`K(* z*AMffG0XG!gI}$fa&SCbw83UVJ81Jf_Gij+v9E4S>*HXKQ|25Kx^bCvY-k;1&T%d7 z5gqFITig#g3oU*(KNghqH>=Oj&)=e5F5@#C6IyedmF4Jj T|FX#9a2YIXSJ$bbJnMe}!K_d? literal 0 HcmV?d00001 diff --git a/unit_cell/unit_cell_param_sweep.py b/unit_cell/unit_cell_param_sweep.py new file mode 100644 index 0000000..8873808 --- /dev/null +++ b/unit_cell/unit_cell_param_sweep.py @@ -0,0 +1,203 @@ +import meep as mp +mp.verbosity(0) + +import numpy as np +import matplotlib.pyplot as plt +from mpi4py import MPI + +comm = MPI.COMM_WORLD +rank = comm.Get_rank() +size = comm.Get_size() + +resolution = 50 # pixels/unit length + +wvl_min = 1.5 +wvl_max = 1.6 +wvl_cen = 0.5*(wvl_min+wvl_max) +fmin = 1/wvl_max +fmax = 1/wvl_min +fcen = 1/wvl_cen +df = fmax-fmin +nfreq = 20 + +dpml = 1.0 #PML thickness +dsub = 1.5 + wvl_cen #substrate thickness +dpad = dsub #padding between cell and PML + +k_point = mp.Vector3(0,0,0) #used for PBC + +n_si = 3.45 #silicon refractive index +si = mp.Medium(index=n_si) +n_sio2 = 1.45 +sio2 = mp.Medium(index=n_sio2) +n_air = 1.0 #air refractive index +air = mp.Medium(index=n_air) + +P = 0.7 #unit cell pitch +h = 0.9 #pillar height + +def unitCell_paramsweep(D): + print("Rank {}: Running simulation for diameter {}".format(rank, D)) + sim_x_size = dpml + dsub + dpad + dpml + sim_z_size = P + sim_y_size = P + + cell_size = mp.Vector3(sim_x_size, sim_y_size, sim_z_size) + + source_position = mp.Vector3(-wvl_cen,0,0) #source position in the middle of the substrate + source = mp.Source( + mp.GaussianSource(fcen, fwidth=df), + component = mp.Ez, + center = source_position, + size = mp.Vector3(0, sim_y_size, sim_z_size), + ) + + pml_layers = [mp.PML(dpml,direction=mp.X)] + + monitor_position = mp.Vector3(h + wvl_cen,0,0) + obs_point = monitor_position + stop_cond = mp.stop_when_fields_decayed(50, mp.Ez, obs_point, 1e-9) + #stop_cond = 10 + sim = mp.Simulation( + resolution = resolution, + cell_size = cell_size, + boundary_layers = pml_layers, + k_point = k_point, + default_material = sio2, + sources = [source], + ) + + flux_mon = sim.add_flux( + fcen, + df, + nfreq, + mp.FluxRegion(center=monitor_position, size=mp.Vector3(0, sim_y_size, sim_z_size)) + ) + + sim.run(until_after_sources=stop_cond) + input_flux = mp.get_fluxes(flux_mon) + + sim.reset_meep() + + pillar_geo = mp.Cylinder( + material = si, + radius = D/2, + height = h, + center = mp.Vector3(h/2,0,0), + axis=mp.Vector3(1,0,0), + ) + + substrate_geo = mp.Block( + material = sio2, + size = mp.Vector3(dsub+dpml, mp.inf, mp.inf), + center = mp.Vector3(-(dsub+dpml)/2,0,0), + ) + + geometries = [substrate_geo, pillar_geo] + + sim = mp.Simulation( + resolution = resolution, + cell_size = cell_size, + boundary_layers = pml_layers, + geometry = geometries, + k_point = k_point, + sources = [source], + ) + + mode_mon = sim.add_flux( + fcen, df, nfreq, mp.FluxRegion(center=monitor_position, size=mp.Vector3(0, sim_y_size, sim_z_size)) + ) + + sim.run(until_after_sources=stop_cond) + + freqs = mp.get_eigenmode_freqs(mode_mon) + res = sim.get_eigenmode_coefficients( + mode_mon, [1], eig_parity=mp.ODD_Z + mp.EVEN_Y + ) + coeffs = res.alpha + + mode_wvl = [1/freqs[nf] for nf in range(nfreq)] + mode_tran = [abs(coeffs[0,nf,0])**2 / input_flux[nf] for nf in range(nfreq)] + mode_phase = [np.angle(coeffs[0,nf,0]) for nf in range(nfreq)] + + return mode_wvl, mode_tran, mode_phase + +if rank == 0: + diameters = np.linspace(0.1, 0.6, 24) +else: + diameters = None +diameters = comm.bcast(diameters, root=0) +local_diameters = np.array_split(diameters, size)[rank] +local_indices = np.array_split(np.arange(len(diameters)), size)[rank] + +local_mode_tran_arr = np.zeros((len(diameters), nfreq)) +local_mode_phase_arr = np.zeros((len(diameters), nfreq)) + +for local_i, (global_i, D) in enumerate(zip(local_indices,local_diameters)): + mode_wvl, mode_tran, mode_phase = unitCell_paramsweep(D) + local_mode_tran_arr[global_i,:] = mode_tran + local_mode_phase_arr[global_i,:] = mode_phase + + +mode_tran_arr = comm.gather(local_mode_tran_arr, root=0) +mode_phase_arr = comm.gather(local_mode_phase_arr, root=0) + +if rank == 0: + final_mode_tran_arr = np.zeros((len(diameters), nfreq)) + final_mode_phase_arr = np.zeros((len(diameters), nfreq)) + num_diameters_per_rank = len(diameters) // size + for i in range(0, size): + for j in range(num_diameters_per_rank): + global_index = i * num_diameters_per_rank + j + final_mode_tran_arr[global_index, :] = mode_tran_arr[i][global_index, :] + final_mode_phase_arr[global_index, :] = mode_phase_arr[i][global_index, :] + + + np.save("unit_cell_param_sweep_diameters.npy", diameters) + np.save("unit_cell_param_sweep_mode_wvl.npy", mode_wvl) + np.save("unit_cell_param_sweep_mode_tran.npy", final_mode_tran_arr) + np.save("unit_cell_param_sweep_mode_phase.npy", final_mode_phase_arr) + + plt.figure() + plt.subplot(1,2,1) + plt.pcolormesh( + mode_wvl, + diameters, + final_mode_tran_arr, + cmap='hot_r', + shading='gouraud', + vmin=0, + vmax=final_mode_tran_arr.max(), + ) + + plt.axis([wvl_min, wvl_max, diameters[0], diameters[-1]]) + plt.xlabel("Wavelength (um)") + plt.xticks([t for t in np.linspace(wvl_min, wvl_max, 3)]) + plt.ylabel("Pillar Diameter (um)") + plt.title("Transmission") + cbar=plt.colorbar() + cbar.set_ticks([t for t in np.arange(0, 1.2, 0.2)]) + + plt.subplot(1,2,2) + plt.pcolormesh( + mode_wvl, + diameters, + final_mode_phase_arr, + cmap='RdBu', + shading='gouraud', + vmin=final_mode_phase_arr.min(), + vmax=final_mode_phase_arr.max(), + ) + + plt.axis([wvl_min, wvl_max, diameters[0], diameters[-1]]) + plt.xlabel("Wavelength (um)") + plt.xticks([t for t in np.linspace(wvl_min, wvl_max, 3)]) + plt.ylabel("Pillar Diameter (um)") + plt.title("phase") + cbar=plt.colorbar() + cbar.set_ticks([t for t in np.arange(-3, 4, 1)]) + + + plt.subplots_adjust(wspace=0.5) + plt.savefig("unit_cell_param_sweep.pdf") + diff --git a/unit_cell/unit_cell_param_sweep_diameters.npy b/unit_cell/unit_cell_param_sweep_diameters.npy new file mode 100644 index 0000000000000000000000000000000000000000..b46d9806ba6a992787b3c72565a5e36e1f70c729 GIT binary patch literal 320 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh= zXCxM+0{I$7COVor3bhL41Fl&!L13r-=fpefy|wn+t4-y9To-%DKIo^xyX&e)?dxMf ztmF1;4i>X?Z98RuLlvaptUZ%!H19m+3-%JyZl zs(sOZkj!=ajzSRQhCNvH!A<)uS07k1x!$%v)ee%lV-Mm*^WL-n1~RVjzCHVSL5*8k f5A212JwgC(L$aQzI)WTrFMnJbR zasVqBZ!8$1=mfAJ_2XCi$FCyT$x5cx z6M&9PmXWxAOpAFYxF{HwEf7b=+B&+rw^!NVu$)@rIFt|oGcz=SX{5%-ldnJ3>=z5z zFTMN0^xF2D&$}t5H*8ji_xD<2#5dP&3`PO}?hM92c6Zv$ynRJNtnyn+1q)2Y+PT)P zl1;@)+jiP2)w~x({UKqunfBmx>E5sD5TWXQA%Z+DAzcl}F;~rfE*38>ZeNu+J8!Y{ zy1Ic=8kG3pm3E_ZB~vI@O?24Oxz(Lny>G9hyuv)kIsNWB62U5cEfgMW{`@)X~!lT#N7#zADoE#bWo}i;wQMnj|O{STgM!?%pu0 zNm{~eNmdMpC&z$8O-(X_Qx)qKZ(l0rk zltj+(Pd<>Fdx&_`q#{;J^K9)xH_L`jKR=$x#IxG0^Uw5-Q*=-y6up!C{1u<2(cZ%d zA4EAxq4eH4+vM99EcVKCbpIW-tG>^L(*s^%GY{6W<6_^~xtyY331`ZQ=Do2*q;r^7 zKk`?c$b2<@ep!TX;i2LOa`d-}%=NoI4|evFAOc_Z5)1HOW>M*-f!98{Jbq+$EfRiX zA5*&Qsa}3z1NF z-5FAuo~s*usTMa;^37i*nJ*9LRJujRm{q=Yldk8_t=z3u;=(LB?AtP}I{9MbufHYX zHV851Iw>XB8WuztG#a~;=9|O*FujTFtrfX@-Z>(J=Trf_&9@C>UUo8qe9m=MPiF|Q zYmdfhX$v`$Cgxa7Nb+9W##S*o1$`EyHw~HcN+I_g2lM)Gifn z*?u=;5MY97HTEQmsb(=XEANllw0x=d&U>v^0c^qM0awpTRk~n{D&j`$UnAqE>*8de z<^^$Jmj<)RHCL?|GGtZyhH5_K(2F~a=PQJO!rFueD&z2a_N*%iWOUij?RL4@Ji=hQ z6d@vO1r+;HS8D3Gvif?O`&0(?8K!*JjbcBM8x}6qHG#jRHF15_UneycEESw{@mRM^ zpK$7l62HS7d~w^+JpV&ckJa7NCBkBRLs56XRg`I+BXGV*m+``ks%53d|J}N~FB#R@ zHKn(kixw7JOkOD+p+S|51fP8N+ux5~=SrSjVkpL( zr3WbqN@1AxyRp+JoLa3U7L;h>Jae6RbWR?g%P2^WHdy$kojTpv{oIi6-H9}Q#9|D7 zL1e6#&eLt`gjF;0d>*(?h#)L*auD0`WLpE4Z+MGf@nw-hN=;~ZH|43 z8K@sgPzZF3<{{4q1H?h(LKXtey)AHxW<$e@|3JqWTTO9pmZoaL|~G_T6m&dehK)Q&eJr8ChCfb+eHova?>Of zL|vCyLl~Z@&L4-9hvHj}I6WIRoJW16rI;OwJ4NGJjcjL*Xpvw@D1YEGJIt~Bds`o{j#a%2qc`%&>@_uFt24y9|q__kCW!iIB-O`YaWTGQ&3_oYKwEwc%=#dROGH*xfI z%*VKD(G4Y&k3&|!JUDM=z&>_bN~9kyfsR%ch$dyab-6A}eE)Mv{l{h5;29l1UCAka z&14ZqZsykv&91fjhB8%dIllJg87lu9lMR2Wqx7^$Pzy!vy| zR&>tt`)ABfj?%LQ;VqAJiDvG@2`vCZl(1s_;K6$moB1412`_(fqJ| z@cGp>;nyPa*;;e4aS23RJhOw|6gNAq`T7+TC&f1<$qJ}btgpPfZ6Ol>ETv(>FVC~4 ze`PU9{$gs*;O$4vkE>Z1)0SKa#?OdQsPAZ|q|%a{S&7;ISV;6%Sfje2g4&$zlgtdo zrEEh%DzrTFE?>yFbmh~$TA8)W><35uAxT{ z`}y#=Qk5x%wjz|rMGdV!jw$9BZjRFnb>bg#24r^|iR?;^h(vC*-N18xNl@3O_=>;m zL$Jw1-5tBhrz+!N@;vEI%r$aD+0VP}o5N1|Y;V8l^a#zf9Q7wLp+qYN%aq)5>U@Yt zhmxCf(=%#{ptPpmPZ)d`<6fNv>n!>8g>h0WrB0UI=Z>m{N6TU~C z$xhZ6A3cWct+$@^vA`^xZ4azxpm;ekd zaVV5q#xoVt#R0E^H0w`#Wa1ewEp+M^h_@L4w;;T-T zzFjBuE*}eu4`GCKp~%~`3xq4%9>ym>V0C$c?GLLB)g&3{3g4W_ zRf;MZf*L+dkky+u(GQ-uFUQ&;r)Qxvt3hdZrmnl_BZe83|Nib;4EJdWtAot&bP~08 z;o?N!mC{^S7N_tuB`@<{WBbDT56LEM+ctul39o-E%W!rR?lC%$(;WD4!0qR_r~g zXV@c9mpxb%5GB=O)g;=YJCedFStbk9Q#xHu3m!Bcu6lV~oO9r#;$06al3`@kP5xTB zxJzCv`$TuQ0>hnh1WS{fE!=EP$ux?4HhEgK79|^=M`un7I%itSfDML`EDj*CK4>o0 zSgl>!;wu&W=uDl68Xx*6bW7GL65-Odc7|ruS{6L{7YEc8bm>8uV_lsB@0oo%SzcaL z!Cy*jceKr+EDRleJ9EX}O!3~RgxYYGK=DJT%a7&LF$3PiobQjBeQkM1RBcA|_8Z}4 zW;XP(+3GEC)Ie$4rQQoI=QlKYR;I$3E0@^sd5S739Oy}>O3bPMO#EiUFZp)=+0jrr)9#}%9wXlS~X>a_7zXWc8c zi|NG{yU-eOy||d)fY#MRcEhkyCMM?7_@$?oNP;-s7LeN{eI_60 zDMJ%q=Du~>r=xdPm9+`VXS*eMFIVLxcY2FEMxgGdyrei)b3*Nyd!t72rW`|{YVMxS z2TIxL&z}gQLxd=|O__eL(q z=*^(Mvq|2E6!tU+fr3!`p3R<(criv@~_8c0Tggn{*N=W zXq*TIJ`=IDpK2D%_zc+{j#TBy-u)^4SBmt!Wj+*Epc6J@35CztyjKPiiKcd|CL9EVD-t|Ijzl#NqbCCQpS|P-L4_2kOm*feUh;Tqe3mWwNKul#gc>2JwDth~PtpAx? z^k<3tnOAJJKi-e4tbW2EV^OHa2qZ3&~RtfOl{UL>+!yW@@TU`l^a>KU_^T)B<& zeu1t?H$+8=KBS7Z#pqT;ZQk?kC*EtMbu<3QHX!F0j^GjJxcY@hn7|*%PE;LD(FQ)V z>JyrI%nxtPjXg2w>*HQ~QB`wReeBHJ$e2g&>*Kd4*XfL9PtNg5si~fz=sjxuWf=vTL|$!`B^*g+?YQYqgJg*i7%v;>8{7e z_ghXh;JmP#N~0)5QhJllmCTjR-?^P=T=k`%e(Q7cTkc>uY-f~pfiP2VUg#1*zs~!& zP&IYMhbf@>Dc-wMp79*Ri{q*j?2 zP5(N-)mX9dHhr!Rxpkr3=Zd?4d_nGg;cAyp8KZl5-*9EW7JqA?M|*_8es-1q@D2(> z|9cuE5J+vP`T&ZM!a9r{wAe`q;tDgCXs(aXlj5h*Cw}%Eud+wE=Q^*YmKJL=UYr0N z6+C=~?jD{iq0BWxJ^EwnPE<`taC3yO5&;$ZZ+1jAl$C+oDB+V zkNF-TE1}Og#IYAn-c@dO0b=$ef2`*9uD4 zMXA6&%Z82>`3GX&QtiZ88XJv>^QSkc=gD;YX2G{3(_*ifZ82MSUG0(bYT3AKwyRo8 ze~FtlZODfMVffL{S3U4!x6&15Mlp(J%kW1eatk-6r6wi6unSs|T^vW$RtZH)r~3)I zmweWnJ-azTHn*%+92BJl<7V4-QjO7Y?m*^`4HR-0rs$=LuB^oL>2(wfsc|uShHhx~ z4kfyT_uSTN=s0haF{&HWubF8)Fb2Utz!Ls0d<4#M{M*_Q{BPl-)Vt#Gna=D{hEljv zt6YGtM>#lLSpE8#Cvb^v@J4R_#liT96gK$rO5sct3Nt9b|%Z% zmMR!(zI@hE{ljiPQ*V$^%rVgFh zbUse+F1e<&LWrNS9^On~#nKYJ#|iWLM#u5B;!)1#ec5=1z-0x*zmI^cZ}b^ z#_&AR5e8-=QPEU7PD$<^*ABlLbb+p%s)6R4^V!I;Gp|k05>8MifZT%nu0Lx>jeraW zKiAw`*3u`Csos9&x>V=eitMe1H8%7t#SQy!^U6`6o^$DM=tnL^~>d&!>$T%Zg z-;q(5Ao#3GJ6wIrB&XX|Iyo9Y?(FUW!R)3tTjDXoTJxRLu|@3R%#S?poXyDMBByFe ziywP1Otjgh8ebWe3Js{(=GW(u0j53}zYT50*x z*EVtW8(xHQcJdJlJpxa{QON(cB`O*s&K*$VDG$))8NAW44IfRa^Ox5iCuQ+rzT$nY|6lwpVRKiOQ<`mYM0w9WzE*4H2h2na+MtVwq(d zcohDQ^FkA(n|6_E+jw{l>28{ZXL-#C?xWtj7or?wJ>lk`E2ZzXa+Wtz<6N3e+{z~- zgY%kTk&l-c!$nS?8Lbcmr+3Jun(;6>NJd$nP6)N9LqB5i)liqB%Xq@NqEa+vYrk#m zptr@_*Ly<#MX*|odB@kU2V+lB57EKe5(iaTk0~;$Rpl`u3ny$H5W; zH5Vv)z{DhffsMB06i1@(lZtHeZagkJ&sPE@M$7Dc$JbKjCY&@gw%C<&FO+vzqz`@B zd5?)9u@HN?LYFaye|eMzc?6$FfMJ-xAE;Yay`H=cO1z39kehjE!?@(mb!UI-`UaYx zRpP`c&3BAvcxRK!p$#*6#rGR$8dGN+UIB&6zMx zH{eDP_MDr|h%Av&kSb_xG~s+mRk70aO`k^2)LiqE|tekIPbBD{aZV2?o6f`4G}5%5}dsCpR^ zk0R`BK4{VU2V{B>Q7^@xl0^`bgl|oE%uY#&fV}@0OmY$&DajjdfuKZ$3k%^hr?eyy z9%Gk2`}_!gj=;RgKLGhh^CK3X)Qz#lGheuF0G^~!m25}U4IQh`B;h=w<0G&mT;LBR zA*+E-IKf2R?Gt);hQV9S#YqLacYwtVBQn5A1j%JPSc}Oh%h@QZFY)Bh=S#K(TciEC zIzl)vbxD~^jVVZ@Zo_FZVJh*R4n2((uHxsy8C&vr>wLHz5-U+@)yF%x7(K*og~c7S z#NX6ArfJ3Xc5_fot81)9)mM~Pyw?i6u_WWmnXKO0jsAE8tR(sjaB;g|x%QFOHtfWU zEQ-1wk-+yaP7x1tiMLs)9yb|O?uI04$5WGz#gT+&SSfPYzB|$J+~*oo)7OXJpik3x z)aAsnnkt1(b|ZHyGS2%+o<$p#-QKzEcwte&Abs;mbY`=qJJUFH&Sv4t+5vs(`X1Sq znf3e;;yeQU!v4hLAPaaL+Q7szeyZj211bGJce-+2KfJaMed>m{Hvm!`5C!YE@>w!9 zzR4>dF!?ymDr7XD<2;vNex`>A@IhZW#}`s?nctJ1$HW(z?vR)Vw|vq={%ZN{r4Uvu zchgi}#d9W|rR|ic71Q+_L0)~<7DT%6bOn=lh`7e)-WcUR3x8W0@4n(YqgXAhTu{k~ zrtos}>ga_L$DPu{Rm7Gra{TwvB9B0u2>2iPS5TTt?j2B~sU<4odx%AsyBry+H1GFo zlN1`f$)ZA8U7sr@-|D)NZij2=?1UZ3^N*=sJxOxXfe@f+__tjdi8Ri*l4+A_&BUfF0?F5R%=YCfTVdUE6e-Sdx0 z(7C>dvhCM78Ex4T^G#uogcBuGiG`Ddp#Gfzcc{&L_TOV_1XA>+25YYwJ_SSLl_X-`jyjTjGn+--lq z&k<}LRLkTzTQ7uvl?;0FrI(&jNpyooq7>9-`pkSu^m->dFHOEAw@zj^;)}`4>gWUu zyt_OWf|4VeAtpvr88V?C19p1L*(m6!3!<)NUmSCJ`W(%jJSJlzv+zkqJ1%%7#c#he zs9!5KUwEbDOzw`#6wRqqZq43QIy}y3&8KLBdlP8=ho@Fvr$ZheFV;7rJM%?MN}y7AW!vl=>%5GTErcLzj=A+Pd6~jl~+JEN+)Orul|YRhU%oWT0K= zJ8Pj^uI|7Awqw|8D@12TIQaG{)lY*}74er0qgPdu7GGcB6q#sjpDv*}oJSP-^wS=y33YhUi{@4Kg+u+<`d|NPs+0xE*V&M94x*1`Az(LLX6}E{fWV zixPUrybwwU@@(TRK|p)T6t!j(dDhDG1tU|2luZUmxvtq%gz{_1d5LvOn4&2K_MaNF z+Ow5VaO9ftjl5Jz(7@>QoT|R)eOH8+;{8W=tkw96vZt~+c8U1u+O_moTVIa9=7o+O zp}iyE)1QE4gg_O5D4<9veu7+>L~{HCsYs@jtUz>}U3OeI0Xb6R8D1svb)*8RoU|1j z*v4lINp2s}`w@5scqR7V(v?*gA(dt#o(-omR zi!E2QCZ^{Yw}U^*4sD8WYQr<$wX*8SbMDB1dHby^)Ng;%8t0k=PnE2NM3a5Utu*9f zv7A$s=iQnLju=~``LIskPbQxvICaEGk3dQY#2?^kR6@Pe2^g`ghfWJ2ad=9MZzuHJ z__k$>syw zJDmKl*7LrUXh0to-nj%*&P?+0b2Vmp(lzQRS*|!+b@Pi-t{o=1v4c*~x01!bGn0=g z{)wl~@|98P+ZA0O`)X2rzbroAS()lf;?`@An{51J zOU0Ej{^3U`<_IJy_y@-aKzUjKq=p`+$j=Av1sv&xXqX?9{p|MmDW)TpP6qdo)N%G85!cY^eQKa`uFZu*-YUnG@n*qFA<_kDU@bKdf1(kkqPB`bAhXL zX^G24`Sf$mQO4Z>d)l;NbVmdtcN`Uxo0;IR6Ee<`P&yr1Yt2ERIuheSCC`TGmngZu zE%nZis$xNI$B^Qb!6 z^G(-CG(~6}P8P@ytol*QI`@mEy2TEvJU7skVoE-p_$iO8z1;qJbhL_giQv#%A?7-p z7+EdCYs~SRls>$lxRvxm=epMJTQc)h+_=c0wplK(`{)4AlHEV*cYUCafQzs{a10-s z$m?xX(j~}a=NIuck7)S_=z{=uod0~$$M-|i`=9{8fj6}_0ahsge2qb~vP0rr2Gznw z0-wo79#<;gdH&Z87Ad=)pN-=AuNuE`UWY5nw%0c4GQapBX%MjxP(;rnL2opqqpEN& zQDL{u#J<7bF1$J>fKU2^_oq-D4GN1uIijn}^vRtc$_n|`@}_3xiJkhXixYB+8uaL2 z3cZhfk>$fzLHN-*pTK&VQPuS9@#fZ_OXB!T{62R`w*&7QfCRIJF6O^6 z{P-Y^$=gfAkGm~Fr#@VXjZnE-Bgfiza0xsoB|?*WJ*~0LgN}50v8$fttCFQXS1xNI z*9vbOYuDC``^(YvsQvAQxW4;^i72^8#>-1Lm8D9}OU`kM1|-0G_f8}w?4J~!8$Ibn zoKikK^RkJ832#U7@&53{1D88_h{}$a>lxNXvE}%4qQxptxVUWYcZ?}q`RdcT^Vz5Q z0Pk{!(Gl-aQAa^U!9O@~sl+Kc%K|U7mx@{$>L~ii^!tkoNiXk+7taQa#FLSp;8+NZ zp{u+7%+-qA)DdCM`_=gwtv52}z{(h6p`u=EdwMi)+D)*sowiI_`-!eoefm@M9yTfz zEP3Vg>}g`=p%=U{X|X2cqK`$F#YXZPKk1562*s}qoz`_F5SkWc=g&aBkOdejL@y%Fr^}8@2|l00q@na&rFk1Bny|dF^R7tF&s4Fo zB1Q2Ik7?{^{1Z7QFQ^a+B}5{f0lF0>7fs*bDyfyCn)PIZVprJT({Gv4awgB|8Jb*% zm@Sy?eySnTeDIUMwBhEIO_IwU%#OguOU(vW-FuSPsY#z!ei8ht-QUGJSt1pdK+p@5 zYdz7o_4O^Q`x|3_NPKeqy&Ab=M)X!S7dlTIgT`lZN_^7JzQbnsbZCw4;_}jkJcxGE z@4|uoyv=}13_U5{tTIS9kILWztV5M!%-NCq5tIFm- z9$+v*5a3#N#>!wkEZyx~u`cdlz((Vu5g=Ng7TCkc4O}O{#R64xCtMiEuWf&?FZDg< zfe!|S!@zt<0R$K(h(-cY3kU=*@B!C_iaFNZ?l3Td9}2};aNzInlnDre<4zrj*8pNe zaCty5g{(6W-hhkd_$P|wFA@kMV`puR0b((5!52n{qw;Vyw*=xsfFrOgE=a=>V~zch zbhopy#R9h@+W{2>fa)?}h%7KLupodiRKP+&uj*hSK(?A-Az&xc2BQJVnggL7 zK!z1i8%P0UTQI~9NCDLy!4M}f1ZdX<41)nVK%YP$2n=Y<6S&a=M07y(?5wc10PRQQ zi?)l7v)y612O!~phl?C)?622=|4fOWbb6?O|Dix&5RT?H9)Jk{NR^Tvmbd^IG#v0( z;(EY6@WFAAT-w}K;k%TeAc+27Ng%ui2>h|LlytUn!~l##w6GW_U0g!N-1~5hfT)l^ zhu-`J`2np0BNG55_`kjgG!%^c{+j=HJ%B9%ga8Qm;(q`4_y5wz&ut3|LV>p6U=(14 zf2RMO76Rt}KhHnX=pW^Qu%dtRfVTdmU)WFWA%J#ZKlP13g5iRPGCFJvm&OTG0H6y1 z2J^4>kwCvFVA>o$19|^knfe){$$$@UPytI{7JPWS2uR6+;W%Rqq*TCwtM4$S4u<2- zp8!o0=>2<28w@`*8@M?S%#6caEAZj_X`mbi(Bya8_xS;c@UX<;Dg!LU_Xpg<1DFFI zK)^x)hI&Z00|SaXOgVxPKYDWlBYqeuXF&hoa}O8hp^Fme&;^V<>=1Ogz7M@X0QDcU zfc(+tA?1g4_%Y7EEX0rQ-En5)kR=Z=@<;dIEdpRe4%xyP3BYgjJ$-0`4v{a=-=ReT z1745AJe-jM-ZCGy2~?6PBCK- zHn&lRE2-IA`f3$N^_TGrA;QLY_VN)EvQ7!ps&w?bbp@|PcWgB6I4^sJZM)T%)W7Rr zzv`L05oZHW$RzG98JqX>%g&wtSgZB^$!mTlhi}UUY?!A^Up(r_#p-XZ-+bVwKBd59 zK28{VPhFcxJMas0nsjU2srlof5?t!1Mvg~cM{-FeQgl`%Qz|{+&WkjlB_?L4C%zAg ze&~IZ`(?J2SBS(yV5>P%Ah+ZTT53Eh7S6gAZq1>)TAZ|%FM-Tha>nC)>tn4l~Vo+?sW5^%g6$@4@JG zFwxWrc&Z>}2PLi7eRFSV%N4uXNM)@=>5oDA$X3>CX=h^5?d5UDL$-A)(`BDpeBFeP zN?rur_gpsUga-Xw6g$TLGZ(1{4idy zEf(wQAq;^y{hi0};%>tQf&h5j%F_~vH~zVrtCcm_!ramU=;p5`fTD*T*ts~%;DWO` zWrX2SxByfLIJUu|NF;&>3g>`AIk&p$P)_kbk8K0j>W=!|lNTstY^-Z+3sB0hcO& zr6J(}zW=q{5i}%B=y(2vgnl0{6o+MSpPys^2N312G$amp{JR_i@p~TtIR33YxWEy0 z1#!sb-|Zn`s9)s=MZ$5=?_cErjSKub4k%LKcbOsufomAQ)m)5cYdrK_m{b{+mBR0o;l97n%?br2axf14qPP zXhJaPZ|w;of3GVfAov^4Xefa8e`ya8->>ruibmm}{l9qu4hp~5MFU6Mf0qMJu)sn7 zR~q6sox`Alzppo_&~G||!Qe;GaNzcrJ^%;8Z}S$rgrcn?y6?y#Lqri}NC`zEV-cZDC6OkCjL8r(M)$We zB_~2;p66-XX4|%Zo`2z4@7DEVy;*BrpDX6fc>^PJJkAq$SNf{m?JI85s^ZdW*N#cc ziA!H|adUUOeAC*+?W*1X>Hl7Kx^1_m-*&is%WkWeQKHF-D;?S|C+;Wy|17S4nz{BD zyJHb}+I}OsCgdI4WiEz)# zg2Y}sHVW`k3IvB6(fE%c$@4E)Kyb~?^G>ByERmruARA z8D+^{>}c8c8-xoVwKTkpgK*oodHP-=_P|ACap=(ma%Y9rak$NZZ_fr-jti4uj?jhO z8;zUD#AI=S8B~w{6iv3{53T}na51RRiHgm6-pm~}UqTNo(+PK<{{jE(GNmn1C($y8 zFUnn+2Bmi9O?;a~c=6Kip!?ixq|lE2a5is1W?9$l4=k+!v4Lqhr7$XXU##gPK6(iO zxu+wq%G!a}EqA742@clQaX6SX5U_6VcND8b<0z%uxGBDA7MPq3%y4igK>_{YW6wl4 zQKs)tf-1cV@tUvvdNjBSc8&xu+455{g-gfxE2J(Wl_n*F?QO% zlv(gB`0P-+BMI)+zM0XZw24ZC-(B(xtwv!7-rcKhS_S-H3LH`dshGiY2?6WEMbsKn z#-*d$1=1dMi~K8thebBJTzM7*EO%cuHTBvkIuYTT`t1xI9Jz7Y#&w(oCG-0m25Br5 z|6FlQ?7$~P6FPq{ym1Yv2JmO|#8I&HJhImB!Ua^ZPi(0{r4Q6w#Uk5dqexg-Ikd)} z4oXODddCt8-gpwq4Fy=}tu-!1QL_vwdx!Q;%&mc-o*(?vVH7N@-q9i9$pSk0!u)q} zb06S*mhKx|frpWb3h8Rq92oy=1;q{SVH6k1BuCzx2kowuhq3p`FlpT4>@M5}s=z0U zQ`_F43_aI5Z%!tl)W{UoZjiC)>vhr@?hG`?;^hq&9|Vhud9v{qLug1WDORL&9_-#5 zn)K}h8Ok(i6py`OBEbfA1(J>z^4YGc>2Zq*F4Tw)G}e-_Bljd(Sri5;`L&Zb<|89_0Xv) z0I7!k%dSNmpy;DRM(Vz;x!uCj6pW6VZ{>+aM~{HcJk08mWgmLl8*)vzXaStOXgg0z zra)sxKj(7k8d5e{`Z%Fc3BnD6}h>&o49n zS;MJKz~ym+`#70|@wedg)mEp`Wc2z-rOYI_H*H#xU-uiy9@`6(YgZ=ojG?7 ztl^-D(x($0cL~_?M-!{SqEYmyr(Wfa1|8^%SY(aSzM!R&>SXA*1|m*rb?_BZu*@Bs zJpW+}=)7{o&zrG*!1UKQ!oeCmY!7_d-~NsRi*qdh)_#5n1@-b6S8dLNlf~v8g{@@x zS*I#^T#1R21g+HroZkTN3vy3u68*>QlvXvcWzS@%|?9>m{$)fG@(wlR5}Z-fKoTs3n@P;rmbh> z;V-*{%u0+Rl;3m$slVSBEvw_<=&OTGae@R4|94i-K*=!D-Vt|X%6}eY_#}p@3 zervjd4ignB1VtD4r+}kuV{Pf64e(*Nm*-cNAt1WL*Zz1l2h^}l(L4DNkGY2@MEcou zqnWu#Tut91_{vRFPIsol!6tmWdCm%IyQR#OHTwnhRu2@pE8t*Is%hdV9})XnwZje* zA4gaCIX~DC=-|zv;J!K@GHmUX%t|oaK&P7@T4xD{pz7y2t?D>i4fts=W?lj#{w9tu*>qLQlSIWdBR$C6+LhHV2lYjf$rC9 z>Jw5p7!nvRe0YU`-Jzt}>m46M3jVbg;-z%3+N!_HHH-}H7V9c4cQMiHGwYs4jE`Vg z_^6_7DGSsXMr3zv?PK1he*I_f6v8qIbYtEbki5To-$(!nI%~$*y`!*D-zDv#%7P@+ zWqw(+l+6Suv>)bi;mDYcSh~7TJPv9IdCnAECt!E~m`!8(!{~Zs5W!H20n)#=Xe^sh z;4>9go$jSIRCDfGo73szb9_bq%9_=OI#CW3tS zrjcQcK~kDS(K?!|(^Ksetpvu3CwvyoHo+act4%GxiCEZ#((ymWljwBYali8HSztT7 zb|!0r1ictfzx-0&Kn1-dp4C?rC_&7y``)&S5|c!)6iWUFoLTi*d-!m0!SJD$!5RUJ zs1({Gx-p8Prb>;|7w90k?>mQNJ{j6SZe(efucMjipE(~o$^p%1WVCz7CYX4xq4b@L zgpC~L{59q}h14T+MoS+`Dc zZ!+c#aTOPG=_s{@9#o`E!u&>I%e@2B$YaiULBMAkB#V=S%e6=_#IHk=Kc9sh?W|G_ zqKlD1x6M_nu5}RY>YQ6pM8?!{rh9zk=aGH+RqjuMqhO%W%JuzQ0vwCuoqq8UhiROO zcrj+*g02e||BN+T0l#U+8twjDoE_`9*8X!5-7ER2Ap4{T$ntAq7yrXUhN`H0t|D`4w(SKWi)x^?anxZ zKR_&>pZHZA4=qing>ASASgQ)fNW69kb+wR7E}Akx?UdU+gEJH;(=fil8@Yy}`aUtH zr)t5gD9=E-B{s-WGSbzOAYi+)VycR#oIs_Glud|E4Y0PRZ9lz?hXa9ijEV?6MlzC3 zRD9Wi^8NG<9TV9=Yhh-z*i>&|{A?q%!E6Jd+Yg@V(@9uN{L?GAky+$fr*pL4eG>G? z%4AGJ5j-I<_3Kq&Sz>YCN^2ry|Kwi&z@HKPmUs)E8lBd>?bIJV}y~TPk zFlT}g}e=t6XVEuyZ)YM>T{s2>=A2sD+vm*)?ylZHc;%t)uumwqF^F?iPFcM7ZA?linpZP03eO& z$1zTGz=+LwQ!^GGyQz*F+RxRG4)w*Wn$9eMgs9g~BW5X3y_gR>T)B$Mza(Yt@ca&* zmMVAkRxTqhv29N|ue5=W!PkBZOyJ=4X`!A_Z2}f@ai%_Pau{LhQ^||3=Rq>M^NRh8 z47Kh!+b8y~qkT=LXHNF#gO@giN!0=@V5BZH)b@#low(uEh$qjX3Stx%-a8JAJW6q& z4-%owNH1ek7>Ch3(MeI-FC@QW7g73V1-#^n6F9temp;)WUrhg6LOHW>&O_(A!0pFO zuV@+ri75{}I9fCW`15D~TA$~Dcd|EXKN#XM+oRi`{JztP?(V*HHq&YewDoMUEt?8! zcvIbjj;x@9OA2aZo6W!`&ET=31rB~%&Nt+u6R^?j#}%EOV@P~7{$m%B4xUuW`=&6r z^7C9%xqre2;Xn@xGOi zl#EzF^_J`(_N;!eSM+dbzYYgH<#c*a*aRNqE1&i2*6BryysPC!(~BTH8HJoAbHaF| zKRajY8K~5I@P76=94rFL?w*8Ike%wkrpZsibW+t*`}OCMZjeaNto$e#i+aEmUKmC@ z>Q3#+X`TnOkNA$t&?r#eC%M1x?JC;;iRZki=O{?HKF=OsTL3{h3`V**6-%(FOwmzT zKJK5NUbeM-Xt2mdyyUxV2mjP6m;}i7f zVay5|xRr~3G#Dl|N#dZA-r*mLTR!hXt%`zQ#*tvkLbIaK9B_ObpPQsJh%zEG z%p`^wVA`;;|CK8RN)2evXl1P-g__pS|5|+o)i*>0yiRO_7hOZV@8w8XiK)PZRNFKX z?Pyo%B1{4Q$KPGVN{H}5l(&zUJR5!Ezo7GoRF8IYe@wfbvkK0=B#+V;DcD8x5DUYl z1vKC?d|jTqAE=f(L>ZGfpoJ6v?ux@4m>k~ibC5$nstmu7PE=b2&GUDi&CJ`-rk8oS zi^Vc%BDP*h6y?Oa3x=ee43^PQ*r2URCAX9S6HN>}>xDAYkgx39oz3ci72)7$j&HlvlRbbrI1* q20M!l27vV`eBk2%4yajM>-n|-j|oeKb#>nCLH;`r_`fyxg`)@sAF-faOGGJLN?aX75Otrnly-_y7&`TVn=zMy!B<1R<>C}B zypO)(G>^h@S=+8M2?`>a{TD@0J}}y;!;hHex)Y`lyN}!xv_?JpZ+J)Ge8irimtdmK2hf#IyuhqyEb;0hXKbymI4WMQ2d)d*}GL+}8YvWj>oVuibtU~#8i_cpRl$$&Gl?=F2 zDD5D}9o(@l_2h#aJt?$SuBkWnrtr6&kAJBzg+~?US|a`wjOX@P1qD*rzRaZcQ83o$ z^152_z3A6o|8raRQ{cbP8%Q}wK}zvXL31dDAGNOoK8Im^*JPb}8ctzDo|(BqB%)Bw zhu~-mQI40T15y69+Tq3%l)t?4UH?AH@8lhRJ&5u<>%CSD#8C(p)1G*qKw+X{&;IL4 z6qGW={gaa^?Ejw2Y&%6k)V6lB*cl30$8+=Vq*9n#yua8vgMwwusp|eL3ULi?>f3TC zEOG4DtvyenMxOut{5%SQ^Gho?6i|5Q|G7Z7kixw1>_1UO6x!Lx?uL|5ICrJONE7$n z3@whUU#2j;^r~~x6@2ga#v-`aD7>^gjGraFZ-kbZ+Da`-w zY;e38{U|W7pZ9`7$8d~;V=KPD>q}*-UQ*yY8x*&vox-v!?cTdyQ%IFlGEeNLV4`v? zJL^3Kj~9Qr#`IF?ojP2m+)p9GS4$#mkV4ta(aX<1Q)q0;m#FzdffXsJ>@|Wsu`5{m z0ULS1aoDn*gS_?U)X!8d=D*6OZ3^-BH1iu-*#zHReoHXeq8&T!6?GV)*iLe_?Caqx4^x z$JqCr(MhcPT8A%1Kd=ritK{GKj(+Q3R+;xN1)geWaP=4kf3K#|G(1P-x1Z!yyielQ z+dQTR@6qlhOyklqzTOj}77myP!})OyGS?|=y1YY(!F&(gUuQI! zPT{9R7W*Cg8z|n=@dx82(`T&EhITBv``l_V{^Ii$N}(RlYk2P;f%#Qj_gF9l>q9Oh zNlH1B!n!-n)v1Zd``7k1-Uy-apy)@}c`pj_2c3#{Y^D$`X?KOIi*cAU6cvKJ^}JAh zjrbf2_uKV+zc2&>x^o`!cN3_Q*sPy;pTOgm#8BBH0=uoA{_UGWz(S>2qYU}UY{ARy zzJ(N|^AfIiU?2PXaaoS(CzL;!dsm=|z(9eja7;OYPLp08p$r0}p{(5TBLtQ|jP}gg zMc}xN!!|oB0>95HN*q~B;AvsMQPdIwza&yT%Z1Qh65Sgx!9Zb%7v&pe;EQW8t7(9N z)!B^GtM?3CvJ|V*?`B}iKvVt-;@yjqxko-QV7}p@`a?wmp~)4kiu2HJjo~WKKMV{d z4>c8jWneAs{wl{|AbOvT*!<57+!(s4>OR=`sR4*Rw_& zvh**2>cAjzO6na7%saA`MsN1aA*&rKA1L=Kt#-aRD_D?qo8VPHLR1jIid*CfV=apBLP28jq~inuQ8NvPeyr1_TE9vdc^*{0UCrlr z;SKbsOJ08&Gz=SDmEgD}1LL!n{_yc3ut8Td$FiQl-rvb5WCAG&8|X31_aL8iDmOg| zKtFa{-PU%&{yK4j&98(0N3OMh*hk=@=%xoQl!0YUUoyrTxu6iKYvvyr2}i0ciWl)m2kDpJN^V_mvhXLn+a^UtX#3@8`hz~&y}VE$iIVXf{R5dtXU)CPK1%y z)jtHp-xP0@LqFZ~bOq*Np5w=4zFrikGoD35_@@f48(MciKo~!6cRZ0z(UL-33=jmWAOPv0s{N^M|^Dv z#6DNODT#dH3JVtGvHRIw}e&4fa&?$xBUtkU;JHaQ^kAY+Zbv|1k27ay5D%1~RV3YfSlkrgua85n;ciw{a z;q6(~s7k=wT{0{BfqwatiEv0oP)3iVCTVGwXSU#)M7PwX#s3V4(jpGUCx06rXE>~@ZNGf_Qu481i11K zZ*GbqptbS6cg+C`lS-#N)3kUHF>6@5J+2oH#hf*~AT|PO(xJ}ApV;vClw@VZSq_{{ zK4fFm!vPCztCNo>Il!u~ZQeYK3)un+%dXCFAmy%uh2Rekl*9^(ICD5K93p(g^*sj) z)*thndCmc^;Do<UM~bH> z*%18U>5_!|Z19yiH_DA=!$Jv7g(s$Lc;o(cx8KAF(6_ovyDyJGmZQRCsNM(+7kMMQ+2Fq9lnDE?s#$E7$mSy_Douhf^hteD7Qi`i0!)) YpY6bf=@;LE%vN*ZW^1YITNy6=4>KIs$p8QV literal 0 HcmV?d00001 diff --git a/unit_cell/unit_cell_param_sweep_mode_wvl.npy b/unit_cell/unit_cell_param_sweep_mode_wvl.npy new file mode 100644 index 0000000000000000000000000000000000000000..740c4035b6389d7716757c269b72bea3c963c050 GIT binary patch literal 288 zcmbR27wQ`j$;eQ~P_3SlTAW;@Zl$1ZlV+i=qoAIaUsO_*m=~X4l#&V(cT3DEP6dh= zXCxM+0{I$720EHL3bhL411?Q|&FqRgS24W$XX67Ya(6Wh?p9{xipm8D`8s?N9nnk`H_P!=7{fypyx9|FEBU zbA>Vcp&#~=YnS|