function [is_face_hit, hit_position, hit_distance, new_beams] =  face_beam_interaction(beam, faces)
	%% calculates refracted and reflected beam after interaction with a face
	% beam - structure defining the light beam
	%   beam.origin - array [x,y] origin/soutce of the light beam
	%   beam.k - k vector i.e. direction [kx,ky]
	%   beam.intensity - intensity of the beam
	%   beam.face - if ibeam starts from face then its index is here
	% faces cell array of face structures
	% face - structure definiong the beam
	%  face.vertex1 - [x,y] of the 1st point/vertex of the face
	%  face.vertex2 - [x,y] of the 2nd point/vertex of the face
	%  face.n_left  - index of refraction on the left hand side 
	%                 with respect to 1st -> 2nd vertex direction
	%  face.n_right - index of refraction on the right hand side 


	k=beam.k;

	%% we go over all faces to find the closest which beam hits
	Nfaces=size(faces)(2);
	hit_distance=Inf;
	is_face_hit = false;
	hit_position = [NA, NA];
	closest_face_index=NA;
	for i=1:Nfaces
		if ( beam.face == i) continue; end
		face=faces{i};
		[hit_distance_tmp, hit_position_tmp, is_face_hit_tmp] = beam2face_distance(beam,face);
		if ( hit_distance_tmp < hit_distance ) 
			% this is the closer face
			is_face_hit=is_face_hit_tmp;
			hit_position=hit_position_tmp;
			hit_distance=hit_distance_tmp;
			closest_face_index=i;
		end
	end

	if (!is_face_hit) 
		new_beams={};
		return;
	end

	%% closest face
	face=faces{closest_face_index};
	kf=face.vertex2 - face.vertex1; % not a unit vector

	hold on; 
	% draw face
	t=linspace(0,1);
	x=face.vertex1(1)+kf(1)*t;
	y=face.vertex1(2)+kf(2)*t;
	plot(x,y,'k-');
	t=linspace(0,hit_distance);
	% draw beam
	x=beam.origin(1)+k(1)*t;
	y=beam.origin(2)+k(2)*t;
	plot(x,y,'r-');

	% find is beam arriving from left or right. I will use vectors cross product property.
	% if z component of the product 'k x kf' is positive then beam arrives from the left
	if ( ( k(1)*kf(2)-k(2)*kf(1) ) > 0 )
	       % beam coming from the left
	       n1=face.n_left;
	       n2=face.n_right;
	else	       
	       % beam coming from the right
	       n1=face.n_right;
	       n2=face.n_left;
	end

	% normal vector to the face, looks to the left of it 
	nf=[ kf(2), -kf(1) ] / norm(kf);
	% incidence angle calculation
	cos_theta_i = dot(k, nf) / (norm(k)*norm(nf));
	sin_theta_i = - ( k(1)*nf(2)-k(2)*nf(1) ) / (norm(k)*norm(nf));
	% positive angle to the right from normal before incidence to the face
	theta_i = atan2(sin_theta_i, cos_theta_i);
	
	% reflected beam direction
	theta_normal = atan2(nf(2), nf(1));
	theta_reflected = theta_normal + pi - theta_i;
	
	beam_reflected.origin = hit_position;
	beam_reflected.k = [cos(theta_reflected), sin(theta_reflected)];
	beam_reflected.face=closest_face_index;
	new_beams{1} = beam_reflected;


	% refracted beam direction
	% refracted angle with respect to normal
	sin_theta_refracted_rel2normal = n1/n2*sin(theta_i);
	if ( abs(sin_theta_refracted_rel2normal) >=1 ) 
		% total internal reflection
	else
		% beam refracts
		theta_refracted_rel2normal = asin( sin_theta_refracted_rel2normal );
		theta_refracted = theta_normal + theta_refracted_rel2normal;

		beam_refracted.origin = hit_position;
		beam_refracted.k = [cos(theta_refracted), sin(theta_refracted)];
		beam_refracted.face=closest_face_index;
		new_beams{2} = beam_refracted;
	end
end