for (i = 1; i <= model->numvertices; i++) {
    if (maxx < model->vertices[3 * i + 0])
      maxx = model->vertices[3 * i + 0];
    if (minx > model->vertices[3 * i + 0])
      minx = model->vertices[3 * i + 0];

    if (maxy < model->vertices[3 * i + 1])
      maxy = model->vertices[3 * i + 1];
    if (miny > model->vertices[3 * i + 1])
      miny = model->vertices[3 * i + 1];

    if (maxz < model->vertices[3 * i + 2])
      maxz = model->vertices[3 * i + 2];
    if (minz > model->vertices[3 * i + 2])
      minz = model->vertices[3 * i + 2];
  }

  /* calculate model width, height, and depth */
  dimensions[0] = glmAbs(maxx) + glmAbs(minx);
  dimensions[1] = glmAbs(maxy) + glmAbs(miny);
  dimensions[2] = glmAbs(maxz) + glmAbs(minz);
}

/* glmScale: Scales a model by a given amount.
 * 
 * model - properly initialized GLMmodel structure
 * scale - scalefactor (0.5 = half as large, 2.0 = twice as large)
 */
GLvoid glmScale(GLMmodel* model, GLfloat scale)
{
  GLuint i;

  for (i = 1; i <= model->numvertices; i++) {
    model->vertices[3 * i + 0] *= scale;
    model->vertices[3 * i + 1] *= scale;
    model->vertices[3 * i + 2] *= scale;
  }
}

/* glmReverseWinding: Reverse the polygon winding for all polygons in
 * this model.  Default winding is counter-clockwise.  Also changes
 * the direction of the normals.
 * 
 * model - properly initialized GLMmodel structure 
 */
GLvoid glmReverseWinding(GLMmodel* model)
{
  GLuint i, swap;

  assert(model);

  for (i = 0; i < model->numtriangles; i++) {
    swap = T(i).vindices[0];
    T(i).vindices[0] = T(i).vindices[2];
    T(i).vindices[2] = swap;

    if (model->numnormals) {
      swap = T(i).nindices[0];
      T(i).nindices[0] = T(i).nindices[2];
      T(i).nindices[2] = swap;
    }

    if (model->numtexcoords) {
      swap = T(i).tindices[0];
      T(i).tindices[0] = T(i).tindices[2];
      T(i).tindices[2] = swap;
    }
  }

  /* reverse facet normals */
  for (i = 1; i <= model->numfacetnorms; i++) {
    model->facetnorms[3 * i + 0] = -model->facetnorms[3 * i + 0];
    model->facetnorms[3 * i + 1] = -model->facetnorms[3 * i + 1];
    model->facetnorms[3 * i + 2] = -model->facetnorms[3 * i + 2];
  }

  /* reverse vertex normals */
  for (i = 1; i <= model->numnormals; i++) {
    model->normals[3 * i + 0] = -model->normals[3 * i + 0];
    model->normals[3 * i + 1] = -model->normals[3 * i + 1];
    model->normals[3 * i + 2] = -model->normals[3 * i + 2];
  }
}

/* glmFacetNormals: Generates facet normals for a model (by taking the
 * cross product of the two vectors derived from the sides of each
 * triangle).  Assumes a counter-clockwise winding.
 *
 * model - initialized GLMmodel structure
 */
GLvoid glmFacetNormals(GLMmodel* model)
{
  GLuint  i;
  GLfloat u[3];
  GLfloat v[3];
  
  assert(model);
  assert(model->vertices);

  /* clobber any old facetnormals */
  if (model->facetnorms)
    free(model->facetnorms);

  /* allocate memory for the new facet normals */
  model->numfacetnorms = model->numtriangles;
  model->facetnorms = (GLfloat*)malloc(sizeof(GLfloat) *
				       3 * (model->numfacetnorms + 1));

  for (i = 0; i < model->numtriangles; i++) {
    model->triangles[i].findex = i+1;

    u[0] = model->vertices[3 * T(i).vindices[1] + 0] -
           model->vertices[3 * T(i).vindices[0] + 0];
    u[1] = model->vertices[3 * T(i).vindices[1] + 1] -
           model->vertices[3 * T(i).vindices[0] + 1];
    u[2] = model->vertices[3 * T(i).vindices[1] + 2] -
           model->vertices[3 * T(i).vindices[0] + 2];

    v[0] = model->vertices[3 * T(i).vindices[2] + 0] -
           model->vertices[3 * T(i).vindices[0] + 0];
    v[1] = model->vertices[3 * T(i).vindices[2] + 1] -
           model->vertices[3 * T(i).vindices[0] + 1];
    v[2] = model->vertices[3 * T(i).vindices[2] + 2] -
           model->vertices[3 * T(i).vindices[0] + 2];

    glmCross(u, v, &model->facetnorms[3 * (i+1)]);
    glmNormalize(&model->facetnorms[3 * (i+1)]);
  }
}

