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