addition of exact diffraction calculation

1 files changed, 35 insertions, 0 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);
+