/* glmVertexNormals: Generates smooth vertex normals for a model.
 * First builds a list of all the triangles each vertex is in.  Then
 * loops through each vertex in the the list averaging all the facet
 * normals of the triangles each vertex is in.  Finally, sets the
 * normal index in the triangle for the vertex to the generated smooth
 * normal.  If the dot product of a facet normal and the facet normal
 * associated with the first triangle in the list of triangles the
 * current vertex is in is greater than the cosine of the angle
 * parameter to the function, that facet normal is not added into the
 * average normal calculation and the corresponding vertex is given
 * the facet normal.  This tends to preserve hard edges.  The angle to
 * use depends on the model, but 90 degrees is usually a good start.
 *
 * model - initialized GLMmodel structure
 * angle - maximum angle (in degrees) to smooth across
 */
GLvoid glmVertexNormals(GLMmodel* model, GLfloat angle)
{
  GLMnode*  node;
  GLMnode*  tail;
  GLMnode** members;
  GLfloat*  normals;
  GLuint    numnormals;
  GLfloat   average[3];
  GLfloat   dot, cos_angle;
  GLuint    i, avg;

  assert(model);
  assert(model->facetnorms);

  /* calculate the cosine of the angle (in degrees) */
  cos_angle = fcos(angle * M_PI / 180.0f);

  /* nuke any previous normals */
  if (model->normals)
    free(model->normals);

  /* allocate space for new normals */
  model->numnormals = model->numtriangles * 3; /* 3 normals per triangle */
  model->normals = (GLfloat*)malloc(sizeof(GLfloat)* 3* (model->numnormals+1));

  /* allocate a structure that will hold a linked list of triangle
     indices for each vertex */
  members = (GLMnode**)malloc(sizeof(GLMnode*) * (model->numvertices + 1));
  for (i = 1; i <= model->numvertices; i++)
    members[i] = NULL;
  
  /* for every triangle, create a node for each vertex in it */
  for (i = 0; i < model->numtriangles; i++) {
    node = (GLMnode*)malloc(sizeof(GLMnode));
    node->index = i;
    node->next  = members[T(i).vindices[0]];
    members[T(i).vindices[0]] = node;

    node = (GLMnode*)malloc(sizeof(GLMnode));
    node->index = i;
    node->next  = members[T(i).vindices[1]];
    members[T(i).vindices[1]] = node;

    node = (GLMnode*)malloc(sizeof(GLMnode));
    node->index = i;
    node->next  = members[T(i).vindices[2]];
    members[T(i).vindices[2]] = node;
  }

  /* calculate the average normal for each vertex */
  numnormals = 1;
  for (i = 1; i <= model->numvertices; i++) {
    /* calculate an average normal for this vertex by averaging the
       facet normal of every triangle this vertex is in */
    node = members[i];
    if (!node)
      fprintf(stderr, "glmVertexNormals(): vertex w/o a triangle\n");
    average[0] = 0.0; average[1] = 0.0; average[2] = 0.0;
    avg = 0;
    while (node) {
      /* only average if the dot product of the angle between the two
         facet normals is greater than the cosine of the threshold
         angle -- or, said another way, the angle between the two
         facet normals is less than (or equal to) the threshold angle */
      dot = glmDot(&model->facetnorms[3 * T(node->index).findex],
 		    &model->facetnorms[3 * T(members[i]->index).findex]);
      if (dot > cos_angle) {
	node->averaged = GL_TRUE;
	average[0] += model->facetnorms[3 * T(node->index).findex + 0];
	average[1] += model->facetnorms[3 * T(node->index).findex + 1];
	average[2] += model->facetnorms[3 * T(node->index).findex + 2];
	avg = 1;			/* we averaged at least one normal! */
      } else {
	node->averaged = GL_FALSE;
      }
      node = node->next;
    }

    if (avg) {
      /* normalize the averaged normal */
      glmNormalize(average);

      /* add the normal to the vertex normals list */
      model->normals[3 * numnormals + 0] = average[0];
      model->normals[3 * numnormals + 1] = average[1];
      model->normals[3 * numnormals + 2] = average[2];
      avg = numnormals;
      numnormals++;
    }

    /* set the normal of this vertex in each triangle it is in */
    node = members[i];
    while (node) {
      if (node->averaged) {
	/* if this node was averaged, use the average normal */
	if (T(node->index).vindices[0] == i)
	  T(node->index).nindices[0] = avg;
	else if (T(node->index).vindices[1] == i)
	  T(node->index).nindices[1] = avg;
	else if (T(node->index).vindices[2] == i)
	  T(node->index).nindices[2] = avg;
      } else {
	/* if this node wasn't averaged, use the facet normal */
	model->normals[3 * numnormals + 0] = 
	  model->facetnorms[3 * T(node->index).findex + 0];
	model->normals[3 * numnormals + 1] = 
	  model->facetnorms[3 * T(node->index).findex + 1];
	model->normals[3 * numnormals + 2] = 
	  model->facetnorms[3 * T(node->index).findex + 2];
	if (T(node->index).vindices[0] == i)
	  T(node->index).nindices[0] = numnormals;
	else if (T(node->index).vindices[1] == i)
	  T(node->index).nindices[1] = numnormals;
	else if (T(node->index).vindices[2] == i)
	  T(node->index).nindices[2] = numnormals;
	numnormals++;
      }
      node = node->next;
    }
  }
  
  model->numnormals = numnormals - 1;

  /* free the member information */
  for (i = 1; i <= model->numvertices; i++) {
    node = members[i];
    while (node) {
      tail = node;
      node = node->next;
      free(tail);
    }
  }
  free(members);

  /* pack the normals array (we previously allocated the maximum
     number of normals that could possibly be created (numtriangles *
     3), so get rid of some of them (usually alot unless none of the
     facet normals were averaged)) */
  normals = model->normals;
  model->normals = (GLfloat*)malloc(sizeof(GLfloat)* 3* (model->numnormals+1));
  for (i = 1; i <= model->numnormals; i++) {
    model->normals[3 * i + 0] = normals[3 * i + 0];
    model->normals[3 * i + 1] = normals[3 * i + 1];
    model->normals[3 * i + 2] = normals[3 * i + 2];
  }
  free(normals);
}


