function [N,Nz,q]=isboundary(ye,s,aa1,bb1,aa2,bb2,aa3,bb3,aa4,bb4);
global  b out
q=0;N=zeros(3,1);Nz=zeros(2,1);
%%%%%%%%%%%%% Check whether we are not the boundary
c1=ye(1,:)>=bb1;
c2=ye(1,:)<=aa1;
c3=ye(2,:)>=bb2;
c4=ye(2,:)<=aa2;
c5=ye(3,:)>=bb3;
c6=ye(3,:)<=aa3;
c7=s(2)>=bb4;
c8=s(2)<=aa4;
%%%%%%%%%%%%%
lim1=[c1,c2,c3,c4,c5,c6,c7,c8]; % which constraints are enabled
lim2=[c1,c2,c3,c4,c5,c6];
p1=find(lim1);p2=find(lim2);
a=sum(lim1);
if a>0  % if at least a boundary is enabled
  out=1;
  N0=[c1 -c2  0  0  0  0  0   0;
	  0   0  c3 -c4 0  0  0   0;
	  0   0  0   0 c5 -c6 c7 -c8];
  if (c5|c6)
	  N=N0(:,p2);
  else
	  N=N0(:,p1);
  end
    if (sum(lim1(7:8))>0) % if one of the two s_2 boundaries is enabled
		NZ=[0 0 1/b^2 -1/b^2  0 0 0 0; % define appropriate normal vectors
	    	0 0  0      0     0 0 0 0];
		if (c5|c6)
			Nz=NZ(:,p2);
		else
			Nz=NZ(:,p1);
		end
	else		
    	NZ=[0 0 1/b^2 -1/b^2  0 0;
			0 0  0 0  1/b^2 -1/b^2];
    	Nz=NZ(:,p2);
    end
q=size(Nz,2);
end