#include <tk.h>
#if defined(__WIN32__) || defined(_WIN32)
#   define WIN32_LEAN_AND_MEAN
#   include <windows.h>
#   undef WIN32_LEAN_AND_MEAN
#endif
#include <GL/gl.h>
#include <math.h>
#include <memory.h>
#include <assert.h>
#include <string.h>
#include "gencyl.h"
#include "tkogl.h"

#ifndef PI
/* The ubiquitous PI constant */
#define PI 3.14159265358979323846264338327950288
#endif

/* Constants that define the various flags for rendering a generic
 * cylinder */

#define STITCH_ENDS 0x01 /* Stitch together the ends of the cyl */
#define STITCH_LOOPS 0x02 /* Close the cross-sections of the cyl */
#define SHADE_SMOOTH 0x08 /* Render with one normal per vertex */
#define TEXTURE_GEN 0x04 /* Generate texture coordinates */

/*---------------------------------------------------------------------------
 *
 *  General operations with vectors and transformation matrices
 *
 *---------------------------------------------------------------------------*/

typedef GLfloat Vector [3];
typedef GLfloat Matrix [4][4];

static void 
MatrixMult (const Matrix m1, const Matrix m2, Matrix result) 
{
   /* Multiplies m1 by m2 storing the result in result */
   int i, j, k;
   Matrix tmp;

   for (i = 0; i < 4; i++) {
      for (j = 0; j < 4; j++) { 
	 tmp [i][j] = 0.0;
	 for (k = 0; k < 4; k++) {
	    tmp [i][j] += m1 [i][k] * m2 [k][j];
	 }
      }
   }
   for (i = 0; i < 4; i++) {
      for (j = 0; j < 4; j++) { 
	 result [i][j] = tmp [i][j];
      }
   }
}

static void 
TransformVector (const Vector v, const Matrix m, Vector result)
{
   /* Applies affine linear transformation m to v storing the result in 
    * result */
   int i, j;
   Vector tmp;

   for (i = 0; i < 3; i++) {
      tmp [i] = m [i][3];
      for (j = 0; j < 3; j++) {
	 tmp [i] += m [i][j] * v [j];
      }
   }
   for (i = 0; i < 3; i++) {
      result [i] = tmp [i];
   }
}

static void 
NormalizeVector (Vector v)
{
   /* Makes v a unit vector */
   GLfloat hypot = sqrt (v [0] * v [0] +
			 v [1] * v [1] +
			 v [2] * v [2]);
   if (hypot == 0.0) {
      return;
   }
   v [0] /= hypot;
   v [1] /= hypot;
   v [2] /= hypot;
}

static void 
AddVector (const Vector v1, const Vector v2, Vector result)
{
   /* Adds v1 to v2 and stores in result */
   result [0] = v1 [0] + v2 [0];
   result [1] = v1 [1] + v2 [1];
   result [2] = v1 [2] + v2 [2];
}

#define InitVector(v,a,b,c) {v[0]=a; v[1]=b; v[2]=c;}

static void
ComputeTriangleNormal (const Vector v1, const Vector v2, const Vector v3,
		       Vector normal)
{
   /* Computes the normal of the triangle given by the three vertices v1, v2, v3
    * and stores it in the vector given by normal */
   Vector e1, e2; 
   GLfloat hypot;
   InitVector (e1, v2[0] - v1[0], v2[1] - v1[1], v2[2] - v1[2]);
   InitVector (e2, v3[0] - v2[0], v3[1] - v2[1], v3[2] - v2[2]);
   normal [0] = e1[1]*e2[2]-e1[2]*e2[1];
   normal [1] = e2[0]*e1[2]-e2[2]*e1[0];
   normal [2] = e1[0]*e2[1]-e1[1]*e2[0];
   NormalizeVector (normal);
}

