/*
 * Copyright (C) 1997-2009, R3vis Corporation.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * Library General Public License for more details.
 *
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
 * USA, or visit http://www.gnu.org/copyleft/lgpl.html.
 *
 * Original Contributor:
 *   Wes Bethel, R3vis Corporation, Marin County, California
 *   http://www.r3vis.com/
 * Additional Contributor(s):
 *
 * The OpenRM project is located at http://openrm.sourceforge.net/.
 */
/*
 * $Id: rmqdrix.c,v 1.7 2005/06/06 02:04:29 wes Exp $
 * Version: $Name: v180-alpha-02 $
 * $Revision: 1.7 $
 * $Log: rmqdrix.c,v $
 * Revision 1.7  2005/06/06 02:04:29  wes
 * Lots of small additions to clean up compiler warnings.
 *
 * Revision 1.6  2005/02/19 16:40:20  wes
 * Distro sync and consolidation.
 * Repairs to fix memory leak associated with repeated calls to rmPipeNew,
 * rmPipeMakeCurrent, rmPipeClose.
 *
 * Revision 1.5  2005/01/23 17:00:22  wes
 * Copyright updated to 2005.
 *
 * Revision 1.4  2004/03/18 15:49:00  wes
 * Minor tweaking on cone normals in an attempt to eliminate some
 * lighting errors visible in the lightModes gordo demo that show up
 * on nvidia hardware. The changes didn't eliminate the pops, but the
 * normals are now more correct.
 *
 * Revision 1.3  2004/01/16 16:48:35  wes
 * Updated copyright line for 2004.
 *
 * Revision 1.2  2003/02/02 02:07:15  wes
 * Updated copyright to 2003.
 *
 * Revision 1.1.1.1  2003/01/28 02:15:23  wes
 * Manual rebuild of rm150 repository.
 *
 * Revision 1.10  2003/01/16 22:21:17  wes
 * Updated all source files to reflect new organization of header files:
 * all header files formerly located in include/rmaux, include/rmi, include/rmv
 * are now located in include/rm.
 *
 * Revision 1.9  2002/08/29 22:20:32  wes
 *
 * Massive upgrade to accommodate dynamic object reallocation within
 * the component manager, and within the context cache. Use the
 * debug #define DEBUG_LEVEL DEBUG_REALLOC_TRACE to get a printf
 * whenever a realloc occurs. With this upgrade, there are no
 * OpenRM limits on the size of the scene graph. There will be external
 * limits, such as the amount of RAM and the amount of space available
 * to your OpenGL implementation.
 *
 * Revision 1.8  2002/04/30 19:33:26  wes
 * Updated copyright dates.
 *
 * Revision 1.7  2001/07/15 17:13:46  wes
 * Added private_rmCacheDeleteQuadrics to delete display lists associated
 * with all quadrics objects.
 *
 * Revision 1.6  2001/03/31 17:12:39  wes
 * v1.4.0-alpha-2 checkin.
 *
 * Revision 1.5  2000/12/03 22:35:38  wes
 * Mods for thread safety.
 *
 * Revision 1.4  2000/08/31 02:05:28  wes
 * Added attribute push/pop to display-listed code for spheres, cones
 * and cylinders. This fixes a bug whereby glFrontFace settings inside
 * those routines could adversely and unexpectedly affect subsequent
 * draw calls.
 *
 * Revision 1.3  2000/08/23 23:28:16  wes
 * Minor edits for readability.
 *
 * Revision 1.2  2000/04/20 16:29:47  wes
 * Documentation additions/enhancements, some code rearragement.
 *
 * Revision 1.1.1.1  2000/02/28 21:29:40  wes
 * OpenRM 1.2 Checkin
 *
 * Revision 1.1.1.1  2000/02/28 17:18:48  wes
 * Initial entry - pre-RM120 release, source base for OpenRM 1.2.
 *
 */


#include <rm/rm.h>
#include "rmprivat.h"

