/*--------------nanotube1.C on February 9, 2000  -------
This is "nanotube1.C" on February 9,2000. This calculates
the molecular coordinates of a (n,m) nanotube. 

Reference:"Physical Properties of Carbon Nanotubes"
Imperial College Press,
by R. Saito, G. Dresselhaus, and M.S. Dresselhaus

 */
#include <stdio.h>
#include <iostream.h>
#include <math.h>
//#include <complex.h>
//#define pi 4*atan(1)
#define true 1
#define false 0
#define acc 1.421
#define nk 5000
#define eps 1.0e-5
//#define eps 0
//#define a 1.421*sqrt(3.)
int GCD(int,int);
int Gen11(int,int);
int p,q,dr,d,r,s;
double dt;
double pi = 4*atan(1);
double a = 1.421*sqrt(3.);
main()
{	
double sq3, x[nk],y[nk],z[nk],r,c,rs,q1,q2,q3,q4,q5,h1,h2,x1,y1,z1,z3
	,x2,y2,z2,l2,l,t,t2;		
int m,n,N,ii,i,ichk,n60,k,kk,kk2,N_u,jj,iii;  
// double_complex ;
//char *line;
//char out_file[]= "/home/keayb/Cprograms/Nanotubes/Molecules/output/tube.xyz";
// const double epo = 8.854e-12;
/*	printf("  time = %6.4g nsec\n", time);      */
FILE *fp1;
fp1 = fopen("/home/keayb/Cprograms/Nanotubes/Molecules/output/tube.xyz","a+");
// if ((fp1 =fopen(out_file,"a+")) == NULL)
//    {
//	printf("Cannot open output file.\n");
//	}
// I = sqrt(double_complex(-1));
printf(" Enter n,m for C_h\n");
scanf("%d %d", &n, &m);	
printf(" C_h = n,m %d %d \n",n,m);
cout << "Enter the number of unit cells N_u" << endl;
cin >> N_u;
//
//
//
Gen11(n,m);	
printf("Made it to here 1!\n");	
cout << "q =  " << q << "  p =  " << p << "  dr =  " << dr << endl;
cout << "a =  " << a << "  acc =  " << acc << "  eps =  " << eps  << "  pi =  " << pi << endl;
sq3=sqrt(3.);
//a = sq3*acc;
//eps = 1.0e-5;
//
// r; |R| ,  c;|C_h|  ,   t;|T|
//
r = a*sqrt(p*p + q*q + p*q);
c = a*sqrt(n*n + m*m + n*m);
t = sqrt(3.)*c/dr;
cout << "r; |R| =  " << r << "  c;|C_h| =  " << c << "  t;|T| =  " << t  << endl << endl << endl;
printf("t = %6.4g    (n*n + m*m + n*m) = %d\n",t,(n*n + m*m + n*m));
//
// N: the number of hexagon in the unit cell N
//
N = 2*(n*n + m*m + m*n)/dr;
if (2*N > nk ) printf("parameter nk is too small");
//
// rs: radiiius of the tube
rs = c/(2.*pi);
cout << "the number of hexagon in the unit cell N =  " << N << " tube radius rs =  " << rs << "  dr =  " << dr << endl;
cout << endl << endl << endl << endl;
cout << "d =  " << d << " dr =  " << dr << "  dt =  " << dt << endl;
cout << "L/a =  " << sqrt(n*n+m*m+m*n) << " T =  " << r <<","<< s << "  T/a =  " << t/a << endl;
cout << "N =  " << N << " R =  " << p <<","<< q << "  M =  " << (m*p - n*q) << endl;
cout << endl << endl << endl << endl;
//
// q1: the chiral angle for C_h
// q2: the 'chiral' angle for R
// q3: the angle between C_h and R
//
q1 = atan( (sq3*m/(2*n + m) ));
q2 = atan( (sq3*q/(2*p + q) ));
q3 = q1 - q2;
cout << "q1 =  " << q1 << "  q2 =  " << q2 << "  q3 =  " << q3 << endl;
//
// q4: a period of an angle for the A atom
// q5: the difference of the angle between the A and B atoms
//
q4  =2.*pi/N;
q5 = acc*cos((pi/6.0) - q1)/c*2.*pi;
//
// h1:
// h2: Delta z between the A and B atoms
//
h1 = fabs(t)/fabs(sin(q3));
h2 = acc*sin(( pi/6.) - q1);
cout << "q4 =  " << q4 << "  q5 =  " << q5 << "  h1 =  " << h1  << "  h2 =  " << h2 << endl;
ii = 0;
	for(i=0;i<= N-1;i++){
	x1 = 0;
	y1 = 0;
	z1 = 0;
	k = i*fabs(r)/h1;
	x1 = rs*cos(i*q4);
	y1 = rs*sin(i*q4);
	z1 = (i*fabs(r) - k*h1)*sin(q3);
	kk2 = fabs(z1/t) + 1;
// 
// Check the A atom is in the unit cell 0 < z1 < t
// 	
	if (z1 > t -.02 ) z1 = z1 -t*kk2;
	if (z1 < -.02 ) z1 = z1 + t*kk2;
	ii = ii + 1;
	x[ii] = x1;
	y[ii] = y1;
	z[ii] = z1;
// 
// The B atoms
//	
    z3 = ( i*fabs(r) - k*h1  )*sin(q3) - h2;
    ii = ii + 1;
// 
// Check the B atom is in the unit cell 0 < z3 < t
//     
	if (( z3 >= -eps ) && (z3 <= t- eps ) ){
		x2 = rs*cos(i*q4 + q5 );
		y2 = rs*sin(i*q4 + q5 );
		z2 = ( i*fabs(r) - k*h1 )*sin(q3) - h2;
		x[ii] = x2;
		y[ii] = y2;
		z[ii] = z2;
	}else{
		x2 = rs*cos(i*q4 + q5 );
		y2 = rs*sin(i*q4 + q5 );
		z2 = ( i*fabs(r) - (k + 1)*h1 )*sin(q3) - h2;
		kk = fabs(z2/t) + 1;
	  if (z2 > t - eps) z2 = z2 - t*kk;
	  if (z2 < - eps) z2 = z2 + t*kk;
	  	x[ii] = x2;
	  	y[ii] = y2;
		z[ii] = z2;
	  }
	  }
cout << "Total number of atoms =  " << N_u*2*N << endl;	  
if (N_u*N > nk ) cout << "nk is too small!"  << endl;	
ii = 2*N;
for(i=1;i<= 2*N;i++){
	for(jj=0;jj<= N_u-1;jj++){
	iii = ii + i + jj*ii;
	z[iii] = z[i] + (jj + 1)*t;
	x[iii] = x[i];
	y[iii] = y[i];
}}
	FILE *fp2;
fp2 = fopen("/home/keayb/Cprograms/Nanotubes/Molecules/output/en.xyz","a+");
fprintf(fp1,"%d \n",N_u*2*N); 	
fprintf(fp1,"%d x %d nanotube \n",n,m); 
	for(i=1;i<= 2*N*N_u;i++){
	if (fabs(x[i]) < eps ) x[i] = 0;
	if (fabs(y[i]) < eps ) y[i] = 0;
	fprintf(fp1,"C  %10.4g  %10.4g  %10.4g\n",x[i],y[i],z[i]); 
   }
fprintf(fp2,"%d \n",2*N); 
fprintf(fp2,"%6.4g  %6.4g \n",t,acc); 	
	for(i=1;i<= N*2;i++){
	fprintf(fp2," %d  %9.4g  %9.4g %9.4g \n",i,x[i],y[i],z[i] );
	printf(" %d  %9.4g  %9.4g %9.4g %6.4g\n",i,x[i],y[i],z[i] );  
   }   
fclose(fp1);	       
fclose(fp2);
 return 0;
}
//
//
//
int Gen11(int n,int m){	
int np[100],nq[100],n60,j1,j2,l2,ichk;
double l,t2,t,N;
d = GCD(n,m);
	printf(" Gen11 d %d \n",d);
if ((n-m)%(3*d) == 0 ) dr = 3*d;
else dr = d;
	printf(" Gen11 dr %d \n",dr);
l2 = n*n + m*m + n*m;
if (l2 <= 0) printf("ERROR l2 <= 0");
l = sqrt(l2) + eps;
dt = a*sqrt(l2)/pi;
cout << "a =  " << a << "  l2 =  " << l2 << "  sqrt(l2) =  " << sqrt(l2)  
<< " pi =  " << pi << " a*sqrt(l2)/pi =  " << sqrt(l2)/pi << "  dt =  " << dt << endl;
// T
	r = (2*m + n)/dr;
	s = -(2*n + m)/dr;
	t2 = 3*l2/dr/dr;
	t = sqrt(t2) + eps;
// N
 	N = 2*l2/dr;
// R
	ichk = 0;
	if(r == 0)	n60 = 1;
	else n60 = r;
	printf(" Gen11 n60 %d \n",n60);
for(p=-abs(n60);p<=abs(n60);p++){
for(q=-abs(n60);q<=abs(n60);q++){
j2 = r*q - s*p;
//		printf(" Gen11 r,q,s,p,j2 %d %d %d %d %d\n",r,q,s,p,j2);
if (j2 == 1){ 
	j1 = m*p - n*q;
	if ( j1 > 0 && j1 < N){
	ichk = ichk + 1;
		printf(" Gen11 ichk %d \n",ichk);
	np[ichk] = p;
	nq[ichk] = q;
	}}
//printf(" Gen11 p, q, ichk %d %d %d\n",p,q,ichk);
	}}
printf(" Gen11 Made it to here 2!\n");
if (ichk == 0) printf("ERROR - ichk == 0, not found p,q strange!!");
	if (ichk >= 2) printf("ERROR - ichk >= 2, not found p,q strange!!");		
	p = np[1];
	q = nq[1]; 	
cout << "q =  " << q << "  p =  " << p << "  dr =  " << dr << endl;	
//
//
//	
 int GCD(int ii,int jj){
 int i,j,iw,gcd,ir;
 	i = abs(ii);
 	j = abs(jj);
 	if (j > i){
 		iw = j;
 		j = i;
 		i = iw;}
 	if (j == 0){
 		gcd = i;
 		}
 	if (j > 0 && i > 0){
 	do {
 	 ir = i%j;
 	 if (ir == 0)
 	    gcd = j;
 	 else {
// 	 printf("GCD ir %d \n",ir);
 	 i = j;
 	 j = ir;
 	 }
 	 } while (ir > 0);
 	}
//gcd = j;
printf("GCD gcd %d \n",gcd);
	return gcd;
}