added field calculation for diffracted image

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)');