#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <GL/glut.h>
#include "../jlib/jlib.h"
#include <SWIFT.h>
#include <iostream.h>

// CONSTANTS

const double DRAWOUTPUT = true;		// Draw the output.
const int NSTEPS = 0;			// Number of time steps to run.  0 means indefinately.
const bool SAMESIZE = true;			// Whether or not all balls have the same size.
const bool PLANAR = false;			// Put everything in the xz plane, and run 2-d (NOT OPTIMIZED FOR THIS!)
const int NBALLS = 50;
const double BOXSIZE = 10.0;
const double BALLSIZE = 1.0;		// Max ball radius when variable, or fixed ball radius when samesize is true.
const double MAXV = 0.2;
const double MAXANGV = 10.0;
//const double PI = 3.14159265;
const int SEGS_PER_Q = 2;

// View globals

bool usebins = true;			// Use binning to reduce tests.
bool do_collisions = true;		// Do sphere-to-sphere collision detection.
int mode;
double beginx, beginy;
double dis = BOXSIZE*2.5, azim = 0.0, elev = 0.0;
double ddis = 0.0, dazim = 0.0, delev = 0.0;
long timestep = 0;				// How many time steps we've completed.

// Object descriptions.

double size[NBALLS];
double X[NBALLS];
double Y[NBALLS];
double Z[NBALLS];
double dX[NBALLS];
double dY[NBALLS];
double dZ[NBALLS];
double the[NBALLS];
double phi[NBALLS];
double dthe[NBALLS];
double dphi[NBALLS];

int ballList;
int boxid=-1;		// Id of containing box.
int firstid;	// First sphere.
SWIFT_Scene *scene;


bool **collide, **collidenew;	// Whether or not a pair of spheres are currently penetrating.


void Swap(double &a, double &b)
{
	double t = a;
	a = b;
	b = t;
}

void
init_viewer_window()
{
  GLfloat Ambient[] = { 0.2f, 0.2f, 0.2f, 1.0f };  
  GLfloat Diffuse[] = { 0.8f, 0.8f, 0.8f, 1.0f };  
  GLfloat Specular[] = { 0.1f, 0.1f, 0.1f, 1.0f };   
  GLfloat SpecularExp[] = { 50 };              
  GLfloat Emission[] = { 0.1f, 0.1f, 0.1f, 1.0f };   

  glMaterialfv(GL_FRONT, GL_AMBIENT, Ambient);
  glMaterialfv(GL_FRONT, GL_DIFFUSE, Diffuse);
  glMaterialfv(GL_FRONT, GL_SPECULAR, Specular);
  glMaterialfv(GL_FRONT, GL_SHININESS, SpecularExp);
  glMaterialfv(GL_FRONT, GL_EMISSION, Emission);

  glMaterialfv(GL_BACK, GL_AMBIENT, Ambient);
  glMaterialfv(GL_BACK, GL_DIFFUSE, Diffuse);
  glMaterialfv(GL_BACK, GL_SPECULAR, Specular);
  glMaterialfv(GL_BACK, GL_SHININESS, SpecularExp);
  glMaterialfv(GL_BACK, GL_EMISSION, Emission);

  glColorMaterial(GL_FRONT_AND_BACK, GL_DIFFUSE);

  glEnable(GL_COLOR_MATERIAL);

  GLfloat light_position[] = { 1.0, 1.0, 1.0, 0.0 };
  glLightfv(GL_LIGHT0, GL_POSITION, light_position);
  glEnable(GL_LIGHT0);
  glEnable(GL_LIGHTING);
  glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, GL_TRUE);

  glDepthFunc(GL_LEQUAL);
  glEnable(GL_DEPTH_TEST);
  glEnable(GL_CULL_FACE);
  glCullFace(GL_BACK);

  glShadeModel(GL_FLAT);
  glClearColor(0.0, 0.0, 0.0, 0.0);

  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  glFrustum(-0.004,0.004,-0.004,0.004,.01,100.0);

  glMatrixMode(GL_MODELVIEW);
}