static void
ComputeQuadNormal (const Vector v1, const Vector v2, const Vector v3, 
		   const Vector v4, Vector normal)
{
   /* Computes the normal at vertex v1 of the quadrilateral given
    * by v1-v2-v3-v4 and stores it in the vector given by normal */
   Vector e1, e2, e3, e4, n1, n2;
   GLfloat h1, h2;
   InitVector (e1, v2[0] - v1[0], v2[1] - v1[1], v2[2] - v1[2]);
   InitVector (e2, v3[0] - v2[0], v3[1] - v2[1], v3[2] - v2[2]);
   InitVector (e3, v4[0] - v3[0], v4[1] - v3[1], v4[2] - v3[2]);
   InitVector (e4, v1[0] - v4[0], v1[1] - v4[1], v1[2] - v4[2]);
   InitVector (n1, e1[1]*e2[2]-e1[2]*e2[1],
		        e2[0]*e1[2]-e2[2]*e1[0],
		        e1[0]*e2[1]-e1[1]*e2[0]);
   InitVector (n2, e3[1]*e4[2]-e3[2]*e4[1],
		        e4[0]*e3[2]-e4[2]*e3[0],
		        e3[0]*e4[1]-e3[1]*e4[0]);
   h1 = sqrt(n1[0]*n1[0]+n1[1]*n1[1]+n1[2]*n1[2]);
   h2 = sqrt(n2[0]*n2[0]+n2[1]*n2[1]+n2[2]*n2[2]);
   if (h1 * 100 > h2) {
      normal [0] = n1 [0] / h1;
      normal [1] = n1 [1] / h1;
      normal [2] = n1 [2] / h1;
   } 
   else {
      if (h2 == 0.0) return;
      normal [0] = n2 [0] / h2;
      normal [1] = n2 [1] / h2;
      normal [2] = n2 [2] / h2;
   }
}


/*---------------------------------------------------------------------------
 *
 *  Management of cross sections
 *
 *---------------------------------------------------------------------------*/

typedef struct {
   Vector * vtx;
   Vector * normal;
   int nvtx, vtxsize;
} CrossSection;

static CrossSection* 
NewCrossSection ()
{
   /* Allocates a new cross section structure */
   CrossSection * result = (CrossSection*) malloc (sizeof (CrossSection));
   assert (result != NULL);
   result->nvtx = 0;
   result->vtxsize = 8;
   result->vtx = (Vector*) malloc (sizeof (Vector)*result->vtxsize);
   result->normal = NULL;
   assert (result->vtx != NULL);
   return result;
}

static void 
FreeCrossSection (CrossSection * s)
{
   /* Deallocates the memory associated with cross-section s */
   assert (s != NULL);
   assert (s->vtx != NULL);
   if (s->normal != NULL) free(s->normal);
   free (s->vtx);
   free (s);
}

static void
AddCrossSectionVtx (CrossSection * s, GLfloat x, GLfloat y, GLfloat z) 
{
   /* Stores another vertex with coords (x,y,0) as the next vertex in
    * cross section 's' */
   assert (s != NULL);
   if (s->vtxsize == s->nvtx) {
      s->vtxsize *= 2;
      s->vtx = (Vector*) realloc (s->vtx, sizeof(Vector)*s->vtxsize);
      assert (s->vtx != NULL);
   }
   s->vtx [s->nvtx][0] = x;
   s->vtx [s->nvtx][1] = y;
   s->vtx [s->nvtx][2] = z;
   s->nvtx++;
}

static CrossSection* 
PolygonCrossSection (GLfloat radius, int nsides)
{
   /* Returns a cross section which is a regular polygon with 'nsides' sides
    * and radius 'radius' */
   GLfloat x, y, ang, incAng;
   int i;
   CrossSection *cross = NewCrossSection ();
   incAng = PI * 2 / nsides;
   ang = 0.0;
   for (i = 0; i < nsides; i++) {
      x = radius * sin (ang);
      y = radius * cos (ang);
      AddCrossSectionVtx (cross, x, y, 0.0);
      ang += incAng;
   }
   return cross;
}

