% nanotube drawing script

%% generators for lattice of hexagon centres
hd = 1;  % distance between hexagon centers
U1 = hd*[1;0];
U2 = hd*[cos(pi/3);sin(pi/3)];
%% offsets for atoms
bl  = (hd/2)*sec(pi/6); % bondlength
atomup = bl*[cos(pi/6);sin(pi/6)];
atomdown = bl*[cos(pi/6);-sin(pi/6)];
%% vectors representing bonds
bonddir1 = bl*[cos(pi/6),sin(pi/6)];
bonddir2 = bl*[cos(5*pi/6),sin(5*pi/6)];
bonddir3 = bl*[cos(3*pi/2),sin(3*pi/2)];

%%set parameters for nanotube
n = 5; 
m = 7;
ch = 1; ch =ch/abs(ch); % set to ch = -1 for mirror image nanotube

%% compute chiral and unit cell vectors and lengths and tube length
B = n*U1+m*U2;          % chiral vector 
lB = sqrt(B'*B);        % chiral vector length = tube circumference


%% compute unit cell vector
d = gcd(2*m+n,2*n+m);   
t1 =(n+2*m)/d;          
t2 =-(2*n+m)/d;
T = t1*U1+t2*U2;        % unit cell vector 
lT = sqrt(T'*T);        % unit cell length
L1 = 5;                % distance from base point to top of tube
L2 = 5;                % distance from base point to bottom of tube
%L1 = max([L1,1.1*lT]);  % ensure that unit cells can be viewed
%L2 = max([L2,1.1*lT]);  % ensure that unit cells can be viewed


%% set up orthogonal coordinates adapted to tube 
V = B/lB;
W = [V(2); -V(1)];
if det([V,W]) < 0, W = -W; end; 
a = U1; 
b = U2;
if B'*a < 0, a = -a; b = -b;end;
if T'*b == 0, c = a; a = b; b = c; end; 

%% compute lattice coefficients and lattice points
eps = 1/1000; %fudge factor
cHLpts = [];
s2 = floor(lB/(V'*a)+eps)+1;
for s = 0:s2
    t1 = ceil((-L2-s*W'*a)/(W'*b));
    t2 = floor((L1-s*W'*a)/(W'*b));
    if t1 > t2, t3 = t1; t1 = t2; t2 = t3; end; % ensure that t2 >= t1 
    cHLpts = [cHLpts,[s*ones(1,t2-t1+1);t1:t2]];
end;
HLpts = [a,b]*cHLpts; % lattice of hexagon centers

%% compute lattices of attoms points
LAup = [HLpts(1,:)+atomup(1);HLpts(2,:)+atomup(2)];
LAdown = [HLpts(1,:)+atomdown(1);HLpts(2,:)+atomdown(2)];

%% compute bond enpoints offsets from LAup
LBpts1 = [LAup(1,:)+bonddir1(1);LAup(2,:)+bonddir1(2)];
LBpts2 = [LAup(1,:)+bonddir2(1);LAup(2,:)+bonddir2(2)];
LBpts3 = [LAup(1,:)+bonddir3(1);LAup(2,:)+bonddir3(2)];


%% rotate and compute tube points
A = [V,ch*W]';             % create rotation matrix
% transform points into X-Y coordinates  X -> tube lattitude Y -> tube meridian
RLAup = A*LAup;    
RLAdown = A*LAdown;    
RLBpts1 = A*LBpts1;      
RLBpts2 = A*LBpts2;      
RLBpts3 = A*LBpts3;      
RB = A*B;
RT = A*T;
Ra = A*a;
Rb = A*b;
RU1 = A*U1;
RU2 = A*U2;



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

%%% 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);

%%% origin
X0 = R;
Y0 = 0;
Z0 = 0;

%%% lattice points (atoms)
Xu = R*cos(RLAup(1,:)/R);
Yu = R*sin(RLAup(1,:)/R);
Zu = RLAup(2,:);
Xd = R*cos(RLAdown(1,:)/R);
Yd = R*sin(RLAdown(1,:)/R);
Zd = RLAdown(2,:);

%%% bond endpoints 
X1 = R*cos(RLBpts1(1,:)/R);
Y1 = R*sin(RLBpts1(1,:)/R);
Z1 = RLBpts1(2,:);
X2 = R*cos(RLBpts2(1,:)/R);
Y2 = R*sin(RLBpts2(1,:)/R);
Z2 = RLBpts2(2,:);
X3 = R*cos(RLBpts3(1,:)/R);
Y3 = R*sin(RLBpts3(1,:)/R);
Z3 = RLBpts3(2,:);

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


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

figure(2)
plot(HLpts(1,:),HLpts(2,:),'r.')
axis equal
title('lattice of hexagon centres');

figure(3)
clf
hold on
plot(RLAup(1,:),RLAup(2,:),'r.')
plot(RLAdown(1,:),RLAdown(2,:),'g.')
plot(0,0,'k.');
plot(RU1(1),RU1(2),'b.',[RU1(1),0],[RU1(2),0],'b-');
plot(RU2(1),RU2(2),'b.',[RU2(1),0],[RU2(2),0],'b-');
plot(RB(1),RB(2),'k.',[RB(1),0],[RB(2),0],'k-');
plot(RT(1),RT(2),'k.',[RT(1),0],[RT(2),0],'k-');
title('atom lattice points rotated to tube coordinatesshowing basis, chiral, and unit cell vectors')
axis equal
hold off


% set camera parameters
 cdist = 20;
 cvangle = (180/pi)*max([2*atan((L1+L2)/(2*cdist)),atan(3*R/cdist)]);

figure(4)
clf
hold on
hc = surf(Xc,Yc,Zc);
set(hc,'EdgeColor','none');
colormap(ones(256,1)*[1,1,0]);
hb1 = plot3([Xu;X1],[Yu;Y1],[Zu;Z1], 'g');
hb2 = plot3([Xu;X2],[Yu;Y2],[Zu;Z2], 'b');
hb3 = plot3([Xu;X3],[Yu;Y3],[Zu;Z3], 'r');
set([hb1,hb2,hb3],'Linewidth',2);
hau = plot3(Xu,Yu,Zu,'.k'); 
had = plot3(Xd,Yd,Zd,'.k'); 
set([hau,had],'MarkerSize',20);
axis equal
axis vis3d
set(gca,'CameraPosition', [cdist,0,0],'CameraTarget', [0,0,0],'CameraUpVector',[0,0,1], 'CameraViewAngle',cvangle);
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]);
hb1 = plot3([Xu;X1],[Yu;Y1],[Zu;Z1], 'g');
hb2 = plot3([Xu;X2],[Yu;Y2],[Zu;Z2], 'b');
hb3 = plot3([Xu;X3],[Yu;Y3],[Zu;Z3], 'r');
set([hb1,hb2,hb3],'Linewidth',2);
h0 = plot3(X0,Y0,Z0,'*k'); 
hau = plot3(Xu,Yu,Zu,'.k'); 
had = plot3(Xd,Yd,Zd,'.k'); 
set([h0,hau,had],'MarkerSize',20);
huc1 = plot3(Xuc,Yuc, lT*ZucO,'-k'); 
huc2 = plot3(Xuc,Yuc, 0*ZucO,'-k'); 
huc3 = plot3(Xuc,Yuc,-lT*ZucO,'-k'); 
set([huc1,huc2,huc3],'Linewidth',3);
axis equal
axis vis3d
set(gca,'CameraPosition', [cdist,0,0],'CameraTarget', [0,0,0],'CameraUpVector',[0,0,1], 'CameraViewAngle',cvangle);
title(['nanotube: n = ',num2str(n),' m = ',num2str(m)]); 
hold off