void
cb_mouse(int _b, int _s, int _x, int _y)
{
  if (_s == GLUT_UP)
  {
    dis += ddis;
    if (dis < .1) dis = .1;
    azim += dazim;
    elev += delev;
    ddis = 0.0;
    dazim = 0.0;
    delev = 0.0;
    return;
  }

  if (_b == GLUT_RIGHT_BUTTON)
  {
    mode = 0;
    beginy = _y;
    return;
  }
  else
  {
    mode = 1;
    beginx = _x;
    beginy = _y;
  }
}

void
cb_motion(int _x, int _y)
{
  if (mode == 0)
  {
    ddis = dis * (double)(_y - beginy)/200.0;
  }
  else
  {
    dazim = (_x - beginx)/5;
    delev = (_y - beginy)/5;      
  }
  
  glutPostRedisplay();
}

//
// BINNING FUNCTIONS
//

struct BinId
{	int i,j,k;
	BinId(int a, int b, int c) : i(a), j(b), k(c) {}
	BinId() {}
	bool operator==(const BinId& b) const { return i==b.i && j==b.j && k==b.k; }
	bool Neighbor(const BinId& b) const
	{
		return !( i-b.i > 1 || i-b.i < -1 ||
				  j-b.j > 1 || j-b.j < -1 ||
				  k-b.k > 1 || k-b.k < -1 );
	}
};

struct Bin
{	int ball;
	Bin* next;
};

int nbins;
double binwidth;
Bin* ***bins;
BinId *inBin;

void CreateBins()
{
	// Calculate bin size.
	nbins = 2*BOXSIZE / (2*BALLSIZE);	// total box length over ball diameter, rounding down.
	binwidth = 2*BOXSIZE/double(nbins);
	cout << nbins << " bins of length " << binwidth << " in each dimension." << endl;

	// Create bins object.
	bins = new Bin* **[nbins];
	for (int i = 0; i < nbins; ++i)
	{
		bins[i] = new Bin* *[nbins];
		for (int j = 0; j < nbins; ++j)
		{
			bins[i][j] = new Bin*[nbins];
			for (int k = 0; k < nbins; ++k)
				bins[i][j][k] = NULL;
		}
	}

	// Create array that shows which bin each ball is in.
	// The reverse lookup list.
	inBin = new BinId[NBALLS];
}

BinId WhichBin(double x, double y, double z)
{
	return BinId( (x+BOXSIZE)/binwidth,
	              (y+BOXSIZE)/binwidth,
				  (z+BOXSIZE)/binwidth );
}

void AddToBin(int ballId, BinId oldbin = BinId(-10,-10,-10))
{
	BinId b = WhichBin(X[ballId], Y[ballId], Z[ballId]);
	
	// Update reverse lookup list.
	inBin[ballId] = b;

	// Add to that bin and surrounding bins.
	for (int i = -1; i <= +1; ++i)
		for (int j = -1; j <= +1; ++j)
			for (int k = -1; k <= +1; ++k)
			{
				// Test for out of bounds.
				if (b.i+i < 0 || b.i+i >= nbins)
					continue;
				if (b.j+j < 0 || b.j+j >= nbins)
					continue;
				if (b.k+k < 0 || b.k+k >= nbins)
					continue;

				// Test if it was already there from the previous location.
				if (oldbin.Neighbor(BinId(b.i+i,b.j+j,b.k+k)))
					continue;

				// Add to this bin.
				Bin*& head = bins[b.i+i][b.j+j][b.k+k];
				Bin* temp = new Bin;
				temp->ball = ballId;
				temp->next = head;
				head = temp;

				// Active pair testing for this ball and all others in this list.
				for (Bin* cur = head->next; cur != NULL; cur = cur->next)
				{
					scene->Activate(ballId+firstid, cur->ball+firstid);
					//cout << "ON " << ballId << ", " << cur->ball;
					//cout << "  in (" << b.i+i << "," << b.j+j << "," << b.k+k << ")" << endl;
				}
			}
}

