% 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 and length of tube
m = 10; 
n = 1;
W = m*a1+n*a2;         % roll-up point on cylinder
C = sqrt(W'*W);        % tube circumfrence
L1 = 5;                % distance from base point to top of tube
L2 = 5;                % 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; 


% compute lattice coefficients, lattice points
eps = 1/1000; %fudge factor
cLpts = [];
s2 = floor(C/(U'*a)+eps);
for s = 0:s2
    t1 = ceil((-L2-s*V'*a)/(V'*b));
    t2 = floor((L1-s*V'*a)/(V'*b));
    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;
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);

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

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,:);


% 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([RW(1),0],[RW(2),0],'k-',[Ra1(1),0],[Ra1(2),0],'k-',[Ra2(1),0],[Ra2(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
