function n = get_normal_unit_vec(v1,v2)

%Find a unit vector normal to span(v1,v2).
%When v1 and v2 are independent, non-zero vectors,
%the function returns cross(v1,v2)/|cross(v1,v2)|. 
%If they are both zero vectors, it returns [0;0;1].
%If they are proportional vectors, it returns the
%normalized crossproduct of v1 and one of the 
%standard unit vectors.


eps=1e-10; %eps = epsilon

norm_v1 = norm(v1);
norm_v2 = norm(v2);
cross12 = cross(v1,v2);
norm_cross12 = norm(cross12);
if(norm_cross12<eps)
	if(norm_v2>eps)
		v1=[0;0;0];
		norm_v1=0;
	endif
	%if norm_v2<eps, then v2 is already almost zero
endif
if(norm_v1>eps)
	if(norm_v2>eps)
		n = cross12/norm_cross12;
	else%(norm_v2<eps)
		%choose as v2 substitute a
		%vector along x, y or z axis.
		%Best of 3 possiblities is one 
		%along which v1 has smallest projection
		pvec=abs(v1);
		min=1;
		if(pvec(2)<pvec(min))
			min=2;
		endif
		if(pvec(3)<pvec(min))
			min=3;
		endif
		n = get_unit_vec(cross(v1,eye(3)(:,min)));
	endif
else %(norm_v1<eps)
	if(norm_v2>eps)
		pvec=abs(v2);
		min=1;
		if(pvec(2)<pvec(min))
			min=2;
		endif
		if(pvec(3)<pvec(min))
			min=3;
		endif
		n = get_unit_vec(cross(eye(3)(:,min),v2));
	else %(norm_v2<eps)
		n = [0;0;1];
	endif
endif

%tests
%zvec=[0;0;0];
%n=get_normal_unit_vec([1;0;0],[1;0;0])
%n=get_normal_unit_vec(zvec,zvec)
%n=get_normal_unit_vec(zvec,[1;0;1])
%n=get_normal_unit_vec([1;0;0],[0;0;1])