// ===============================================================================
// 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

//
	
	enum {
		TGWriteAsWIENInput,
		TGWriteAsGaussianInput
	};
	
	enum {
		TGBohr,
		TGAngstrom
	};
	
	enum {
		TGHexagonalLattice,
		TGCubicLattice,
		TGFlatLattice
	};
		
//

	//  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");

//

  int							gMValue							= 3;
  int							gNValue							= 3;
  int							gLatticeShape				= TGHexagonalLattice;
  int							gFileFormat					= TGWriteAsWIENInput;
  int							gUnits							= TGAngstrom;
  double					gCC									= kStandardCarbonCarbon;
  
  //  The program will generate the lattice such that the tubule translation
  //  vector is ALWAYS aligned with the z-axis.  Hence, in both cases we typically
  //  don't want to pad in the z direction:
  TVector3D				gGutter							= { kStandardChiralPadding , kStandardChiralPadding , 0 };
  int							gCellCounts[3]			= { 1 , 1 , 1 };
	
	//  We can resize the unit cell in the case of the rolled lattice in order to
	//  keep the C-C bond length correct; this flag is ON by default:
	int							gUseExpansionRatios	= 1;

  TVector3D				a1,a2,Ch,T;
  double					LenCh,LenT;
  int							nPrime,mPrime;

//

  // ----------------------------------------------------------------------
  // * 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,int);
  
  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;
  }

//

  // ----------------------------------------------------------------------
  // * lcm
  // ----------------------------------------------------------------------
  // Updated:	Jeff Frey, 03.19.2001
  // Purpose:	Least common multiple of two integers.
  //
  // Last Mod:	n/a
  
  int lcm(int,int);
  
  int lcm(
    int i,
    int j
  )
  {
    int m,n,r;
    
    (i < j) ? (m = j,n = r = i) : (m = i,n = r = j);
    
    while (n % m)
      n += r;
    
    return n;
  }

//

  // ----------------------------------------------------------------------
  // * setupVectors
  // ----------------------------------------------------------------------
  // Updated:	Jeff Frey, 07.26.2001
  // Purpose:	Based on chirality indices (n,m) calculate the chiral
  //					vector, Ch, and the tubule translate vector, T.  Prints
  //					a lot of the information generated regarding the tubule
  //					as well.
  //
  // Last Mod:	18.NOV.2002: No longer returns alpha,beta values for loop
	//												 loop now uses nPrime,mPrime.
  
  int setupVectors();

  int setupVectors(
  )
  {
    int				TwoMplusN,TwoNplusM;
    double		conv,tmp;
    int				d,dR,N;
    
    //  Set up the basic symmetry vectors, a1 and a2:
    a1.x = sqrt(3.) * gCC * cos(kPi / 6.);
    a1.y = sqrt(3.) * gCC * sin(kPi / 6.);
    a2.x = a1.x;
    a2.y = -a1.y;
    a1.z = a2.z = 0.;
    Vector3D_Rezero(&a1,FLT_EPSILON);
    Vector3D_Rezero(&a2,FLT_EPSILON);
    
    //  Set up Ch:
    Ch.x = (double)gNValue * a1.x + (double)gMValue * a2.x;
    Ch.y = (double)gNValue * a1.y + (double)gMValue * a2.y;
    Ch.z = 0.;
    Vector3D_Rezero(&Ch,FLT_EPSILON);
    LenCh = Vector3D_Magnitude(&Ch);
    
    //  Calculate dR
    d = gcd(gNValue,gMValue);
		if (((gNValue - gMValue) % (3 * d)) == 0)
      dR = 3 * d;
    else
      dR = d;
    
    //  Calculate the tubule translation vector, T:
    TwoMplusN = 2 * gMValue + gNValue;
    TwoNplusM = 2 * gNValue + gMValue;
    nPrime = TwoMplusN / dR;
    mPrime = TwoNplusM / dR;
    T.x = ((double)nPrime * a1.x - (double)mPrime * a2.x);
    T.y = ((double)nPrime * a1.y - (double)mPrime * a2.y);
    T.z = 0.;
		
    Vector3D_Rezero(&T,FLT_EPSILON);
    LenT = Vector3D_Magnitude(&T);
    
    if (gUnits == TGBohr) {
      conv = 1.0 / ANGSTROM_PER_BOHR;
      printf("  All units are in bohr.\n");
    } else {
      conv = 1.0;
      printf("  All units are in angstrom.\n");
    }
    
    //  Calculate the number of hexagonal sub-cells:
    N = 2 * (gMValue * gMValue + gNValue * gNValue + gNValue * gMValue) / 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",gCC * conv);
    printf("  Basis vectors constructed using n = %d and m = %d:\n",gNValue,gMValue);
    printf("    a1 = < %lg , %lg >  a2 = < %lg , %lg >\n",a1.x * conv,a1.y * conv,a2.x * conv,a2.y * conv);
    printf("  Chiral vector Ch constructed as %d(a1) + %d(a2):\n",gNValue,gMValue);
    printf("    Ch = < %lg , %lg >, |Ch| = %lg\n",Ch.x * conv,Ch.y * conv,LenCh * conv);
    printf("  Tubule diameter is: %f\n",LenCh / kPi * conv);
		printf("  Tubule translation vector T constructed as %d(a1) - %d(a2):\n",nPrime,mPrime);
    printf("    T = < %lg , %lg >, |T| = %lg\n",T.x * conv,T.y * conv,LenT * conv);
    printf("  Check for orthogonality:  T.Ch = %lg\n",(tmp = Vector3D_Dot(&T,&Ch) * conv));
    
    return 0;
  }
  
