minor changes to plots and parameters

2026-10-14 23:06:39 +00:00 · 2026-07-16 08:49:54 +02:00
parent 6260bbac4b
commit 9da2beaf75
1 changed files with 144 additions and 101 deletions
@@ -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,16 +140,14 @@ 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())

            b_solutions.append(b[intersection_index])
@@ -154,7 +156,9 @@ 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
@@ -176,17 +180,19 @@ def get_fiber_mode_pattern(n_mode_sel):
                    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