void UpdateBin(int ballId)
{
	BinId newb = WhichBin(X[ballId], Y[ballId], Z[ballId]);
	BinId oldb = inBin[ballId];
	
	// See if the ball has moved bins.
	if (newb == oldb)
		return;		// Same as before; no update.

	//cout << ballId << "->(" << newb.i << "," << newb.j << "," << newb.k << ")" << endl;

	// Otherwise, we need to unbin the ball from its old bin and neighbors.
	for (int i = -1; i <= +1; ++i)
		for (int j = -1; j <= +1; ++j)
			for (int k = -1; k <= +1; ++k)
			{
				BinId b(oldb.i+i, oldb.j+j, oldb.k+k);

				// Test for out of bounds.
				if (b.i < 0 || b.i >= nbins)
					continue;
				if (b.j < 0 || b.j >= nbins)
					continue;
				if (b.k < 0 || b.k >= nbins)
					continue;

				if (b.Neighbor(newb))
					continue;

				// If it's in bounds, then remove from that bin.

				Bin*& head = bins[b.i][b.j][b.k];
				if (head == NULL)
					continue;
				if (head->ball == ballId)
				{
					Bin* temp = head;
					head = head->next;
					delete temp;
					if (head == NULL)
						continue;
				}
				Bin* cur = head;
				scene->Deactivate(ballId+firstid, cur->ball+firstid);
				while (cur->next != NULL)
				{
					if (cur->next->ball == ballId)
					{
						Bin* temp = cur->next;
						cur->next = cur->next->next;
						delete temp;
						break;
					}
					cur = cur->next;
					scene->Deactivate(ballId+firstid, cur->ball+firstid);
				}
			}

	// Add to its new bins.
	AddToBin(ballId, oldb);
}

void StartBinning()
{
	// Clear the pair test matrix.
	scene->Deactivate();
	for (int i = 0; i < NBALLS; ++i)
		AddToBin(i);
}

void StopBinning()
{
	// Empty linked lists from all bins..
	for (int i = 0; i < nbins; ++i)
		for (int j = 0; j < nbins; ++j)
			for (int k = 0; k < nbins; ++k)
			{
				Bin*& head = bins[i][j][k];
				Bin* temp;
				while (head != NULL)
				{
					temp = head;
					head = head->next;
					delete temp;
				}
			}

	scene->Activate();
}


//
//
//

void DoWallCollide(int i)
{
	// Handle collision detection/resolution of ball i with the walls.

	if (X[i]+size[i] > BOXSIZE)
		dX[i] = -fabs(dX[i]);
	if (X[i]-size[i] < -BOXSIZE)
		dX[i] = +fabs(dX[i]);
	if (Y[i]+size[i] > BOXSIZE)
		dY[i] = -fabs(dY[i]);
	if (Y[i]-size[i] < -BOXSIZE)
		dY[i] = +fabs(dY[i]);
	if (Z[i]+size[i] > BOXSIZE)
		dZ[i] = -fabs(dZ[i]);
	if (Z[i]-size[i] < -BOXSIZE)
		dZ[i] = +fabs(dZ[i]);
}

//
// UPDATE FUNCTION
//