/* 
 * 5/3/99
 * the following #define is used to turn on or off normalization of
 * surface normals for quadrics (spheres, cones, cylinders). for the
 * most part, this is needed because we build canonical representations
 * of these objects, then use what is effectively a scaling operation
 * to get them the right size and in the right place. 
 *
 * one problem with how we're doing this is that if GL_NORMALIZE is
 * enabled prior to calling the canonical build routines, it will be
 * turned off (only when QUADRICS_NORMALIZE == 1). however, this may
 * not be a big deal if we call private_rmInitQUadrics() prior to calling
 * any user code (where they might monkey with GL_NORMALIZE)
 */
#define QUADRICS_NORMALIZE  1

/*
 * 10/2000
 * define max number of tesselated models for each style of quadrics
 * object. the model flags for cones and cylinders (RM_CONES_4, etc.)
 * are defined in rmdefs.h - and are interpreted as both the number
 * of radial samples as well as an index into the model table for
 * display list caching. space is allocated for up to 360*2 models
 * for each of cones & cylinders (720*2 total) - which is somewhat
 * wasteful of space, but only a handful of models are actually
 * built and used. those that are built and available for use are listed
 * in coneModels[] and cylinderModels[] below.
 */

#define MAX_SPHERE_MODELS   5
#define MAX_CONE_MODELS     360
#define MAX_CYLINDER_MODELS 360

/* cone models */
static int coneModels[] =
{
    RM_CONES_4,
    RM_CONES_8,
    RM_CONES_12,
    RM_CONES_16,
    RM_CONES_32,
    RM_CONES_64,
    RM_CONES_128
};

static int numConeModels = sizeof(coneModels) / sizeof(int);

/* cylinder models */
static int cylinderModels[] =
{
    RM_CYLINDERS_4,
    RM_CYLINDERS_8,
    RM_CYLINDERS_12,
    RM_CYLINDERS_16,
    RM_CYLINDERS_32,
    RM_CYLINDERS_64,
    RM_CYLINDERS_128
};
static int    numCylinderModels = sizeof(cylinderModels) / sizeof(int);

/* sphere models */
typedef struct triangle
{
    RMvertex3D pt[3];
} triangle;

typedef struct sphereobj
{
    int       ntriangles;
    triangle *list;
} sphereobj;

typedef RMvertex3D point; /* PRIVATE */

/* Six equidistant points lying on the unit sphere */
#define XPLUS {  1,  0,  0 }    /*  X */
#define XMIN  { -1,  0,  0 }    /* -X */
#define YPLUS {  0,  1,  0 }    /*  Y */
#define YMIN  {  0, -1,  0 }    /* -Y */
#define ZPLUS {  0,  0,  1 }    /*  Z */
#define ZMIN  {  0,  0, -1 }    /* -Z */

static RMvertex3D *sphere_verts, *sphere_norms;
static int         sphere_facet_count;
static int         triangles_per_sphere = 0;
static int         tri_ptr;

triangle octahedron[] = {
    { { XPLUS, ZPLUS, YPLUS } },
    { { YPLUS, ZPLUS, XMIN  } },
    { { XMIN , ZPLUS, YMIN  } },
    { { YMIN , ZPLUS, XPLUS } },
    { { XPLUS, YPLUS, ZMIN  } },
    { { YPLUS, XMIN , ZMIN  } },
    { { XMIN , YMIN , ZMIN  } },
    { { YMIN , XPLUS, ZMIN  } }
};

static sphereobj oct = {sizeof(octahedron) / sizeof(octahedron[0]), &octahedron[0]};


/* PRIVATE declarations: build and cache display lists for a set of predefined, 
 * semi-procedural quadric surfaces: cones, cylinders and spheres
 */
GLuint        private_rmBuildCone (int subdivisions, int flipped);
GLuint        private_rmBuildCylinder (int subdivisions, int flipped);
GLuint        private_rmBuildSphere (void);
int           private_initSphereTables (int n);
int           private_freeSphereTables (void);
static point *normalize ( point *);
static point *midpoint (point *, point *);
void          private_yaxis_to_dir (RMmatrix *m, RMvertex3D *dir);


/* precomputed trig tables */
#define TRIG_RES 360
extern double cos_table[TRIG_RES], sin_table[TRIG_RES];


