clear; close all; clc;


addheader=1;

if(~addheader)
	cfile='c.dat';
	hfile='h.dat';
else
	cfile='H2_C.C.ic';
	hfile='H2_C.H2.ic';
end

M=20; % also try 16,20
n=3;


 
  a0 = 1.418;
  a1 = 2.456;
  % a0 is carbon-carbon distance, a1 is next nearest neighbor (a0*sqrt(3)) */

  theta = 2*3.14159265/M;
  r0 = (0.5*a1)/sin(theta/2)
  % r0 is distance from nanotube center to atom */

  densfac=0.0; % 0.5, 2.0 & 0.0 for perfect, overstuff, and understuff
  Lxy=5*r0;
  
  Lz=4.254*n;
  
  Vcnt=pi*r0^2*Lz;
  ncarb=4*M*n;
  
  hdens=densfac*ncarb/Vcnt;
  Vbox=Lxy^2*Lz;
  nhyd=.5*ncarb; % we reduced the density 1/8
  
  nhydin=round(densfac*ncarb);
  nhydout=nhyd-nhydin;
%  nhydout=1.5*nhydin;
  
%  nhyd=nhydin+nhydout;
  
  x=zeros(2*M,1);
  y=zeros(2*M,1);

  
  
   x0 = 0;
  y0 = 0;
  z0 = -Lz/2-4.254;
  % x0, y0, z0 are coordinates of the bottom-center of nanotube */

  
  
	fidc=fopen(cfile,'w');
	fidh=fopen(hfile,'w');
	

cheader=sprintf('RANK %4i %4i %6i\nSIZE %8.4f %8.4f %8.4f\nBEGIN coordinates\n',2,3,ncarb,Lxy,Lxy,Lz)
hheader=sprintf('RANK %4i %4i %6i\nSIZE %8.4f %8.4f %8.4f\nBEGIN coordinates\n',2,3,nhyd,Lxy,Lxy,Lz)
	
if(addheader)
	fprintf(fidc,cheader);
	fprintf(fidh,hheader);
end
	
  
  for i=1:M
  
    x(i) = x0 + r0*cos(i*theta);
    y(i) = y0 + r0*sin(i*theta);
    x(i+M) = x0 + r0*cos((i+0.5)*theta);
    y(i+M) = y0 + r0*sin((i+0.5)*theta);

    %pretty bad way to calculate...*/

	
	

	
    for j=1:n
      fprintf(fidc,'%12.6f %12.6f %12.6f\n', x(i), y(i), z0+3*j*a0);
      fprintf(fidc,'%12.6f %12.6f %12.6f\n', x(i), y(i), z0+(3*j+1)*a0);
      fprintf(fidc,'%12.6f %12.6f %12.6f\n', x(i+M), y(i+M), z0+(3*j+1.5)*a0);
      fprintf(fidc,'%12.6f %12.6f %12.6f\n', x(i+M), y(i+M), z0+(3*j+2.5)*a0);
    end
end



%initialize hydrogen
x=zeros(nhydin,1);
y=zeros(nhydin,1);

r=x;
th=x;

r=r0*rand(nhydin,1);
th=2*pi*rand(nhydin,1);
x=r.*cos(th);
y=r.*sin(th);
z=Lz*rand(nhydin,1)-Lz/2;



for j=1:nhydin
	fprintf(fidh,'%12.6f %12.6f %12.6f\n',x(j),y(j),z(j));
end



% outsize
clear x y z;



i=0;
while(i<nhydout)
	x=Lxy*rand(1)-Lxy/2;
	y=Lxy*rand(1)-Lxy/2;
	z=Lz*rand(1)-Lz/2;
	
	rtmp=sqrt(x^2+y^2);
	if(rtmp>r0)
		i=i+1;
		fprintf(fidh,'%12.6f %12.6f %12.6f\n',x,y,z);
	end
end
	
	
	fclose(fidc);
	fclose(fidh);
	
if(~addheader)
	
	load('h.dat');
	load('c.dat');
	
	figure
	plot3(h(:,1),h(:,2),h(:,3),'r.')
	hold on
	plot3(c(:,1),c(:,2),c(:,3),'b.')
	hold off
	%axis equal
end