static void 
TransformCrossSection (CrossSection* s, const Matrix m)
{
   int i;
   for (i = 0; i < s->nvtx; i++) {
      TransformVector (s->vtx [i], m, s->vtx [i]);
   }
}

static void 
DupCrossSection (const CrossSection* src, CrossSection* dst)
{
   if (dst->vtxsize < src->vtxsize) {
      dst->vtxsize = src->vtxsize;
      dst->vtx = (Vector*) realloc (dst->vtx, sizeof(Vector)*dst->vtxsize);
      assert (dst->vtx != NULL);
   }
   dst->nvtx = src->nvtx;
   memcpy (dst->vtx, src->vtx, sizeof (Vector)*dst->nvtx);
}

/*---------------------------------------------------------------------------
 *
 * Model management
 *
 *---------------------------------------------------------------------------*/

typedef struct {
   double smin, smax, tmin, tmax;
   int flags;
   int ncross, sizecross;
   CrossSection** cross;
} Model;

static Model * 
NewModel () 
{
   /* Allocates a new model and returns a pointer to it */
   Model* result;
   
   result = (Model*) malloc (sizeof (Model));
   result->ncross = 0;
   result->sizecross = 8;
   result->cross = (CrossSection**) malloc (sizeof (CrossSection*) * 
					    result->sizecross);
   assert (result->cross != NULL);
   return result;
}

static void 
FreeModel (Model* model)
{      
   /* Deallocates all memory associated with model model */
   int i;
   for (i = 0; i < model->ncross; i++) {
      FreeCrossSection (model->cross [i]);
   }
   free (model->cross);
}

static void
AddModelCrossSection (Model* model, CrossSection* cross)
{
   /* Adds another cross section to the model */
   if (model->sizecross == model->ncross) {
      model->sizecross *= 2;
      model->cross = (CrossSection**) realloc (model->cross,
				sizeof (CrossSection*) * model->sizecross);
   }
   model->cross [model->ncross++] = cross;
}

static void
ComputeModelNormals (Model* model)
{
   /* Computes normals for each vertex of the model */
   int icross, ivtx;
   int flags = model->flags;
   for (icross = 0; icross < model->ncross; icross++) {
      CrossSection *thisCross = model->cross [icross];
      CrossSection *nextCross = 
	 (icross+1<model->ncross || (flags&STITCH_ENDS) ? 
	  model->cross [(icross+1)%model->ncross] : NULL);
      CrossSection *prevCross = 
	 (icross-1 >= 0 || (flags&STITCH_ENDS) ?
	  model->cross [(icross+model->ncross-1)%model->ncross] : NULL);
      int n = thisCross->nvtx;
      thisCross->normal = (Vector*) malloc (sizeof (Vector) * n);
      for (ivtx = 0; ivtx < n; ivtx++) {
	 Vector* this = &(thisCross->vtx [ivtx]);
	 Vector* next = (((ivtx < n-1) || (flags&STITCH_LOOPS)) ?
			 &(thisCross->vtx [(ivtx+1)%n]) : NULL);
	 Vector* prev = (((ivtx != 0) || (flags&STITCH_LOOPS)) ?
			 &(thisCross->vtx [(ivtx+n-1)%n]) : NULL);
	 Vector* left = (nextCross == NULL ? NULL : 
			 &(nextCross->vtx [ivtx * nextCross->nvtx / n]));
	 Vector* right = (prevCross == NULL ? NULL :
			  &(prevCross->vtx [ivtx * prevCross->nvtx / n]));
	 Vector normal = { 0.0, 0.0, 0.0 };	 
	 if (left != NULL) {
	    Vector tmp;
	    Vector* dPrev =&(nextCross->vtx[((ivtx+n-1)%n)*nextCross->nvtx/n]);
	    Vector* dNext =&(nextCross->vtx[((ivtx+1)%n)*nextCross->nvtx/n]);
	    if (prev != NULL) {
	       ComputeQuadNormal (*this, *prev, *dPrev, *left, tmp);
	       AddVector (normal, tmp, normal);
	    }
	    if (next != NULL) {
	       ComputeQuadNormal (*this, *left, *dNext, *next, tmp);
	       AddVector (normal, tmp, normal);
	    }
	 }
	 if (right != NULL) {
	    Vector tmp;
	    Vector* dPrev =&(prevCross->vtx[((ivtx+n-1)%n)*prevCross->nvtx/n]);
	    Vector* dNext =&(prevCross->vtx[((ivtx+1)%n)*prevCross->nvtx/n]);
	    if (next != NULL) {
	       ComputeQuadNormal (*this, *next, *dNext, *right, tmp);
	       AddVector (normal, tmp, normal);
	    }
	    if (prev != NULL) {
	       ComputeQuadNormal (*this, *right, *dPrev, *prev, tmp);
	       AddVector (normal, tmp, normal);
	    }
	 }
	 NormalizeVector (normal);
	 memcpy (&(thisCross->normal [ivtx]), normal, sizeof (Vector));
      }
   }
}