/* PRIVATE */
void
private_rmInitQuadrics (RMcontextCache *cache)
{
    /*
     * this routine gets called in order to build cached models
     * for cones, spheres and cylinders. 
     * 10/2000 - updated to 1) malloc space inside a context cache
     * for display lists indices and 2) to deposit precomputed
     * display lists into the context cache.
     */
    
    int        i;

    cache->sphereIDs = (GLuint *)malloc(sizeof(GLuint)*MAX_SPHERE_MODELS);
    cache->coneIDs = (GLuint *)malloc(sizeof(GLuint)*MAX_CONE_MODELS);
    cache->flipConeIDs = (GLuint *)malloc(sizeof(GLuint)*MAX_CONE_MODELS);
    cache->cylinderIDs = (GLuint *)malloc(sizeof(GLuint)*MAX_CYLINDER_MODELS);
    cache->flipCylinderIDs = (GLuint *)malloc(sizeof(GLuint)*MAX_CYLINDER_MODELS);

    /*
     * mark the tables containing OpenGL display list indices as
     * uninitialized.
     */
    for (i = 0; i < TRIG_RES; i++)
    {
	cache->coneIDs[i] = cache->cylinderIDs[i] =  -1;
	cache->flipConeIDs[i] = cache->flipCylinderIDs[i] =  -1;
    }

    /*
     * build canonical cones.
     * range of cones from 4 to 16 subdivisions around a circle
     */
    for (i = 0; i < numConeModels; i++)
    {
	int indx = coneModels[i];	/* RM_CONES_4, RM_CONES_8, etc. */

	cache->coneIDs[indx] = private_rmBuildCone(indx, 0); 
	cache->flipConeIDs[indx] = private_rmBuildCone(indx, 1); 
    }
    
    /*
     * build canonical spheres
     *  value    number of triangles
     *    1             8
     *    2             32
     *    3             128
     *    4             512
     *    n             2*4^n
     */

    for (i = 1; i < MAX_SPHERE_MODELS; i++)
    {
	private_initSphereTables(i);
	cache->sphereIDs[i] = private_rmBuildSphere();
	private_freeSphereTables();
    }

    /* build canonical cylinders */
    for (i = 0; i < numCylinderModels; i++)
    {
	int indx = cylinderModels[i];	/* RM_CYLINDERS_4, RM_CYLINDERS_8, etc. */

	cache->cylinderIDs[indx] = private_rmBuildCylinder(indx, 0);
	cache->flipCylinderIDs[indx] = private_rmBuildCylinder(indx, 1);
    }
}

/* PRIVATE */
void
private_rmCacheDeleteQuadrics (RMcontextCache *cache)
{
    /*
     * 7/4/01 w.bethel
     * this routine will release all OpenGL display lists associated
     * with a context cache. it is assumed that the caller owns the
     * OpenGL rendering context.
     */

    GLuint *lists;
    int i;

    lists = cache->sphereIDs;
    if (lists != NULL)
    {
	for (i=0;i<MAX_SPHERE_MODELS;i++)
	    if (glIsList(lists[i]) == GL_TRUE)
		glDeleteLists(lists[i], 1);
    }

    lists = cache->coneIDs;
    if (lists != NULL)
    {
	for (i=0;i<MAX_CONE_MODELS;i++)
	    if (glIsList(lists[i]) == GL_TRUE)
		glDeleteLists(lists[i],1);
    }

    lists = cache->flipConeIDs;
    if (lists != NULL)
    {
	for (i=0;i<MAX_CONE_MODELS;i++)
	    if (glIsList(lists[i]) == GL_TRUE)
		glDeleteLists(lists[i], 1);
    }

    lists = cache->cylinderIDs;
    if (lists != NULL)
    {
	for (i=0;i<MAX_CYLINDER_MODELS;i++)
	    if (glIsList(lists[i]) == GL_TRUE)
		glDeleteLists(lists[i], 1);
    }

    lists = cache->flipCylinderIDs;
    if (lists != NULL)
    {
	for (i=0;i<MAX_CYLINDER_MODELS;i++)
	    if (glIsList(lists[i]) == GL_TRUE)
		glDeleteLists(lists[i], 1);
    }

    free(cache->sphereIDs);
    cache->sphereIDs = 0;
    
    free(cache->coneIDs);
    cache->coneIDs = 0;
    
    free(cache->flipConeIDs);
    cache->flipConeIDs = 0;
    
    free(cache->cylinderIDs);
    cache->cylinderIDs = 0;
    
    free(cache->flipCylinderIDs);
    cache->flipCylinderIDs = 0;
}

