function ln = insdbnd(bndx,bndy,x,y)
%ln = insdbnd(bndx,bndy,x,y)
%ln = insdbnd([meshnod(bndele(:,2),2) meshnod(bndele(:,3),2)], ...
%   [meshnod(bndele(:,2),3) meshnod(bndele(:,3),3)],x,y)
% or
%ln = insdbnd[meshnod(bndele(:,2),2) meshnod(bndele(:,2),3),x,y)
%        In this mode it auto connects last to first
% if ln ==0
% outside
%text(x,y,'outside')
%else
% inside
%text(x,y,'inside')
%end
[m,n]=size(bndx);
if n~=2
bndx = [ [bndx], [bndx(2:m) ; bndx(1)] ];
bndy = [ [bndy], [bndy(2:m) ; bndy(1)] ];
end

xbnd = bndx(:,1)-x +1e-8;
ybnd = bndy(:,1)-y +1e-8;
xbnd2 = bndx(:,2)-x +1e-8;
ybnd2 = bndy(:,2)-y +1e-8;
%plot( xbnd2+x,ybnd2+y,'r',x,y,'ob')
% The origin of the boundary is moved over the point

% Draw a line from center to the right
% Find all of the line segments which straddle it.
% l = find( (xbnd>0 | xbnd2>0 ) & sign(ybnd)~=sign(ybnd2) );

%[xbnd xbnd2 ybnd ybnd2]

lp = find( (sign(ybnd)+sign(ybnd2))==0 ) ;
lz = find(sign(ybnd)==0);
lx = find( xbnd(lp)>=0 & xbnd2(lp)>=0 ) ;  % both x's positive
lc = find( (sign(xbnd(lp))+sign(xbnd2(lp)))==0 ) ; % just one x positive

if ~isempty(lc)
  dy = (ybnd2(lp(lc))-ybnd(lp(lc)))/(xbnd2(lp(lc))-xbnd(lp(lc)))...
      *(0-xbnd(lp(lc))) + ybnd(lp(lc)) - 0   ;
  for i = 1:length(lc)
   if  ( (xbnd(lp(lc(i)))<=0 & ybnd(lp(lc(i)))>=0 & dy(i)>=0) |...
         (xbnd(lp(lc(i)))<=0 & ybnd(lp(lc(i)))<=0 & dy(i)<=0) |...  
         (xbnd(lp(lc(i)))>=0 & ybnd(lp(lc(i)))<=0 & dy(i)>=0) |...     
         (xbnd(lp(lc(i)))>=0 & ybnd(lp(lc(i)))>=0 & dy(i)<=0)  )
    lx = [lx ;lc(i)];
   end
  end
end

l = lp(lx) ;



if ~isempty(lz)
disp([ lz'])
end

ln =  mod(length(l) ,2);
%if(ln==0)
%disp(lp')
%disp([x y])
%end

%if ln ==0
% outside
%text(x,y,'outside')
%else
% inside
%text(x,y,'inside')
%end