static void
ComputeModelBBox (Tcl_Interp* interp, Model* model)
{
   /* Computes the model's bounding box and stores it as two elements
      in the Interp's result: one for the minimum corner and another for
      the maximum corner */
   Vector min = { 9e99, 9e99, 9e99 };
   Vector max = { -9e99, -9e99, -9e99 };
   char buf [200];
   int icross, ivtx, icoord;
   for (icross = 0; icross < model->ncross; icross++) {
      CrossSection *cross = model->cross [icross];
      for (ivtx = 0; ivtx < cross->nvtx; ivtx++) {
	 for (icoord = 0; icoord < 3; icoord++) {
	    if (cross->vtx [ivtx][icoord] < min[icoord]) {
	       min [icoord] = cross->vtx [ivtx][icoord];
	    }
	    if (cross->vtx [ivtx][icoord] > max[icoord]) {
	       max [icoord] = cross->vtx [ivtx][icoord];
	    }
	 }
      }
   }
   sprintf (buf, "min %f %f %f", min[0], min[1], min[2]);
   Tcl_AppendElement (interp, buf);
   sprintf (buf, "max %f %f %f", max[0], max[1], max[2]);
   Tcl_AppendElement (interp, buf);
}

static void 
GenTexCoords (Model* model, int icross, int ivtx) 
{
   /* Generates texture coordinates for vertex ivtx of cross section icross
      of model */
   GLfloat texcoords [2];
   if (!(model->flags & TEXTURE_GEN)) return;
   if (model->flags & STITCH_ENDS) {
      texcoords [0] = (model->smin + icross * (model->smax-model->smin) /
		       model->ncross);
   }
   else {
      texcoords [0] = (model->smin + icross * (model->smax-model->smin) /
		       (model->ncross-1));
   }
   if (model->flags & STITCH_LOOPS) {
      texcoords [1] = (model->tmin + ivtx * (model->tmax-model->tmin) /
		       (model->cross[icross%(model->ncross)]->nvtx));
   } 
   else {
      texcoords [1] = (model->tmin + ivtx * (model->tmax-model->tmin) /
		       (model->cross[icross%(model->ncross)]->nvtx-1));
   }
   glTexCoord2fv (texcoords);
}
   