/* PRIVATE
 *
 * rmuCone is a private, internal routine. it's purpose is to
 * draw a cone given a cone description (input parms). this routine
 * is not intended to be called by applications. it is intended to
 * be called from the "cone primitive" draw func, is similar in function
 * gluCylinder.
 */
void
rmuCone (void (*colorfunc) (const float *),
	 float *color,
	 RMvertex3D *p1,
	 RMvertex3D *p2,
	 float radius1,
	 int subdivisions,
	 RMpipe *p)
{
    float      fmag;
    float      t, dt;
    float      radius2 = 0.0;
    double     mag;
    RMvertex3D dir;
    RMmatrix   rot, trans;
    GLuint     *list;

    if (subdivisions > TRIG_RES)
	subdivisions = TRIG_RES;
    
    t = 0.0;
    dt = (float)(TRIG_RES) / (float)(subdivisions);

    if (colorfunc)
	(*colorfunc)(color);
    
    rmVertex3DDiff(p2, p1, &dir);
    rmVertex3DMagNormalize(&dir, &mag);
    fmag = mag;

    if (mag == 0.0)
	mag = 1.0e5;

    if ((radius1 == 0.0F) && (radius2 == 0.0F))
    {
	glBegin(GL_LINES);
	glVertex3fv((GLfloat *)p1);
	glVertex3fv((GLfloat *)p2);
	glEnd();
	return;
    }

    private_yaxis_to_dir(&rot, &dir);
    rmMatrixTranspose(&rot, &rot);
    glPushMatrix();
    rmMatrixIdentity(&trans);
    
    trans.m[3][0] = p1->x;
    trans.m[3][1] = p1->y;
    trans.m[3][2] = p1->z;
    
    {
	RMmatrix m;

	rmMatrixIdentity(&m);
	m.m[0][0] = radius1;
	m.m[1][1] = fmag;
	m.m[2][2] = radius1;
	rmMatrixMultiply(&m, &rot, &m);		/* scale, then rotate.. */
	rmMatrixMultiply(&m, &trans, &m);	/* then translate */

	glMultMatrixf(&(m.m[0][0]));
    }
    
    if (rot.m[1][1] == -1.0F)
	list = p->contextCache->flipConeIDs; /* flipped */
    else
	list = p->contextCache->coneIDs; /* not flipped */

    if (list[subdivisions] != (unsigned int)(-1))
    {
	glCallList(list[subdivisions]);
	glPopMatrix();
	return;
    }
    else
	fprintf(stderr, " rmuCone(): no display list for this cone. \n");

    glPopMatrix();
}