void cb_idle()
{
	// Do one step of the movement.

	SWIFT_Real trans[3];
	SWIFT_Real rotation[9];

	for (int i = 0; i < NBALLS; ++i)
	{
		X[i] += dX[i];
		Y[i] += dY[i];
		Z[i] += dZ[i];
		the[i] += dthe[i];
		phi[i] += dphi[i];

		trans[0] = X[i];
		trans[1] = Y[i];
		trans[2] = Z[i];
		rotation[0] = cos(phi[i]) - 0;
		rotation[1] = sin(phi[i]) + 0;
		rotation[2] = 0;
		rotation[3] = 0 - cos(the[i])*sin(phi[i]);
		rotation[4] = 0 + cos(the[i])*cos(phi[i]);
		rotation[5] = sin(the[i]);
		rotation[6] = sin(the[i])*sin(phi[i]);
		rotation[7] = -sin(the[i])*cos(phi[i]);
		rotation[8] = cos(the[i]);
		scene->Set_Object_Transformation(firstid+i, rotation, trans);

		if (usebins)
			UpdateBin(i);
	}

	//
	// Collision detection between spheres.
	//

	// Clear the new collision matrix.
	for (i = 0; i < NBALLS; ++i)
		for (int j = 0; j < NBALLS; ++j)
			collidenew[i][j] = false;

	int n;
	int *oids;
	if (do_collisions)	// If collision detection is on.
		scene->Query_Intersection(false, n, &oids);
	else
		n = 0;
	for (i = 0; i < n; ++i)
	{
		if (oids[2*i] == boxid)
			DoWallCollide(oids[2*i+1] - firstid);
		else if (oids[2*i+1] == boxid)
			DoWallCollide(oids[2*i] - firstid);
		else
		{
			int p1 = oids[2*i] - firstid;
			int p2 = oids[2*i+1] - firstid;

			collidenew[p1][p2] = collidenew[p2][p1] = true;

			if (collide[p1][p2])
			{
				//cout << "Re-collision." << endl;
				continue;	// Saw this collision last time step, don't do the collision again.
			}

			Swap(dX[p1],dX[p2]);
			Swap(dY[p1],dY[p2]);
			Swap(dZ[p1],dZ[p2]);
			Swap(dthe[p1],dthe[p2]);
			Swap(dphi[p1],dphi[p2]);
			//cout << "Bump " << p1 << ", " << p2 << endl;
		}
	}

	for (i = 0; i < NBALLS; ++i)  DoWallCollide(i);

	// Swap the collide matricies.
	bool **temp = collide;
	collide = collidenew;
	collidenew = temp;
	
	glutPostRedisplay();

	timestep++;
	if (timestep == NSTEPS)
	{
		double elapsed = ReadTimer();
		cout << elapsed << " sec per " << NSTEPS << " steps.  " << elapsed/NSTEPS << " sec/step."  << endl;
		cin.get();
		exit(0);
	}
}

void cb_keyboard(unsigned char key, int Mx, int My) 
{
	switch(key) 
	{
	case 'q': 
		exit(0);
	case 'b':
		if (!SAMESIZE)
			break;
		usebins = !usebins;
		cout << "Binning is " << (usebins ? "ON" : "OFF") << endl;
		if (usebins)
			StartBinning();
		else
			StopBinning();
		break;
	case 'c':
		do_collisions = !do_collisions;
		cout << "Collision detection is " << (do_collisions ? "ON" : "OFF") << endl;
		break;
	case ' ':	// Pause
		cin.get();
		break;
	default:
		cb_idle();
		break;
	}

	glutPostRedisplay();
}

void
cb_display()
{
	if (!DRAWOUTPUT)
		return;

	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

	glLoadIdentity(); 
	glTranslatef(0.0, 0.0, -(dis+ddis));
	glRotated(elev+delev, 1.0, 0.0, 0.0);
	glRotated(azim+dazim, 0.0, 1.0, 0.0);
	glRotated(90.0,-1.0,0.0,0.0);
	glColor3f(1,1,1);

	for (int i = 0; i < NBALLS; ++i)
	{
		glPushMatrix();
		glTranslated(X[i],Y[i],Z[i]);
		glScalef(size[i], size[i], size[i]);
		glRotated(-the[i], 1.0, 0.0, 0.0);
		glRotated(phi[i], 0.0, 0.0, 1.0);
		glCallList(ballList);
		glPopMatrix();
	}

	glFlush();
	glutSwapBuffers();
}

inline SWIFT_Real Trimzero(double f)
{
	return (fabs(f) < 1e-14) ? 0.0 : f;
}

#define XCOORD(x,y)  Trimzero(cos(x)*cos(y))
#define YCOORD(x,y)  Trimzero(sin(y))
#define ZCOORD(x,y)  Trimzero(sin(x)*cos(y))

