// ===============================================================================
// TubeGen.cpp																						   ©2001-2002, J.T. Frey
// ===============================================================================
// Written:		J.T. Frey, 07.26.2001
// Purpose:		Main program file for the TubeGen app.  This file contains the
//						routines which generate the atomic basis for a unit cell and the
//						k-point list for a unit cell.  Also contains the main program
//						which processes user commands.
//
// Last Mod:	20.FEB.2002: Switch to hexagonal unit cell; k-point generator is
//												 no longer necessary.
//						18.NOV.2002: Modifications to position generation loop [3.0.5]
//						03.DEC.2002: Cell expansion/contraction to preserve rolled C-C
//												 bond lengths added as an option [3.0.6]

//  Include system headers:
#include <fstream>
#include <iostream>
#include <iomanip>

#include <stdio.h>
#include <stdlib.h>
#include <ctype.h>
#include <math.h>

//  Include some of our own headers:
#include "CrystalCell.h"
#include "CoreConst.h"
#include "Constants.h"

#include "ANSR.h"
#include "eprintf.h"

//

  //  The commands can only be this long:
  #define CMD_BUFFER_LEN		256
  
  //  For boolean stuff:
  #ifndef YES
    #define YES 1
  #endif
  #ifndef NO
    #define NO 0
  #endif
	
	#define			ANGSTROM_PER_BOHR					0.529

//
	
	enum {
	  TVector_a1 = 0,
		TVector_a2 = 1,
		TVector_aB = 2
	};
	
//

	typedef struct STubuleGlobals {
		double							aCC;
		TVector3D						A[3];
		double							lenA[3];
		TVector3D						gutter;
		int									replicate[3];
	} STubuleGlobals;
	
//

	typedef struct STubuleOptions {
		int									oUnits				: 1;
		int									oFormat				: 4;
		int									oLattice			: 2;
		int									oRelaxTubule	: 1;
		int									oVerbose			: 1;
		double							oRadiusConv;
		double							oErrorConv;
		double							oGammaConv;
	} STubuleOptions;
	
	enum {
	  EUnits_Angstrom			= 0,
		EUnits_Bohr					= 1
	};
	
	enum {
	  EFormat_Gaussian		= 0,
		EFormat_WIEN				= 1
	};
	
	enum {
	  ELattice_Hexagonal	= 0,
		ELattice_Cubic 			= 1,
		ELattice_Planar			= 2
	};
	
	#define RADIUS_CONVERGENCE			1e-8
	#define ERROR_CONVERGENCE				1e-8
	#define GAMMA_CONVERGENCE				1e-12
	
	#define SET_DEFAULT_OPTIONS(S)	(S).oUnits				= EUnits_Angstrom;		\
																	(S).oFormat				= EFormat_Gaussian;		\
																	(S).oLattice			= ELattice_Hexagonal;	\
																	(S).oRelaxTubule	= 1;									\
																	(S).oVerbose			= 1;									\
																	(S).oRadiusConv		= RADIUS_CONVERGENCE;	\
																	(S).oErrorConv		= ERROR_CONVERGENCE;	\
																	(S).oGammaConv		= GAMMA_CONVERGENCE;

//

	typedef struct STubuleParameters {
		STubuleGlobals*			globals;
		STubuleOptions			options;
		
		int									n,m;
		int									nprime,mprime;
		int									d,dR;
		
		double							gamma[3];
		double							lenCh,lenT;
		TVector3D						Ch,T;
		double							r;
	} STubuleParameters;

//

	//  We use an ANSRDB initialized from /usr/share/periodic.table
	//  which, if the file isn't there, will just be the default
	//  anyway:
	ANSRDB					gPeriodicTable("/usr/share/periodic.table");

//

  // ----------------------------------------------------------------------
  // * mungeWhitespace
  // ----------------------------------------------------------------------
  // Updated:		Jeff Frey, 07.26.2001
  // Purpose:		Removes white space (space, tab, newline characters)
  //						from the head and tail of a string.
  //
  // Last Mod:	n/a
  
  void mungeWhitespace(char*);
  
  void
  mungeWhitespace(
    char*		buffer
  )
  {
    char*		p1 = buffer;
    char*		p2 = buffer + strlen(buffer);
		size_t	l;
		
    while (isblank(*p1)) p1++;
		while (isblank(*(--p2)));
		l = p2 - p1 + 1;
    memmove(buffer,p1,l);
		buffer[l] = '\0';
  }
  
//

  // ----------------------------------------------------------------------
  // * gcd
  // ----------------------------------------------------------------------
  // Updated:	Jeff Frey, 03.19.2001
  // Purpose:	Greatest common divisor of two integers.
  //
  // Last Mod:	n/a
  
  int gcd(
    int i,
    int j
  )
  {
    int m,n,r;
    
    (i < j) ? (m = j,n = i) : (m = i,n = j);
    
    if (n == 0)
      return m;
      
    while (r = m % n) {
      m = n;
      n = r;
    }
    
    return n;
  }
	