static void 
RenderModelQuadStrip (Model* model) 
{
   /* Renders the surface between each pair of cross sections as a 
      quad strip */

   int icross, previcross, ivtx;
   int ncross = model->ncross;
   int flags = model->flags;
   GLfloat texcoords [2] = { 0.0, 0.0 };
   if (flags&STITCH_ENDS) {
      icross = 0;
      previcross = ncross;
   }
   else {
      icross = 1;
      previcross = 0;
   }
   while (icross < ncross) {
      CrossSection *a = model->cross [(icross+ncross-1)%ncross];
      CrossSection *b = model->cross [icross];     
      int na = a->nvtx;
      int nb = b->nvtx;
      int n = (nb > na ? nb : na);
      glBegin (GL_TRIANGLE_STRIP);
      texcoords [1] = icross * 1.0 / (ncross-1);
      if (flags & SHADE_SMOOTH) {
	 /* One normal per vertex */
	 for (ivtx = 0; ivtx < n; ivtx++) {
	    int ia = ivtx * na / n;
	    int ib = ivtx * nb / n;
	    GenTexCoords (model, previcross, ia);
	    glNormal3fv (a->normal [ia]);
	    glVertex3fv (a->vtx [ia]);
	    GenTexCoords (model, icross, ib);
	    glNormal3fv (b->normal [ib]);
	    glVertex3fv (b->vtx [ib]);
	 }
	 if (flags & STITCH_LOOPS) {
	    GenTexCoords (model, previcross, na);
	    glNormal3fv (a->normal [0]);
	    glVertex3fv (a->vtx [0]);
	    GenTexCoords (model, icross, nb);
	    glNormal3fv (b->normal [0]);
	    glVertex3fv (b->vtx [0]);
	 }
      }
      else {
	 /* One normal per face */
	 for (ivtx = 0; ivtx < n; ivtx++) {
	    int ia = ivtx * na / n;
	    int ib = ivtx * nb / n;
	    int ianext = (ivtx+1) * na / n % na;
	    int ibnext = (ivtx+1) * nb / n % nb;
	    Vector normal;
	    AddVector (a->normal [ia], b->normal [ib], normal);
	    AddVector (a->normal [ianext], normal, normal);
	    AddVector (b->normal [ibnext], normal, normal);
	    NormalizeVector (normal);
	    GenTexCoords (model, previcross, ia);
	    glVertex3fv (a->vtx [ia]);
	    GenTexCoords (model, icross, ib);
	    glVertex3fv (b->vtx [ib]);
	    glNormal3fv (normal);
	 }
	 if (flags & STITCH_LOOPS) {
	    Vector normal;
	    AddVector (a->normal [na-1], b ->normal [nb-1], normal);
	    AddVector (a->normal [0], normal, normal);
	    AddVector (b->normal [0], normal, normal);
	    NormalizeVector (normal);
	    glVertex3fv (a->vtx [0]);
	    GenTexCoords (model, previcross, na);
	    glVertex3fv (b->vtx [0]);
	    GenTexCoords (model, icross, nb);
	    glNormal3fv (normal);
	 }
      }
      glEnd ();
      icross++;
      previcross = (previcross+1) % ncross;
   }
}

/*--------------------------------------------------------------------------
 *
 *  Main Procedure for generating generic cylinders
 *
 *--------------------------------------------------------------------------*/

