% nanotube drawing script

% lattice generators
r = 1;
a1 = r*cot(pi/6)*[1;0];
a2 = r*cot(pi/6)*[cos(pi/3);sin(pi/3)];
arcdir1 = r*[0;1];
arcdir2 = r*[cos(7*pi/6);sin(7*pi/6)];
arcdir3 = r*[cos(11*pi/6);sin(11*pi/6)];

% set circumference, length of tube and unit cell length
m = 2; 
n = 3;
W = m*a1+n*a2;          % roll-up point on cylinder
C = sqrt(W'*W);         % tube circumference
d = gcd(2*m+n,2*n+m);   % compute unit cell vector
t1 =(m+2*n)/d;          
t2 =-(2*m+n)/d;
T = t1*a1+t2*a2;
ucL = sqrt(T'*T)       % unit cell distance
L1 = max([10,ucL]);     % distance from base point to top of tube
L2 = max([10,ucL]);     % distance from base point to bottom of tube



% set up orthogonal coordinates adapted to tube cylinder
U = W/C;
V = [U(2); -U(1)];
if det([U,V]) < 0, V = -V; end;
a = a1; 
b = a2;
if U'*a < 0, a = -a; b = -b;end;
if V'*b == 0, c = a; a = b; b = c; end; 

% compute lattice coefficients, lattice points
eps = 1/1000; %fudge factor
cLpts = [];
s2 = floor(C/(U'*a)+eps)+1;
for s = 0:s2
    t1 = ceil((-L2-s*V'*a)/(V'*b));
    t2 = floor((L1-s*V'*a)/(V'*b));
    if t1 > t2, t3 = t1; t1 = t2; t2 = t3; end; % ensure that t2 >= t1 
    cLpts = [cLpts,[s*ones(1,t2-t1+1);t1:t2]];
end;
Lpts = [a,b]*cLpts; % lattice points

% compute arc endpoints
Lpts1 = [Lpts(1,:)+arcdir1(1);Lpts(2,:)+arcdir1(2)];
Lpts2 = [Lpts(1,:)+arcdir2(1);Lpts(2,:)+arcdir2(2)];
Lpts3 = [Lpts(1,:)+arcdir3(1);Lpts(2,:)+arcdir3(2)];


% rotate and compute tube points

A = [U,V]';           % create rotation matrix
RLpts = A*Lpts;       % transform points into X-Y coordinates 
RLpts1 = A*Lpts1;     % transform points into X-Y coordinates 
RLpts2 = A*Lpts2;     % transform points into X-Y coordinates 
RLpts3 = A*Lpts3;     % transform points into X-Y coordinates 
RW = A*W;
RT = A*T;
Ra = A*a;
Rb = A*b;
Ra1 = A*a1;
Ra2 = A*a2;



% transform to 3-D tube coordinates
R = C/(2*pi);

%% background cylinder
eps = 0.97;
ndiv = 50;
Rc = eps*R;
Lu = 2*pi*Rc;
Lv = L1+L2;
du = Lu/ndiv;
Uc = 0:du:Lu;
Vc = ((-L1):du:L2)';
mc = length(Vc);
nc = length(Uc);
Xc = ones(mc,1)*(Rc*cos(Uc/Rc));
Yc = ones(mc,1)*(Rc*sin(Uc/Rc));
Zc = Vc*ones(1,nc);

%% lattice points (atoms)
X = R*cos(RLpts(1,:)/R);
Y = R*sin(RLpts(1,:)/R);
Z = RLpts(2,:);

%% arcendpoints (bonds)
X1 = R*cos(RLpts1(1,:)/R);
Y1 = R*sin(RLpts1(1,:)/R);
Z1 = RLpts1(2,:);
X2 = R*cos(RLpts2(1,:)/R);
Y2 = R*sin(RLpts2(1,:)/R);
Z2 = RLpts2(2,:);
X3 = R*cos(RLpts3(1,:)/R);
Y3 = R*sin(RLpts3(1,:)/R);
Z3 = RLpts3(2,:);

%% unit cell boundaries
lgirth = 2*pi*linspace(0,C,100)';
Xuc = R*cos(lgirth/R);
Yuc = R*sin(lgirth/R);
ZucO = ones(size(lgirth));


% draw pictures
figure(1)
plot(cLpts(1,:),cLpts(2,:),'r.')
title('lattice coefficients');

figure(2)
plot(Lpts(1,:),Lpts(2,:),'r.')
axis equal
title('lattice points');

figure(3)
clf
hold on
plot(RLpts(1,:),RLpts(2,:),'r.')
plot(RLpts1(1,:),RLpts1(2,:),'g.')
plot(Ra1(1),Ra1(2),'k.');
plot(Ra2(1),Ra2(2),'k.');
plot(RW(1),RW(2),'k.');
plot(RT(1),RT(2),'k.');
plot([RW(1),0],[RW(2),0],'k-',[Ra1(1),0],[Ra1(2),0],'k-',[Ra2(1),0],[Ra2(2),0],'k-',[RT(1),0],[RT(2),0],'k-');
plot([Ra(1),0],[Ra(2),0],'b-',[Rb(1),0],[Rb(2),0],'b-');


title('lattice points rotated to tube coordinates')
axis equal
hold off

figure(4)
clf
hold on
hc = surf(Xc,Yc,Zc);
set(hc,'EdgeColor','none');
colormap(ones(256,1)*[1,1,0]);
hl1 = plot3([X;X1],[Y;Y1],[Z;Z1], 'r');
hl2 = plot3([X;X2],[Y;Y2],[Z;Z2], 'g');
hl3 = plot3([X;X3],[Y;Y3],[Z;Z3], 'b');
set([hl1,hl2,hl3],'Linewidth',2);
ha = plot3(X,Y,Z,'.k'); 
ha1 = plot3(X1,Y1,Z1,'.k'); 
set([ha,ha1],'MarkerSize',20);
axis equal
axis vis3d
view(0,0)
title(['nanotube: m = ',num2str(m),' n = ',num2str(n)]); 
hold off

figure(5)
clf
hold on
hc = surf(Xc,Yc,Zc);
set(hc,'EdgeColor','none');
colormap(ones(256,1)*[1,1,0]);
hl1 = plot3([X;X1],[Y;Y1],[Z;Z1], 'r');
hl2 = plot3([X;X2],[Y;Y2],[Z;Z2], 'g');
hl3 = plot3([X;X3],[Y;Y3],[Z;Z3], 'b');
set([hl1,hl2,hl3],'Linewidth',2);
ha = plot3(X,Y,Z,'.k'); 
ha1 = plot3(X1,Y1,Z1,'.k'); 
set([ha,ha1],'MarkerSize',20);
hauc1 = plot3(Xuc,Yuc, ucL*ZucO,'-k'); 
hauc2 = plot3(Xuc,Yuc, 0*ZucO,'-k'); 
hauc3 = plot3(Xuc,Yuc,-ucL*ZucO,'-k'); 
set([hauc1,hauc2,hauc3],'Linewidth',3);
axis equal
axis vis3d
view(0,0)
title(['nanotube: m = ',num2str(m),' n = ',num2str(n)]); 
hold off