1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
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);
|
