From 45492233d4a5e0ae3cdad387d3403d1eb2f87d33 Mon Sep 17 00:00:00 2001 From: Edoardo Vergerio Date: Tue, 8 Sep 2026 16:56:56 +0200 Subject: [PATCH] implemented and did some testing about applying a conic filter on phase masks. This should smooth out outliers at the cost of a possible reduction of total overlap efficiency. Several example devices have been added for testing purposes. --- filters.py | 86 ++++----- mpid.py | 172 ++++++++++++------ ...g_converter.py => single_converter_symm.py | 13 +- testing_filtering.py => sorter_to_gauss.py | 1 - sorter_to_ident.py | 86 +++++++++ 5 files changed, 250 insertions(+), 108 deletions(-) rename testing_converter.py => single_converter_symm.py (84%) rename testing_filtering.py => sorter_to_gauss.py (99%) create mode 100644 sorter_to_ident.py diff --git a/filters.py b/filters.py index 4083713..1da675f 100644 --- a/filters.py +++ b/filters.py @@ -1,9 +1,9 @@ -import autograd.numpy as np +import autograd.numpy as npa from numpy.typing import ArrayLike as ArrayLikeType import matplotlib.pyplot as plt filter_radius_um = 1 -def _centered(arr: np.ndarray, newshape: ArrayLikeType) -> np.ndarray: +def _centered(arr: npa.ndarray, newshape: ArrayLikeType) -> npa.ndarray: """Formats the output of an FFT to center the zero-frequency component. A helper function borrowed from SciPy: @@ -17,15 +17,15 @@ def _centered(arr: np.ndarray, newshape: ArrayLikeType) -> np.ndarray: Returns: The input array with the zero-frequency component as the central element. """ - newshape = np.asarray(newshape) - currshape = np.array(arr.shape) + 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))] return arr[tuple(myslice)] -def _edge_pad(arr: np.ndarray, pad: np.ndarray) -> np.ndarray: +def _edge_pad(arr: npa.ndarray, pad: npa.ndarray) -> npa.ndarray: """Border-pads the edges of an array. Used to preprocess the design weights prior to convolution with the filter. @@ -42,36 +42,36 @@ def _edge_pad(arr: np.ndarray, pad: np.ndarray) -> np.ndarray: A 2d array with border padding. """ # fill sides - left = np.tile(arr[0, :], (pad[0][0], 1)) - right = np.tile(arr[-1, :], (pad[0][1], 1)) - top = np.tile(arr[:, 0], (pad[1][0], 1)).transpose() - bottom = np.tile(arr[:, -1], (pad[1][1], 1)).transpose() + 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 = np.tile(arr[0, 0], (pad[0][0], pad[1][0])) - top_right = np.tile(arr[-1, 0], (pad[0][1], pad[1][0])) - bottom_left = np.tile(arr[0, -1], (pad[0][0], pad[1][1])) - bottom_right = np.tile(arr[-1, -1], (pad[0][1], pad[1][1])) + 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 np.concatenate( + return npa.concatenate( ( - np.concatenate((top_left, top, top_right)), - np.concatenate((left, arr, right)), - np.concatenate((bottom_left, bottom, bottom_right)), + 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 np.concatenate((top, arr, bottom), axis=1) + 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 np.concatenate((left, arr, right), axis=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: np.ndarray, pad_to: np.ndarray) -> np.ndarray: +def _quarter_to_full_kernel(arr: npa.ndarray, pad_to: npa.ndarray) -> npa.ndarray: """Constructs the full kernel from its nonnegative quadrant. Args: @@ -83,24 +83,24 @@ def _quarter_to_full_kernel(arr: np.ndarray, pad_to: np.ndarray) -> np.ndarray: Returns: The complete kernel. """ - pad_size = pad_to - 2 * np.array(arr.shape) + 1 + pad_size = pad_to - 2 * npa.array(arr.shape) + 1 - top = np.zeros((pad_size[0], arr.shape[1])) - bottom = np.zeros((pad_size[0], arr.shape[1] - 1)) - middle = np.zeros((pad_to[0], pad_size[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 = np.flipud(arr[1:, :]) - bottom_left = np.fliplr(arr[:, 1:]) - bottom_right = np.flipud( - np.fliplr(arr[1:, 1:]) + 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 - return np.concatenate( + return npa.concatenate( ( - np.concatenate((top_left, top, top_right)), + npa.concatenate((top_left, top, top_right)), middle, - np.concatenate((bottom_left, bottom, bottom_right)), + npa.concatenate((bottom_left, bottom, bottom_right)), ), axis=1, ) @@ -113,34 +113,36 @@ def conic_filter(phase_mask,radius,space_grid): 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() - xv = np.arange(0, Nx//2-1) * (space_grid.get_resolution()) - yv = np.arange(0, Ny//2-1) * (space_grid.get_resolution()) - Xv, Yv = np.meshgrid(xv, yv) - - h = np.where(Xv**2 + Yv**2 < radius**2, (1 - np.sqrt(abs(Xv**2 + Yv**2)) / radius),0) + #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 = _quarter_to_full_kernel(h, 3 * np.array([Nx, Ny])) + h = npa.where(Xv**2 + Yv**2 < radius**2, (1 - npa.sqrt(abs(Xv**2 + Yv**2)) / radius),0) - h = h / np.sum(h) + h = _quarter_to_full_kernel(h, 3 * npa.array([Nx, Ny])) + + h = h / npa.sum(h) #make a plot of the filter kernel - H = np.fft.fft2(h) + 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 = np.fft.fft2(phase_mask_padded) + phase_mask_fft = npa.fft.fft2(phase_mask_padded) Y = phase_mask_fft * H y = _centered( - np.real(np.fft.ifft2(Y)), + npa.real(npa.fft.ifft2(Y)), (Nx, Ny) ) diff --git a/mpid.py b/mpid.py index af7970b..a74291a 100644 --- a/mpid.py +++ b/mpid.py @@ -5,6 +5,7 @@ import os import pyfftw import nlopt from autograd import tensor_jacobian_product, grad as gd +import autograd.numpy as npa from filters import conic_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: @@ -62,6 +63,8 @@ def make_1Dplot_of(x_axis,y_axis,plot_zoom_x,save_name="plot_1D",) -> None: plt.close() class SpaceGrid: + ''' + ''' def __init__(self, side_pixels: int, unit_cell_pitch: float, prop_distances: np.ndarray, materials_refr_index: np.ndarray, xy_symmetry=False): self.Nx = side_pixels#for the moment we consider only square grids self.Ny = side_pixels @@ -346,7 +349,7 @@ class Simulation: pyfftw.interfaces.cache.disable() pyfftw.config.PLANNER_EFFORT = planner_effort - def set_optimization(self, opt_param_min, opt_param_max, ms_max_avail_diameter, phase_mask_init_size, opt_max_eval, conic_filter_radius, results_while_solving=True, custom_mask = None): + def set_optimization(self, opt_param_min, opt_param_max, ms_max_avail_diameter, phase_mask_init_size, opt_max_eval, conic_filter_radius = None, results_while_solving=True, custom_mask = None): ''' Sets the optimization parameters for the simulation. Args: @@ -362,24 +365,33 @@ class Simulation: self.opt_max_eval = opt_max_eval self.results_while_solving = results_while_solving self.custom_mask = custom_mask #future: one might be able to directly pass a custom mask - self.filter_radius = conic_filter_radius - 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():]) - self.filtered_phase_mask = np.concatenate((filtered_phase_mask_1, filtered_phase_mask_2)) + if conic_filter_radius is None: + self.use_filtering = False + else: + self.use_filtering = True + self.filter_radius = None + + if self.use_filtering: + self.filter_radius = conic_filter_radius + 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():]) + self.filtered_phase_mask = np.concatenate((filtered_phase_mask_1, filtered_phase_mask_2)) + else: + self.filtered_phase_mask = self.phase_mask mask_upper_bounds = self.norm_phase_max * np.where( self.space_grid.get_X().flatten()**2 + self.space_grid.get_Y().flatten()**2 < (self.ms_max_avail_diameter/2)**2, 1, 0) - mask_upper_bounds = np.concatenate((mask_upper_bounds, mask_upper_bounds),axis=0) - mask_lower_bounds = np.zeros_like(mask_upper_bounds) + mask_upper_bounds_full = np.concatenate((mask_upper_bounds, mask_upper_bounds),axis=0) + mask_lower_bounds_full = np.zeros_like(mask_upper_bounds_full) if self.space_grid.xy_symmetry: mask_upper_bounds_symm = mask_upper_bounds.reshape((self.space_grid.get_Nx(),self.space_grid.get_Ny()))[:self.space_grid.get_nx()+1, :self.space_grid.get_ny()+1].flatten() - mask_upper_bounds_symm = np.concatenate((mask_upper_bounds_symm, mask_upper_bounds_symm),axis=0) - mask_lower_bounds_symm = np.zeros_like(mask_upper_bounds_symm) - self.chosen_mask_upper_bounds = mask_upper_bounds_symm - self.chosen_mask_lower_bounds = mask_lower_bounds_symm + mask_upper_bounds_symm_full = np.concatenate((mask_upper_bounds_symm, mask_upper_bounds_symm),axis=0) + mask_lower_bounds_symm_full = np.zeros_like(mask_upper_bounds_symm_full) + self.chosen_mask_upper_bounds = mask_upper_bounds_symm_full + self.chosen_mask_lower_bounds = mask_lower_bounds_symm_full else: - self.chosen_mask_upper_bounds = mask_upper_bounds - self.chosen_mask_lower_bounds = mask_lower_bounds + self.chosen_mask_upper_bounds = mask_upper_bounds_full + self.chosen_mask_lower_bounds = mask_lower_bounds_full def get_sources(self): @@ -391,6 +403,9 @@ class Simulation: def get_space_grid(self): return self.space_grid + def get_ms_max_avail_diameter(self): + return self.ms_max_avail_diameter + def get_central_wavelength(self): return self.central_wavelength @@ -453,26 +468,37 @@ class Simulation: ''' Initializes the phase masks for the simulation. This method should be called once, before iterating over the targets. - If multiple targets are present in this simulation instance, the first initialization uses the first target + If multiple source/targets are present in this simulation instance, the mask + is locally initialized around each source and target. ''' - source_to_initialize = self.sources[0] - target_to_initialize = self.targets[0] + + source_Efields = [source.get_field(flatten=True) for source in self.sources] + source_centers = [source.get_center() for source in self.sources] - source_Efield = source_to_initialize.get_field(flatten=True) - source_center = source_to_initialize.get_center() - target_Efield = target_to_initialize.get_field(flatten=True) - target_center = target_to_initialize.get_center() + target_Efields = [target.get_field(flatten=True) for target in self.targets] + target_centers = [target.get_center() for target in self.targets] self.w_norm = np.zeros((2*self.space_grid.get_S())) - circular_mask_target = self._build_circular_mask(target_center, self.phase_mask_init_size) - circular_mask_source = self._build_circular_mask(source_center, self.phase_mask_init_size) + circular_mask_source = [self._build_circular_mask(source_center, self.phase_mask_init_size) for source_center in source_centers] + circular_mask_target = [self._build_circular_mask(target_center, self.phase_mask_init_size) for target_center in target_centers] - source_pattern_phase = (np.angle(source_Efield) + 2 * np.pi) % (2 * np.pi) * circular_mask_source #make sure the phase is between 0 and 2pi - target_pattern_phase = (np.angle(target_Efield) + 2 * np.pi) % (2 * np.pi) * circular_mask_target #make sure the phase is between 0 and 2pi - - self.w_norm[:self.space_grid.get_S()] += source_pattern_phase / (2 * np.pi) #normalize the source phase to be between 0 and 1 - self.w_norm[self.space_grid.get_S():] += target_pattern_phase / (2 * np.pi) #normalize the target phase to be between 0 and 1 + source_pattern_phase = [(np.angle(source_Efields[i]) + 2 * np.pi) % (2 * np.pi) * circular_mask_source[i] for i in range(len(source_Efields))] #make sure the phase is between 0 and 2pi + target_pattern_phase = [(np.angle(target_Efields[i]) + 2 * np.pi) % (2 * np.pi) * circular_mask_target[i] for i in range(len(target_Efields))] #make sure the phase is between 0 and 2pi + + #initialize mask on source side + cumulated_source_phase = np.zeros_like(self.w_norm[:self.space_grid.get_S()]) + for _, source_phase in enumerate(source_pattern_phase): + cumulated_source_phase += source_phase / (2 * np.pi) #normalize the source phase to be between 0 and 1 + + self.w_norm[:self.space_grid.get_S()] = ((cumulated_source_phase + 2*np.pi) % (2*np.pi)) / (2*np.pi) + + #initialize mask on target side + cumulated_target_phase = np.zeros_like(self.w_norm[self.space_grid.get_S():]) + for _, target_phase in enumerate(target_pattern_phase): + cumulated_target_phase += target_phase / (2 * np.pi) #normalize the target phase to be between 0 and 1 + + self.w_norm[self.space_grid.get_S():] = ((cumulated_target_phase + 2*np.pi) % (2*np.pi)) / (2*np.pi) if self.space_grid.xy_symmetry: w_norm_symm_1 = self._symmetrize(self.w_norm[:self.space_grid.get_S()]) @@ -481,25 +507,36 @@ class Simulation: self.w_norm_symm[:] = np.concatenate((w_norm_symm_1, w_norm_symm_2), axis=0) self.phase_mask[:] = self._phase_given_w(self.w_norm) - - #make a plot of the initialized phase mask - plt.imshow(self.phase_mask[:self.space_grid.get_S()].reshape((self.space_grid.get_Nx(), self.space_grid.get_Ny())), cmap='gray') - plt.title("Initialized Phase Mask") - plt.colorbar() - plt.savefig("results/initialized_phase_mask_1.pdf") - plt.imshow(self.phase_mask[self.space_grid.get_S():].reshape((self.space_grid.get_Nx(), self.space_grid.get_Ny())), cmap='gray') - plt.title("Initialized Phase Mask") - plt.colorbar() - plt.savefig("results/initialized_phase_mask_2.pdf") - def _phase_filter(self, phase_mask): - '''Applies a conic filter on the phase mask''' - phase_mask_filtered = conic_filter( - phase_mask, + 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. + Note that due to the circular behavior of phase, a direct averaging of phase values can lead to incorrect results. + The following code averages the sine and cosine components separately and then reconstructs the phase mask through the arctangent. + A practical example might help: if the phase mask contains values close to 0 and 2π, a direct averaging would incorrectly yield a value around π. + By averaging the sine and cosine components separately, the correct average near 0 (or 2π) is obtained. + + Args: + phase_mask (ndarray): The input phase mask to be filtered. + + Returns: + ndarray: The filtered phase mask. + ''' + + 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, ) - return phase_mask_filtered.flatten() + + phase_mask_filtered_sin = conic_filter( + phase_mask_sin, + self.filter_radius, + self.space_grid, + ) + return npa.arctan2(phase_mask_filtered_sin, phase_mask_filtered_cos) def _forward_propagate(self, source_object: Source): @@ -533,6 +570,11 @@ class Simulation: def _reduce_symmetrized_gradient(self, g_full): + ''' + Utility function that takes a spatial full size gradient and reduces it + to a smaller, symmetrized representation. In other words, it combines + the contributions from the other quadrants into the selected quadrant. + ''' g_2d = g_full.reshape((self.space_grid.get_Nx(), self.space_grid.get_Ny())) third = g_2d[self.space_grid.get_nx()+1:,:self.space_grid.get_ny()+1] @@ -612,12 +654,15 @@ class Simulation: self.w_norm[:] = all_ms self.phase_mask[:] = self._phase_given_w(self.w_norm) - + phase_mask_1 = self.phase_mask[:self.space_grid.get_S()] + phase_mask_2 = self.phase_mask[self.space_grid.get_S():] + if self.use_filtering: + filtered_phase_mask_1 = self._phase_filter(phase_mask_1) + filtered_phase_mask_2 = self._phase_filter(phase_mask_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():]) - - self.filtered_phase_mask[:] = np.concatenate((filtered_phase_mask_1, filtered_phase_mask_2)) + self.filtered_phase_mask[:] = npa.concatenate((filtered_phase_mask_1, filtered_phase_mask_2)) + else: + self.filtered_phase_mask[:] = self.phase_mask #after the phase mask is set, we test the same mask across all targets C_all_wavelengths = [] @@ -625,20 +670,19 @@ class Simulation: for target in self.targets: - C, full_grad1, full_grad2 = self._value_and_grad(target) + C, full_grad_1, full_grad_2 = self._value_and_grad(target) C_all_wavelengths.append(C) - if self.space_grid.xy_symmetry: - grad_partial_1 = self._reduce_symmetrized_gradient(full_grad1) - grad_partial_2 = self._reduce_symmetrized_gradient(full_grad2) - else: - grad_partial_1 = full_grad1 - grad_partial_2 = full_grad2 + jac_grad_1 = tensor_jacobian_product(self._phase_filter)(phase_mask_1, full_grad_1) + jac_grad_2 = tensor_jacobian_product(self._phase_filter)(phase_mask_2, full_grad_2) - jac_grad_1 = tensor_jacobian_product(self._phase_filter)(self.phase_mask[:self.space_grid.get_S()], grad_partial_1) - jac_grad_2 = tensor_jacobian_product(self._phase_filter)(self.phase_mask[self.space_grid.get_S():], grad_partial_2) - jac_grad = np.concatenate((jac_grad_1, jac_grad_2),axis=0) + if self.space_grid.xy_symmetry: + jac_grad_1_symm = self._reduce_symmetrized_gradient(jac_grad_1) + jac_grad_2_symm = self._reduce_symmetrized_gradient(jac_grad_2) + jac_grad = np.concatenate((jac_grad_1_symm, jac_grad_2_symm), axis=0) + else: + jac_grad = np.concatenate((jac_grad_1, jac_grad_2),axis=0) grad_partial[:] += jac_grad @@ -683,15 +727,23 @@ class Simulation: output_field.append(self.intermediate_fields[-1]) np.save(f"{self.results_folder}/output_fields_{i}.npy", output_field[-1]) - make_2Dplot_of(output_field[-1], self.space_grid.get_X(), self.space_grid.get_Y(), plot_zoom_x=60, plot_zoom_y=60, choose_quantity="amplitude", save_name=f"{self.results_folder}/output_fields_{i}_plot") + make_2Dplot_of(output_field[-1], self.space_grid.get_X(), self.space_grid.get_Y(), plot_zoom_x=self.get_ms_max_avail_diameter()/2, plot_zoom_y=self.get_ms_max_avail_diameter()/2, choose_quantity="amplitude", save_name=f"{self.results_folder}/output_fields_{i}_plot") + make_2Dplot_of(output_field[-1], self.space_grid.get_X(), self.space_grid.get_Y(), plot_zoom_x=self.get_ms_max_avail_diameter()/2, plot_zoom_y=self.get_ms_max_avail_diameter()/2, choose_quantity="phase", save_name=f"{self.results_folder}/output_fields_{i}_phase_plot") #phase masks phase_mask_1 = self._phase_given_w(self.w_norm[:self.space_grid.get_S()]).reshape((self.space_grid.get_Nx(), self.space_grid.get_Ny())) phase_mask_2 = self._phase_given_w(self.w_norm[self.space_grid.get_S():]).reshape((self.space_grid.get_Nx(), self.space_grid.get_Ny())) - make_2Dplot_of(phase_mask_1, self.space_grid.get_X(), self.space_grid.get_Y(), plot_zoom_x=60, plot_zoom_y=60, choose_quantity="phase_mask", save_name=f"{self.results_folder}/phase_mask_1_plot") - make_2Dplot_of(phase_mask_2, self.space_grid.get_X(), self.space_grid.get_Y(), plot_zoom_x=60, plot_zoom_y=60, choose_quantity="phase_mask", save_name=f"{self.results_folder}/phase_mask_2_plot") + filtered_phase_mask_1 = self._phase_filter(phase_mask_1) + filtered_phase_mask_2 = self._phase_filter(phase_mask_2) + make_2Dplot_of(phase_mask_1, self.space_grid.get_X(), self.space_grid.get_Y(), plot_zoom_x=self.get_ms_max_avail_diameter()/2, plot_zoom_y=self.get_ms_max_avail_diameter()/2, choose_quantity="phase_mask", save_name=f"{self.results_folder}/phase_mask_1_plot") + make_2Dplot_of(phase_mask_2, self.space_grid.get_X(), self.space_grid.get_Y(), plot_zoom_x=self.get_ms_max_avail_diameter()/2, plot_zoom_y=self.get_ms_max_avail_diameter()/2, choose_quantity="phase_mask", save_name=f"{self.results_folder}/phase_mask_2_plot") + make_2Dplot_of(filtered_phase_mask_1, self.space_grid.get_X(), self.space_grid.get_Y(), plot_zoom_x=self.get_ms_max_avail_diameter()/2, plot_zoom_y=self.get_ms_max_avail_diameter()/2, choose_quantity="phase_mask", save_name=f"{self.results_folder}/filtered_phase_mask_1_plot") + make_2Dplot_of(filtered_phase_mask_2, self.space_grid.get_X(), self.space_grid.get_Y(), plot_zoom_x=self.get_ms_max_avail_diameter()/2, plot_zoom_y=self.get_ms_max_avail_diameter()/2, choose_quantity="phase_mask", save_name=f"{self.results_folder}/filtered_phase_mask_2_plot") + np.save(f"{self.results_folder}/optimized_phase_mask_1.npy", phase_mask_1) np.save(f"{self.results_folder}/optimized_phase_mask_2.npy", phase_mask_2) + np.save(f"{self.results_folder}/filtered_phase_mask_1.npy", filtered_phase_mask_1) + np.save(f"{self.results_folder}/filtered_phase_mask_2.npy", filtered_phase_mask_2) def get_best_opt_val(self): if self.is_optimized_flag: diff --git a/testing_converter.py b/single_converter_symm.py similarity index 84% rename from testing_converter.py rename to single_converter_symm.py index 95b2c5f..adcb758 100644 --- a/testing_converter.py +++ b/single_converter_symm.py @@ -1,6 +1,8 @@ import mpid as pid import numpy as np - +import os +if not os.path.exists("results"): + os.makedirs("results") #Space grid definition N = 1025 unit_cell_pitch = 700e-9 @@ -42,16 +44,16 @@ sources = [pid.Source(input_field=gaussian_field_normalized, sources[0].plot2D(plot_zoom_x=38,plot_zoom_y=38,choose_quantity='amplitude', save_name='results/input_field') #define targets #import field data from files and create target objects -target_files = ["fiber_modes/LP21_field_distribution_1025.npy"] +target_files = ["fiber_modes/LP21_field_distribution_lda_1550_N_1025.npy",] target_fields = [np.load(file) for file in target_files] -norm_target_fields = [field / np.sqrt(np.sum(np.abs(field)**2)) for field in target_fields] -targets = [pid.Target(input_field=norm_target_fields[i], +targets = [pid.Target(input_field=target_fields[i], associated_source=sources[i], wavelength=wvl, space_grid=space_grid, - wvl_unit_measure=space_grid.nm) for i, (file, wvl) in enumerate(zip(target_files, wavelengths_to_test_nm))] + wvl_unit_measure=space_grid.nm, + field_norm = True) for i, (file, wvl) in enumerate(zip(target_files, wavelengths_to_test_nm))] targets[0].plot2D(plot_zoom_x=38,plot_zoom_y=38,choose_quantity='amplitude', save_name='results/target_field') @@ -70,6 +72,7 @@ simulation.set_optimization( 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.5e-6, ) simulation.run() diff --git a/testing_filtering.py b/sorter_to_gauss.py similarity index 99% rename from testing_filtering.py rename to sorter_to_gauss.py index ac4b205..1a92055 100644 --- a/testing_filtering.py +++ b/sorter_to_gauss.py @@ -56,7 +56,6 @@ target2 = pid.Target(input_field='fiber_modes/LP01_field_distribution_lda_1550_N field_norm=True) #make a plot of sources and targets - target.plot2D(plot_zoom_x=100, plot_zoom_y=100, choose_quantity="amplitude", save_name="results/target_plot") target2.plot2D(plot_zoom_x=100, plot_zoom_y=100, choose_quantity="amplitude", save_name="results/target2_plot") source1.plot2D(plot_zoom_x=100, plot_zoom_y=100, choose_quantity="amplitude", save_name="results/source1_plot") diff --git a/sorter_to_ident.py b/sorter_to_ident.py new file mode 100644 index 0000000..93b2246 --- /dev/null +++ b/sorter_to_ident.py @@ -0,0 +1,86 @@ +import mpid as pid +import numpy as np +import os +if not os.path.exists("results"): + os.makedirs("results") + +#Space grid definition +N = 1025 +unit_cell_pitch = 700e-9 + +n0 = 1.45 +n1 = 1.0 +n2 = 1.45 +mat_refr_index = np.array([n0, n1, n2]) + +d0 = 476e-6 #propagation distance: source - MS1 [m] +d1 = 510e-6 #propagation distance: MS1 - MS2 [m] +d2 = d0 #propagation distance: MS2 - target [m] +d = np.array([d0, d1, d2]) + +ms_max_available_diameter = 122e-6 #[um] +output_fiber_core_diameter = 19e-6 #[um] +opt_max_eval = 140 + +space_grid = pid.SpaceGrid( + side_pixels=N, + unit_cell_pitch = unit_cell_pitch, + prop_distances=d, + materials_refr_index=mat_refr_index, +) + +wavelength = 1.55e-6 + +#define the sources +source1 = pid.Source(input_field='fiber_modes/LP11_field_distribution_lda_1550_N_1025.npy', + wavelength=wavelength, + space_grid=space_grid, + field_norm=True) +source2 = pid.Source(input_field='fiber_modes/LP21_field_distribution_lda_1550_N_1025.npy', + wavelength=wavelength, + space_grid=space_grid, + field_norm=True) + + +target = pid.Target(input_field='fiber_modes/LP11_field_distribution_lda_1550_N_1025.npy', + center=(20e-6,0), + associated_source=source1, + wavelength=wavelength, + space_grid=space_grid, + field_norm=True) +target2 = pid.Target(input_field='fiber_modes/LP21_field_distribution_lda_1550_N_1025.npy', + center=(-20e-6,0), + associated_source=source2, + wavelength=wavelength, + space_grid=space_grid, + field_norm=True) + +#make a plot of sources and targets +target.plot2D(plot_zoom_x=100, plot_zoom_y=100, choose_quantity="amplitude", save_name="results/target_plot") +target2.plot2D(plot_zoom_x=100, plot_zoom_y=100, choose_quantity="amplitude", save_name="results/target2_plot") +source1.plot2D(plot_zoom_x=100, plot_zoom_y=100, choose_quantity="amplitude", save_name="results/source1_plot") +source2.plot2D(plot_zoom_x=100, plot_zoom_y=100, choose_quantity="amplitude", save_name="results/source2_plot") + +sim = pid.Simulation( + sources=[source1, source2], + targets=[target, target2], + central_wavelength=wavelength, + space_grid=space_grid, + n_cpu=1 +) + +opt_param_min = 0 +opt_param_max = 1 + +sim.set_optimization(opt_param_min=opt_param_min, + opt_param_max=opt_param_max, + ms_max_avail_diameter=ms_max_available_diameter, + phase_mask_init_size=output_fiber_core_diameter, + opt_max_eval=opt_max_eval, + conic_filter_radius=1.5e-6) + +sim.run() + +best_opt_val, distances = sim.get_best_opt_and_distances() +print("Best optimization value:", best_opt_val) +print("Distances:", distances) \ No newline at end of file