//

	typedef struct SGammaCalcEntity {
		double		l;
		double		h;
		double		theta;
	} SGammaCalcEntity;
	
	typedef struct SGammaCalculation {
	  double						gammaCh,gammaT;
		double						r;
		double						aCC;
		SGammaCalcEntity	params[4];
	} SGammaCalculation;
	
	double
	calculateGammaCh(
		SGammaCalculation*	g,
		int									i,
		int									j
	)
	{
		double					denominator,numerator,t1,t2,rsqrd;
		
		rsqrd = g->r * g->r;
		numerator = g->aCC * g->aCC;
		
		//  Calculate and subtract gammaT term:
		t2 = g->gammaT * g->gammaT;
		t1 = g->params[i].h - g->params[j].h;
		numerator -= t2 * t1 * t1;
		
		//  Calculate denominator:
		t2 = 2.0 * rsqrd;
		t1 = 1.0 - cos( g->params[i].theta - g->params[j].theta );
		denominator = t2 * t1;
		
		return sqrt( numerator / denominator );
	}
	
	double
	calculateGammaT(
		SGammaCalculation*	g,
		int									i,
		int									j
	)
	{
		double					denominator,numerator,t1,t2,rsqrd;
		
		rsqrd = g->r * g->r;
		numerator = g->aCC * g->aCC;
		
		//  Calculate angular part of numerator:
		t2 = 2.0 * rsqrd * g->gammaCh * g->gammaCh;
		t1 = 1.0 - cos( g->params[i].theta - g->params[j].theta );
		numerator -= t2 * t1;
		
		//  Calculate denominator:
		denominator = g->params[i].h - g->params[j].h;
		denominator *= denominator;
		
		return sqrt( numerator / denominator );
	}
	
	double
	calculateDistance(
		SGammaCalculation*	g,
		int									i,
		int									j
	)
	{
		double					accum,t1,t2,t3,tmp,rsqrd;
		
		rsqrd = g->r * g->r;
		accum = 0.0;
		
		//  X-term:
		t3 = 2.0;
		t1 = cos(g->params[i].theta); t2 = t1 * t1; t3 *= t1;
		t1 = cos(g->params[j].theta); t2 += t1 * t1; t3 *= t1;
		t2 -= t3;
		tmp = g->gammaCh * g->gammaCh * rsqrd * t2;
		accum += tmp;
		
		//  Y-term:
		t3 = 2.0;
		t1 = sin(g->params[i].theta); t2 = t1 * t1; t3 *= t1;
		t1 = sin(g->params[j].theta); t2 += t1 * t1; t3 *= t1;
		t2 -= t3;
		tmp = g->gammaCh * g->gammaCh * rsqrd * t2;
		accum += tmp;
		
		//  Z-term:
		t1 = g->gammaT * g->gammaT;
		t2 = g->params[i].h - g->params[j].h;
		accum += t1 * t2 * t2;
		
		return sqrt( accum );
		
	}
	
	void
	calculateGammas(
		double			r,
		double			*gammaCh,
		double			*gammaT
	)
	{
		SGammaCalculation		gCalc;
		double							tmp;
		int									iters,iterdone = 0;
		
		//  Initial values:
		gCalc.gammaCh = gCalc.gammaT = 1.0;
		gCalc.r = r;
		gCalc.aCC = gCC;
		
		gCalc.params[0].l				= 0;
		gCalc.params[0].h				= 0;
		gCalc.params[0].theta		= 0;
		
		gCalc.params[1].l				= tmp = gCC * Ch.x / LenCh; // <-- < aCC * X , Ch > / ||Ch||
		gCalc.params[1].h				= gCC * T.x / LenT;   // <-- < aCC * X , T > / ||T||
		gCalc.params[1].theta		= tmp / r;
		
		gCalc.params[2].l				= tmp = Vector3D_Dot(&a1,&Ch) / LenCh;
		gCalc.params[2].h				= Vector3D_Dot(&a1,&T) / LenT;
		gCalc.params[2].theta		= tmp / r;
		
		gCalc.params[3].l				= tmp = Vector3D_Dot(&a2,&Ch) / LenCh;
		gCalc.params[3].h				= Vector3D_Dot(&a2,&T) / LenT;
		gCalc.params[3].theta		= tmp / r;
		
		printf("  %-2s | %-12s | %-12s | %-12s\n","i","   l","   h","   theta");
		printf("  ------------------------------------------------------\n");
		for ( iters = 0 ; iters < 4 ; iters++ )
		  printf("  %-2d | %12.8f | %12.8f | %12.8f\n",iters,gCalc.params[iters].l,
				gCalc.params[iters].h,gCalc.params[iters].theta);
		printf("  ------------------------------------------------------\n");
		
		iters = 0;
		gCalc.gammaCh = calculateGammaCh(&gCalc,0,1);
		while (!iterdone) {
			tmp = calculateGammaT(&gCalc,1,3);
			if (fabs(tmp - gCalc.gammaT) < DBL_EPSILON)
				iterdone = 1;
			gCalc.gammaT = tmp;
			iters++;
			tmp = calculateGammaCh(&gCalc,0,1);
			if ((!iterdone) && (fabs(tmp - gCalc.gammaCh) < DBL_EPSILON))
				iterdone = 1;
			gCalc.gammaCh = tmp;
		}
		
		printf("  Expansion/Compression ratios converged to %1.0g",DBL_EPSILON);
		if (iters > 0) {
			printf(" in %d cycle",iters);
			if (iters > 1)
				printf("s\n");
			else
				printf("\n");
		} else {
			printf(" immediately.\n");
		}
		printf("  gammaCh = %lf [cylindrical expansion ratio]\n",gCalc.gammaCh);
		printf("  gammaT  = %lf [axial expansion ratio]\n",gCalc.gammaT);
	
		if (gammaCh)
			*gammaCh = gCalc.gammaCh;
		if (gammaT)
			*gammaT = gCalc.gammaT;
		
		printf("  [0,1]   = %lf\n",calculateDistance(&gCalc,0,1));
		printf("  [1,2]   = %lf\n",calculateDistance(&gCalc,1,2));
		printf("  [1,3]   = %lf\n",calculateDistance(&gCalc,1,3));
	
	}

