% Section 1.8 Brin & Struck, the horseshoe map

  close all;
  clear x;
  clear y;
  clear a;
  clear b;
  clear t;
  N=100000;
  x=3*rand(N,1)-1;
  y=rand(N,1);
  figure; hold on; axis equal; 
  scatter(x,y,'.');
   for k=1:N
	[v1,v2,t0]=hs_map(x(k),y(k));
	if t0>0
   	   a(k)=v1;
	   b(k)=v2;
	   t(k)=t0;
        end 
     end % for k

     figure; hold on; axis equal;
     t1=find(t==1);
     scatter(a(t1),b(t1),'k.');
     t1=find(t==2);
     scatter(a(t1),b(t1),'g.');
     t1=find(t==3);
     scatter(a(t1),b(t1),'b.');
     t1=find(t==4);
     scatter(a(t1),b(t1),'r.');
     t1=find(t==5);
     scatter(a(t1),b(t1),'c.');
     t1=find(t>0);
     figure; hold on; axis equal;
     t1=find(t==1);
     scatter(x(t1),y(t1),'k.');
     t1=find(t==2);
     scatter(x(t1),y(t1),'g.');
     t1=find(t==3);
     scatter(x(t1),y(t1),'b.');
     t1=find(t==4);
     scatter(x(t1),y(t1),'r.');
     t1=find(t==5);
     scatter(x(t1),y(t1),'c.');

     t1=find(t>0);
     x=a(t1);
     y=b(t1);
     t=t(t1);

   for iter=1:3
      for k=1:length(x)
	[v1,v2,t0]=hs_map(x(k),y(k));
        x(k)=v1;
        y(k)=v2;
     end % for k
     figure; hold on; axis equal;
     scatter(x,y,'.');

     figure; hold on; axis equal;
     t1=find(t==1);
     scatter(x(t1),y(t1),'k.');
     t1=find(t==3);
     scatter(x(t1),y(t1),'g.');
     t1=find(t==3);
     scatter(x(t1),y(t1),'b.');
     t1=find(t==4);
     scatter(x(t1),y(t1),'r.');
     t1=find(t==5);
     scatter(x(t1),y(t1),'c.');
  end % for iter
