6 files changed, 149 insertions, 13 deletions
diff --git a/beam_intensity_at_point_from_image.m b/beam_intensity_at_point_from_image.m
new file mode 100644
index 0000000..f305d16
--- /dev/null
+++ b/beam_intensity_at_point_from_image.m
@@ -0,0 +1,35 @@
+function intensity=beam_intensity_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
new file mode 100644
index 0000000..49c2aff
--- /dev/null
+++ b/demo_diffraction_on_a_disk.m
@@ -0,0 +1,43 @@
+% light wavelength
+lambda=795e-9; % Rb D1 line
+
+% set up mask (source image)
+Nx_s=301;
+Ny_s=301;
+img_source=zeros(Ny_s,Nx_s); 
+
+src_size = 1.0*12.5e-3;
+xpos_s=linspace( -src_size, src_size, Nx_s);
+ypos_s=linspace( -src_size, src_size, Ny_s);
+
+[xmesh,ymesh]=meshgrid(xpos_s,ypos_s);
+w=0.0014*3/2;
+R=1;
+z_before_mask=LaguerreGaussianE([0,0,q_(w,R,lambda),lambda],xmesh,ymesh,'cart');
+
+% disk mask
+mask_R = (w*0.33);
+mask = mask_transparency(xmesh,ymesh, mask_R);
+
+% resulting field after mask
+img_source= z_before_mask .* mask;
+%img_source= z_before_mask;
+
+
+% set up target 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;
+
+
+% generating image of the mask accounting the diffraction
+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
+
+% plot source and target images
+show_diffraction(img_source, xpos_s, ypos_s, img_t, xpos_t, ypos_t);
diff --git a/diffracted_image_at_target.m b/diffracted_image_at_target.m
new file mode 100644
index 0000000..afc92fd
--- /dev/null
+++ b/diffracted_image_at_target.m
@@ -0,0 +1,20 @@
+function img_target = diffracted_image_at_target(img_source, xpos_s, ypos_s, image_target, xpos_t, ypos_t, z_t, lambda)
+
+[Ny,Nx]=size(image_target);
+
+xim_t=linspace(xpos_t(1), xpos_t(end), Nx);
+yim_t=linspace(ypos_t(1), ypos_t(end), Ny);
+
+img_target=zeros(Ny, Nx);
+
+for 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;
+    end
+end
+
+
+
+
+
diff --git a/gauss_diffraction_on_aperture.m b/gauss_diffraction_on_aperture.m
index 6f72e85..bd962f6 100644
--- a/gauss_diffraction_on_aperture.m
+++ b/gauss_diffraction_on_aperture.m
@@ -7,7 +7,6 @@
 % the magnitude of the coefficients of the various modes in the
 % decomposition.
 
-
 ploton=[1 0];
 overlaponly=0; showfigure=0;
 
@@ -16,14 +15,16 @@ clear domain;
 screensize=0.0125;
 nptsr=500;
 nptstheta=100;
-accuracy=0.001;
-n=30;
+accuracy=0.0001;
+n=50;
+
+
 
 [rmesh,thetamesh,xmesh,ymesh]=polarmesh([0,screensize,nptsr],[0 2*pi nptstheta],'lin');
 domain(:,:,1)=rmesh; domain(:,:,2)=thetamesh;
 
-w=0.002;
-R=-1e3;
+w=0.0014*3/2;
+R=1;
 lambda=0.795e-6;
 q=q_(w,R,lambda);
 
@@ -33,8 +34,9 @@ deltax=0;
 deltay=0;
 wfactor=1;
 
-mask_R = (w/2);
-mask = 1.0*((xmesh.^2+ymesh.^2) > mask_R^2);
+mask_R = (w*0.33);
+mask = mask_transparency( xmesh, ymesh, mask_R );
+
 %mask = 1.0;
 
 if overlaponly
@@ -51,7 +53,7 @@ if overlaponly
     return
 end
 
-z_before_mask=LaguerreGaussianE([0,0,q_(w*wfactor,R,lambda),lambda],xmesh+deltax,ymesh+deltay,'cart');
+z_before_mask=LaguerreGaussianE([0,0,q_(w,R,lambda),lambda],xmesh+deltax,ymesh+deltay,'cart');
 zin=mask.*z_before_mask;
 clear tmat tmat_in
 tmat_in(1:n,1,1)=1;