/* glmLinearTexture: Generates texture coordinates according to a
 * linear projection of the texture map.  It generates these by
 * linearly mapping the vertices onto a square.
 *
 * model - pointer to initialized GLMmodel structure
 */
GLvoid glmLinearTexture(GLMmodel* model)
{
  GLMgroup *group;
  GLfloat dimensions[3];
  GLfloat x, y, scalefactor;
  GLuint i;
  
  assert(model);

  if (model->texcoords)
    free(model->texcoords);
  model->numtexcoords = model->numvertices;
  model->texcoords=(GLfloat*)malloc(sizeof(GLfloat)*2*(model->numtexcoords+1));
  
  glmDimensions(model, dimensions);
  scalefactor = 2.0f / 
    glmAbs(glmMax(glmMax(dimensions[0], dimensions[1]), dimensions[2]));

  /* do the calculations */
  for(i = 1; i <= model->numvertices; i++) {
    x = model->vertices[3 * i + 0] * scalefactor;
    y = model->vertices[3 * i + 2] * scalefactor;
    model->texcoords[2 * i + 0] = (x + 1.0f) / 2.0f;
    model->texcoords[2 * i + 1] = (y + 1.0f) / 2.0f;
  }
  
  /* go through and put texture coordinate indices in all the triangles */
  group = model->groups;
  while(group) {
    for(i = 0; i < group->numtriangles; i++) {
      T(group->triangles[i]).tindices[0] = T(group->triangles[i]).vindices[0];
      T(group->triangles[i]).tindices[1] = T(group->triangles[i]).vindices[1];
      T(group->triangles[i]).tindices[2] = T(group->triangles[i]).vindices[2];
    }    
    group = group->next;
  }

#if 0
  printf("glmLinearTexture(): generated %d linear texture coordinates\n",
	  model->numtexcoords);
#endif
}

/* glmSpheremapTexture: Generates texture coordinates according to a
 * spherical projection of the texture map.  Sometimes referred to as
 * spheremap, or reflection map texture coordinates.  It generates
 * these by using the normal to calculate where that vertex would map
 * onto a sphere.  Since it is impossible to map something flat
 * perfectly onto something spherical, there is distortion at the
 * poles.  This particular implementation causes the poles along the X
 * axis to be distorted.
 *
 * model - pointer to initialized GLMmodel structure
 */
GLvoid glmSpheremapTexture(GLMmodel* model)
{
  GLMgroup* group;
  GLfloat theta, phi, rho, x, y, z, r;
  GLuint i;
  
  assert(model);
  assert(model->normals);

  if (model->texcoords)
    free(model->texcoords);
  model->numtexcoords = model->numnormals;
  model->texcoords=(GLfloat*)malloc(sizeof(GLfloat)*2*(model->numtexcoords+1));
     
  for (i = 1; i <= model->numnormals; i++) {
    z = model->normals[3 * i + 0];	/* re-arrange for pole distortion */
    y = model->normals[3 * i + 1];
    x = model->normals[3 * i + 2];
    r = fsqrt((x * x) + (y * y));
    rho = fsqrt((r * r) + (z * z));
      
    if(r == 0.0f) {
	theta = 0.0f;
	phi = 0.0f;
    } else {
      if(z == 0.0f)
	phi = 3.14159265f / 2.0f;
      else
	phi = facos(z / rho);