From 9da2beaf75f307bf3b3e4368efd44e5babb8b63c Mon Sep 17 00:00:00 2001 From: Edoardo Vergerio Date: Thu, 16 Jul 2026 08:49:54 +0200 Subject: [PATCH] minor changes to plots and parameters --- metaprop_fibermode_doublet_multi.py | 245 ++++++++++++++++------------ 1 file changed, 144 insertions(+), 101 deletions(-) diff --git a/metaprop_fibermode_doublet_multi.py b/metaprop_fibermode_doublet_multi.py index 5e8d923..795db08 100644 --- a/metaprop_fibermode_doublet_multi.py +++ b/metaprop_fibermode_doublet_multi.py @@ -22,7 +22,7 @@ from scipy.special import jv, kv wavelength = 1.55e-6 k0 = 2 * np.pi / wavelength -opt_max_eval = 250 +opt_max_eval = 300 unit_cell_pitch = 400e-9 #equivalent to spatial sampling @@ -43,25 +43,21 @@ sampling_period = xs[1] - xs[0] d0 = 345e-6 #propagation distance: source - MS1 d1 = 500e-6 #propagation distance: MS1 - MS2 -d2 = 500e-6 #propagation distance: MS2 - target +d2 = 300e-6 #propagation distance: MS2 - target d = [d1, d2] #d = [d1,d2] d1:distance MS1-MS2, d2: distance MS2-target phase_min = 0 phase_max = 2 * np.pi ############################## -############################## - -############################## - ############################## #Source -beam_waist = 9e-6 -source_spacing = 150e-6 +beam_waist = 4e-6 +source_spacing = 300e-6 ############################## ############################## -#plan pyfftw objects for forward and inverse FFTs +#PYFFTW PLANNING SETUP fft_input_array = pyfftw.empty_aligned((Nx, Ny), dtype='complex128', n=16) fft_output_array = pyfftw.empty_aligned((Nx, Ny), dtype='complex128', n=16) ifft_input_array = pyfftw.empty_aligned((Nx, Ny), dtype='complex128', n=16) @@ -97,8 +93,10 @@ P_2_dagger_nat = np.fft.ifftshift(P_2).T.conj() ############################## #Plot zoom parameters [um] -zoom_x = 60 -zoom_y = 60 +zoom_x = 350 +zoom_x_outmode = 50 +zoom_y = 350 +zoom_y_outmode = 50 full_view_x = size_x * 1e6 / 2 full_view_y = size_y * 1e6 / 2 ############################## @@ -119,9 +117,15 @@ def get_pattern() -> np.ndarray: return pattern.flatten() -def get_fiber_mode_pattern(n_mode_sel): - '''Returns the target fiber mode pattern''' - print("Computing target fiber mode pattern for LP %d1" % n_mode_sel) +def get_fiber_mode_pattern(mode_list:list)->list: + '''Returns the fiber mode pattern for the selected n index. + Args: + mode_list: the list of the selected fiber mode indices + Returns: + The fiber mode patterns corresponding to the selected indices returned + as a list of 2D arrays with the same order as the input mode_list + ''' + print("Computing fiber patterns for mode indices: ", mode_list) lda = 1.55e-6 n_co = 1.4630 @@ -136,15 +140,13 @@ def get_fiber_mode_pattern(n_mode_sel): V_list = np.linspace(2, 5.1, 20) V_obs = v_1550 #b, V = np.meshgrid(b_list, V_list) - n_list = [0,1,2] b_solutions_at_V_obs = [] - for n in n_list: + for n in mode_list: b_solutions = [] for V in V_list: LHS = np.sqrt(1-b) * jv(n+1, V * np.sqrt(1-b)) / jv(n, V * np.sqrt(1-b)) RHS = np.sqrt(b) * kv(n+1, V * np.sqrt(b)) / kv(n, V * np.sqrt(b)) - intersection_index = np.argmax(np.where(np.isclose(LHS, RHS, atol=0.000001),1,0).flatten()) @@ -154,39 +156,43 @@ def get_fiber_mode_pattern(n_mode_sel): RHS = np.sqrt(b) * kv(n+1, V_obs * np.sqrt(b)) / kv(n, V_obs * np.sqrt(b)) b_solutions_at_V_obs.append(b[np.argmax(np.where(np.isclose(LHS, RHS, atol=0.00001),1,0).flatten())]) - b_mode_sel = b_solutions_at_V_obs[n_mode_sel] + output_patterns = [] + for enum_index,n_mode_sel in enumerate(mode_list): + b_mode_sel = b_solutions_at_V_obs[enum_index] - k_co = 2 * np.pi / lda * n_co - k_cl = 2 * np.pi / lda * n_cl + k_co = 2 * np.pi / lda * n_co + k_cl = 2 * np.pi / lda * n_cl - beta = np.sqrt(b_mode_sel * (k_co**2 - k_cl**2) + k_cl**2) + beta = np.sqrt(b_mode_sel * (k_co**2 - k_cl**2) + k_cl**2) - chi_co = v_1550 * np.sqrt(1 - b_mode_sel) / a - chi_cl = v_1550 * np.sqrt(b_mode_sel) / a + chi_co = v_1550 * np.sqrt(1 - b_mode_sel) / a + chi_cl = v_1550 * np.sqrt(b_mode_sel) / a - E_field = np.zeros_like(X, dtype=complex) + E_field = np.zeros_like(X, dtype=complex) - for i in range(Nx): - for j in range(Ny): - r = np.sqrt(X[i,j]**2 + Y[i,j]**2) - if r < 10*a: - phi = np.arctan2(Y[i,j], X[i,j]) - if r < a:#inside core - E_field[i,j] = jv(n_mode_sel,chi_co * r) / jv(n_mode_sel, chi_co * a) * np.cos(n_mode_sel * phi) - else:#outside core - E_field[i,j] = kv(n_mode_sel, chi_cl * r) / kv(n_mode_sel, chi_cl * a) * np.cos(n_mode_sel * phi) + for i in range(Nx): + for j in range(Ny): + r = np.sqrt(X[i,j]**2 + Y[i,j]**2) + if r < 10*a: + phi = np.arctan2(Y[i,j], X[i,j]) + if r < a:#inside core + E_field[i,j] = jv(n_mode_sel,chi_co * r) / jv(n_mode_sel, chi_co * a) * np.cos(n_mode_sel * phi) + else:#outside core + E_field[i,j] = kv(n_mode_sel, chi_cl * r) / kv(n_mode_sel, chi_cl * a) * np.cos(n_mode_sel * phi) - return E_field + output_patterns.append(E_field) + return output_patterns def get_source_field(beam_waist,center_x=0,center_y=0,): - '''Returns a normalized Gaussian source field for a given beam waist. Normalization consists in scaling - the field such that its integrated intensity is 1. + '''Returns a 2D normalized Gaussian source field for a given beam waist and centered at the specified (x,y) coordinates. + Normalization consists in scaling the field such that its integrated intensity is 1. + Args: beam_waist: the waist of the Gaussian beam (in meters) center_x: the x-coordinate of the center of the Gaussian beam (in meters) center_y: the y-coordinate of the center of the Gaussian beam (in meters) Returns: - The normalized Gaussian source field + The 2D normalized Gaussian source field. ''' if center_x != 0 or center_y != 0: @@ -197,16 +203,24 @@ def get_source_field(beam_waist,center_x=0,center_y=0,): gaussian_field = np.exp(-(rho)**2 / (beam_waist)**2) #Gaussian input field, flattened gaussian_field_normalized = (gaussian_field / np.sqrt(np.sum(np.abs(gaussian_field)**2))) #normalize the input field to have power = 1 + return gaussian_field_normalized def propagate_source(source_to_propagate): - '''Returns the propagated source field (flattened) just before metasurface 1 (MS1). This involves propagating the - source for a distance d0: from its plane of definition up to MS1.''' - fft_input_array[:,:] = source_to_propagate.reshape((Nx, Ny)) + '''Returns the propagated source field (flattened) just before metasurface 1 (MS1). + This involves propagating the source for a distance d0: from its plane of definition up to MS1. + Args: + source_to_propagate: the source field to be propagated (flattened) + Returns: + The propagated source field (flattened) just before MS1. + ''' + + fft_input_array[:,:] = np.fft.ifftshift(source_to_propagate.reshape((Nx, Ny))) fft_operator() ifft_input_array[:,:] = fft_operator.output_array * P_0_nat ifft_operator() - propagated_source_field = ifft_operator.output_array.flatten() + propagated_source_field = np.fft.fftshift(ifft_operator.output_array).flatten() + return propagated_source_field def make_1Dplot_of(x_axis,y_axis,plot_zoom_x=zoom_x,save_name="plot_1D",) -> None: @@ -263,7 +277,7 @@ def make_2Dplot_of(given_field,choose_quantity="amplitude",save_name="plot",plot plt.close() def shift_spatial_grid(shift_amount_x,shift_amount_y): - '''Shifts the spatial grid by a given amount. This is useful for centering the source field on the metasurface. + '''Shifts the simulation spatial grid by a given physical (true space) amount. Args: shift_amount_x: amount to shift the grid in x direction (in meters) shift_amount_y: amount to shift the grid in y direction (in meters) @@ -273,14 +287,10 @@ def shift_spatial_grid(shift_amount_x,shift_amount_y): Yp = Y + shift_amount_y return Xp, Yp - ########################### #initialize source field -source_field_1 = get_source_field(beam_waist,center_x=-source_spacing/2,center_y=0) -source_field_1_2d = source_field_1.reshape((Nx, Ny)) #reshape to 2D for plotting - -source_field_2 = get_source_field(beam_waist,center_x=source_spacing/2,center_y=0) -source_field_2_2d = source_field_2.reshape((Nx, Ny)) #reshape to 2D for plotting +source_field_1 = get_source_field(beam_waist, center_x=-source_spacing/2, center_y=0) +source_field_2 = get_source_field(beam_waist, center_x=source_spacing/2, center_y=0) #source field integrated intensity source_1_power = np.sum(np.abs(source_field_1)**2) @@ -288,10 +298,10 @@ source_2_power = np.sum(np.abs(source_field_2)**2) print("Source 1 field integrated intensity: %.4e" % source_1_power) print("Source 2 field integrated intensity: %.4e" % source_2_power) #make a plot and save the source field -make_2Dplot_of(source_field_1_2d, choose_quantity="amplitude", save_name="source_1_field_amplitude", plot_zoom_x=3*zoom_x, plot_zoom_y=3*zoom_y) -make_2Dplot_of(source_field_1_2d, choose_quantity="phase", save_name="source_1_field_phase", plot_zoom_x=3*zoom_x, plot_zoom_y=3*zoom_y) -make_2Dplot_of(source_field_2_2d, choose_quantity="amplitude", save_name="source_2_field_amplitude", plot_zoom_x=3*zoom_x, plot_zoom_y=3*zoom_y) -make_2Dplot_of(source_field_2_2d, choose_quantity="phase", save_name="source_2_field_phase", plot_zoom_x=3*zoom_x, plot_zoom_y=3*zoom_y) +make_2Dplot_of(source_field_1, choose_quantity="amplitude", save_name="source_1_field_amplitude", plot_zoom_x=full_view_x, plot_zoom_y=full_view_y) +make_2Dplot_of(source_field_1, choose_quantity="phase", save_name="source_1_field_phase", plot_zoom_x=full_view_x, plot_zoom_y=full_view_y) +make_2Dplot_of(source_field_2, choose_quantity="amplitude", save_name="source_2_field_amplitude", plot_zoom_x=full_view_x, plot_zoom_y=full_view_y) +make_2Dplot_of(source_field_2, choose_quantity="phase", save_name="source_2_field_phase", plot_zoom_x=full_view_x, plot_zoom_y=full_view_y) ########################### ########################### @@ -312,21 +322,19 @@ input_field_2_intensity = np.sum(np.abs(input_field_2)**2) print("Input field integrated intensity: %.4e" % input_field_2_intensity) #make a plot and save the input field -make_2Dplot_of(input_field_1_2d, choose_quantity="amplitude", save_name="input_field_1_amplitude", plot_zoom_x=zoom_x*3, plot_zoom_y=zoom_y*3) -make_2Dplot_of(input_field_1_2d, choose_quantity="phase", save_name="input_field_1_phase", plot_zoom_x=zoom_x*3, plot_zoom_y=zoom_y*3) -make_2Dplot_of(input_field_2_2d, choose_quantity="amplitude", save_name="input_field_2_amplitude", plot_zoom_x=zoom_x*3, plot_zoom_y=zoom_y*3) -make_2Dplot_of(input_field_2_2d, choose_quantity="phase", save_name="input_field_2_phase", plot_zoom_x=zoom_x*3, plot_zoom_y=zoom_y*3) +make_2Dplot_of(input_field_1_2d, choose_quantity="amplitude", save_name="input_field_1_amplitude", plot_zoom_x=zoom_x, plot_zoom_y=zoom_y) +make_2Dplot_of(input_field_1_2d, choose_quantity="phase", save_name="input_field_1_phase", plot_zoom_x=zoom_x, plot_zoom_y=zoom_y) +make_2Dplot_of(input_field_2_2d, choose_quantity="amplitude", save_name="input_field_2_amplitude", plot_zoom_x=zoom_x, plot_zoom_y=zoom_y) +make_2Dplot_of(input_field_2_2d, choose_quantity="phase", save_name="input_field_2_phase", plot_zoom_x=zoom_x, plot_zoom_y=zoom_y) #list of input fields to be propagated within the optimization loop input_field_list = [input_field_1, input_field_2] -#input_field_list = [input_field_1] ########################### ########################### #initialize target field target_mode_indices = [1,2] -fiber_mode_pattern_1 = get_fiber_mode_pattern(target_mode_indices[0]) -fiber_mode_pattern_2 = get_fiber_mode_pattern(target_mode_indices[1]) +fiber_mode_pattern_1, fiber_mode_pattern_2 = get_fiber_mode_pattern(target_mode_indices) target_Efield_1_2d = fiber_mode_pattern_1 / np.sqrt(np.sum(np.abs(fiber_mode_pattern_1)**2)) target_Efield_1 = target_Efield_1_2d.flatten() #normalized target field (intensity = 1) @@ -341,7 +349,6 @@ make_2Dplot_of(target_Efield_2_2d, choose_quantity="amplitude", save_name="targe make_2Dplot_of(target_Efield_2_2d, choose_quantity="phase", save_name="target_field_2_phase", plot_zoom_x=zoom_x, plot_zoom_y=zoom_y) target_field_list = [target_Efield_1, target_Efield_2] -#target_field_list = [target_Efield_1] ########################### def forward_propagate(input_field:np.ndarray,weights:np.ndarray) -> np.ndarray: @@ -355,26 +362,26 @@ def forward_propagate(input_field:np.ndarray,weights:np.ndarray) -> np.ndarray: phase_mask = phase_given_w(weights) fft_input_array[:,:] = np.fft.ifftshift((input_field * np.exp(1j * phase_mask[:S])).reshape((Nx, Ny))) fft_operator() - print("BEFORE") + #print("BEFORE") fft_power = np.sum(np.abs(fft_operator.output_array)**2)/S non_propagating_mask = np.where((wavelength * nu_parallel)**2 > 1, 0, 1) non_propagating_mask = np.fft.ifftshift(non_propagating_mask) propagating_power = np.sum(np.abs(fft_operator.output_array * non_propagating_mask)**2)/S - print("FFT power: %.4e, Propagating power: %.4e, Non-propagating power: %.4e" % (fft_power, propagating_power, fft_power - propagating_power)) + #print("FFT power: %.4e, Propagating power: %.4e, Non-propagating power: %.4e" % (fft_power, propagating_power, fft_power - propagating_power)) #FIRST PROPAGATION ifft_input_array[:,:] = fft_operator.output_array * P_1_nat #element wise product - print("AFTER") + #print("AFTER") fft_power = np.sum(np.abs(ifft_operator.input_array)**2)/S non_propagating_mask = np.where((wavelength * nu_parallel)**2 > 1, 0, 1) non_propagating_mask = np.fft.ifftshift(non_propagating_mask) propagating_power = np.sum(np.abs(ifft_operator.input_array * non_propagating_mask)**2)/S - print("FFT power: %.4e, Propagating power: %.4e, Non-propagating power: %.4e" % (fft_power, propagating_power, fft_power - propagating_power)) + #print("FFT power: %.4e, Propagating power: %.4e, Non-propagating power: %.4e" % (fft_power, propagating_power, fft_power - propagating_power)) ifft_operator() intermediate_fields[0] = np.fft.fftshift(ifft_operator.output_array.copy()) - print("Field intensity after first phase mask and propagation: %.4e" % np.sum(np.abs(intermediate_fields[0])**2)) + #print("Field intensity after first phase mask and propagation: %.4e" % np.sum(np.abs(intermediate_fields[0])**2)) fft_input_array[:,:] = np.fft.ifftshift((np.fft.fftshift(ifft_operator.output_array) * np.exp(1j * phase_mask[S:].reshape((Nx, Ny))))) fft_operator() @@ -392,14 +399,15 @@ def forward_propagate(input_field:np.ndarray,weights:np.ndarray) -> np.ndarray: return intermediate_fields -def adjoint_propagate(input_field_2d,target_Efield_2d,intermediate_fields): +def adjoint_propagate(input_field_2d,target_Efield_2d,intermediate_fields,input_weights): '''Backpropagates the output field using the adjoint of the propagation matrix P. Returns: - the gradient of the cost function with respect to the phase mask (flattened) + the gradients of the cost function with respect to the phase mask parameters (flattened) ''' print("Computing adjoint...") start_time = time.time() + phase_mask = phase_given_w(input_weights) Phi_1_dagger = np.exp(-1j * phase_mask[:S]).reshape((Nx, Ny)) Phi_2_dagger = np.exp(-1j * phase_mask[S:]).reshape((Nx, Ny)) @@ -424,6 +432,16 @@ def adjoint_propagate(input_field_2d,target_Efield_2d,intermediate_fields): def value_and_grad(weights): + ''' + Returns the value of the cost function and its gradient with respect to the phase mask parameters for + every available set of input/target fields. + Args: + weights: the current phase mask parameters (flattened) + Returns: + values: the value of the cost function for each input/target field pair + gradients: the gradient of the cost function with respect to the phase mask parameters for each input/target field pair + ''' + values = [] gradients = [] @@ -432,8 +450,8 @@ def value_and_grad(weights): intermediate_fields = forward_propagate(input_field, weights) - circular_mask = np.where(rho.flatten() < 10*6e-6, 1, 0) - C_s = np.abs(np.sum(intermediate_fields[-1].flatten() * circular_mask * np.conj(target_field_list[i]))) ** 2 + #circular_mask = np.where(rho.flatten() < 10*6e-6, 1, 0) + C_s = np.abs(np.sum(intermediate_fields[-1].flatten() * np.conj(target_field_list[i]))) ** 2 C = np.real(np.sum(C_s)) values.append(C) @@ -441,7 +459,9 @@ def value_and_grad(weights): adjoint_propagate( input_field_2d = input_field.reshape((Nx, Ny)), target_Efield_2d = target_field_list[i].reshape((Nx, Ny)), - intermediate_fields = intermediate_fields) + intermediate_fields = intermediate_fields, + input_weights = weights + ) ) return np.array(values), np.array(gradients) @@ -449,6 +469,14 @@ def value_and_grad(weights): opt_history = [] def obj_fun(weights, grad): ''' + Objective function for the optimization. Computes the overall cost function by taking the mean value + of the individual cost functions for each input/target field pair. Also computes the gradient of the overall cost function + with respect to the weights parameters. + Args: + weights: the current normalized phase mask parameters (flattened) + grad: the gradient of the overall cost function with respect to the weights parameters (flattened) + Returns: + average_obj_val: the average value of the cost function across all input/target field pairs ''' w_norm[:] = weights @@ -469,31 +497,19 @@ def obj_fun(weights, grad): make_2Dplot_of(np.abs(grad[S:].reshape((Nx, Ny))), choose_quantity="amplitude", save_name="gradient_magnitude", plot_zoom_x=full_view_x, plot_zoom_y=full_view_y) opt_history.append(average_obj_val) - - plt.figure() - plt.plot(opt_history) - plt.xlabel('Iteration') - plt.ylabel('Cost function') - plt.title('Optimization History') - plt.tight_layout() - plt.savefig("results/modeconv1_optimization_history.pdf") - plt.close() + make_1Dplot_of([np.arange(len(opt_history))], [opt_history], plot_zoom_x=len(opt_history), save_name="optimization_history") return average_obj_val -if __name__ == "__main__": - ########################### - #initialize the normalized phase masks and global phase mask - w_norm = np.zeros((2*S)) #concatenated normalized parameters for both masks - +def initialize_masks(): X1,Y1 = shift_spatial_grid(-source_spacing/2,0) rho1 = np.sqrt(X1**2 + Y1**2) - circular_mask_1 = np.where(rho1.flatten() < 10*6e-6, 1, 0) + circular_mask_1 = np.where(rho1.flatten() < 3*beam_waist, 1, 0) X2,Y2 = shift_spatial_grid(source_spacing/2,0) rho2 = np.sqrt(X2**2 + Y2**2) - circular_mask_2 = np.where(rho2.flatten() < 10*6e-6, 1, 0) + circular_mask_2 = np.where(rho2.flatten() < 3*beam_waist, 1, 0) target_pattern_1_norm_phase = (np.angle(target_Efield_1) + 2 * np.pi) % (2 * np.pi) / (2*np.pi) target_pattern_1_norm_phase = np.roll(target_pattern_1_norm_phase, shift=int(source_spacing/sampling_period/2), axis=0) @@ -502,10 +518,16 @@ if __name__ == "__main__": target_pattern_2_norm_phase = np.roll(target_pattern_2_norm_phase, shift=-int(source_spacing/sampling_period/2), axis=0) #on metasurface 1, set initial guess on normalized phase elements w_norm[:S] += circular_mask_1 * target_pattern_1_norm_phase + circular_mask_2 * target_pattern_2_norm_phase - #on metasurface 2, generate random initial guess for normalized phase elements - #w_norm[S:] += np.random.rand(S) * circular_mask_1 + np.random.rand(S) * circular_mask_2 + + +if __name__ == "__main__": + ########################### + #initialize the normalized phase masks and global phase mask + w_norm = np.zeros((2*S)) #concatenated normalized parameters for both masks + initialize_masks() phase_mask = phase_given_w(w_norm) + #make an initial plot of the phase masks make_2Dplot_of(phase_mask[:S].reshape((Nx, Ny)), choose_quantity="phase_mask", save_name="initial_phase_mask_1", plot_zoom_x=full_view_x, plot_zoom_y=full_view_y) make_2Dplot_of(phase_mask[S:].reshape((Nx, Ny)), choose_quantity="phase_mask", save_name="initial_phase_mask_2", plot_zoom_x=full_view_x, plot_zoom_y=full_view_y) @@ -519,9 +541,6 @@ if __name__ == "__main__": mask_upper_bounds = np.concatenate((mask_upper_bounds, mask_upper_bounds),axis=0) mask_lower_bounds = np.zeros_like(mask_upper_bounds) - - #epigraph_and_weights = np.insert(w_norm,0,0) #initialize epigraph variable to 0 - #print("#####Shape of epigraph_and_weights: ", epigraph_and_weights.shape) ########################### #initialize nlopt solver @@ -550,17 +569,42 @@ if __name__ == "__main__": all_source_on_field = input_field_list[0] + input_field_list[1] verify_out_field_all = forward_propagate(all_source_on_field, opt_weights)[-1] - make_2Dplot_of(verify_out_field_1, choose_quantity="amplitude", save_name="modeconv1_optimized_output_field_source1_amplitude", plot_zoom_x=zoom_x, plot_zoom_y=zoom_y) - make_2Dplot_of(verify_out_field_1, choose_quantity="phase", save_name="modeconv1_optimized_output_field_source1_phase", plot_zoom_x=zoom_x, plot_zoom_y=zoom_y) - make_2Dplot_of(verify_out_field_2, choose_quantity="amplitude", save_name="modeconv1_optimized_output_field_source2_amplitude", plot_zoom_x=zoom_x, plot_zoom_y=zoom_y) - make_2Dplot_of(verify_out_field_2, choose_quantity="phase", save_name="modeconv1_optimized_output_field_source2_phase", plot_zoom_x=zoom_x, plot_zoom_y=zoom_y) - make_2Dplot_of(verify_out_field_all, choose_quantity="amplitude", save_name="modeconv1_optimized_output_field_allSource_amplitude", plot_zoom_x=zoom_x, plot_zoom_y=zoom_y) - make_2Dplot_of(verify_out_field_all, choose_quantity="phase", save_name="modeconv1_optimized_output_field_allSource_phase", plot_zoom_x=zoom_x, plot_zoom_y=zoom_y) + make_2Dplot_of(verify_out_field_1, choose_quantity="amplitude", save_name="modeconv1_optimized_output_field_source1_amplitude", plot_zoom_x=zoom_x_outmode, plot_zoom_y=zoom_y_outmode) + make_2Dplot_of(verify_out_field_1, choose_quantity="phase", save_name="modeconv1_optimized_output_field_source1_phase", plot_zoom_x=zoom_x_outmode, plot_zoom_y=zoom_y_outmode) + make_2Dplot_of(verify_out_field_2, choose_quantity="amplitude", save_name="modeconv1_optimized_output_field_source2_amplitude", plot_zoom_x=zoom_x_outmode, plot_zoom_y=zoom_y_outmode) + make_2Dplot_of(verify_out_field_2, choose_quantity="phase", save_name="modeconv1_optimized_output_field_source2_phase", plot_zoom_x=zoom_x_outmode, plot_zoom_y=zoom_y_outmode) + make_2Dplot_of(verify_out_field_all, choose_quantity="amplitude", save_name="modeconv1_optimized_output_field_allSource_amplitude", plot_zoom_x=zoom_x_outmode, plot_zoom_y=zoom_y_outmode) + make_2Dplot_of(verify_out_field_all, choose_quantity="phase", save_name="modeconv1_optimized_output_field_allSource_phase", plot_zoom_x=zoom_x_outmode, plot_zoom_y=zoom_y_outmode) + + #make 1D slices of target/outpu fields to compare + slice_y_index = Ny//2 + slice_x = xs * 1e6 + slice_target_1 = target_Efield_1.reshape((Nx, Ny))[slice_y_index,:] + slice_target_2 = target_Efield_2.reshape((Nx, Ny))[slice_y_index,:] + slice_output_1 = verify_out_field_1.reshape((Nx, Ny))[slice_y_index,:] + slice_output_2 = verify_out_field_2.reshape((Nx, Ny))[slice_y_index,:] + slice_output_all = verify_out_field_all.reshape((Nx, Ny))[slice_y_index,:] + #amplitude and phase target 1 + make_1Dplot_of([slice_x, slice_x], [np.abs(slice_target_1), np.abs(slice_output_1)], plot_zoom_x=zoom_x_outmode, save_name="modeconv1_output_vs_target_source1_amplitude") + make_1Dplot_of([slice_x, slice_x], [np.angle(slice_target_1), np.angle(slice_output_1)], plot_zoom_x=zoom_x_outmode, save_name="modeconv1_output_vs_target_source1_phase") + #amplitude and phase target 2 + make_1Dplot_of([slice_x, slice_x], [np.abs(slice_target_2), np.abs(slice_output_2)], plot_zoom_x=zoom_x_outmode, save_name="modeconv1_output_vs_target_source2_amplitude") + make_1Dplot_of([slice_x, slice_x], [np.angle(slice_target_2), np.angle(slice_output_2)], plot_zoom_x=zoom_x_outmode, save_name="modeconv1_output_vs_target_source2_phase") + else: verify_out_field_1 = forward_propagate(input_field_list[0], opt_weights)[-1] make_2Dplot_of(verify_out_field_1, choose_quantity="amplitude", save_name="modeconv1_optimized_output_field_source1_amplitude", plot_zoom_x=zoom_x, plot_zoom_y=zoom_y) make_2Dplot_of(verify_out_field_1, choose_quantity="phase", save_name="modeconv1_optimized_output_field_source1_phase", plot_zoom_x=zoom_x, plot_zoom_y=zoom_y) - + + #make a slice of the target/output fields to compare + slice_y_index = Ny//2 + slice_x = xs * 1e6 + slice_target_1 = target_Efield_1.reshape((Nx, Ny))[slice_y_index,:] + slice_output_1 = verify_out_field_1.reshape((Nx, Ny))[slice_y_index,:] + #amplitude and phase target 1 + make_1Dplot_of([slice_x, slice_x], [np.abs(slice_target_1), np.abs(slice_output_1)], plot_zoom_x=zoom_x_outmode, save_name="modeconv1_output_vs_target_source1_amplitude") + make_1Dplot_of([slice_x, slice_x], [np.angle(slice_target_1), np.angle(slice_output_1)], plot_zoom_x=zoom_x_outmode, save_name="modeconv1_output_vs_target_source1_phase") + #save the output fields for later use saveFields = False if saveFields: @@ -571,10 +615,9 @@ if __name__ == "__main__": # Reshape fields to 2D for visualization phase_mask_1 = phase_given_w(opt_weights[:S]).reshape((Nx, Ny)) phase_mask_2 = phase_given_w(opt_weights[S:]).reshape((Nx, Ny)) - phase_mask_zoom = 150 # Plot: Phase masks - make_2Dplot_of(phase_mask_1, choose_quantity="phase_mask", save_name="phase_mask_1", plot_zoom_x=phase_mask_zoom, plot_zoom_y=phase_mask_zoom) - make_2Dplot_of(phase_mask_2, choose_quantity="phase_mask", save_name="phase_mask_2", plot_zoom_x=phase_mask_zoom, plot_zoom_y=phase_mask_zoom) + make_2Dplot_of(phase_mask_1, choose_quantity="phase_mask", save_name="phase_mask_1", plot_zoom_x=full_view_x, plot_zoom_y=full_view_y) + make_2Dplot_of(phase_mask_2, choose_quantity="phase_mask", save_name="phase_mask_2", plot_zoom_x=full_view_x, plot_zoom_y=full_view_y) #save the phase masks for later use saveMasks = True