/* PRIVATE */
GLuint
private_rmBuildCone (int subdivisions,
		     int flipped)
{
    int        i;
    float      t, dt;
    float      yconst = 0.5F;
    RMvertex3D v, n;
    GLuint     list;

    list = glGenLists(1);
    glNewList(list, GL_COMPILE);

    glPushAttrib(GL_POLYGON_BIT);

    /* get vertex winding */
    if (flipped)
	glFrontFace(GL_CW);
    else
	glFrontFace(GL_CCW);

#if QUADRICS_NORMALIZE
    glEnable(GL_NORMALIZE);
#endif

    t = 0.0;
    dt = (float)(TRIG_RES) / (float)(subdivisions);

    glBegin(GL_TRIANGLE_STRIP);

    for (i = 0; i <= subdivisions; i++)
    {
	int indx  = (int)t;

	indx = indx % TRIG_RES;
	
	n.x = v.x = cos_table[indx];
	n.z = v.z = sin_table[indx];
	v.y = 0.0;
	n.y = yconst;
	rmVertex3DNormalize(&n);

	glNormal3fv((GLfloat *)&n);
	glVertex3fv((GLfloat *)&v);

	indx = t + (dt * 0.5);
	indx = indx % TRIG_RES;

	v.x = v.z = 0.0;
	v.y = 1.0;

	n.x = cos_table[indx];
	n.y = yconst;
	n.z = sin_table[indx];

	rmVertex3DNormalize(&n);

	glNormal3fv((GLfloat *)&n);
	glVertex3fv((GLfloat *)&v);

	t += dt;
    }

    glEnd();

    /* do the cap */
    glBegin(GL_TRIANGLE_FAN);

    n.x = n.z = 0.0F;
    n.y = -1.0F;
    glNormal3fv((GLfloat *)&n);

    v.x = v.y = v.z = 0.0F;
    glVertex3fv((GLfloat *)&v);
    
    for (i = 0; i <= subdivisions; i++)
    {
	int indx = (int)t;

	indx = indx % TRIG_RES;

	v.x = cos_table[indx];
	v.z = sin_table[indx];
	glVertex3fv((GLfloat *)&v);

	t += dt;
    }

    glEnd();
#if QUADRICS_NORMALIZE
    glDisable(GL_NORMALIZE);
#endif

#if 0
    /* this hunk of code shows the normals */
    glDisable(GL_LIGHTING);
    glBegin(GL_LINES);

    t = 0.0;

    for (i = 0; i < subdivisions; i++, t += dt)
    {
	int indx = (int)t;

	indx = indx % TRIG_RES;
	
	n.x = v.x = cos_table[indx];
	n.z = v.z = sin_table[indx];
	v.y = 0.0;
	n.y = yconst;
	rmVertex3DNormalize(&n);

	glVertex3fv((GLfloat *)&v);

	v.x += n.x;
	v.y += n.y;
	v.z += n.z;
	glVertex3fv((GLfloat *)&v);

	indx = t + (dt * 0.5);
	indx = indx % TRIG_RES;

	v.x = v.z = 0.0;
	v.y = 1.0;

	n.x = cos_table[indx];
	n.z = sin_table[indx];
	n.y = yconst;
        rmVertex3DNormalize(&n);

	glVertex3fv((GLfloat *)&v);

	v.x += n.x;
	v.y += n.y;
	v.z += n.z;
	glVertex3fv((GLfloat *)&v);
	
	t += dt;
    }

    glEnd();
#endif

    glPopAttrib();
    
    glEndList();

    return(list);
}


/*
 * rmuCylinder is a replacement for gluCylinder, but uses one rather
 * than two radius values. it's purpose in life is to draw cylinders
 * given position & radius information. it is intended to be called
 * from the cylinder primitive draw routine.
 */
void
rmuCylinder (void (*colorfunc) (const float *),
	     float *color,
	     RMvertex3D *p1,
	     RMvertex3D *p2,
	     float radius1,
	     int subdivisions,
	     RMpipe *p)
{
    float      fmag;
    double     mag;
    RMvertex3D dir;
    RMmatrix   rot, trans, scale, m;
    GLuint     *list;

    if (colorfunc)
	(*colorfunc)(color);
    
    /* check for zero radius cylinders, draw them as a line segment */
    rmVertex3DDiff(p2, p1, &dir);
    rmVertex3DMagNormalize(&dir, &mag);
    fmag = mag;

    if (radius1 == 0.0F)
    {
	glBegin(GL_LINES);
	glVertex3fv((GLfloat *)p1);
	glVertex3fv((GLfloat *)p2);
	glEnd();
	return;
    }

    private_yaxis_to_dir(&rot, &dir);
    
    if (rot.m[1][1] == -1.0F)
	list = p->contextCache->flipCylinderIDs;  /* flipped */
    else
	list = p->contextCache->cylinderIDs; /* not flipped */
    
    rmMatrixTranspose(&rot, &rot);

    glPushMatrix();
    rmMatrixIdentity(&trans);
    
    trans.m[3][0] = p1->x;
    trans.m[3][1] = p1->y;
    trans.m[3][2] = p1->z;

    rmMatrixIdentity(&scale);
    scale.m[0][0] = scale.m[2][2] = radius1;
    scale.m[1][1] = fmag;

    rmMatrixMultiply(&scale, &rot, &m);	/* scale, then rotate.. */
    rmMatrixMultiply(&m, &trans, &m);	/* then translate */
     
    glMultMatrixf(&(m.m[0][0]));

    if (list[subdivisions] != (unsigned int)(-1))
	glCallList(list[subdivisions]);
    else
	rmWarning("rmuCylinder() warning: no display list for this cylinder. \n");

    glPopMatrix();
}