//

	double
	calculateError(
		STubuleParameters*	params,
		int									vec,
		double							theta
	)
	{
		double			accum;
		double			tmp1,tmp2;
		
		tmp1 = params->globals->lenA[vec];
		accum = tmp1 * tmp1;			//  lenA[i]^2
		
		tmp1 = params->r;
		tmp1 *= tmp1;													//  r^2
		tmp1 *= 2.0;													//  2 r^2
		tmp2 = 1.0 - cos(theta);							//  1 - cos(theta)
		accum -= tmp1 * tmp2;			//  lenA[i]^2 - 2 r^2 (1 - cos(theta))
		
		tmp1 = params->gamma[vec];
		tmp1 *= tmp1;													//  gamma[i]^2
		tmp2 = Vector3D_Dot(&(params->globals->A[vec]),&(params->T)) / params->lenT;
		tmp2 *= tmp2;
		accum -= tmp1 * tmp2;			//  lenA[i]^2 - 2 r^2 (1 - cos(theta)) - gamma[i]^2 (<a[i],T>/|T|)^2
		
		return accum;
	}
	
//

	double
	calculateErrorDerivative(
		STubuleParameters*	params,
		int									vec,
		double							theta
	)
	{
		double			accum;
		double			tmp1,tmp2,tmp3;
		
		tmp1 = params->r;
		tmp1 *= tmp1;													//  r^2
		tmp1 *= 8.0 * M_PI;										//  8 pi r^2
		tmp1 *= sin(theta);										//  8 pi r^2 sin(theta)
		tmp3 = 1.0 / params->lenCh;
		tmp3 *= tmp3;													//  1 / |Ch|^2
		tmp2 = Vector3D_Dot(&(params->globals->A[vec]),&(params->Ch)) * tmp3;
		accum = -tmp1 * tmp2;
		
		tmp1 = 4.0 * params->gamma[vec];
		tmp2 = Vector3D_Dot(&(params->globals->A[vec]),&(params->T)) / params->lenT;
		tmp2 *= tmp2;
		accum -= tmp1 * tmp2;
		
		accum *= calculateError(params,vec,theta);
		
		return accum;
	}

//

	double
	calculateTheta(
		STubuleParameters*	params,
		int									vec
	)
	{
		double				tmp1 = 2.0 * M_PI * params->gamma[vec];
		double				tmp2 = Vector3D_Dot(&(params->globals->A[vec]),&(params->Ch));
		double				tmp3 = 1.0 / params->lenCh;
		
		tmp3 *= tmp3;
		return tmp1 * tmp2 * tmp3;
	}

//

	void
	setupGlobals(
		STubuleGlobals*		globals,
		double						ccBond,
		double						gutterx,
		double						guttery,
		double						gutterz,
		int								replicatex,
		int								replicatey,
		int								replicatez
	)
	{
		//  Set the bond length:
		if (ccBond > 0.0) {
			globals->aCC = ccBond;
			
			//  Construct a1:
			Vector3D_Set(&(globals->A[TVector_a1]),1.5 * ccBond,0.5 * sqrt(3.0) * ccBond,0.0);
			globals->lenA[TVector_a1] = Vector3D_Magnitude(&(globals->A[TVector_a1]));
			
			//  Construct a2:
			Vector3D_Set(&(globals->A[TVector_a2]),1.5 * ccBond,-0.5 * sqrt(3.0) * ccBond,0.0);
			globals->lenA[TVector_a2] = Vector3D_Magnitude(&(globals->A[TVector_a2]));
			
			//  Construct a3 (bond vector):
			Vector3D_Set(&(globals->A[TVector_aB]),ccBond,0.0,0.0);
			globals->lenA[TVector_aB] = ccBond; //Vector3D_Magnitude(&(globals->A[TVector_aB]));
    
			Vector3D_Rezero(&(globals->A[TVector_a1]),FLT_EPSILON);
			Vector3D_Rezero(&(globals->A[TVector_a2]),FLT_EPSILON);
			Vector3D_Rezero(&(globals->A[TVector_aB]),FLT_EPSILON);
		}
		
		//  Set the gutters:
		if (gutterx >= 0.0)
		  globals->gutter.x = gutterx;
		if (guttery >= 0.0)
		  globals->gutter.y = guttery;
		if (gutterz >= 0.0)
		  globals->gutter.z = gutterz;
		
		//  Set the replication:
		globals->replicate[0] = (replicatex > 0)?(replicatex):(1);
		globals->replicate[1] = (replicatey > 0)?(replicatey):(1);
		globals->replicate[2] = (replicatez > 0)?(replicatez):(1);
	}
	
	#define SET_DEFAULT_GLOBALS(S)		setupGlobals(&S,1.421,1.6735,1.6735,0.0,1,1,1)

//

	void
	setupVectors(
		STubuleParameters*	params
	)
	{
		double			tmp;
		
		//  Calculate Ch:
		tmp = (double)params->n * params->gamma[TVector_a1];
		Vector3D_Scalar(&(params->globals->A[TVector_a1]),tmp,&(params->Ch));
		tmp = (double)params->m * params->gamma[TVector_a2];
		Vector3D_ScaledSum(&(params->Ch),tmp,&(params->globals->A[TVector_a2]),&(params->Ch));
		Vector3D_Rezero(&(params->Ch),FLT_EPSILON);
		//  Store the length of Ch:
		params->lenCh = Vector3D_Magnitude(&(params->Ch));
		
		//  Calculate T:
		tmp = (double)params->nprime * params->gamma[TVector_a1];
		Vector3D_Scalar(&(params->globals->A[TVector_a1]),tmp,&(params->T));
		tmp = -( (double)params->mprime * params->gamma[TVector_a2] );
		Vector3D_ScaledSum(&(params->T),tmp,&(params->globals->A[TVector_a2]),&(params->T));
		Vector3D_Rezero(&(params->T),FLT_EPSILON);
		//  Store the length of T:
		params->lenT = Vector3D_Magnitude(&(params->T));
		
		//  Calculate r:
		params->r = 0.5 * M_1_PI * params->lenCh;
	}
	