@@ -111,7 +113,7 @@ if length(ploton)>=2 & ploton(2)==1
 end
 
 figure(3)
-z_test=LaguerreGaussianE([30,0,q_(w*wfactor,R,lambda),lambda],xmesh+deltax,ymesh+deltay,'cart');
+z_test=LaguerreGaussianE([100,0,q_(w,R,lambda),lambda],xmesh+deltax,ymesh+deltay,'cart');
 h=pcolor(xmesh,ymesh,abs(z_test).^2); set(h,'edgecolor','none'); axis square; colorbar; drawnow; shg;
 title('high order LG');
 
@@ -119,8 +121,8 @@ title('high order LG');
 figure(4)
 [rmesh,thetamesh,xmesh,ymesh]=polarmesh([0,screensize/10,nptsr],[0 2*pi nptstheta],'lin');
 domain(:,:,1)=rmesh; domain(:,:,2)=thetamesh;
-mask = 1.0*((xmesh.^2+ymesh.^2) > mask_R^2);
-z_before_mask_magn=LaguerreGaussianE([0,0,q_(w*wfactor,R,lambda),lambda],xmesh+deltax,ymesh+deltay,'cart');
+mask = mask_transparency( xmesh, ymesh, mask_R );
+z_before_mask_magn=LaguerreGaussianE([0,0,q_(w,R,lambda),lambda],xmesh+deltax,ymesh+deltay,'cart');
 zin_magn=mask.*z_before_mask_magn;
 zout_magn=recompose(domain,'lg',coeffs,[q,lambda,accuracy]);
 subplot(221);
@@ -142,8 +144,8 @@ figure(5)
 clear mask domain;
 [rmesh,thetamesh,xmesh,ymesh]=polarmesh([0,screensize/5,nptsr],[0 2*pi 2],'lin');
 domain(:,:,1)=rmesh; domain(:,:,2)=thetamesh;
-mask = 1.0*((xmesh.^2+ymesh.^2) > mask_R.^2);
-z_before_mask_magn=LaguerreGaussianE([0,0,q_(w*wfactor,R,lambda),lambda],xmesh,ymesh,'cart');
+mask = mask_transparency( xmesh, ymesh, mask_R );
+z_before_mask_magn=LaguerreGaussianE([0,0,q_(w,R,lambda),lambda],xmesh,ymesh,'cart');
 zin_magn=mask.*z_before_mask_magn;
 zout_magn=recompose(domain,'lg',coeffs,[q,lambda,accuracy]);
 
diff --git a/mask_transparency.m b/mask_transparency.m
new file mode 100644
index 0000000..0a858bd
--- /dev/null
+++ b/mask_transparency.m
@@ -0,0 +1,7 @@
+function mask = mask_transparency( xmesh, ymesh, mask_R )
+% mask_R is the radius of the aperture
+
+% solid disk mask
+mask = 1.0*((xmesh.^2+ymesh.^2) > mask_R^2);
+
+end
diff --git a/show_diffraction.m b/show_diffraction.m
new file mode 100644
index 0000000..3421b0f
--- /dev/null
+++ b/show_diffraction.m
@@ -0,0 +1,29 @@
+function show_diffraction(img_source, xpos_s, ypos_s, img_target, xpos_t, ypos_t)
+
+[Ny_s,Nx_s] = size(img_source);
+[Ny_t,Nx_t] = size(img_target);
+
+% plot mask image
+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).^2); 
+colorbar; colormap('gray');
+xlabel('x (m)');
+ylabel('y (m)');
+title('mask image');
+
+% plot diffraction image
+figure(2);
+x_t=linspace(xpos_t(1), xpos_t(end), Nx_t);
+y_t=linspace(ypos_t(1), ypos_t(end), Ny_t);
+% I use gamma correction due to luck of the display dynamic range
+% pixel brightness correction, display dependent
+%gamma=1/4;
+gamma=1;
+imagesc(x_t, y_t, (abs(img_target).^2).^(gamma));  
+colorbar; colormap('gray');
+xlabel('x (m)');
+ylabel('y (m)');
+title('target image');
+