/* PRIVATE */
GLuint
private_rmBuildCylinder (int subdivisions,
			 int flipped)
{
    int        i;
    float      flip;
    float      t, dt;
    RMvertex3D v, n;
    GLuint     list;

    list = glGenLists(1);
    glNewList(list, GL_COMPILE);

    glPushAttrib(GL_POLYGON_BIT | GL_ENABLE_BIT);

    if (flipped == 1)
	flip = -1.0;
    else
	flip = 1.0;

    t = 0.;
    dt = (float)(TRIG_RES)/(float)(subdivisions);

#if QUADRICS_NORMALIZE
    glEnable(GL_NORMALIZE);
#endif

    /* get vertex winding */
    if (flipped)
	glFrontFace(GL_CW);
    else
	glFrontFace(GL_CCW);

    glBegin(GL_TRIANGLE_STRIP);

    for (i = 0; i <= subdivisions; i++, t+=dt)
    {
	int indx = (int)t;

	indx = indx % TRIG_RES;

	v.x = cos_table[indx] ;
	n.x = cos_table[indx] * flip;
	n.y = v.y = 0.0;
	
	v.z = sin_table[indx] ;
	n.z = sin_table[indx] ;

	glNormal3fv((float *)&n);
	glVertex3fv((float *)&v);

	v.x = cos_table[indx];
	v.y = 1.0F;
	v.z = sin_table[indx];

	glVertex3fv((float *)&v);
    }
    
    glEnd();
#if QUADRICS_NORMALIZE
    glDisable(GL_NORMALIZE);
#endif

    glPopAttrib();
    glEndList();

    return(list);
}


/*
 * rmuSphere is a replacement for gluSphere. it draws spheres based upon
 * position, size and color input information. it is not intended to be
 * called except from the sphere primitive draw routine.
 */
void
rmuSphere (void (*colorfunc) (const float *),
	   float *color,
	   RMvertex3D *p1,
	   float radius,
	   int model_switch,
	   RMpipe *p)
{
    RMmatrix m;
    GLuint *sphereLists;

    if (colorfunc)
	(*colorfunc)(color);

    rmMatrixIdentity(&m);
    m.m[0][0] = m.m[1][1] = m.m[2][2] = radius;

    m.m[3][0] = p1->x;
    m.m[3][1] = p1->y;
    m.m[3][2] = p1->z;

    glPushMatrix();
    glMultMatrixf(&(m.m[0][0]));

    /* want better range checking on "model switch" */

    sphereLists = p->contextCache->sphereIDs;

    if (sphereLists[model_switch] != (unsigned int)(-1))
	glCallList(sphereLists[model_switch]);
    else
	fprintf(stderr, " rmuSphere: no display list for this model_switch. \n");
    glPopMatrix();
}


/* PRIVATE */
GLuint
private_rmBuildSphere (void)
{
    /*
     * this routine is part of a sequence - it uses variables
     * computed by other routines (ugly).
     */
    int i;
    GLuint n;

    n = glGenLists(1);
    glNewList(n, GL_COMPILE); /* ?? */

    glPushAttrib(GL_POLYGON_BIT);

#if QUADRICS_NORMALIZE
    glEnable(GL_NORMALIZE);
#endif

    if (n & 0x1)
	glFrontFace(GL_CCW);
    else
	glFrontFace(GL_CW);
    
    glBegin(GL_TRIANGLES);
    for (i = 0; i < sphere_facet_count * 3; i++)
    {
	glNormal3fv((GLfloat *)(sphere_norms + i));
	glVertex3fv((GLfloat *)(sphere_verts + i));
    }
    glEnd();

#if QUADRICS_NORMALIZE
    glDisable(GL_NORMALIZE);
#endif
    
    glPopAttrib();

    glEndList();
    return(n);
}