//

	void
	optimizeGammas(
		STubuleParameters*	params
	)
	{
		int									rIters = 0;
		double							delta_r;
		
		do {
			int								vec;
			
			//  Optimize gamma[i]:
			for ( vec = 0 ; vec < 3 ; vec++ ) {
				double					theta,error,dgamma = 1.0;
				int							gIters = 0;
				
			  theta = calculateTheta(params,vec);
				error = calculateError(params,vec,theta);
				error *= error;
				while ((error > params->options.oErrorConv) && (fabs(dgamma) > params->options.oGammaConv)) {
					double		denom = 1.0 / calculateErrorDerivative(params,vec,theta);
					
					dgamma = 0.5 * error * denom;
					params->gamma[vec] -= dgamma;
					theta = calculateTheta(params,vec);
					error = calculateError(params,vec,theta);
					error *= error;
					gIters++;
				}
			}
			
			//  Calculate new Ch, T, and r:
			delta_r = params->r;
			setupVectors(params);
			delta_r = params->r - delta_r;
			rIters++;
		} while (fabs(delta_r) > params->options.oRadiusConv);
		
		//  Print information if necessary:
		if (params->options.oVerbose) {
			printf("  -- Relaxing tubule to appropriate bond lengths ------------------\n");
			printf("  Iterative relaxation of rolled structure converged.\n");
			printf("    Convergence of radius to %1.0g and expansion ratios to %1.0g",	\
					params->options.oRadiusConv,params->options.oGammaConv \
			);
			if (rIters > 0) {
				printf(" in %d cycle",rIters);
				if (rIters > 1)
					printf("s\n");
				else
					printf("\n");
			} else {
				printf(" immediately.\n");
			}
			printf("    gamma1  = %lf [graphite basis a1 scaling]\n",params->gamma[0]);
			printf("    gamma2  = %lf [graphite basis a2 scaling]\n",params->gamma[1]);
			printf("    gamma3  = %lf [graphite basis bond scaling]\n\n",params->gamma[2]);
			printf("  NOTE:  Tubule radius and length will be modified to adjust for\n");
			printf("         carbon-carbon bond compression during rolling.\n");
			printf("  -----------------------------------------------------------------\n\n");
		}
	}
	
//

	void
	setupParameters(
		STubuleParameters*	params,
		int									nval,
		int									mval
	)
	{
		//  Set chirality indices:
		params->n			= nval;
		params->m			= mval;
		
		//  Calculate d and dR:
		params->d = gcd(nval,mval);
		if (((nval - mval) % (3 * params->d)) == 0)
      params->dR = 3 * params->d;
    else
      params->dR = params->d;
		
		//  Calculate nprime and mprime:
		params->nprime = (2 * mval + nval) / params->dR;
		params->mprime = (2 * nval + mval) / params->dR;
		
		//  Gammas to 1.0:
		params->gamma[0] = 1.0;
		params->gamma[1] = 1.0;
		params->gamma[2] = 1.0;
		
		//  Setup the vectors:
		setupVectors(params);
		
		//  Optimize the gammas if necessary:
		if (params->options.oRelaxTubule)
		  optimizeGammas(params);
		
		//  Display information about the cell now, if necessary:
		if (params->options.oVerbose) {
			double			conv;
			int					N;
			TVector3D		v;
			
			printf("  -- Summary of Graphitic Basis -----------------------------------\n");
			switch (params->options.oUnits) {
			  case EUnits_Bohr:
				  conv = 1.0 / ANGSTROM_PER_BOHR;
					printf("  All units are in bohr.\n");
					break;
				
				default:
				  conv = 1.0;
					printf("  All units are in angstrom.\n");
			}
			
			//  Calculate the number of hexagonal sub-cells:
			N = 2 * (mval * mval + nval * nval + nval * mval) / params->dR;
			
			if (N > 1)
				printf("  Lattice consists of %d hexagonal sub-cells.\n",N);
			else
				printf("  Lattice consists of a single hexagonal sub-cell.\n");
			printf("  Nearest neighbor carbon-carbon distance as:  %lg\n",params->globals->aCC * conv);
			printf("  Basis vectors constructed using n = %d and m = %d:\n",nval,mval);
			v = params->globals->A[TVector_a1];
			printf("    a1 = %lg < %lg , %lg >  ",params->gamma[TVector_a1],v.x * conv,v.y * conv);
			v = params->globals->A[TVector_a2];
			printf("a2 = %lg < %lg , %lg >\n",params->gamma[TVector_a2],v.x * conv,v.y * conv);
			printf("  Chiral vector Ch constructed as %d(a1) + %d(a2):\n",nval,mval);
			v = params->Ch;
			printf("    Ch = < %lg , %lg >, |Ch| = %lg\n",v.x * conv,v.y * conv,params->lenCh * conv);
			printf("  Tubule radius is: %f\n",params->r * conv);
			printf("  Tubule translation vector T constructed as %d(a1) - %d(a2):\n",params->nprime,params->mprime);
			v = params->T;
			printf("    T = < %lg , %lg >, |T| = %lg\n",v.x * conv,v.y * conv,params->lenT * conv);
			printf("  Check for orthogonality:  T.Ch = %lg\n",Vector3D_Dot(&v,&(params->Ch)) * conv * conv);
		
			printf("  -----------------------------------------------------------------\n\n");
		}
	}
	
