diff options
| -rw-r--r-- | beam_field_at_point_from_image.m | 35 | ||||
| -rw-r--r-- | demo_diffraction_on_a_disk.m | 36 | ||||
| -rw-r--r-- | diffracted_image_at_target.m | 4 | ||||
| -rw-r--r-- | show_diffraction.m | 4 |
4 files changed, 69 insertions, 10 deletions
diff --git a/beam_field_at_point_from_image.m b/beam_field_at_point_from_image.m new file mode 100644 index 0000000..8cd4f50 --- /dev/null +++ b/beam_field_at_point_from_image.m @@ -0,0 +1,35 @@ +function E_tot=beam_field_at_point_from_image(x,y,z, img, xpos, ypos, lambda) + +[Ny,Nx] = size(img); + +xim=linspace(xpos(1), xpos(end), Nx); +yim=linspace(ypos(1), ypos(end), Ny); + +[ind_x, ind_y] = meshgrid(1:Nx, 1:Ny); +[rvec_x, rvec_y] = meshgrid(xim-x, yim-y); % rvec_x and rvec_y are matrices of Ny x Nx +rvec_z=z+rvec_x*0; % to get same matrix size as rvec_x and rvec_y + +dist=sqrt(rvec_x.^2+rvec_y.^2+rvec_z.^2); + + +E_ampl = img ./dist .* exp(1i*2*pi.*dist./lambda); + +% simple case: far field - Fraunhofer limit + +E_tot= (sum(E_ampl(:))); +%intensity = E_tot*conj(E_tot); + +%%%% Below is exact case: Fresnel diffraction +%kvec_x= rvec_x./dist; +%kvec_y= rvec_y./dist; +%kvec_z= rvec_z./dist; +% +%Ex= E_ampl * kvec_x; +%Ey= E_ampl * kvec_y; +%Ez= E_ampl * kvec_z; +%% summing all E field components +%Ex_tot = sum (Ex(:)); +%Ey_tot = sum (Ey(:)); +%Ez_tot = sum (Ez(:)); +%intensity = Ex_tot*conj(Ex_tot) + Ey_tot*conj(Ey_tot) + Ez_tot*conj(Ez_tot); + diff --git a/demo_diffraction_on_a_disk.m b/demo_diffraction_on_a_disk.m index 49c2aff..a066da6 100644 --- a/demo_diffraction_on_a_disk.m +++ b/demo_diffraction_on_a_disk.m @@ -12,7 +12,7 @@ ypos_s=linspace( -src_size, src_size, Ny_s); [xmesh,ymesh]=meshgrid(xpos_s,ypos_s); w=0.0014*3/2; -R=1; +R=.3; % positive R is focusing along propagation z_before_mask=LaguerreGaussianE([0,0,q_(w,R,lambda),lambda],xmesh,ymesh,'cart'); % disk mask @@ -24,14 +24,33 @@ img_source= z_before_mask .* mask; %img_source= z_before_mask; -% set up target image +%% set up target image +% square image Nx_t=101; Ny_t=101; img_t=zeros(Ny_t,Nx_t); -tgt_size = 2*12.5e-3; -xpos_t=linspace( -tgt_size, tgt_size, Nx_t); -ypos_t=linspace( -tgt_size, tgt_size, Ny_t); -z_t=10; +tgt_size_x = 12.5e-3; +tgt_size_y = tgt_size_x; +xpos_t=linspace( -tgt_size_x, tgt_size_x, Nx_t); +ypos_t=linspace( -tgt_size_y, tgt_size_y, Ny_t); + +% one quadrant, suitable for symmetrical images +%xpos_t=linspace( 0, tgt_size, Nx_t); +%ypos_t=linspace( 0, tgt_size, Ny_t); + +% slice along x axis for radial simmetry +%ypos_t=linspace( 0, 0, Ny_t); + +%% full image which mimics camera +Nx_t=101; % use Nx_t = 680 to mimic camera +Ny_t=int32(Nx_t*(480/640)); +img_t=zeros(Ny_t,Nx_t); +tgt_size_x = (640*7e-6)/2; % each pixel is 7 um +tgt_size_y = (480*7e-6)/2; +xpos_t=linspace( -tgt_size_x, tgt_size_x, Nx_t); +ypos_t=linspace( -tgt_size_y, tgt_size_y, Ny_t); + +z_t=.50; % generating image of the mask accounting the diffraction @@ -39,5 +58,10 @@ tic; % start timer img_t=diffracted_image_at_target(img_source, xpos_s, ypos_s, img_t, xpos_t, ypos_t, z_t, lambda); toc % show evaluation time +imwrite(img_t/max(img_t(:)),'diffracted_image.png') + % plot source and target images show_diffraction(img_source, xpos_s, ypos_s, img_t, xpos_t, ypos_t); + +% radial slice of the target image +figure(3); plot( xpos_t, abs( img_t(1,:) ) ) diff --git a/diffracted_image_at_target.m b/diffracted_image_at_target.m index f53ac12..90de657 100644 --- a/diffracted_image_at_target.m +++ b/diffracted_image_at_target.m @@ -9,8 +9,8 @@ img_target=zeros(Ny, Nx); parfor i=1:Nx for k=1:Ny - intensity = beam_intensity_at_point_from_image(xim_t(i), yim_t(k), z_t, img_source, xpos_s, ypos_s, lambda); - img_target(k,i) = intensity; + field = beam_field_at_point_from_image(xim_t(i), yim_t(k), z_t, img_source, xpos_s, ypos_s, lambda); + img_target(k,i) = field; end end diff --git a/show_diffraction.m b/show_diffraction.m index 84a44d0..3421b0f 100644 --- a/show_diffraction.m +++ b/show_diffraction.m @@ -7,7 +7,7 @@ function show_diffraction(img_source, xpos_s, ypos_s, img_target, xpos_t, ypos_t figure(1); x_s=linspace(xpos_s(1), xpos_s(end), Nx_s); y_s=linspace(ypos_s(1), ypos_s(end), Ny_s); -imagesc(x_s,y_s, abs(img_source)); +imagesc(x_s,y_s, abs(img_source).^2); colorbar; colormap('gray'); xlabel('x (m)'); ylabel('y (m)'); @@ -21,7 +21,7 @@ y_t=linspace(ypos_t(1), ypos_t(end), Ny_t); % pixel brightness correction, display dependent %gamma=1/4; gamma=1; -imagesc(x_t, y_t, (abs(img_target)).^(gamma)); +imagesc(x_t, y_t, (abs(img_target).^2).^(gamma)); colorbar; colormap('gray'); xlabel('x (m)'); ylabel('y (m)'); |