/* PRIVATE */
static void
build_return_arrays (sphereobj *obj,
		     RMvertex3D **vlist,
		     RMvertex3D **nlist,
		     int *ntriangles_return)
{
    int         i, ntriangles;
    RMvertex3D *v;
    RMvertex3D *n;
    triangle   *t;

    ntriangles = obj->ntriangles;

    v = rmVertex3DNew(ntriangles * 3);
    n = rmVertex3DNew(ntriangles * 3);

    for (i = 0; i < ntriangles; i++)
    {
	t = &(obj->list[i]);
	VCOPY(t->pt, v +(i * 3));
	VCOPY(&(t->pt[2]), v + (i * 3) + 2);
	VCOPY(&(t->pt[1]), v + (i * 3) + 1);
    }
    
    for (i = 0; i < (ntriangles * 3); i++)
    {
	VCOPY(v + i, n + i);
	rmVertex3DNormalize(n + i);
    }

    *vlist = v;
    *nlist = n;
    *ntriangles_return = ntriangles;
}


/* PRIVATE */
int
private_rmGetTrianglesPerSphere (void)
{
    return(triangles_per_sphere);
}


/* PRIVATE */
int
private_initSphereTables (int n)
{
    /*
     * the sphere tesselator builds internal triangle and normal tables
     * which will be dispatched whenever a sphere is "drawn."  the
     * first parameter is an integer which indicates the granularity
     * of subdivision.  the tesselator will create equal-area triangles.
     *
     *  value    number of triangles
     *    1             8
     *    2             32
     *    3             128
     *    4             512
     *    n             2*4^n
     */
    private_rmSphereTesselator(n, &sphere_verts, &sphere_norms, &sphere_facet_count);
    triangles_per_sphere = 2 * (int)(pow(4.0, (double)n));

    return(0);
}


/* PRIVATE */
int
private_freeSphereTables ()
{
    if (sphere_verts != NULL)
	rmVertex3DDelete(sphere_verts);
    if (sphere_norms != NULL)
	rmVertex3DDelete(sphere_norms);

    return(0);
}


/* PRIVATE */
void
private_rmSphereTesselator (int subdivisions,
			    RMvertex3D **vlist,
			    RMvertex3D **nlist,
			    int *ntriangles_return)
{
    int        ccwflag = 1;	/* Reverse vertex order if true */
    int        i;
    int        level;       	/* Current subdivision level */
    int	       maxlevel = 1;    /* Maximum subdivision level */
    sphereobj *old, *new1;
 
    tri_ptr = 0;
    maxlevel = subdivisions;
   
    if (ccwflag)
    {
        /* Reverse order of points in each triangle */
        for (i = 0; i < oct.ntriangles; i++)
	{
            point tmp = oct.list[i].pt[0];
            oct.list[i].pt[0] = oct.list[i].pt[2];
            oct.list[i].pt[2] = tmp;
        }
    }
    old = &oct;

    /* Subdivide each starting triangle (maxlevel - 1) times */
    for (level = 1; level < maxlevel; level++)
    {
	/* Allocate a new object */
	new1 = (sphereobj *)malloc(sizeof(sphereobj));
	if (new1 == NULL)
	{
	    rmError("malloc error in sphere tesselator");
	}
	new1->ntriangles = old->ntriangles * 4;

	/* Allocate 4* the number of points in the current approximation */
	new1->list  = (triangle *)malloc(new1->ntriangles * sizeof(triangle));
	if (new1->list == NULL)
	{
	    new1 = (sphereobj *)malloc(sizeof(sphereobj));
	}

	/* Subdivide each triangle in the old approximation and normalize
	 *  the new points thus generated to lie on the surface of the unit
	 *  sphere.
	 * Each input triangle with vertices labelled [0,1,2] as shown
	 *  below will be turned into four new triangles:
	 *
	 *			Make new points
	 *			    a = (0+2)/2
	 *			    b = (0+1)/2
	 *			    c = (1+2)/2
	 *	  1
	 *	 /\		Normalize a, b, c
	 *	/  \
	 *    b/____\ c		Construct new triangles
	 *    /\    /\		    [0,b,a]
	 *   /	\  /  \		    [b,1,c]
	 *  /____\/____\	    [a,b,c]
	 * 0	  a	2	    [a,c,2]
	 */
	for (i = 0; i < old->ntriangles; i++)
	{
	    triangle *oldt = &old->list[i];
	    triangle *newt = &new1->list[i*4];
	    point     a, b, c;

	    a = *normalize(midpoint(&oldt->pt[0], &oldt->pt[2]));
	    b = *normalize(midpoint(&oldt->pt[0], &oldt->pt[1]));
	    c = *normalize(midpoint(&oldt->pt[1], &oldt->pt[2]));

	    newt->pt[0] = oldt->pt[0];
	    newt->pt[1] = b;
	    newt->pt[2] = a;
	    newt++;

	    newt->pt[0] = b;
	    newt->pt[1] = oldt->pt[1];
	    newt->pt[2] = c;
	    newt++;

	    newt->pt[0] = a;
	    newt->pt[1] = b;
	    newt->pt[2] = c;
	    newt++;

	    newt->pt[0] = a;
	    newt->pt[1] = c;
	    newt->pt[2] = oldt->pt[2];
	}

	if (level > 1)
	{
	    free(old->list);
	    free(old);
	}

	/* Continue subdividing new triangles */
	old = new1;
    }

    build_return_arrays(old, vlist, nlist, ntriangles_return);

    if (old != &oct)
    {
	free(old->list);
	free(old);
    }
}