//

  // ----------------------------------------------------------------------
  // * generateUnitCell
  // ----------------------------------------------------------------------
  // Updated:	Jeff Frey, 12.06.2002
  // Purpose:	Creates a CrystalCell object containing the nanotube lattice.
	//					Note that this function expects the params structure to have
	//					globals and options already set-up.
  //
  // Last Mod:	n/a
  
  CrystalCell*
  generateUnitCell(
		STubuleParameters*	params,
		int									nval,
		int									mval
	)
  {
    CrystalCell*			cell = NULL;
    double						a,b,c;
    double						r,lenV;
    int								i,j,k;
		int								iMax,jMax,iMin,jMin;
    TVector3D					v;
		TVector3D					center;
    TPoint3D					p,p1;
		
		//  Display lattice-type information:
		if (params->options.oVerbose) {
			switch (params->options.oLattice) {
				case ELattice_Cubic:
					printf("  Producing rolled, cubic nanotube lattice.\n\n");
					break;
					
				case ELattice_Planar:
					printf("  Producing planar nanotube lattice.\n\n");
					break;
					
				default:
					printf("  Producing rolled, hexagonal nanotube lattice.\n\n");
					break;
			}
		}
		
		//  Setup the parameters for the cell:
		setupParameters(params,nval,mval);
		
		//  Create the crystal cell:
		switch (params->options.oLattice) {
			case ELattice_Cubic:
			case ELattice_Hexagonal: {
				a = 2.0 * (params->r + params->globals->gutter.x);
				b = 2.0 * (params->r + params->globals->gutter.y);
				c = params->lenT + 2.0 * params->globals->gutter.z;
				if (params->options.oLattice == ELattice_Hexagonal)
				  cell = new CrystalCell(a,b,c,90.0,90.0,120.0);
				else
				  cell = new CrystalCell(a,b,c,90.0,90.0,90.0);
				center = cell->GetRealBasisVector1();
				v = cell->GetRealBasisVector2();
				Vector3D_Scalar(&center,0.5,&center);
				Vector3D_ScaledSum(&center,0.5,&v,&center);
				break;
			}
			
			case ELattice_Planar: {
				a = params->lenCh + 2.0 * params->globals->gutter.x;
				b = params->lenCh + 2.0 * params->globals->gutter.y;
				c = params->lenT + 2.0 * params->globals->gutter.z;
				cell = new CrystalCell(a,b,c,90.0,90.0,90.0);
				break;
			}
		}
		
		//  Begin generating coordinates:
		iMin = (params->nprime < 0)?(params->nprime):(0);
		iMax = ((params->n + params->nprime) > params->n)?(params->n + params->nprime):(params->n);
		jMin = (-params->mprime > 0)?(-params->mprime):(0);
		jMax = ((params->m + params->mprime) > params->m)?(params->m + params->mprime):(params->m);
		for ( i = iMin ; i <= iMax ; i++ ) {
			for ( j = jMin ; j <= jMax ; j++) {
        //  And finally, we loop over the two atoms in the
        //  hexagonal graphite basis, giving us
        //  i(a1) + j(a2)   and   i(a1) + j(a2) + <C-C,0,0>
        for ( k = 0 ; k < 2 ; k++ ) {
        
          //  Construct i(a1) + j(a2):
          v.x = (double)i * params->gamma[TVector_a1] * params->globals->A[TVector_a1].x + \
								(double)j * params->gamma[TVector_a2] * params->globals->A[TVector_a2].x;
          v.y = (double)i * params->gamma[TVector_a1] * params->globals->A[TVector_a1].y + \
								(double)j * params->gamma[TVector_a2] * params->globals->A[TVector_a2].y;
          v.z = 0.0;
          
          //  Second time through we add a C-C bond displacement:
          if (k == 1)
						Vector3D_ScaledSum(&v,params->gamma[TVector_aB],&(params->globals->A[TVector_aB]),&v);
					Vector3D_Rezero(&v,FLT_EPSILON);
          
          lenV = Vector3D_Magnitude(&v);
          
          //  Check v; if it's a zero vector that's really easy; otherwise
          //  we need to project onto Ch and T to get fractional coordinates
          //  along those axes.
          p.y = 0.5;
          if (lenV < FLT_EPSILON) {
            p.x = p.z = 0.0;
          } else {
            p.x = Vector3D_Dot(&v,&(params->Ch)) / (params->lenCh * params->lenCh);
            p.z = Vector3D_Dot(&v,&(params->T)) / (params->lenT * params->lenT);
            if (fabs(p.x) < FLT_EPSILON) p.x = 0.0;
            if (fabs(p.z) < FLT_EPSILON) p.z = 0.0;
          }
          //  If point "p" is within [0,1) in x and z, we have a point:
          if ((p.x < 1.0) && (p.x >= 0.0) && (p.z < 1.0) && (p.z >= 0.0)) {
            //  Check if we're too close to 1.0:
            if ((1.0 - p.x > FLT_EPSILON) && (1.0 - p.z > FLT_EPSILON)) {
              //  This is the rolled- vs. flat-specific stuff:
							if ((params->options.oLattice == ELattice_Cubic) ||
									(params->options.oLattice == ELattice_Hexagonal)) {
                //  theta = 2(pi) times displacement along chiral vector:
                double		theta = k2Pi * p.x;
								
								//  Redefine the point as a polar coordinate in xy-plane:
								p.x = params->r * cos(theta) + center.x;
								p.y = params->r * sin(theta) + center.y;
								p.z *= params->lenT;
								cell->DidAddAtomAtCartesianPoint(6,p);
              } else {
                p.x = (p.x * params->lenCh + params->globals->gutter.x) / a;
                p.z = (p.z * params->lenT + params->globals->gutter.z) / c;
								cell->DidAddAtomAtFractionalPoint(6,p);
              }
            }
          }
        }
      }
    }
		
    if (params->options.oVerbose) {
		  if (cell)
        printf("  Cell generation complete.  %d basis points defined.\n",cell->GetBasisCount());
			else
			  printf("  Cell generate failed!\n");
		}
    
    return cell;
	}
	