//

  // ----------------------------------------------------------------------
  // * generateUnitCell
  // ----------------------------------------------------------------------
  // Updated:	Jeff Frey, 07.26.2001
  // Purpose:	Creates a CrystalCell object using the dimensions generated
  //					by the setupVectors function (coupled with the user's
  //					gutter dimensions for the cell).  Fills the cell with
  //					unique carbon atom positions either as a planar
  //					lattice or as a rolled, tubule structure.  The routine
  //					returns either NULL or the resulting CrystalCell object.
  //
  // Last Mod:	n/a

  CrystalCell* generateUnitCell();
  
  CrystalCell*
  generateUnitCell()
  {
    CrystalCell*			cell = NULL;
    double						a,b,c;
    double						r;
    int								i,j,k;
		int								iMax,jMax,iMin,jMin;
		char							didCalcGamma = 0;
    TVector3D					v;
		TVector3D					center;
    TPoint3D					p,p1;
    double						LenV;
		double						gammaCh = 1.0,gammaT = 1.0;
    
		switch (gLatticeShape) {
      case TGHexagonalLattice: {
        printf("  Producing rolled, hexagonal nanotube lattice.\n\n");
				break;
			}
			case TGCubicLattice: {
        printf("  Producing rolled, cubic nanotube lattice.\n\n");
				break;
			}
			default: {
				printf("  Producing flat nanotube lattice.\n\n");
				break;
			}
		}
      
    //  Setup at a1 and a2, the two basis vectors,
    //  Ch (chiral vector), and T (tubule translation vector):
		printf("  -- Generating Graphitic Basis -----------------------\n");
    setupVectors();
		printf("  -----------------------------------------------------\n\n");
    
    //  Calculate the radius of the tube (used by the rolling routine):
    r = 0.5 * LenCh * M_1_PI;
		
		//  Create the crystal cell:
    if ((gLatticeShape == TGHexagonalLattice) || (gLatticeShape == TGCubicLattice)) {
			a = 2.0 * (r + gGutter.x);
			b = 2.0 * (r + gGutter.y);
			c = LenT + 2.0 * gGutter.z;
			if (gLatticeShape == TGHexagonalLattice) {
			  //b *= 2.0 / sqrt(3.0);
		    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);
		} else {
		  a = LenCh + 2.0 * gGutter.x;
			b = LenCh + 2.0 * gGutter.y;
			c = LenT + 2.0 * gGutter.z;
		  cell = new CrystalCell(a,b,c,90.0,90.0,90.0);
		}
		
		if (gUseExpansionRatios) {
			double		tmp;
			//  If requested, figure out the radial and axial expansion/compression
			//  scalars:
			printf("  -- Radial/Axial Expansion/Compression ---------------\n");
			calculateGammas(r,&gammaCh,&gammaT);
			
			//  Now that we have the gammaT, we need to reset the c-dimension
			//  of the cell:
			cell->SetDimensionC( cell->GetDimensionC() * gammaT );
			//  Now that we have the gammaCh, we need to reset the a,b-dimensions
			//  of the cell:
			tmp = 2.0 * (gammaCh - 1.0) * r;
			cell->SetDimensionA( cell->GetDimensionA() + tmp );
			cell->SetDimensionB( cell->GetDimensionB() + tmp );
			printf("  Cell redimensioned to match expansions/compressions.\n");
			printf("  -----------------------------------------------------\n\n");
		}
    
    //  Begin generating coordinates:
		iMin = (nPrime < 0)?(nPrime):(0);
		iMax = ((gNValue + nPrime) > gNValue)?(gNValue + nPrime):(gNValue);
		jMin = (-mPrime > 0)?(-mPrime):(0);
		jMax = ((gMValue + mPrime) > gMValue)?(gMValue + mPrime):(gMValue);
		printf("[%d %d | %d %d]\n",iMin,iMax,jMin,jMax);
		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 * a1.x + (double)j * a2.x;
          v.y = (double)i * a1.y + (double)j * a2.y;
          v.z = 0.0;
          
          //  Second time through we add a C-C displacement along the x-axis:
          if (k == 1)
            v.x += gCC;
          
          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,&Ch) / (LenCh * LenCh);
            p.z = Vector3D_Dot(&v,&T) / (LenT * 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 ((gLatticeShape == TGHexagonalLattice) || (gLatticeShape == TGCubicLattice)) {
                //  theta = 2(pi) times displacement along chiral vector:
                double		theta = k2Pi * p.x;
								
								//  Redefine the point as a polar coordinate in xy-plane:
								//  Are we doing expansion ratio fixing of rolled C-C bond lengths?
								if (gUseExpansionRatios) {
									
									p.x = gammaCh * r * cos(theta) + center.x;
									p.y = gammaCh * r * sin(theta) + center.y;
									p.z *= LenT * gammaT;
								} else {
									p.x = r * cos(theta) + center.x;
									p.y = r * sin(theta) + center.y;
									p.z *= LenT;
								}
								cell->DidAddAtomAtCartesianPoint(6,p);
              } else {
                p.x = (p.x * LenCh + gGutter.x) / a;
                p.z = (p.z * LenT + gGutter.z) / c;
								cell->DidAddAtomAtFractionalPoint(6,p);
              }
            }
          }
        }
      }
    }
    if (cell)
      printf("  Cell generation complete.  %d basis points defined.\n",cell->GetBasisCount());
    
    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;
		
		//  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.0.6 [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,flat}\n");
        printf("    chirality           #,#\n");
        printf("    cell_count          #,#,#\n");
        printf("    fix_roll_lengths    {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            => ");
        if (gFileFormat == TGWriteAsWIENInput)
          printf("wien\n");
        else
          printf("gaussian\n");
        printf("  units             => ");
        if (gUnits == TGBohr)
          printf("bohr\n");
        else
          printf("angstrom\n");
        printf("  cc_bond           => ");
        if (gUnits == TGBohr)
          printf("%lg (bohr)\n",gCC / ANGSTROM_PER_BOHR);
        else
          printf("%lg (angstrom)\n",gCC);
        printf("  gutter            => ");
        if (gUnits == TGBohr)
          printf("%lg,%lg,%lg (bohr)\n",gGutter.x / ANGSTROM_PER_BOHR,gGutter.y / ANGSTROM_PER_BOHR,gGutter.z / ANGSTROM_PER_BOHR);
        else
          printf("%lg,%lg,%lg (angstrom)\n",gGutter.x,gGutter.y,gGutter.z);
        printf("  shape             => ");
        if (gLatticeShape == TGHexagonalLattice)
          printf("hexagonal (rolled)\n");
				else if (gLatticeShape == TGCubicLattice)
          printf("cubic (rolled)\n");
        else
          printf("flat\n");
        printf("  chirality         => %d,%d\n",gNValue,gMValue);
        printf("  cell_count        => %d,%d,%d\n",gCellCounts[0],gCellCounts[1],gCellCounts[2]);
        printf("  fix_roll_lengths  => %s\n",((gUseExpansionRatios)?("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)
            gFileFormat = TGWriteAsWIENInput;
          else if (strcmp(varName,"gaussian") == 0)
            gFileFormat = TGWriteAsGaussianInput;
          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)
            gUnits = TGBohr;
          else if (strcmp(varName,"angstrom") == 0)
            gUnits = TGAngstrom;
          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)
            if (gUnits == TGBohr)
              gCC = tmp * ANGSTROM_PER_BOHR;
            else
              gCC = tmp;
          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) {
            if (gUnits == TGBohr)
              gGutter.x = x * ANGSTROM_PER_BOHR;
            else
              gGutter.x = x;
          }
          if (c >= 2) {
            if (gUnits == TGBohr)
              gGutter.y = y * ANGSTROM_PER_BOHR;
            else
              gGutter.y = y;
          }
          if (c >= 3) {
            if (gUnits == TGBohr)
              gGutter.z = z * ANGSTROM_PER_BOHR;
            else
              gGutter.z = z;
          }
          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)
            gLatticeShape = TGHexagonalLattice;
          else if (strcmp(varName,"cubic") == 0)
            gLatticeShape = TGCubicLattice;
          else if (strcmp(varName,"flat") == 0)
            gLatticeShape = TGFlatLattice;
          else
            eprintf("  %@yellow;WARNING%@reset;:  Invalid parameter value; variable not reset.\n");
        } else if (strcmp(varName,"chirality") == 0) {
          //  Chirality parameters (n,m):
          int		n,m;
          int		c;
          
          *end = saved;
          varName = end;
          while ((*varName != '\0') && (!isdigit(*varName))) varName++;
          c = sscanf(varName,"%u,%u",&n,&m);
          if (c >= 1) {
            gNValue = n;
          }
          if (c >= 2) {
            gMValue = m;
          }
          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++;
          c = sscanf(varName,"%u,%u,%u",&x,&y,&z);
          if (c >= 1) {
            gCellCounts[0] = x;
          }
          if (c >= 2) {
            gCellCounts[1] = y;
          }
          if (c >= 3) {
            gCellCounts[2] = z;
          }
          if (c < 1)
            eprintf("  %@yellow;WARNING%@reset;:  Invalid parameter value; variable not reset.\n");
        } else if (strcmp(varName,"fix_roll_lengths") == 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)
            gUseExpansionRatios = 1;
          else if (strcmp(varName,"no") == 0)
            gUseExpansionRatios = 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();
      } 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];
						
						if (gFileFormat == TGWriteAsWIENInput) {
							if ((gLatticeShape == TGHexagonalLattice) || (gLatticeShape == TGCubicLattice))
								sprintf(filename,"%02d%02dr.struct\0",gNValue,gMValue);
							else
								sprintf(filename,"%02d%02df.struct\0",gNValue,gMValue);
						} else {
							if ((gLatticeShape == TGHexagonalLattice) || (gLatticeShape == TGCubicLattice))
								sprintf(filename,"%02d%02dr.com\0",gNValue,gMValue);
							else
								sprintf(filename,"%02d%02df.com\0",gNValue,gMValue);
						}
						file = new ofstream(filename);
					} else {
						end = fName + 1;
						while ((*end != '\0') && (!isblank(*end))) end++;
						*end = '\0';
						file = new ofstream(fName);
					}
					
					if (file) {
						if (gFileFormat == TGWriteAsWIENInput) {
							file->form("(%d,%d) Carbon Nanotube\n",gNValue,gMValue);
							if (gLatticeShape == TGHexagonalLattice) {
								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);
							} else {
								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);
							}
							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";
						} else {
							file->form("# uhf/3-21g\n\n (%d,%d) Carbon Nanotube",gNValue,gMValue);
							if ((gCellCounts[0] > 1) || (gCellCounts[1] > 1) || (gCellCounts[2] > 1))
								file->form(" (%d,%d,%d) replication",gCellCounts[0],gCellCounts[1],gCellCounts[2]);
							*file << "\n\n 0    1\n";
							theCell->Propogate(gCellCounts[0],gCellCounts[1],gCellCounts[2],file,kCrystalCellPropogateCentered);
							*file << "\n";
						}
						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;
}