/* Normalize a point p */
static point *
normalize (point *p)
{
    double       mag;
    static point r;

    r = *p;
    mag = r.x * r.x + r.y * r.y + r.z * r.z;
    if (mag != 0.0) 
    {
	mag = 1.0 / sqrt(mag);
	r.x *= mag;
	r.y *= mag;
	r.z *= mag;
    }
    return &r;
}


/* Return the midpoint on the line between two points */
static point *
midpoint (point *a,
	  point *b)
{
    static point r;

    r.x = (a->x + b->x) * 0.5;
    r.y = (a->y + b->y) * 0.5;
    r.z = (a->z + b->z) * 0.5;

    return &r;
}


/* PRIVATE */
void
private_yaxis_to_dir (RMmatrix *m,
		      RMvertex3D *dir)
{
    /*
     * compute a rotation matrix that will transform from the y-axis
     * (0,1,0) to the vector in "dir". the cross product of (0,1,0)
     * and "dir" (let's call it a,b,c) is (-c,0,a). the angle of
     * rotation is acos(b). the transpose of resulting matrix is:
     *
      |  A^2 * B + C^2                 A*B*C - A*C   |
      |  -------------        A        -----------   |
      |    A^2 + C^2                    A^2 + C^2    |
      |                                              |
      |       -A              B            -C        |
      |                                              |
      |   A*B*C - A*C                 B * C^2 + A^2  |
      |   -----------         C       -------------  |
      |    A^2 + C^2                    A^2 + C^2    |
      |                                              |
     */
    
    double a2c2;
    double m11, m12, m13, m31, m32, m33;
    
    rmMatrixIdentity(m);

    a2c2 = (dir->x * dir->x) + (dir->z * dir->z);

    if (a2c2 == 0.0)  /* a degenerate condition; make an identity matrix */
    {
        m11 = 1.0;
        m12 = 0.0;
        m13 = 0.0;
 
        m31 = 0.0;
        m32 = 0.0;
        m33 = 1.0;
 
    }
    else
    {
        a2c2 = 1.0 / a2c2;
 
	m11 = (dir->x * dir->x * dir->y) + (dir->z * dir->z);
        m11 *= a2c2;
 
	m12 = -1.0 * dir->x;
 
	m13 = (dir->x * dir->y * dir->z) - (dir->x * dir->z);
	
        m13 *= a2c2;
 
        m31 = m13;
 
	m32 = -1.0 * dir->z;
 
	m33 = (dir->y * dir->z * dir->z) + (dir->x * dir->x);
        m33 *= a2c2;
    }
    m->m[0][0] = m11;
    m->m[1][0] = m12;
    m->m[2][0] = m13;
    m->m[0][2] = m31;
    m->m[1][2] = m32;
    m->m[2][2] = m33;
    m->m[0][1] = dir->x;
    m->m[1][1] = dir->y;
    m->m[2][1] = dir->z;
}
/* EOF */