//

  // ----------------------------------------------------------------------
  // * main
  // ----------------------------------------------------------------------
  // Updated:	Jeff Frey, 07.26.2001
  // Purpose:	The entry point of the program.  This is one big long
  //					function which may be slightly bloated, but factoring it
  //					out into "chunks" and making the design a bit more
  //					modular would be of no true benefit.  Essentially, the
  //					program processes commands coming from a "comFile" while
  //					may be an actual file or may be the standard input stream.
  //
  // Last Mod:	n/a

int main (int argc, const char * argv[])
{
    char							command[CMD_BUFFER_LEN];
    FILE*							comFile = stdin;
    CrystalCell*			theCell = NULL;
    int								printPrompt = YES;
		int								echoCommands = NO;
		unsigned					i = 1;
		TVector3D					v;
		STubuleGlobals		globals;
		STubuleParameters	params;
		int								n = 3,m = 3;
		
		//  Set up the options, etc:
		SET_DEFAULT_GLOBALS(globals);
		SET_DEFAULT_OPTIONS(params.options);
		params.globals = &globals;
		
		//  If any argument is '--quiet' then we redirect stdout to
		//  /dev/null:
		while (i < argc) {
		  if (strcmp(argv[i],"--quiet") == 0) {
				freopen("/dev/null","w",stdout);
				freopen("/dev/null","w",stderr);
			}
			i++;
		}
		
		//  Are we interactive or not?
		if ((fseek(stdin,0,SEEK_END) == 0) && (ftell(stdin) > 0)) {
		  rewind(stdin);
		  echoCommands = YES;
		}
    
    //  Print out a header which shows version info, etc:
    eprintf("%@bold;TubeGen%@reset; - Carbon Nanotube Stucture Generator\n");
    printf("Version 3.1 [Build %d, %s %s]\n",__BUILD__,__TIME__,__DATE__);
    printf("Copyright (c) 2000-2002, J. T. Frey\n\n");
    
    //  If the user specified a "script" file, open it; otherwise we stick
    //  to using stdin:
    if (argc > 1) {
      if ((comFile = fopen(argv[1],"r")) == NULL)
        comFile = stdin;
    }
    
    //  This loop runs as long as the quit/exit command is not given.  When
    //  that happens, a simple "break" statement pops us out!
    while(1) {
      
      //  Read the next command, drop the terminating newline
      //  character, drop any leading whitespace:
      if (printPrompt) {
        printf("> ");
        fflush(stdout);
      }
			
      fgets(command,CMD_BUFFER_LEN,comFile);
			if (feof(comFile)) break;
			
      if (command[strlen(command) - 1] == '\n')
        command[strlen(command) - 1] = '\0';
      mungeWhitespace(command);
        
      //  Watch for comments; the '#' sign at the head of a line
      //  indicates this!
      if (*command == '#') {
        printPrompt = NO;
        continue;
      }
      printPrompt = YES;
        
      //  Print the command if we're reading a script:
      if ((comFile != stdin) || (echoCommands))
        printf("%s\n",command);
      
      //  Scan the command:
      if ((*command == '?') || (strcmp(command,"help") == 0)) {
        //  Display help page
        printf("  Command Summary:\n");
        printf("    help                print this information\n");
        printf("    status              display state variables\n");
        printf("    set VARIABLE VALUE  reset a state variable\n");
        printf("    generate            create the cell and basis\n");
        printf("    ftranslate          apply a fractional translation to the atomic basis\n");
        printf("    ctranslate          apply a Cartesian translation to the atomic basis\n");
        printf("    print               print the cell to the display\n");
        printf("    save                save the cell to a file\n");
        printf("    exit                exit the program\n");
        printf("  Variables:\n");
        printf("    format              {wien,gaussian}\n");
        printf("    units               {bohr,angstrom}\n");
        printf("    cc_bond             #\n");
        printf("    gutter              #,#,#\n");
        printf("    shape               {hexagonal,cubic,planar}\n");
        printf("    chirality           #,#\n");
        printf("    cell_count          #,#,#\n");
        printf("    relax_tube          {yes,no}\n");
      } else if ((strcmp(command,"exit") == 0) || (strcmp(command,"quit") == 0)) {
        //  Exit the program
        break;
      } else if (strcmp(command,"status") == 0) {
        //  Display state variables
        printf("  format            => ");
				switch(params.options.oFormat) {
				  case EFormat_WIEN:
            printf("wien\n");
						break;
					case EFormat_Gaussian:
					  printf("gaussian\n");
						break;
					default:
					  printf("\n");
				}
				
        printf("  units             => ");
				switch(params.options.oUnits) {
				  case EUnits_Bohr:
					  printf("bohr\n");
						break;
					case EUnits_Angstrom:
					  printf("angstrom\n");
						break;
					default:
					  printf("\n");
				}
				
        printf("  cc_bond           => ");
				switch(params.options.oUnits) {
				  case EUnits_Bohr:
					  printf("%lg (bohr)\n",globals.aCC / ANGSTROM_PER_BOHR);
						break;
					case EUnits_Angstrom:
					  printf("%lg (angstrom)\n",globals.aCC);
						break;
					default:
					  printf("\n");
				}
				
        printf("  gutter            => ");
				v = globals.gutter;
				switch(params.options.oUnits) {
				  case EUnits_Bohr:
					  printf("%lg,%lg,%lg (bohr)\n",v.x / ANGSTROM_PER_BOHR,v.y / ANGSTROM_PER_BOHR,v.z / ANGSTROM_PER_BOHR);
						break;
					case EUnits_Angstrom:
					  printf("%lg,%lg,%lg (angstrom)\n",v.x,v.y,v.z);
						break;
					default:
					  printf("\n");
				}
        
				printf("  shape             => ");
				switch(params.options.oLattice) {
				  case ELattice_Hexagonal:
						printf("hexagonal (rolled)\n");
						break;
					case ELattice_Cubic:
					  printf("cubic (rolled)\n");
						break;
					case ELattice_Planar:
						printf("planar\n");
						break;
					default:
					  printf("\n");
				}
				
        printf("  chirality         => %d,%d\n",n,m);
				
        printf("  cell_count        => %d,%d,%d\n", \
						globals.replicate[0],globals.replicate[1],globals.replicate[2] \
				);
				
        printf("  relax_tube        => %s\n",((params.options.oRelaxTubule)?("yes"):("no")));
				
      } else if (strstr(command,"set") == command) {
        //  Reset a state variable value
        char*		varName = command + 3;
        char*		end;
        char		saved;
        
        while ((*varName != '\0') && (!isalpha(*varName))) varName++;
        end = varName + 1;
        while ((*end != '\0') && ((isalpha(*end)) || (*end == '_'))) end++;
        saved = *end;
        *end = '\0';
        
        if (strcmp(varName,"format") == 0) {
          //  Exported file type:
          *end = saved;
          varName = end;
          while ((*varName != '\0') && (!isalpha(*varName))) varName++;
          end = varName + 1;
          while ((*end != '\0') && ((isalnum(*end)) || (*end == '_'))) end++;
          if (strcmp(varName,"wien") == 0)
            params.options.oFormat = EFormat_WIEN;
          else if (strcmp(varName,"gaussian") == 0)
            params.options.oFormat = EFormat_Gaussian;
          else
            eprintf("  %@yellow;WARNING%@reset;:  Invalid parameter value; variable not reset.\n");
        }
				
				else if (strcmp(varName,"units") == 0) {
          //  Units:
          *end = saved;
          varName = end;
          while ((*varName != '\0') && (!isalpha(*varName))) varName++;
          end = varName + 1;
          while ((*end != '\0') && ((isalnum(*end)) || (*end == '_'))) end++;
          if (strcmp(varName,"bohr") == 0)
            params.options.oUnits = EUnits_Bohr;
          else if (strcmp(varName,"angstrom") == 0)
            params.options.oUnits = EUnits_Angstrom;
          else
            eprintf("  %@yellow;WARNING%@reset;:  Invalid parameter value; variable not reset.\n");
        }
				
				else if (strcmp(varName,"cc_bond") == 0) {
          //  Carbon-Carbon bond length:
          double	tmp;
          
          *end = saved;
          varName = end;
          while ((*varName != '\0') && (!isdigit(*varName))) varName++;
          if (sscanf(varName,"%lg",&tmp) == 1) {
						switch(params.options.oUnits) {
						  case EUnits_Bohr:
							  tmp *= ANGSTROM_PER_BOHR;
								break;
						}
            setupGlobals(&globals,tmp,-1.0,-1.0,-1.0,0,0,0);
          } else
            eprintf("  %@yellow;WARNING%@reset;:  Invalid parameter value; variable not reset.\n");
        }
				
				else if (strcmp(varName,"gutter") == 0) {
          //  Unit cell padding (gutter):
          double	x,y,z;
          int		c;
          
          *end = saved;
          varName = end;
          while ((*varName != '\0') && (!isdigit(*varName))) varName++;
          c = sscanf(varName,"%lg,%lg,%lg",&x,&y,&z);
          if (c >= 1) {
						switch(params.options.oUnits) {
						  case EUnits_Bohr:
							  x *= ANGSTROM_PER_BOHR;
								break;
						}
						setupGlobals(&globals,0.0,x,-1.0,-1.0,0,0,0);
          }
          if (c >= 2) {
						switch(params.options.oUnits) {
						  case EUnits_Bohr:
							  y *= ANGSTROM_PER_BOHR;
								break;
						}
						setupGlobals(&globals,0.0,-1.0,y,-1.0,0,0,0);
          }
          if (c >= 3) {
						switch(params.options.oUnits) {
						  case EUnits_Bohr:
							  z *= ANGSTROM_PER_BOHR;
								break;
						}
						setupGlobals(&globals,0.0,-1.0,-1.0,z,0,0,0);
          }
          if (c < 1)
            eprintf("  %@yellow;WARNING%@reset;:  Invalid parameter value; variable not reset.\n");
        }
				
				else if (strcmp(varName,"shape") == 0) {
          //  Lattice shape:
          *end = saved;
          varName = end;
          while ((*varName != '\0') && (!isalpha(*varName))) varName++;
          end = varName + 1;
          while ((*end != '\0') && ((isalnum(*end)) || (*end == '_'))) end++;
          if (strcmp(varName,"hexagonal") == 0)
            params.options.oLattice = ELattice_Hexagonal;
          else if (strcmp(varName,"cubic") == 0)
            params.options.oLattice = ELattice_Cubic;
          else if (strcmp(varName,"planar") == 0)
            params.options.oLattice = ELattice_Planar;
          else
            eprintf("  %@yellow;WARNING%@reset;:  Invalid parameter value; variable not reset.\n");
        }
				
				else if (strcmp(varName,"chirality") == 0) {
          //  Chirality parameters (n,m):
          int		tn,tm;
          int		c;
          
          *end = saved;
          varName = end;
          while ((*varName != '\0') && (!isdigit(*varName))) varName++;
          c = sscanf(varName,"%u,%u",&tn,&tm);
          if (c >= 1) {
            n = tn;
          }
          if (c >= 2) {
            m = tm;
          }
          if (c < 1)
            eprintf("  %@yellow;WARNING%@reset;:  Invalid parameter value; variable not reset.\n");
        }
				
				else if (strcmp(varName,"cell_count") == 0) {
          //  How many cells to generate in each direction:
          unsigned	x,y,z;
          int		c;
          
          *end = saved;
          varName = end;
          while ((*varName != '\0') && (!isdigit(*varName))) varName++;
					x = y = z = 0;
          c = sscanf(varName,"%u,%u,%u",&x,&y,&z);
          if (c < 1)
            eprintf("  %@yellow;WARNING%@reset;:  Invalid parameter value; variable not reset.\n");
					else
					  setupGlobals(&globals,0.0,0.0,0.0,0.0,x,y,z);
        }
				
				else if (strcmp(varName,"relax_tube") == 0) {
          //  Use expansion ratios in rolled cell:
          *end = saved;
          varName = end;
          while ((*varName != '\0') && (!isalpha(*varName))) varName++;
          end = varName + 1;
          while ((*end != '\0') && ((isalnum(*end)) || (*end == '_'))) end++;
          if (strcmp(varName,"yes") == 0)
            params.options.oRelaxTubule = 1;
          else if (strcmp(varName,"no") == 0)
            params.options.oRelaxTubule = 0;
          else
            eprintf("  %@yellow;WARNING%@reset;:  Please specify yes or no.\n");
        } else
          eprintf("  %@yellow;WARNING%@reset;:  Unknown variable.\n");
      }
			
			else if (strcmp(command,"generate") == 0) {
        //  Generate a unit cell; if one already exists, we overwrite it:
        if (theCell)
          delete theCell;
        theCell = generateUnitCell(&params,n,m);
      }
			
			else if (strcmp(command,"print") == 0) {
        if (theCell)
          theCell->print(&cout);
        else
          eprintf("  %@red;ERROR%@reset;:  No cell generated yet.\n");
      }
			
			else if (strstr(command,"save") == command) {
        //  Save a unit cell
        ofstream*		file = NULL;
        char*				fName = command + 4;
        char*				end;
        char				saved;
        
				if (theCell == NULL) {
          eprintf("  %@red;ERROR%@reset;:  No cell generated yet.\n");
				} else {
					while ((*fName != '\0') && (isblank(*fName))) fName++;
					if (*fName == '\0') {
						char		filename[32];
						int			charcnt = 0;
						
						switch (params.options.oLattice) {
							case ELattice_Hexagonal:
							case ELattice_Cubic:
								charcnt = sprintf(filename,"%02d%02dr.",n,m);
								break;
							case ELattice_Planar:
								charcnt = sprintf(filename,"%02d%02df.",n,m);
								break;
							default:
								charcnt = sprintf(filename,"tubegen_out.");
								break;
						}
						switch (params.options.oFormat) {
						  case EFormat_Gaussian:
							  sprintf(filename + charcnt,"com\0");
								break;
							case EFormat_WIEN:
							  sprintf(filename + charcnt,"struct\0");
								break;
						}
						file = new ofstream(filename);
					} else {
						end = fName + 1;
						while ((*end != '\0') && (!isblank(*end))) end++;
						*end = '\0';
						file = new ofstream(fName);
					}
					
					if (file) {
						switch (params.options.oFormat) {
						
							case EFormat_WIEN: {
								file->form("(%d,%d) Carbon Nanotube\n",n,m);
								switch (params.options.oLattice) {
								  case ELattice_Hexagonal: {
										file->form("H   LATTICE,NONEQUIV. ATOMS: 1\n");
										file->form("MODE OF CALC=RELA\n");
										file->form("%10.6lf%10.6lf%10.6lf%10.6lf%10.6lf%10.6lf\n",theCell->GetDimensionA() / ANGSTROM_PER_BOHR,
												theCell->GetDimensionB() / ANGSTROM_PER_BOHR,theCell->GetDimensionC() / ANGSTROM_PER_BOHR,90.0,90.0,120.0);
										break;
									}
									
									case ELattice_Cubic: {
										file->form("P   LATTICE,NONEQUIV. ATOMS: 1\n");
										file->form("MODE OF CALC=RELA\n");
										file->form("%10.6lf%10.6lf%10.6lf%10.6lf%10.6lf%10.6lf\n",theCell->GetDimensionA() / ANGSTROM_PER_BOHR,
												theCell->GetDimensionB() / ANGSTROM_PER_BOHR,theCell->GetDimensionC() / ANGSTROM_PER_BOHR,90.0,90.0,90.0);
										break;
									}
								}
								for ( unsigned i = 0 ; i < theCell->GetBasisCount() ; i++ ) {
									TPoint3D		pt = theCell->GetCoordinate(i);
									if (i == 0)
										*file << "ATOM= ";
									else
										*file << "      ";
									file->form("-1: X=%10.8lf Y=%10.8lf Z=%10.8lf\n",pt.x,pt.y,pt.z);
									if (i == 0)
										file->form("          MULT=%2d          ISPLIT= 8\n",theCell->GetBasisCount());
								}
								*file << "C          NPT=  781  R0=0.00010000 RMT=    1.3000   Z:  6.0\n";
								*file << "                     1.0000000 0.0000000 0.0000000\n";
								*file << "                     0.0000000 1.0000000 0.0000000\n";
								*file << "                     0.0000000 0.0000000 1.0000000\n";
								*file << "   0 SYMMETRY OPERATIONS:\n";
								break;
							}
							
							case EFormat_Gaussian: {
								file->form("# uhf/3-21g\n\n (%d,%d) Carbon Nanotube",n,m);
								if ((globals.replicate[0] > 1) || (globals.replicate[1] > 1) || (globals.replicate[2] > 1))
									file->form(" (%d,%d,%d) replication",globals.replicate[0],globals.replicate[1],globals.replicate[2]);
								*file << "\n\n 0    1\n";
								theCell->Propogate(globals.replicate[0],globals.replicate[1],globals.replicate[2],file,kCrystalCellPropogateCentered);
								*file << "\n";
								break;
							}
						}
						delete file;
					}
				}
      }
			
			else if (strstr(command,"ftranslate") == command) {
        //  Apply a fractional translation:
        TVector3D		delta;
        int					c;
        char*				varName = command + 9;
        
        while ((*varName != '\0') && (!isdigit(*varName)) && (*varName != '-') && (*varName != '.') && (*varName != '+'))
          varName++;
        if (theCell) {
          c = sscanf(varName,"%lg,%lg,%lg",&delta.x,&delta.y,&delta.z);
          if (c == 3) {
            if ((delta.x > -1.0) && (delta.x < 1.0) && \
								(delta.y > -1.0) && (delta.y < 1.0) && \
								(delta.z > -1.0) && (delta.z < 1.0))
							theCell->ApplyFractionalTranslation(delta);
						else
							eprintf("  %@yellow;WARNING%@reset;:  Invalid translation vector; all values must be in (-1,1) in fractional coordinates.\n");
          } else
					  eprintf("  %@yellow;WARNING%@reset;:  Invalid translation vector; transformation not done.\n");
        } else
          eprintf("  %@red;ERROR%@reset;:  No cell generated yet.\n");
      }
			
			else if (strstr(command,"ctranslate") == command) {
        //  Apply a cartesian translation:
        TVector3D		delta;
        int					c;
        char*				varName = command + 9;
        
        while ((*varName != '\0') && (!isdigit(*varName)) && (*varName != '-') && (*varName != '.') && (*varName != '+'))
          varName++;
        if (theCell) {
          c = sscanf(varName,"%lg,%lg,%lg",&delta.x,&delta.y,&delta.z);
          if (c == 3) {
						TVector3D		ctv = theCell->GetCellTranslationVector();
						
						//  If projection of the translation onto the cell translation vector
						//  is -1.0 < p < 1.0 then it's a vector within the proper displacement
						//  bounds:
            if ( fabs(Vector3D_Dot(&ctv,&delta) / Vector3D_Dot(&ctv,&ctv)) < 1.0)
							theCell->ApplyCartesianTranslation(delta);
						else
							eprintf("  %@yellow;WARNING%@reset;:  Invalid Cartesian translation vector.\n");
          } else
					  eprintf("  %@yellow;WARNING%@reset;:  Invalid Cartesian translation vector.\n");
        } else
          eprintf("  %@red;ERROR%@reset;:  No cell generated yet.\n");
      }
			
			else if (*command != '\0')
        //  Command not understood
        eprintf("  %@red;ERROR%@reset;:  %s\n",command);
    }
    
    return 0;
}