function plot_beams_and_faces_figure(border_limits, img, faces, fig_handle)
	%% plot faces and beams images
	% border_limits has coordinates of left bottom and right top coners
	xlb=border_limits(1);
	ylb=border_limits(2);
	xrt=border_limits(3);
	yrt=border_limits(4);

	[Ny,Nx]=size(img);
	xc=linspace(xlb,xrt, Nx);
	yc=linspace(ylb,yrt, Ny);

	figure(fig_handle);
	%% dummy plot just to put axis in the proper directions
	plot(xlb,ylb,'.', xrt,yrt, '.');
	hold on;
	
	indx_ne_zero = (img != 0);
	indx_eq_zero = (img == 0);
	%% intensity to log scale
	img( indx_ne_zero ) = log10( img(indx_ne_zero) );
	% make sure that the places with no beams get smallest image value
	smallest_intensity = min(min(img(indx_ne_zero)));
	largest_intensity  = max(max(img(indx_ne_zero)));
	img( indx_eq_zero ) = smallest_intensity-(largest_intensity-smallest_intensity)*.4;

	%% colormap
	Ncolors=64;
	R=linspace(0.3,1.0, Ncolors);
	G=R*0;
	B=R*0;
	% no beams will be white
	R(1)=1; G(1)=1; B(1)=1;
	colormap([R',G',B']);
	%% plot ray images
	imagesc(xc,yc, img); colorbar;
	hold on;

	%% plot all faces
	Nf=size(faces)(2); % number of faces
	for i=1:Nf
		hold on;
		t=linspace(0,1,100);
		xf=faces{i}.vertex1(1) + (faces{i}.vertex2(1)-faces{i}.vertex1(1))*t;
		yf=faces{i}.vertex1(2) + (faces{i}.vertex2(2)-faces{i}.vertex1(2))*t;
		plot(xf,yf, 'k-');
	end

	axis('equal');
	hold off;



end