function ind=icefront_maker(index,x,y,node_on_icefront,element_on_boundary);
%ICEFRONT_MAKER - ???
%
%   This program takes element_on_boundary and node_on_icefront 
%   x and  y and index and makes an ordered icefront with it,
%   turning counter clockwise.
%
%   Usage:
%      ind=icefront_maker(index,x,y,node_on_icefront,element_on_boundary);

pos=find(node_on_icefront);
pos2=find(x(pos)==max(x(pos)));
xa=x(pos(pos2));
ya=y(pos(pos2));

temp=node_on_icefront;
temp(pos(pos2))=0;

pos=find(temp);
pos2=find(x(pos)==max(x(pos)));
xb=x(pos(pos2));
yb=y(pos(pos2));



plot(xa,ya,'r*')
plot(xb,yb,'r*')


n00=find(x==xa & y==ya);
n10=find(x==xb & y==yb);
n1=n10;
n0=n00;
cont=[n0 n1];

   
while (n1~=n00),
   
   els=find(index(:,1)==n1 | index(:,2)==n1 | index(:,3)==n1);
   possible_nodes=zeros(2*length(els),1);
   for n=1:length(els),
   other_pos=find(index(els(n),:)~=n1);
   possible_nodes(2*n-1,:)= index(els(n),other_pos(1));
   possible_nodes(2*n,:)=index(els(n),other_pos(2));
	end
potential_node=possible_nodes(1);

n0n1=[x(n1)-x(n0)
   y(n1)-y(n0)];
n0n1_unit=n0n1/norm(n0n1);
ny=flipud(n0n1_unit);
ny(1)=-ny(1);
possible_angles=zeros(length(possible_nodes),1);

for n=1:length(possible_nodes),
   if (possible_nodes(n)~=n0),
	    u=[x(possible_nodes(n))-x(n1)
  	    y(possible_nodes(n))-y(n1)];
  	    ux=sum(u.*n0n1_unit);
  	    uy=sum((u-ux*n0n1_unit).*ny);
         possible_angles(n)=atan2(uy,ux);
      else
         possible_angles(n)=NaN;
   end
end
pos=find(possible_angles==min(possible_angles));
potential_node=possible_nodes(pos(1));

cont=[cont
   [n1 potential_node]];
n0=n1;
n1=potential_node;
end

ind=cont(2:length(cont),1);