int
GenericCylinder (Tcl_Interp *interp, int argc, char* argv [])
{

#define ERRMSG(msg) \
   { Tcl_AppendResult (interp, (msg), (char*) NULL);\
     result = TCL_ERROR;\
     goto done; }

#define ERRMSG2(msg1, msg2) \
   { Tcl_AppendResult (interp, (msg1), (msg2), (char*) NULL);\
     result = TCL_ERROR;\
     goto done; }

#define I { \
      1, 0, 0, 0, \
      0, 1, 0, 0, \
      0, 0, 1, 0, \
      0, 0, 0, 1  \
   }

   int result = TCL_OK;

   Matrix transf = I;
   CrossSection *currentCross = NewCrossSection ();
   Model* model = NewModel ();
   int iarg;
   int dlist = -1;

   model->flags = STITCH_LOOPS | SHADE_SMOOTH;

   AddCrossSectionVtx (currentCross, -1.0, -1.0, 0.0);
   AddCrossSectionVtx (currentCross, -1.0, 1.0, 0.0);
   AddCrossSectionVtx (currentCross, 1.0, 1.0, 0.0);
   AddCrossSectionVtx (currentCross, 1.0, -1.0, 0.0);

   for (iarg = 2; iarg < argc; iarg++) {	
      int len = strlen (argv [iarg]);
      if (strncmp (argv [iarg], "-plot", len) == 0) {
	 CrossSection* cross = NewCrossSection ();
	 DupCrossSection (currentCross, cross);
	 TransformCrossSection (cross, transf);
	 AddModelCrossSection (model, cross);
      }
      else if (strncmp (argv [iarg], "-displaylist", len) == 0) {
	 iarg++;
         if (strcmp (argv [iarg], "none") == 0) {
            dlist = 0;
         }
	 else {
	    result = Tcl_GetInt (interp, argv [iarg], &dlist);
	    if (result != TCL_OK) goto done;
	 }
      }
      else if (strncmp (argv [iarg], "-cross", len) == 0) {
	 FreeCrossSection (currentCross);
	 currentCross = NewCrossSection ();
	 while (iarg+3 < argc && !isalpha (argv [iarg+1][1])) {
	    double x, y, z;
	    if (Tcl_GetDouble (interp, argv [iarg+1], &x) != TCL_OK ||
		Tcl_GetDouble (interp, argv [iarg+2], &y) != TCL_OK ||
		Tcl_GetDouble (interp, argv [iarg+3], &z) != TCL_OK) {
	       ERRMSG ("\n error in -cross");
	    }
	    AddCrossSectionVtx (currentCross, x, y, z);
	    iarg += 3;
	 }
      }
      else if (strncmp (argv [iarg], "-polygon", len) == 0) {
	 double radius;
	 int nsides;
	 if (iarg+2 >= argc ||
	     Tcl_GetDouble (interp, argv [iarg+1], &radius) != TCL_OK ||
	     Tcl_GetInt (interp, argv [iarg+2], &nsides) != TCL_OK) {
	    ERRMSG ("\nError in -polygon");
	 }
	 iarg += 2;
	 FreeCrossSection (currentCross);
	 currentCross = PolygonCrossSection (radius, nsides);
      }
      else if (strncmp (argv [iarg], "-stitch", len) == 0) {
	 if (iarg+1 >= argc) ERRMSG ("No value for -stitch");
	 iarg++;
	 len = strlen (argv [iarg]);
	 if (strncmp (argv [iarg], "both", len) == 0) {
	    model->flags |= STITCH_LOOPS | STITCH_ENDS;
	 } else if (strncmp (argv [iarg], "none", len) == 0) {
	    model->flags &= ~(STITCH_LOOPS | STITCH_ENDS);
	 } else if (strncmp (argv [iarg], "ends", len) == 0) {
	    model->flags |= STITCH_ENDS;
	    model->flags &= ~STITCH_LOOPS;
	 } else if (strncmp (argv [iarg], "loops", len) == 0) {
	    model->flags &= ~STITCH_ENDS;
	    model->flags |= STITCH_LOOPS;
	 } else {
	    ERRMSG2 ("Should be 'both', 'none', 'loops' or 'ends':",
		     argv [iarg]);
	 }
      }
      else if (strncmp (argv [iarg], "-shade", len) == 0) {
	 if (iarg+1 >= argc) ERRMSG ("No value for -shade");
	 iarg++;
	 len = strlen (argv [iarg]);
	 if (strncmp (argv [iarg], "smooth", len) == 0) {
	    model->flags |= SHADE_SMOOTH;
	 } else if (strncmp (argv [iarg], "flat", len) == 0) {
	    model->flags &= ~SHADE_SMOOTH;
	 } else {
	    ERRMSG2 ("Should be 'flat' or 'smooth':",
		     argv [iarg]);
	 }
      }
      else if (strncmp (argv [iarg], "-texgen", len) == 0) {
	 double smin, smax, tmin, tmax;
	 model->flags |= TEXTURE_GEN;
	 if (iarg+4 >= argc ||
	     Tcl_GetDouble (interp, argv [iarg+1], &smin) != TCL_OK ||
	     Tcl_GetDouble (interp, argv [iarg+2], &smax) != TCL_OK ||
	     Tcl_GetDouble (interp, argv [iarg+3], &tmin) != TCL_OK ||
	     Tcl_GetDouble (interp, argv [iarg+4], &tmax) != TCL_OK)
	    ERRMSG ("\nError in -rotate");
	 iarg += 4;
	 model->smin = smin;
	 model->smax = smax;
	 model->tmin = tmin;
	 model->tmax = tmax;
      }
      else if (strncmp (argv [iarg], "-identity", len) == 0) {
	 Matrix ident = I;
	 memcpy (transf, ident, sizeof (Matrix));
      }
      else if (strncmp (argv [iarg], "-rotate", len) == 0) {
	 double angle, x, y, z, norm, sint, cost;
	 Matrix m = I;
	 if (iarg+4 >= argc ||
	     Tcl_GetDouble (interp, argv [iarg+1], &angle) != TCL_OK ||
	     Tcl_GetDouble (interp, argv [iarg+2], &x) != TCL_OK ||
	     Tcl_GetDouble (interp, argv [iarg+3], &y) != TCL_OK ||
	     Tcl_GetDouble (interp, argv [iarg+4], &z) != TCL_OK)
	    ERRMSG ("\nError in -rotate");
	 iarg += 4;
	 norm = sqrt (x*x+y*y+z*z);
	 if (norm == 0.0) ERRMSG ("Null Vector");
	 x /= norm; y /= norm; z /= norm;
	 angle *= PI/180;
	 sint = sin(angle);
	 cost = cos(angle);
	 m [0][0] = x*x + cost * (1 - x*x) + sint * 0; 
	 m [0][1] = x*y + cost * (0 - x*y) + sint * (-z); 
	 m [0][2] = x*z + cost * (0 - x*z) + sint * (y);
	 m [1][0] = y*x + cost * (0 - y*x) + sint * (z);
	 m [1][1] = y*y + cost * (1 - y*y) + sint * 0;
	 m [1][2] = y*z + cost * (0 - y*z) + sint * (-x);
	 m [2][0] = z*x + cost * (0 - z*x) + sint * (-y);
	 m [2][1] = z*y + cost * (0 - z*y) + sint * (x);
	 m [2][2] = z*z + cost * (1 - z*z) + sint * 0; 
	 MatrixMult (m, transf, transf);
      }
      else if (strncmp (argv [iarg], "-translate", len) == 0) {
	 double x, y, z;
	 Matrix m = I;
	 if (iarg+3 >= argc ||
	     Tcl_GetDouble (interp, argv [iarg+1], &x) != TCL_OK ||
	     Tcl_GetDouble (interp, argv [iarg+2], &y) != TCL_OK ||
	     Tcl_GetDouble (interp, argv [iarg+3], &z) != TCL_OK)
	    ERRMSG ("\nError in -translate");
	 iarg += 3;
	 m [0][3] = x;  m [1][3] = y; m [2][3] = z;
	 MatrixMult (m, transf, transf);
      }
      else if (strncmp (argv [iarg], "-scale", len) == 0) {
	 double x, y, z;
	 Matrix m = I;
	 if (iarg+3 >= argc ||
	     Tcl_GetDouble (interp, argv [iarg+1], &x) != TCL_OK ||
	     Tcl_GetDouble (interp, argv [iarg+2], &y) != TCL_OK ||
	     Tcl_GetDouble (interp, argv [iarg+3], &z) != TCL_OK)
	    ERRMSG ("\nError in -scale");
	 iarg += 3;
	 m [0][0] = x;  m [1][1] = y; m [2][2] = z;
	 MatrixMult (m, transf, transf);
      }
      else {
	 ERRMSG2 ("No such option: ", argv [iarg]);
      }
   }

done:
   if (result == TCL_OK) {
      char buf [100];
      if (dlist == -1) dlist = glGenLists (1);
      if (dlist != 0) glNewList (dlist, GL_COMPILE);
      ComputeModelNormals (model);
      RenderModelQuadStrip (model);
      if (dlist != 0) {
         glEndList ();
         sprintf (buf, "displaylist %d", dlist);
      }
      Tcl_AppendElement (interp, buf);
      ComputeModelBBox (interp, model);
   }
   FreeModel (model);
   return result;
}