void InitScene()
{
	scene = new SWIFT_Scene();

	SWIFT_Real *boxverts = new SWIFT_Real[8*3];
	int *boxfaces = new int[6*4];
	int *boxvalence = new int[6];

	for (int i = 0; i < 6; ++i)
		boxvalence[i] = 4;

	boxverts[3*0+0] = -1;	// -x -y -z
	boxverts[3*0+1] = -1;
	boxverts[3*0+2] = -1;
	boxverts[3*1+0] = -1;	// -x -y +z
	boxverts[3*1+1] = -1;
	boxverts[3*1+2] = +1;
	boxverts[3*2+0] = -1;	// -x +y -z
	boxverts[3*2+1] = +1;
	boxverts[3*2+2] = -1;
	boxverts[3*3+0] = -1;	// -x +y +z
	boxverts[3*3+1] = +1;
	boxverts[3*3+2] = +1;
	boxverts[3*4+0] = +1;	// +x -y -z
	boxverts[3*4+1] = -1;
	boxverts[3*4+2] = -1;
	boxverts[3*5+0] = +1;	// +x -y +z
	boxverts[3*5+1] = -1;
	boxverts[3*5+2] = +1;
	boxverts[3*6+0] = +1;	// +x +y -z
	boxverts[3*6+1] = +1;
	boxverts[3*6+2] = -1;
	boxverts[3*7+0] = +1;	// +x +y +z
	boxverts[3*7+1] = +1;
	boxverts[3*7+2] = +1;

	boxfaces[0*4+0] = 0;	// -x
	boxfaces[0*4+1] = 1;
	boxfaces[0*4+2] = 3;
	boxfaces[0*4+3] = 2;
	boxfaces[1*4+0] = 4;	// +x
	boxfaces[1*4+1] = 5;
	boxfaces[1*4+2] = 7;
	boxfaces[1*4+3] = 6;
	boxfaces[2*4+0] = 0;	// -y
	boxfaces[2*4+1] = 1;
	boxfaces[2*4+2] = 5;
	boxfaces[2*4+3] = 4;
	boxfaces[3*4+0] = 2;	// +y
	boxfaces[3*4+1] = 3;
	boxfaces[3*4+2] = 7;
	boxfaces[3*4+3] = 6;
	boxfaces[4*4+0] = 0;	// -z
	boxfaces[4*4+1] = 2;
	boxfaces[4*4+2] = 6;
	boxfaces[4*4+3] = 4;
	boxfaces[5*4+0] = 1;	// +z
	boxfaces[5*4+1] = 3;
	boxfaces[5*4+2] = 7;
	boxfaces[5*4+3] = 5;

	//scene->Add_Convex_Object(boxverts, boxfaces, 8, 6, boxid, true,
	//	       DEFAULT_ORIENTATION, DEFAULT_TRANSLATION, BOXSIZE,
	//		   DEFAULT_BOX_SETTING, DEFAULT_BOX_ENLARGE_REL, DEFAULT_BOX_ENLARGE_ABS,
	//		   boxvalence);

	int NFACES = 2*(2*(SEGS_PER_Q-1)*(4*SEGS_PER_Q) + (4*SEGS_PER_Q));
	int *faces = new int[3*NFACES];
	SWIFT_Real *verts = new SWIFT_Real[3*3*NFACES];
	int vi = 0;
	int fi = 0;

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

	glBegin(GL_QUADS);
	for (int x = 0; x < 4*SEGS_PER_Q; x++)
		for (int y = -SEGS_PER_Q+1; y < SEGS_PER_Q-1; y++)
	{
		double az0 = x/double(SEGS_PER_Q) * PI/2;
		double az1 = (x+1)/double(SEGS_PER_Q) * PI/2;
		double el0 = y/double(SEGS_PER_Q) * PI/2;
		double el1 = (y+1)/double(SEGS_PER_Q) * PI/2;
		glNormal3f(XCOORD(az0,el0), YCOORD(az0,el0), ZCOORD(az0,el0));
		glVertex3f(XCOORD(az0,el0), YCOORD(az0,el0), ZCOORD(az0,el0));
		glVertex3f(XCOORD(az0,el1), YCOORD(az0,el1), ZCOORD(az0,el1));
		glVertex3f(XCOORD(az1,el1), YCOORD(az1,el1), ZCOORD(az1,el1));
		glVertex3f(XCOORD(az1,el0), YCOORD(az1,el0), ZCOORD(az1,el0));

		faces[fi++] = vi;
		verts[vi*3+0] = XCOORD(az0,el0);
		verts[vi*3+1] = YCOORD(az0,el0);
		verts[vi*3+2] = ZCOORD(az0,el0);
		vi++;
		faces[fi++] = vi;
		verts[vi*3+0] = XCOORD(az0,el1);
		verts[vi*3+1] = YCOORD(az0,el1);
		verts[vi*3+2] = ZCOORD(az0,el1);
		vi++;
		faces[fi++] = vi;
		verts[vi*3+0] = XCOORD(az1,el0);
		verts[vi*3+1] = YCOORD(az1,el0);
		verts[vi*3+2] = ZCOORD(az1,el0);
		vi++;

		faces[fi++] = vi;
		verts[vi*3+0] = XCOORD(az1,el0);
		verts[vi*3+1] = YCOORD(az1,el0);
		verts[vi*3+2] = ZCOORD(az1,el0);
		vi++;
		faces[fi++] = vi;
		verts[vi*3+0] = XCOORD(az0,el1);
		verts[vi*3+1] = YCOORD(az0,el1);
		verts[vi*3+2] = ZCOORD(az0,el1);
		vi++;
		faces[fi++] = vi;
		verts[vi*3+0] = XCOORD(az1,el1);
		verts[vi*3+1] = YCOORD(az1,el1);
		verts[vi*3+2] = ZCOORD(az1,el1);
		vi++;
	}
	glEnd();

	glBegin(GL_TRIANGLES);
	for (x = 0; x < 4*SEGS_PER_Q; x++)
	{
		double az0 = x/double(SEGS_PER_Q) * PI/2;
		double az1 = (x+1)/double(SEGS_PER_Q) * PI/2;

		double el0 = (SEGS_PER_Q-1)/double(SEGS_PER_Q) * PI/2;
		double el1 = SEGS_PER_Q/double(SEGS_PER_Q) * PI/2;
		glNormal3f(XCOORD(0.5*(az0+az1),0.5*(el0+el1)), YCOORD(0.5*(az0+az1),0.5*(el0+el1)), ZCOORD(0.5*(az0+az1),0.5*(el0+el1)));
		glVertex3f(XCOORD(az0,el0), YCOORD(az0,el0), ZCOORD(az0,el0));
		glVertex3f(XCOORD(az0,el1), YCOORD(az0,el1), ZCOORD(az0,el1));
		glVertex3f(XCOORD(az1,el0), YCOORD(az1,el0), ZCOORD(az1,el0));

		faces[fi++] = vi;
		verts[vi*3+0] = XCOORD(az0,el0);
		verts[vi*3+1] = YCOORD(az0,el0);
		verts[vi*3+2] = ZCOORD(az0,el0);
		vi++;
		faces[fi++] = vi;
		verts[vi*3+0] = XCOORD(az0,el1);
		verts[vi*3+1] = YCOORD(az0,el1);
		verts[vi*3+2] = ZCOORD(az0,el1);
		vi++;
		faces[fi++] = vi;
		verts[vi*3+0] = XCOORD(az1,el0);
		verts[vi*3+1] = YCOORD(az1,el0);
		verts[vi*3+2] = ZCOORD(az1,el0);
		vi++;

		el0 = -SEGS_PER_Q/double(SEGS_PER_Q) * PI/2;
		el1 = (-SEGS_PER_Q+1)/double(SEGS_PER_Q) * PI/2;
		glNormal3f(XCOORD(0.5*(az0+az1),0.5*(el0+el1)), YCOORD(0.5*(az0+az1),0.5*(el0+el1)), ZCOORD(0.5*(az0+az1),0.5*(el0+el1)));
		glVertex3f(XCOORD(az0,el0), YCOORD(az0,el0), ZCOORD(az0,el0));
		glVertex3f(XCOORD(az0,el1), YCOORD(az0,el1), ZCOORD(az0,el1));
		glVertex3f(XCOORD(az1,el1), YCOORD(az1,el1), ZCOORD(az1,el1));

		faces[fi++] = vi;
		verts[vi*3+0] = XCOORD(az0,el0);
		verts[vi*3+1] = YCOORD(az0,el0);
		verts[vi*3+2] = ZCOORD(az0,el0);
		vi++;
		faces[fi++] = vi;
		verts[vi*3+0] = XCOORD(az0,el1);
		verts[vi*3+1] = YCOORD(az0,el1);
		verts[vi*3+2] = ZCOORD(az0,el1);
		vi++;
		faces[fi++] = vi;
		verts[vi*3+0] = XCOORD(az1,el1);
		verts[vi*3+1] = YCOORD(az1,el1);
		verts[vi*3+2] = ZCOORD(az1,el1);
		vi++;
	}
	glEnd();
	glEndList();

	// Add unit-sized, origin-centered sphereoid to SWIFT scene.
	// We'll copy it many times (with different sizes and positions) and delete it later.
	//int unitid;
	//scene->Add_Convex_Object(verts, faces, NFACES*3, NFACES, unitid, false);

	cout << NBALLS << " balls." << endl;
	cout << NFACES << " triangles per ball." << endl;

	// Choose initial size, location, and velocity for each sphere.
	for (i = 0; i < NBALLS; ++i)
	{
		if (SAMESIZE)
			size[i] = BALLSIZE;
		else
			size[i] = Random(BALLSIZE/100.0, BALLSIZE);
		X[i] = Random(-BOXSIZE+size[i], +BOXSIZE-size[i]);
		Y[i] = Random(-BOXSIZE+size[i], +BOXSIZE-size[i]);
		if (PLANAR)
			Y[i] = 0;
		Z[i] = Random(-BOXSIZE+size[i], +BOXSIZE-size[i]);
		dX[i] = Random(-MAXV, MAXV);
		dY[i] = Random(-MAXV, MAXV);
		if (PLANAR)
			dY[i] = 0;
		dZ[i] = Random(-MAXV, MAXV);
		the[i] = Random(0, 2*PI);
		phi[i] = Random(-PI, +PI);
		dthe[i] = Random(-MAXANGV, +MAXANGV);
		dphi[i] = Random(-MAXANGV, +MAXANGV);


		int foo;
		//scene->Add_Convex_Object(NULL, NULL, 0, 0,   // Unused geometry input.
		//           foo, false, DEFAULT_ORIENTATION, DEFAULT_TRANSLATION, size[i],
		//		   DEFAULT_BOX_SETTING, DEFAULT_BOX_ENLARGE_REL, DEFAULT_BOX_ENLARGE_ABS, NULL,
		//		   unitid);
		scene->Add_Convex_Object(verts, faces, NFACES*3, NFACES,
		           foo, false, DEFAULT_ORIENTATION, DEFAULT_TRANSLATION, size[i]);
		//scene->Add_Convex_Object("sphere1.tri", foo, false, DEFAULT_ORIENTATION, DEFAULT_TRANSLATION, size[i]);
		if (i == 0)
			firstid = foo;	// Store the id of the first sphere.

	}

	// Delete the unit sphere.
	//scene->Delete_Object(unitid);

	delete[] faces;
	delete[] verts;
	delete[] boxverts;
	delete[] boxfaces;
	delete[] boxvalence;

	// Do a dummy query outside the timing loop to get warmed up.
	int n;
	int *oids;
	scene->Query_Intersection(false, n, &oids);
}

void InitData()
{
	collide = new bool*[NBALLS];
	collidenew = new bool*[NBALLS];
	for (int i = 0; i < NBALLS; ++i)
	{
		collide[i] = new bool[NBALLS];
		for (int j = 0; j < NBALLS; ++j)
			collide[i][j] = false;
		collidenew[i] = new bool[NBALLS];
	}

	CreateBins();

	if (usebins)
		StartBinning();
}

void main(int argc, char **argv)
{
  // init glut

  glutInit(&argc, argv);
  glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH | GLUT_MULTISAMPLE);

  // create the window

  glutInitWindowSize(500, 500);
  glutInitWindowPosition(500, 100);
  glutCreateWindow("259 Hw2");

  // set OpenGL graphics state -- material props, perspective, etc.

  init_viewer_window();

  // set the callbacks

  glutDisplayFunc(cb_display);
  glutMouseFunc(cb_mouse);
  glutMotionFunc(cb_motion);  
  glutKeyboardFunc(cb_keyboard);
  glutIdleFunc(cb_idle);

  // Create models
  SeedRandom(1);
  InitScene();

  InitData();

  if (SAMESIZE)
	  cout << "Binning is " << (usebins?"ON":"OFF") << endl;

  StartTimer();

  glutMainLoop();
}