From f6fd791fffa7572887b7b957004894b8ffafbc54 Mon Sep 17 00:00:00 2001 From: Edoardo Vergerio Date: Fri, 19 Jun 2026 17:17:45 +0200 Subject: [PATCH] general testing of parameters. Regarding metaprop_fibermode_doublet, the spatial sampling has been fixed and the grid size is controlled by the number of pixels. I have seen that the mask-mask and mask-output plane distance is crucial and should not be chosen randomly, at the moment i have tested empirically some dimensions. The parameters in this commit seem to work ok, but we need to understand. Then, intensity conversation is still an issue, since the pattern at the output has 95/100, meaning 5/100 is lost somewhere. We really need to understand if this is an issue or not. With these settings, we might want to increase the window size, but then see the next point. I have introduced boundary boxes on phase masks, this is to avoid the algorithm using pixels that are clearly outside and only see noise. If we extend the simulation domain, we ought to extend these boundaries, otherwise the issue will present itself again. For what pixelation is concerned, honestly i expected worse, though there is still the need to apply some kind of filter to eradicate spurious single pixels. That should be doable using autograd, since the heavylifting is done by the adjoint, using autograd just for the filter part should not require too much computing effort. This needs to be tested. --- metaprop_asm.py | 28 +++++++++++----------- metaprop_fibermode_doublet.py | 44 +++++++++++++++++++++-------------- 2 files changed, 42 insertions(+), 30 deletions(-) diff --git a/metaprop_asm.py b/metaprop_asm.py index 33b0052..ee1557f 100644 --- a/metaprop_asm.py +++ b/metaprop_asm.py @@ -6,14 +6,12 @@ k0 = 2 * np.pi / wavelength unit_cell_pitch = 400e-9 -Nx = 1024*10 +Nx = 1024*2 Ny = Nx S = Nx * Ny -#size_x = unit_cell_pitch * Nx -#size_y = unit_cell_pitch * Ny -size_x = 7968e-6 -size_y = size_x +size_x = unit_cell_pitch * Nx +size_y = unit_cell_pitch * Ny norm_phase_min = 0 norm_phase_max = 1 @@ -63,8 +61,8 @@ def forward_propagate(phase_mask): KX, KY = np.meshgrid(kappas, kappas) K_parallel = np.sqrt(KX**2 + KY**2) nu_parallel = K_parallel / (2 * np.pi) - #phase_factor = k0 * z_prop * np.sqrt(1 - (wavelength * nu_parallel) ** 2 + 0*1j) - phase_factor = k0 * z_prop * (1 - 0.5*((wavelength * nu_parallel) ** 2 + 0*1j)) + phase_factor = k0 * z_prop * np.sqrt(1 - (wavelength * nu_parallel) ** 2 + 0*1j) + #phase_factor = k0 * z_prop * (1 - 0.5*((wavelength * nu_parallel) ** 2 + 0*1j)) P = np.exp(1j * phase_factor) fft_field_propagated = np.multiply(field_fft,P) #element wise product @@ -79,7 +77,7 @@ if __name__ == "__main__": output_field, P = forward_propagate(phase_mask_test) - plt.figure(figsize=(10, 12)) + plt.figure(figsize=(10, 10)) plt.subplot(4, 2, 1) #plt.imshow(np.abs(field_fft), extent=(-ks/2, ks/2, -ks/2, ks/2)) plt.imshow(np.abs(output_field.reshape((Nx,Ny))), extent=(-size_x/2*1e6, size_x/2*1e6, -size_y/2*1e6, size_y/2*1e6)) @@ -103,16 +101,16 @@ if __name__ == "__main__": #make this subplot take two places in the figure plt.subplot(4, 2, 3) - plt.plot(xs*1e6, np.abs(output_field.reshape((Nx,Ny))[Nx//2, :])**2, 'b-', label='output') - plt.plot(xs*1e6, gaussian_profile**2, 'r--', label='theory') + plt.plot(xs*1e6, np.abs(output_field.reshape((Nx,Ny))[Nx//2, :])**2 / np.max(np.abs(output_field.reshape((Nx,Ny))[Nx//2, :])**2), 'b-', label='output') + plt.plot(xs*1e6, gaussian_profile**2 / np.max(gaussian_profile**2), 'r--', label='theory') plt.xlabel('x (um)') plt.ylabel('Intensity (a.u.)') - plt.legend() + plt.legend(loc='center left') plt.subplot(4, 2, 4) plt.plot(xs*1e6, np.abs(input_field.reshape((Nx,Ny))[Nx//2, :])**2, 'g--', label='input') plt.xlabel('x (um)') plt.ylabel('Intensity (a.u.)') - plt.legend() + plt.legend(loc='center left') plt.subplot(4,1,3) plt.plot(xs*1e6, np.angle(output_field.reshape((Nx,Ny))[Nx//2, :])%(2*np.pi), 'b-', label='output phase') @@ -120,7 +118,8 @@ if __name__ == "__main__": plt.plot(xs*1e6, np.angle(input_field.reshape((Nx,Ny))[Nx//2, :]), 'g--', label='input phase') plt.xlabel('x (um)') plt.ylabel('Phase (rad)') - plt.legend() + plt.legend(loc='center left') + plt.subplot(4,1,4) plt.plot(xs*1e6, np.angle(output_field.reshape((Nx,Ny))[Nx//2, :])%(2*np.pi), 'b-', label='output phase') @@ -129,5 +128,8 @@ if __name__ == "__main__": plt.xlabel('x (um)') plt.xlim(np.max(size_x/4*1e6), np.max(size_x/2*1e6)) plt.ylabel('Phase (rad)') + plt.legend(loc='center left') + + plt.tight_layout() plt.savefig("propagation_verification.pdf") \ No newline at end of file diff --git a/metaprop_fibermode_doublet.py b/metaprop_fibermode_doublet.py index d9de014..04b80b4 100644 --- a/metaprop_fibermode_doublet.py +++ b/metaprop_fibermode_doublet.py @@ -18,11 +18,11 @@ from scipy.special import jv, kv wavelength = 1.55e-6 k0 = 2 * np.pi / wavelength -opt_max_eval = 250 +opt_max_eval = 200 unit_cell_pitch = 400e-9 #equivalent to spatial sampling -Nx = 8192 #pixels per dimension +Nx = 2048 #pixels per dimension Ny = Nx S = Nx * Ny @@ -40,9 +40,6 @@ fft_operator = pyfftw.FFTW(fft_input_array, fft_output_array, axes=(0,1), direct ifft_operator = pyfftw.FFTW(ifft_input_array, ifft_output_array, axes=(0,1), direction='FFTW_BACKWARD', threads=n_cpu, flags=('FFTW_MEASURE','FFTW_DESTROY_INPUT',)) ###################### -norm_phase_min = 0 -norm_phase_max = 1 - phase_min = 0 phase_max = 2 * np.pi @@ -54,8 +51,8 @@ rho = np.sqrt(X**2 + Y**2) sampling_period = xs[1] - xs[0] -d1 = 4000e-6 #propagation distance -d2 = 4000e-6 #propagation distance +d1 = 1000e-6 #propagation distance +d2 = 1000e-6 #propagation distance d = [d1, d2] #d = [d1,d2] d1:distance MS1-MS2, d2: distance MS2-target def phase_given_w(w): @@ -136,7 +133,7 @@ def get_fiber_mode_pattern(): ########################### #initialize input field -beam_waist = 300e-6 +beam_waist = 80e-6 gaussian_field = np.exp(-rho**2 / (beam_waist)**2) #Gaussian input field, flattened input_field_2d = (gaussian_field / np.sqrt(np.sum(np.abs(gaussian_field)**2))) #normalize the input field to have power = 1 input_field = input_field_2d.flatten() #normalize the input field to have power = 1 @@ -169,6 +166,15 @@ w_norm[:S] += circular_mask_1 target_pattern_phase = (np.angle(target_Efield) + 2 * np.pi) % (2 * np.pi) #make sure the phase is between 0 and 2pi w_norm[S:] += target_pattern_phase / (2 * np.pi) #normalize the target phase to be between 0 and 1 phase_mask = phase_given_w(w_norm) + +#define bounding box for the optimization. The outer mask is fixed to zero phase, while the inner mask can vary between 0 and 2pi phase. +#the bounding box is the outer shell of a square with thickness 100 um +norm_phase_min = 0 +norm_phase_max = 1 +mask_buffer = 100e-6 +mask_upper_bounds = norm_phase_max * np.where((np.abs(X.flatten()) > size_x/2 - mask_buffer) | (np.abs(Y.flatten()) > size_y/2 - mask_buffer), 0, 1) +mask_upper_bounds = np.concatenate((mask_upper_bounds, mask_upper_bounds),axis=0) +mask_lower_bounds = np.zeros_like(mask_upper_bounds) ########################### ########################### @@ -208,8 +214,9 @@ def forward_propagate() -> None: ifft_input_array[:,:] = fft_operator.output_array * P_1_nat #element wise product ifft_operator() intermediate_fields[0] = ifft_operator.output_array.copy() + print("Field intensity after first phase mask and propagation: %.4e" % np.sum(np.abs(intermediate_fields[0])**2)) - fft_input_array[:,:] = (ifft_operator.output_array.flatten() * np.exp(1j * phase_mask[S:])).reshape((Nx, Ny)) + fft_input_array[:,:] = (ifft_operator.output_array * np.exp(1j * phase_mask[S:].reshape((Nx, Ny)))) fft_operator() ifft_input_array[:,:] = fft_operator.output_array * P_2_nat #element wise product ifft_operator() @@ -244,8 +251,7 @@ def adjoint_propagate(): grad_C_phi_2 = 2 * np.real(-1j * intermediate_fields[0].T.conj() * fft_input_array).flatten() fft_operator() - adjoint_field_2_propagated_fft = fft_operator.output_array - ifft_input_array[:,:] = adjoint_field_2_propagated_fft * P_1_dagger_nat #element wise + ifft_input_array[:,:] = fft_operator.output_array * P_1_dagger_nat #element wise ifft_operator() grad_C_phi_1 = 2 * np.real(-1j * input_field_2d.T.conj() * ifft_operator.output_array * Phi_1_dagger).flatten() @@ -301,8 +307,10 @@ def cost_fun(x, grad): if __name__ == "__main__": #initialize nlopt solver solver = nlopt.opt(nlopt.LD_CCSAQ, 2*S) - solver.set_lower_bounds(norm_phase_min) - solver.set_upper_bounds(norm_phase_max) + print("Type of bounds",type(mask_lower_bounds.flatten()), type(mask_upper_bounds.flatten())) + print("Size of bounds",mask_lower_bounds.flatten().shape, mask_upper_bounds.flatten().shape) + solver.set_lower_bounds(mask_lower_bounds.flatten()) + solver.set_upper_bounds(mask_upper_bounds.flatten()) solver.set_max_objective(cost_fun) solver.set_maxeval(opt_max_eval) solver.set_param("dual_ftol_rel", 1e-7) @@ -394,7 +402,7 @@ if __name__ == "__main__": # Plot 5: Phase mask 1 fig, ax = plt.subplots(figsize=(10, 8)) - im = ax.imshow((phase_mask_1+2*np.pi)%(2*np.pi),extent=(-size_x*1e6/2, size_x*1e6/2, -size_y*1e6/2, size_y*1e6/2), origin='lower', cmap='viridis') + im = ax.imshow((phase_mask_1+2*np.pi)%(2*np.pi),extent=(-size_x*1e6/2, size_x*1e6/2, -size_y*1e6/2, size_y*1e6/2), origin='lower', cmap='hsv', vmin=0, vmax=2*np.pi) ax.set_title('Phase Mask 1') ax.set_xlabel('x') ax.set_ylabel('y') @@ -402,11 +410,12 @@ if __name__ == "__main__": plt.colorbar(im, ax=ax) plt.tight_layout() plt.savefig("results/phase_mask_1.pdf") - plt.close() + plt.show() + #plt.close() # Plot 6: Phase mask 2 fig, ax = plt.subplots(figsize=(10, 8)) - im = ax.imshow((phase_mask_2+2*np.pi)%(2*np.pi),extent=(-size_x*1e6/2, size_x*1e6/2, -size_y*1e6/2, size_y*1e6/2), origin='lower', cmap='viridis') + im = ax.imshow((phase_mask_2+2*np.pi)%(2*np.pi),extent=(-size_x*1e6/2, size_x*1e6/2, -size_y*1e6/2, size_y*1e6/2), origin='lower', cmap='hsv', vmin=0, vmax=2*np.pi) ax.set_title('Phase Mask 2') ax.set_xlabel('x') ax.set_ylabel('y') @@ -414,7 +423,8 @@ if __name__ == "__main__": plt.colorbar(im, ax=ax) plt.tight_layout() plt.savefig("results/phase_mask_2.pdf") - plt.close() + plt.show() + #plt.close() # Plot 7: Cuts along y=0 axis - Amplitude comparison y_center_idx = Ny // 2