function [processed_beams] = beam_trace(beams, faces, borders, border_limits)
	%% trace beam and all it reflections and refractions on faces
	% for beams and faces structure description see face_beam_interaction.m
	% border similar to faces cell array which enclose the beam
	% i.e. it does not leave the polygon consisting of border faces

	intensity_threshold = 1e-3; % intensity of the weakest beam which we still trace
	processed_beams = {};

	Nbeams=size(beams)(2);
	if ( Nbeams == 0 )
		% no more beam to trace
		return;
	end

	% trace each beam in beams cell array
	for i=1:Nbeams
		beam=beams{i};
		[is_face_hit, hit_position, hit_distance, new_beams] =  face_beam_interaction(beam, faces);
		if (is_face_hit)
			beams{i}.hit_position=hit_position;
			%% remove beams which are to weak to trace
			N_new_beams=size(new_beams)(2);
			intense_enough_beams={};
			intense_beams_counter=0;
			for k=1:N_new_beams
				if ( new_beams{k}.intensity >= intensity_threshold )
					intense_beams_counter=intense_beams_counter + 1;
					intense_enough_beams{intense_beams_counter}=new_beams{k};
				end
			end
			[new_processed_beams] = beam_trace(intense_enough_beams, faces, borders, border_limits);
			processed_beams=cat(2,processed_beams, new_processed_beams);
		else
			% beam does not hit face but it should stop and borders
			beam.face=NA;
			[is_face_hit, hit_position, hit_distance, new_beams] =  face_beam_interaction(beam, borders);
			if (!is_face_hit)
				error('borders are badly defined, the beam misses them');
			end
			beams{i}.hit_position=hit_position;
		end
	end
	processed_beams=cat(2,beams, processed_beams);
end