%% set up lattice and generators
a = 1;
s  = sqrt(3)*a;
U1 = s*[cos(0), sin(0)];
U2 = s*[cos(pi/3), sin(pi/3)];
U3 = U2-U1;
U4 = -U1;
U5 = -U2;
U6 = -U3;

%% set up chiral vector C and unit cell vector T and vector K
n0 = 7, m0 = 9
nm = [n0,m0];
u0 = 3; v0 = 4;
uv = [u0,v0];
check = det([uv; nm])

% limits on number of hexagons in  x and y directions
 N1 = -20; N2 = 20;
 M1 = -20; M2 = 20;

C = n0*U1+m0*U2;
d = gcd(m0,n0);
t1 = -(n0+2*m0)/d;
t2 = (2*n0+m0)/d;
T = t1*U1+t2*U2;
l2C = C*C';

%% set up hexagons, plotting window, hexagons and cylinder lines


 
% hexagons 
 t = pi/6+2*pi*(0:6)/6;
 hex = a*[cos(t);sin(t)];
 hX = hex(1,:);
 hY = hex(2,:);
 shex = 0.75*hex;
 shX = shex(1,:);
 shY = shex(2,:);
 
 % plotting window
 fs = 5;
 minx = round(s*N1-fs);   maxx = round(s*N2+fs);
 miny = round(sin(pi/3)*s*M1-fs); maxy = round(sin(pi/3)*s*M2+fs);
 
 % cylinder lines
 TLminx = [minx,-minx*C(1)/C(2)];
 TLmaxx = [maxx,-maxx*C(1)/C(2)];
 TLminy = [-C(2)*miny/C(1),miny];
 TLmaxy = [-C(2)*maxy/C(1),maxy];
 if TLminx(2) < miny 
     T1 = TLminx; 
 else 
     T1 = TLminy;
 end;
 if TLmaxx(2) > maxy 
     T2 = TLmaxx; 
 else 
     T2 = TLmaxy;
 end;

 TLCminx = [minx,(l2C-minx*C(1))/C(2)];
 TLCmaxx = [maxx,(l2C-maxx*C(1))/C(2)];
 TLCminy = [(l2C-C(2)*miny)/C(1),miny];
 TLCmaxy = [(l2C-C(2)*maxy)/C(1),maxy];
 if TLCminx(2) < miny 
     TC1 = TLCminx; 
 else 
     TC1 = TLCminy;
 end;
 if TLCmaxx(2) > maxy 
     TC2 = TLCmaxx; 
 else 
     TC2 = TLCmaxy;
 end;

%% draw hexagons chiral and unit cell vectors
  
 figure(1)
 clf
  hold on
 for m = M1:2:M2
     k = -m/2;
     for n=(k+N1):(k+N2)
         U = n*U1+m*U2;
         plot(hX+U(1),hY+U(2));
     end
     k = -m/2-1;
     for n=(k+N1):(k+N2)
         U = n*U1+(m+1)*U2;
         plot(hX+U(1),hY+U(2));
     end
 end
 axis equal 
 plot([0,C(1)],[0,C(2)],'-ko','LineWidth',2);
 plot([0,T(1)],[0,T(2)],'-ko','LineWidth',2);
 plot([T1(1),T2(1)],[T1(2),T2(2)],'-r');
 plot([TC1(1),TC2(1)],[TC1(2),TC2(2)],'-r');
 
 
 
 %% start plotting hexagons in order
 dt =0.00000001;
 u = 0; v= 0;
 for j = 0:200     
     U = u*U1+v*U2;
     z = floor(U*C'/l2C);
     u=u-z*n0; v=v-z*m0;
     U = U-z*C;
     patch(shX+U(1),shY+U(2),'r')    
     u =u+u0; v=v+v0;
  pause(dt)
 end;
 u = 0; v= 0;
 for j = 0:200     
     U = u*U1+v*U2;
     z = floor(U*C'/l2C);
     u=u-z*n0; v=v-z*m0;
     U = U-z*C; 
     patch(shX+U(1),shY+U(2),'r')    
     u =u-u0; v=v-v0;
  pause(dt)
 end;

 hold off