%  Given the trace of a random walk as  the adjacency matrix
%  G and the vertex locations verts, compute the list of 
%  complementary components of the trace with the unbounded 
%  component listed first and create the adjacency matriix
%  between componets, based on sharing an edg

function [dist,norm_diam]=find_comps(dir_edges,out_list,verts,steps,print_flag)

   N = length(verts);
   DE = size(dir_edges,1);
   % create a list of directed edges

   % for each directed edge, findthe next edge; the
   % starting point of the next edge is the terminal 
   % point of the given edge and is the next edge in 
   % counter clockwise order around this common vertex

   rot=[0,-1;1,0];
   flag=0;
   next=zeros(DE,1);
   for k=1:DE
       v1=dir_edges(k,1);
       v2=dir_edges(k,2);
       xy1=[verts(v1,1);verts(v1,2)];
       xy2=[verts(v2,1);verts(v2,2)];
       delta=xy1-xy2;
       out=out_list(v2,:);
       for j=1:4
          delta=rot*delta;
          for n=1:out(1)
             [k,j,n];
	     w=out(1+n);
             out_edge=out(1+n);
	     w=dir_edges(out_edge,2); 
             xy3=[verts(w,1);verts(w,2)];
	         delta2=xy3-xy2;
	         if delta==delta2
		        next(k)=out_edge;
		        delta=2*delta;
	         end % if 
         end % for n
      end % for j
   end % for k
   
   if length(unique(next))< length(next)
       error=3
       next
       return;
   end % if
   
   % find the loops 
   found=zeros(DE,1);
   counter1=0;
   clear comps;
   clear comp_size;
   for k=1:DE
      if found(k)==0
	 found(k)=1;
	 counter1=counter1+1;
	 comps(k)=counter1;
	 start=k;
	 current=start;
	 comp=k;
	 while next(current)~= start
        comp=[comp,next(current)];
	    current=next(current);
	    comps(current)=counter1;
	    [start,current];
	    found(current)=1;
	 end % while 
     comp_size(counter1)=length(comp);
	 %figure; hold on; axis equal; axis off;
	 %plot_edge_list(dir_edges(comp,:),verts)
     end % if 
   end % for k 

   % plot the trace
     if print_flag > 0
        figure; hold on; axis equal; axis off;
        title({['N = ',num2str(steps),',  The trace'],['  ']})
        plot_edge_list_thick(dir_edges,verts)
     end % if print_flag

   % compute the approximate diameters of the comps
   num_comps=max(comps);
   diam=zeros(num_comps,1);
   for k=1:num_comps
      list=find(comps==k);
      diam(k)=find_comp_diam(list,dir_edges,verts);
   end % for k
   [max_diam,outer_comp]=max(diam);
   norm_diam=diam/diam(outer_comp);

   % find the unbounded component
   if 0==1
      Mx=max(verts(:,1));
      s=find(verts(:,1)==Mx);
      [My,t]=max(verts(s,2));
      NE=s(t);
      out=find(dir_edges(:,1)==NE);
      if length(out)==1
         outer_edge=out;
      else
         out1=out(1);
         v2=dir_edges(out1,2);
         y2=verts(v2,2);
         if y2==My-1
            outer_edge=out1;
         else
            outer_edge=out(2);
         end % if
      end % if
      clear new_list;
      outer_comp=comps(outer_edge);
   end % if 0==1
   
   %  put outer boundary edges in order for using 'fill'
      list =find(comps==outer_comp);
      new_list(1)=list(1);
      for j=2:length(list)
         new_list(j)=next(new_list(j-1));
      end % for j
      
   % plot filled trace    
   if print_flag > 2
      figure; hold on; axis equal; axis off;
      title({['N = ',num2str(steps),',  Filled trace'],['  ']})
      color=[0,0,0];
      plot_comp(new_list,dir_edges,verts,color);
   end % if print_flag

   % plot outer boundary of trace
   if print_flag > 1
      figure; hold on; axis equal; axis off;
      title({['N = ',num2str(steps),',  Outer boundary'],[' ']})
      color=[0,0,0];
      for k=1:length(new_list)
         x(k)=verts(dir_edges(new_list(k),1),1);
         y(k)=verts(dir_edges(new_list(k),1),2);
     end % for k 
     plot(x,y,'k','LineWidth',2);
   end % if print_flag

   num_comps=max(comps);


   % plot the comps
   if print_flag > 4
      setc=flipud(jet(num_comps));
      figure; hold on; axis equal; axis off;
       title({['N = ',num2str(steps),', Number of components = ',num2str(num_comps)],['  ']})
      %plot_edge_list_thick(dir_edges,verts);

      for k=2:num_comps
         list =find(comps==k);
         % sort the componet edges into the correct 
         % order for plotting
         clear new_list;
         new_list(1)=list(1);
         for j=2:length(list)
             new_list(j)=next(new_list(j-1));
         end % for j
         color=setc(k,:);
         plot_comp(new_list,dir_edges,verts,color);
      end % for k
      comp_size;
      plot_edge_list_thick(dir_edges,verts);
   end % if print_flag
   
   % find adjacency between components
   H=sparse(num_comps,num_comps);
   for k=1:2:length(dir_edges)
          H(comps(k),comps(k+1))=1;
          H(comps(k+1),comps(k))=1;
   end % for k 
   %figure
   %spy(H)
   
   % find spanning tree of comps with unbounded
   % component at the root
   outer_comp;
   num_comps;
   [tree,dist]=find_spanning_tree(H,outer_comp);
   depth = max(dist);
   
   % plot components according to distance from 
   % unbounded component
   if print_flag > 0
      setc=flipud(jet(max(dist)));
      figure; hold on; axis equal; axis off;
      title({['N = ',num2str(steps),', Number of components = ',num2str(num_comps), ...
                    ', Depth = ',num2str(max(dist))],['   ']});
      %plot_edge_list_thick(dir_edges,verts);
      for k=1:num_comps
         if k~=outer_comp
            list =find(comps==k);
            % sort the componet edges into the correct 
            % order for plotting
            clear new_list;
            new_list(1)=list(1);
            for j=2:length(list)
               new_list(j)=next(new_list(j-1));
            end % for j
            color=setc(dist(k),:);
            plot_comp(new_list,dir_edges,verts,color);
        end % if 
      end % for k
      %colormap(setc);
      %colorbar
      plot_edge_list_thick(dir_edges,verts);

   end % if print_flag

   % compute the approximate diameters of the comps
   diam=zeros(num_comps,1);
   for k=1:num_comps
      list=find(comps==k);
      diam(k)=find_comp_diam(list,dir_edges,verts);
   end % for k
   norm_diam=diam/diam(outer_comp);
   diam;
   norm_diam;
   max(norm_diam);

   %figure; hold on; grid on;
   %scatter(dist,norm_diam,'filled')

   
return 


