#include <GL/glut.h>
#include <fstream.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <iostream>

// some types for later use
typedef unsigned char byte;
typedef byte rgb[3];			// these are 3 byte pixels in RGB order

// a few colors to use
static rgb red = { 255, 0, 0 };
static rgb blue = { 0, 0, 255 };
static rgb black = {0,0,0};
int num_vert;


/***Now is the global scene information variables.***/

double cop[3];		//The center of projection.
double dov[3];		//The direction of view, lookp - eyep.
double up_dir[3];  //the up-direction vector, picked up during file input.
double foc;			//The focal length.
double background[3];//The background color.
int xres = 512;
int yres = 512;

struct sphere {
	double rad, c[3];			//Spatial properties.
	double amb[3], dif[3], spec[3], refl, transp, index; //Material properties.
	int specpow;
	sphere *next;

	sphere(double r, double x[]) {
		rad = r; c[0] = x[0]; c[1] = x[1]; c[2] = x[2]; next = NULL;
	}
} *sphere_head;

struct plane {
	double p[3], n[3];  //A point and normal of the plane.
	double amb[3], dif[3], spec[3], refl, transp, index; //Material properties.
	int specpow;
	plane *next;

	plane(double x[], double norm[]) {
		p[0] = x[0]; p[1] = x[1]; p[2] = x[2];
		n[0] = norm[0]; n[1] = norm[1]; n[2] = norm[2];
		next = NULL;
	}
} *plane_head;

struct light {
	double p[3];
	double amb[3], dif[3], spec[3];
	light *next;

	light(double x[]) {
		p[0] = x[0]; p[1] = x[1]; p[2] = x[2]; next = NULL;
		//color[0] = c[0]; color[1] = c[1]; color[2] = c[2];
	}
} *light_head;

//Triangles will be filled out later.
struct triangle {
	double p1[3], p2[3], p3[3];
};

struct vertex {
	rgb color;
	double x, y, z;

	void fill_vertex(double gx, double gy, double gz, const rgb p) {
		x = gx; y = gy; z = gz;
		color[0] = p[0]; color[1] = p[1]; color[2] = p[2];
	}
};

// an object to represent the framebuffer where I draw
struct FrameBuffer {
	int width, height;
	rgb *pix;
	double *pz; //pixel z value.

	FrameBuffer(int w, int h) { width=w; height=h; pix = new rgb[width*height]; 
								pz = new double[width*height];}
	void Clear(const rgb p);
	void Test();

} *fb; // global framebuffer object for use in callbacks

struct Projection {
	double mat[3][3];

	Projection(double dov[], double foc);
	//Initializes the world-from-image projection matrix.

	void Convert(double v[]);
	//Applies mat to v, where v is image coords [i j 1], returning 
	//v = mat*v, a vector in the direction of that image plane point,
	//such that cop + v is on the image plane.
} *Pwfi; //global projection matrix.

static vinc = sizeof(vertex);

/*
struct TriList {
	vertex *list;

	TriList(int num) {list = (vertex *) malloc(num*vinc);}
	//TriList(){}
} *tl; //global vertext list, interpreted as triangle list.
*/
	
Projection::Projection(double dov[], double foc)
/***
This function initializes a projection matrix for a pinhole camera model
using the One True Way (see UNC's Bishop).  Meters are the unit here, and 
the 36mm x 36mm film is modeled. 
***/
{

	double dov_norm = sqrt(double(pow(dov[0],2)+pow(dov[1],2)+pow(dov[2],2)));

	double foc_v[3];
		foc_v[0] = foc*dov[0]/dov_norm;
		foc_v[1] = foc*dov[1]/dov_norm;
		foc_v[2] = foc*dov[2]/dov_norm;
	double foc_norm = sqrt(double(pow(foc_v[0],2)+pow(foc_v[1],2)+pow(foc_v[2],2)));

	//Unit vector in foc_v direction.
	double foc_vu[3];
		foc_vu[0] = foc_v[0]/foc_norm;
		foc_vu[1] = foc_v[1]/foc_norm;
		foc_vu[2] = foc_v[2]/foc_norm;

	double up_foc_dot = up_dir[0]*foc_vu[0]+ up_dir[1]*foc_vu[1]+ up_dir[2]*foc_vu[2];
	
	//The b image plane vector.
	double b[3];
		b[0] = up_dir[0]-up_foc_dot*foc_vu[0];
		b[1] = up_dir[1]-up_foc_dot*foc_vu[1];
		b[2] = up_dir[2]-up_foc_dot*foc_vu[2];
	double b_norm = sqrt(double(pow(b[0],2)+pow(b[1],2)+pow(b[2],2)));
	b[0]=0.036*b[0]/b_norm; b[1]=0.036*b[1]/b_norm; b[2]=0.036*b[2]/b_norm;

	//The a image plane vector, the cross of foc_vu and b.
	double a[3];
		a[0] = foc_vu[1]*b[2]-b[1]*foc_vu[2];
		a[1] = b[0]*foc_vu[2]-foc_vu[0]*b[2];
		a[2] = foc_vu[0]*b[1]-b[0]*foc_vu[1];
	double a_norm = sqrt(double(pow(a[0],2)+pow(a[1],2)+pow(a[2],2)));
	a[0]=0.036*a[0]/a_norm; a[1]=0.036*a[1]/a_norm; a[2]=0.036*a[2]/a_norm;

	//o, a vector pointing to the lower right hand corner of the image plane.
	double o[3];
		o[0] = foc_v[0] - 0.5*b[0] - 0.5*a[0];
		o[1] = foc_v[1] - 0.5*b[1] - 0.5*a[1];
		o[2] = foc_v[2] - 0.5*b[2] - 0.5*a[2];

	//Resize for the resolution, and off we go...
	a[0] = a[0]/xres; a[1] = a[1]/xres; a[2] = a[2]/xres;
	b[0] = b[0]/yres; b[1] = b[1]/yres; b[2] = b[2]/yres;

	for (int i = 0; i < 3; i ++) {
		mat[i][0] = a[i];
		mat[i][1] = b[i];
		mat[i][2] = o[i];
	}

	//Now, a ray going out into the world from the cop through the pixel (i,j)
	//in the image plane is cop + mat*(i j 1)'.
}

void Projection::Convert(double v[])
/***
Simply apply the projection matrix to the referenced vector v, and return the 
new v.
***/
{
	int i,j;
	double mat_b[3] = {0.0,0.0,0.0};

	for (i = 0; i < 3; i ++) {
		for (j = 0; j < 3; j ++) {
			mat_b[i] += mat[i][j]*v[j];
		}
	}

	v[0] = mat_b[0]; v[1] = mat_b[1]; v[2] = mat_b[2];
}

void FrameBuffer::Clear(const rgb p)
{
	register byte r = p[0];
	register byte g = p[1];
	register byte b = p[2];
	register byte* p0 = pix[0];
	register byte* p1 = p0 + sizeof(rgb)*width*height;
	while(p0 < p1) {
		p0[0] = r;
		p0[1] = g;
		p0[2] = b;
		p0 += 3;
	}
}

void Parser()
//This function will fill out the global scene information from
//an input file.
{
	background[0] = 0.0; background[1] = 0.0; background[2] = 0.0;
}

void Parser_Test()
//Fills out global scene information with test data.
{
	double spec[3],amb[3],dif[3], x[3], n[3];
	int i;
	sphere *ball, *curs;
	plane *plast, *curp;
	light *llast, *curl;
 
	 
	background[0] = 0.5; background[1] = 0.5; background[2] = 0.75;
	cop[0] = .098; cop[1] = .0; cop[2] = .08; 
	dov[0] = -.5; dov[1] = 0; dov[2] = -.08; 
	up_dir[0] = 0.0; up_dir[1] = 0.0; up_dir[2] = 5.0;
	foc = 0.036;

	x[0] = 0.098; x[1] = 0.0; x[2] = .08;
	amb[0] = .3; amb[1] = .3; amb[2] = .3;
	dif[0] = .5; dif[1] = .5; dif[2] = .5; 
	spec[0] = .9; spec[1] = .9; spec[2] = .9; 
	sphere_head = new sphere(0.04,x);
	for (i = 0; i < 3; i ++) {
		sphere_head->amb[i] = amb[i];
		sphere_head->dif[i] = dif[i];
		sphere_head->spec[i] = spec[i];
	}
	sphere_head->specpow = 64;
	sphere_head->refl = .0;
	sphere_head->transp = .9;
	sphere_head->index = 1.5;

	x[0] = 0.0; x[1]= -0.05; x[2]= 0.0;
	amb[0] = .4; amb[1] = .2; amb[2] = .2;
	dif[0] = .6; dif[1] = .3; dif[2] = .3; 
	spec[0] = .4; spec[1] = .4; spec[2] = .4;
	ball = new sphere(.04,x);
	for (i = 0; i < 3; i ++) {
		ball->amb[i] = amb[i];
		ball->dif[i] = dif[i];
		ball->spec[i] = spec[i];
	}
	ball->specpow = 64;
	ball->refl = .3;
	ball->transp = 0.0;
	ball->index = 1.0;
	sphere_head->next = ball;

	x[0] = 0.0; x[1]= 0.05; x[2]= 0.0;
	amb[0] = .2; amb[1] = .4; amb[2] = .2;
	dif[0] = .3; dif[1] = .6; dif[2] = .3; 
	spec[0] = .4; spec[1] = .4; spec[2] = .4;
	curs = new sphere(.04,x);
	for (i = 0; i < 3; i ++) {
		curs->amb[i] = amb[i];
		curs->dif[i] = dif[i];
		curs->spec[i] = spec[i];
	}
	curs->specpow = 64;
	curs->refl = .3;
	curs->transp = 0.0;
	curs->index = 1.0;
	ball->next = curs;
	ball = curs;
	
	/* This the set of spheres for model1.

	x[0] = 0.272166; x[1]= 0.272166; x[2]= 0.544331;
	amb[0] = .035; amb[1] = .0325; amb[2] = .025;
	dif[0] = .5; dif[1] = .45; dif[2] = .35; 
	spec[0] = .8; spec[1] = .8; spec[2] = .8; 
	sphere_head = new sphere(0.166667,x);
	for (i = 0; i < 3; i ++) {
		sphere_head->amb[i] = amb[i];
		sphere_head->dif[i] = dif[i];
		sphere_head->spec[i] = spec[i];
	}
	sphere_head->specpow = 64;
	sphere_head->refl = 1.0;
	sphere_head->transp = 0.0;
	sphere_head->index = 1.0;

	x[0] = 0.643951; x[1]= 0.172546; x[2]= 0.000000;
	amb[0] = .035; amb[1] = .0325; amb[2] = .025;
	dif[0] = .5; dif[1] = .45; dif[2] = .35; 
	spec[0] = .8; spec[1] = .8; spec[2] = .8;
	ball = new sphere(.166667,x);
	for (i = 0; i < 3; i ++) {
		ball->amb[i] = amb[i];
		ball->dif[i] = dif[i];
		ball->spec[i] = spec[i];
	}
	ball->specpow = 64;
	ball->refl = 1.0;
	ball->transp = 0.0;
	ball->index = 1.0;
	sphere_head->next = ball;

	x[0] = 0.172546; x[1]= 0.643951; x[2]= 0.000000;
	amb[0] = .035; amb[1] = .0325; amb[2] = .025;
	dif[0] = .5; dif[1] = .45; dif[2] = .35; 
	spec[0] = .8; spec[1] = .8; spec[2] = .8;
	curs = new sphere(.166667,x);
	for (i = 0; i < 3; i ++) {
		curs->amb[i] = amb[i];
		curs->dif[i] = dif[i];
		curs->spec[i] = spec[i];
	}
	curs->specpow = 64;
	curs->refl = 1.0;
	curs->transp = 0.0;
	curs->index = 1.0;
	ball->next = curs;
	ball = curs;

	x[0] = -0.371785; x[1]= 0.099619; x[2]= 0.544331;  
	amb[0] = .035; amb[1] = .0325; amb[2] = .025;
	dif[0] = .5; dif[1] = .45; dif[2] = .35; 
	spec[0] = .8; spec[1] = .8; spec[2] = .8;
	curs = new sphere(.166667,x);
	for (i = 0; i < 3; i ++) {
		curs->amb[i] = amb[i];
		curs->dif[i] = dif[i];
		curs->spec[i] = spec[i];
	}
	curs->specpow = 64;
	curs->refl = 1.0;
	curs->transp = 0.0;
	curs->index = 1.0;
	ball->next = curs;
	ball = curs;

	x[0] = -0.471405; x[1]= 0.471405; x[2]= 0.000000;    
	amb[0] = .035; amb[1] = .0325; amb[2] = .025;
	dif[0] = .5; dif[1] = .45; dif[2] = .35; 
	spec[0] = .8; spec[1] = .8; spec[2] = .8;
	curs = new sphere(.166667,x);
	for (i = 0; i < 3; i ++) {
		curs->amb[i] = amb[i];
		curs->dif[i] = dif[i];
		curs->spec[i] = spec[i];
	}
	curs->specpow = 64;
	curs->refl = 1.0;
	curs->transp = 0.0;
	curs->index = 1.0;
	ball->next = curs;
	ball = curs;

	x[0] = -0.643951; x[1]= -0.172546; x[2]= 0.000000;
	amb[0] = .035; amb[1] = .0325; amb[2] = .025;
	dif[0] = .5; dif[1] = .45; dif[2] = .35; 
	spec[0] = .8; spec[1] = .8; spec[2] = .8;
	curs = new sphere(.166667,x);
	for (i = 0; i < 3; i ++) {
		curs->amb[i] = amb[i];
		curs->dif[i] = dif[i];
		curs->spec[i] = spec[i];
	}
	curs->specpow = 64;
	curs->refl = 1.0;
	curs->transp = 0.0;
	curs->index = 1.0;
	ball->next = curs;
	ball = curs;

	x[0] = 0.099619; x[1]= -0.371785; x[2]= 0.544331;  
	amb[0] = .035; amb[1] = .0325; amb[2] = .025;
	dif[0] = .5; dif[1] = .45; dif[2] = .35; 
	spec[0] = .8; spec[1] = .8; spec[2] = .8;
	curs = new sphere(.166667,x);
	for (i = 0; i < 3; i ++) {
		curs->amb[i] = amb[i];
		curs->dif[i] = dif[i];
		curs->spec[i] = spec[i];
	}
	curs->specpow = 64;
	curs->refl = 1.0;
	curs->transp = 0.0;
	curs->index = 1.0;
	ball->next = curs;
	ball = curs;

	x[0] = -0.172546; x[1]= -0.643951; x[2]= 0.000000;    
	amb[0] = .035; amb[1] = .0325; amb[2] = .025;
	dif[0] = .5; dif[1] = .45; dif[2] = .35; 
	spec[0] = .8; spec[1] = .8; spec[2] = .8;
	curs = new sphere(.166667,x);
	for (i = 0; i < 3; i ++) {
		curs->amb[i] = amb[i];
		curs->dif[i] = dif[i];
		curs->spec[i] = spec[i];
	}
	curs->specpow = 64;
	curs->refl = 1.0;
	curs->transp = 0.0;
	curs->index = 1.0;
	ball->next = curs;
	ball = curs;

	x[0] = 0.471405; x[1]= -0.471405; x[2]= 0.000000;    
	amb[0] = .035; amb[1] = .0325; amb[2] = .025;
	dif[0] = .5; dif[1] = .45; dif[2] = .35; 
	spec[0] = .8; spec[1] = .8; spec[2] = .8;
	curs = new sphere(.166667,x);
	for (i = 0; i < 3; i ++) {
		curs->amb[i] = amb[i];
		curs->dif[i] = dif[i];
		curs->spec[i] = spec[i];
	}
	curs->specpow = 64;
	curs->refl = 1.0;
	curs->transp = 0.0;
	curs->index = 1.0;
	ball->next = curs;
	ball = curs;

	x[0] = 0.0; x[1]= 0.0; x[2]= 0.0;    
	amb[0] = .035; amb[1] = .0325; amb[2] = .025;
	dif[0] = .5; dif[1] = .45; dif[2] = .35; 
	spec[0] = .8; spec[1] = .8; spec[2] = .8;
	curs = new sphere(.5,x);
	for (i = 0; i < 3; i ++) {
		curs->amb[i] = amb[i];
		curs->dif[i] = dif[i];
		curs->spec[i] = spec[i];
	}
	curs->specpow = 64;
	curs->refl = 1.0;
	curs->transp = 0.0;
	curs->index = 1.0;
	ball->next = curs;
	ball = curs;

	//This was the end of the model1 set */
	
	
/*


	x[0] = 0.5; x[1] = 0.5; x[2] = .0;
	amb[0] = .3; amb[1] = .3; amb[2] = .3;
	dif[0] = .6; dif[1] = .6; dif[2] = .6; 
	spec[0] = .3; spec[1] = .3; spec[2] = .3; 
	sphere_head = new sphere(0.4,x);
	for (i = 0; i < 3; i ++) {
		sphere_head->amb[i] = amb[i];
		sphere_head->dif[i] = dif[i];
		sphere_head->spec[i] = spec[i];
	}
	sphere_head->specpow = 64;
	sphere_head->refl = 0.0;
*/

/*balls for the room.

	x[0] = 0.6; x[1] = 0.6; x[2] = .0;
	amb[0] = .3; amb[1] = 3; amb[2] = .3;
	dif[0] = .5; dif[1] = .5; dif[2] = .5; 
	spec[0] = .3; spec[1] = .3; spec[2] = .3; 
	sphere_head = new sphere(0.3,x);
	for (i = 0; i < 3; i ++) {
		sphere_head->amb[i] = amb[i];
		sphere_head->dif[i] = dif[i];
		sphere_head->spec[i] = spec[i];
	}
	sphere_head->specpow = 64;
	sphere_head->refl = 1.0;
	sphere_head->transp = 0.0;
	sphere_head->index = 1.0;

	amb[0] = .4; amb[1] = .0; amb[2] = .0;
	dif[0] = .7; dif[1] = .0; dif[2] = .0; 
	spec[0] = .1; spec[1] = .1; spec[2] = .1; 
	x[0] = 0.8; x[1] = 0.6; x[2] = .3;
	ball = new sphere(.05,x);
	for (i = 0; i < 3; i ++) {
		ball->amb[i] = amb[i];
		ball->dif[i] = dif[i];
		ball->spec[i] = spec[i];
	}
	ball->specpow = 64;
	ball->refl = 0.0;
	ball->transp = 0.0;
	ball->index = 1.0;
	sphere_head->next = ball;

	amb[0] = .2; amb[1] = .2; amb[2] = .2;
	dif[0] = .4; dif[1] = .4; dif[2] = .4; 
	spec[0] = .2; spec[1] = .2; spec[2] = .2; 
	x[0] = .0; x[1] = 0.0; x[2] = 0.0;
	curs = new sphere(.2,x);
	for (i = 0; i < 3; i ++) {
		curs->amb[i] = amb[i];
		curs->dif[i] = dif[i];
		curs->spec[i] = spec[i];
	}
	curs->specpow = 64;
	curs->refl = 0.1;
	curs->transp = 0.0;
	curs->index = 1.0;
	ball->next = curs;
	ball = curs;

	amb[0] = .2; amb[1] = .2; amb[2] = .2;
	dif[0] = .4; dif[1] = .4; dif[2] = .4; 
	spec[0] = .2; spec[1] = .2; spec[2] = .2; 
	x[0] = 1.3; x[1] = .8; x[2] = 0.8;
	curs = new sphere(.2,x);
	for (i = 0; i < 3; i ++) {
		curs->amb[i] = amb[i];
		curs->dif[i] = dif[i];
		curs->spec[i] = spec[i];
	}
	curs->specpow = 64;
	curs->refl = 0.0;
	curs->transp = .7;
	curs->index = .95;
	ball->next = curs;
	ball = curs;
*/
	/*
	x[0] = 0.11; x[1] = 0.15; x[2] = 0.07;
	amb[0] = .2; amb[1] = .4; amb[2] = .2;
	dif[0] = .4; dif[1] = .5; dif[2] = .4; 
	spec[0] = .3; spec[1] = .3; spec[2] = .3; 
	sphere_head = new sphere(0.1,x);
	for (i = 0; i < 3; i ++) {
		sphere_head->amb[i] = amb[i];
		sphere_head->dif[i] = dif[i];
		sphere_head->spec[i] = spec[i];
	}
	sphere_head->specpow = 64;
	sphere_head->refl = 0.0;
	sphere_head->transp = 0.0;
	sphere_head->index = 1.0;

	amb[0] = .3; amb[1] = .0; amb[2] = .0;
	dif[0] = .4; dif[1] = .0; dif[2] = .0; 
	spec[0] = .4; spec[1] = .4; spec[2] = .4; 
	x[0] = .35; x[1] = 0.01; x[2] = -0.04;
	curs = new sphere(.1,x);
	for (i = 0; i < 3; i ++) {
		curs->amb[i] = amb[i];
		curs->dif[i] = dif[i];
		curs->spec[i] = spec[i];
	}
	curs->specpow = 64;
	curs->refl = 0.2;
	curs->index = 1.1;
	curs->transp = .8;
	sphere_head->next = curs;
	ball = curs;


 
  //These are balls on the axes. Transparency numbers will need to be 
  //added.

	amb[0] = .2; amb[1] = .0; amb[2] = .0;
	dif[0] = .6; dif[1] = .0; dif[2] = .0; 
	spec[0] = .3; spec[1] = .0; spec[2] = .0; 
	x[0] = 0.04; x[1] = 0.0; x[2] = .21;
	ball = new sphere(.07,x);
	for (i = 0; i < 3; i ++) {
		ball->amb[i] = amb[i];
		ball->dif[i] = dif[i];
		ball->spec[i] = spec[i];
	}
	ball->specpow = 32;
	sphere_head->next = ball;


	amb[0] = .2; amb[1] = .0; amb[2] = .0;
	dif[0] = .5; dif[1] = .0; dif[2] = .0; 
	spec[0] = .3; spec[1] = .0; spec[2] = .0; 
	x[0] = .5; x[1] = 0.0; x[2] = 0.0;
	ball = new sphere(.2,x);
	for (i = 0; i < 3; i ++) {
		ball->amb[i] = amb[i];
		ball->dif[i] = dif[i];
		ball->spec[i] = spec[i];
	}
	ball->specpow = 32;
	sphere_head->next = ball;

	amb[0] = .0; amb[1] = .2; amb[2] = .0;
	dif[0] = .0; dif[1] = .5; dif[2] = .0; 
	spec[0] = .0; spec[1] = .3; spec[2] = .0; 
	x[0] = .0; x[1] = 0.5; x[2] = 0.0;
	curs = new sphere(.2,x);
	for (i = 0; i < 3; i ++) {
		curs->amb[i] = amb[i];
		curs->dif[i] = dif[i];
		curs->spec[i] = spec[i];
	}
	curs->specpow = 32;
	ball->next = curs;
	ball = curs;

	amb[0] = .0; amb[1] = .0; amb[2] = .2;
	dif[0] = .0; dif[1] = .0; dif[2] = .5; 
	spec[0] = .0; spec[1] = .0; spec[2] = .3; 
	x[0] = .0; x[1] = 0.0; x[2] = 0.5;
	curs = new sphere(.2,x);
	for (i = 0; i < 3; i ++) {
		curs->amb[i] = amb[i];
		curs->dif[i] = dif[i];
		curs->spec[i] = spec[i];
	}
	curs->specpow = 32;
	ball->next = curs;
	ball = curs;






*/
	amb[0] = .2; amb[1] = .2; amb[2] = .3;
	dif[0] = .4; dif[1] = .4; dif[2] = .6; 
	spec[0] = .4; spec[1] = .4; spec[2] = .4;
	x[0] = 0; x[1] = 0; x[2] = 0;
	n[0] = .0; n[1] = .0; n[2] = 1.0;
	plane_head = new plane(x,n);
	for (i = 0; i < 3; i ++) {
		plane_head->amb[i] = amb[i];
		plane_head->dif[i] = dif[i];
		plane_head->spec[i] = spec[i];
	}
	plane_head->specpow = 32;
	plane_head->refl = .3;
	plane_head->index = 1.0;
	plane_head->transp = 0.0;

	amb[0] = .1; amb[1] = .1; amb[2] = .1;
	dif[0] = .1; dif[1] = .1; dif[2] = .1; 
	spec[0] = .9; spec[1] = .9; spec[2] = .9;
	x[0] = -0.1; x[1] = 0; x[2] = 0;
	n[0] = 1.0; n[1] = .0; n[2] = 0.0;
	curp = new plane(x,n);
	for (i = 0; i < 3; i ++) {
		curp->amb[i] = amb[i];
		curp->dif[i] = dif[i];
		curp->spec[i] = spec[i];
	}
	plane_head->next = curp;
	curp->specpow = 30;
	curp->refl = 1.0;
	curp->transp = 0.0;
	curp->index = 1.0;
	plast = curp;

	amb[0] = .1; amb[1] = .1; amb[2] = .1;
	dif[0] = .1; dif[1] = .1; dif[2] = .1; 
	spec[0] = .9; spec[1] = .9; spec[2] = .9;
	x[0] = .1; x[1] = 0.0; x[2] = 0;
	n[0] = -1.0; n[1] = 0.0; n[2] = 0.0;
	curp = new plane(x,n);
	for (i = 0; i < 3; i ++) {
		curp->amb[i] = amb[i];
		curp->dif[i] = dif[i];
		curp->spec[i] = spec[i];
	}
	plast->next = curp;
	curp->specpow = 30;
	curp->refl = 1.0;
	curp->index = 1.0;
	curp->transp = .0;
	plast = curp;

/*
	amb[0] = .1; amb[1] = .1; amb[2] = .1;
	dif[0] = .1; dif[1] = .1; dif[2] = .1; 
	spec[0] = .4; spec[1] = .4; spec[2] = .4;
	x[0] = 0.0; x[1] = 0; x[2] = 0;
	n[0] = 1.0; n[1] = .0; n[2] = 0.0;
	curp = new plane(x,n);
	for (i = 0; i < 3; i ++) {
		curp->amb[i] = amb[i];
		curp->dif[i] = dif[i];
		curp->spec[i] = spec[i];
	}
	plane_head->next = curp;
	curp->specpow = 128;
	curp->refl = .9;
	curp->transp = 0.0;
	curp->index = 1.0;
	plast = curp;

	amb[0] = .1; amb[1] = .1; amb[2] = .1;
	dif[0] = .1; dif[1] = .1; dif[2] = .1; 
	spec[0] = .4; spec[1] = .4; spec[2] = .4;
	x[0] = 1.0; x[1] = 0.0; x[2] = 0;
	n[0] = -1.0; n[1] = 0.0; n[2] = 0.0;
	curp = new plane(x,n);
	for (i = 0; i < 3; i ++) {
		curp->amb[i] = amb[i];
		curp->dif[i] = dif[i];
		curp->spec[i] = spec[i];
	}
	plast->next = curp;
	curp->specpow = 128;
	curp->refl = .9;
	plast = curp;


	amb[0] = .1; amb[1] = .1; amb[2] = .1;
	dif[0] = .1; dif[1] = .1; dif[2] = .1; 
	spec[0] = .4; spec[1] = .4; spec[2] = .4;
	x[0] = 0; x[1] = .0; x[2] = 0;
	n[0] = 0.0; n[1] = 1.0; n[2] = 0.0;
	curp = new plane(x,n);
	for (i = 0; i < 3; i ++) {
		curp->amb[i] = amb[i];
		curp->dif[i] = dif[i];
		curp->spec[i] = spec[i];
	}
	plast->next = curp;
	curp->specpow = 128;
	curp->refl = .9;
	curp->transp = 0.0;
	curp->index = 1.0;
	plast = curp;


	amb[0] = .3; amb[1] = .; amb[2] = .3;
	dif[0] = .6; dif[1] = .; dif[2] = .8; 
	spec[0] = .0; spec[1] = .0; spec[2] = .0;
	x[0] = 0; x[1] = 0; x[2] = -0.5;
	n[0] = .0; n[1] = .0; n[2] = 1.0;
	plane_head = new plane(x,n);
	for (i = 0; i < 3; i ++) {
		plane_head->amb[i] = amb[i];
		plane_head->dif[i] = dif[i];
		plane_head->spec[i] = spec[i];
	}
	plane_head->specpow = 32;
	plane_head->refl = .0;
	plane_head->index = 1.0;
	plane_head->transp = 0.0;
*/

/*
The planes for the room.

	amb[0] = .2; amb[1] = .2; amb[2] = .4;
	dif[0] = .4; dif[1] = .4; dif[2] = .6; 
	spec[0] = .4; spec[1] = .4; spec[2] = .4;
	x[0] = 0; x[1] = 0; x[2] = 0;
	n[0] = 1.0; n[1] = .0; n[2] = 0.0;
	plane_head = new plane(x,n);
	for (i = 0; i < 3; i ++) {
		plane_head->amb[i] = amb[i];
		plane_head->dif[i] = dif[i];
		plane_head->spec[i] = spec[i];
	}
	plane_head->specpow = 32;
	plane_head->refl = 1.0;
	plane_head->index = 1.0;
	plane_head->transp = 0.0;


	amb[0] = .2; amb[1] = .2; amb[2] = .4;
	dif[0] = .4; dif[1] = .4; dif[2] = .6; 
	spec[0] = .4; spec[1] = .4; spec[2] = .4;
	x[0] = .0; x[1] = .0; x[2] = .0;
	n[0] = .0; n[1] = 1.0; n[2] = 0.0;
	curp = new plane(x,n);
	for (i = 0; i < 3; i ++) {
		curp->amb[i] = amb[i];
		curp->dif[i] = dif[i];
		curp->spec[i] = spec[i];
	}
	plane_head->next = curp;
	curp->specpow = 32;
	curp->refl = 1.0;
	curp->transp = 0.0;
	curp->index = 1.0;
	plast = curp;

	amb[0] = .2; amb[1] = .2; amb[2] = .2;
	dif[0] = .2; dif[1] = .2; dif[2] = .2; 
	spec[0] = .0; spec[1] = .0; spec[2] = .0;
	x[0] = 0; x[1] = .0; x[2] = .0;
	n[0] = 0.0; n[1] = .0; n[2] = 1.0;
	curp = new plane(x,n);
	for (i = 0; i < 3; i ++) {
		curp->amb[i] = amb[i];
		curp->dif[i] = dif[i];
		curp->spec[i] = spec[i];
	}
	plast->next = curp;
	curp->specpow = 32;
	curp->refl = .7;
	curp->transp = 0.0;
	curp->index = 1.0;
	plast = curp;

	//colored planes.

	amb[0] = .2; amb[1] = .2; amb[2] = .4;
	dif[0] = .2; dif[1] = .2; dif[2] = .6; 
	spec[0] = .0; spec[1] = .0; spec[2] = .0;
	x[0] = 0; x[1] = .0; x[2] = 1.5;
	n[0] = 0.0; n[1] = .0; n[2] = -1.0;
	curp = new plane(x,n);
	for (i = 0; i < 3; i ++) {
		curp->amb[i] = amb[i];
		curp->dif[i] = dif[i];
		curp->spec[i] = spec[i];
	}
	plast->next = curp;
	curp->specpow = 32;
	curp->refl = 0.0;
	curp->transp = 0.0;
	curp->index = 1.0;
	plast = curp;

	amb[0] = .3; amb[1] = .2; amb[2] = .2;
	dif[0] = .4; dif[1] = .2; dif[2] = .2; 
	spec[0] = .0; spec[1] = .0; spec[2] = .0;
	x[0] = 1.5; x[1] = .0; x[2] = 0.0;
	n[0] = -1.0; n[1] = .0; n[2] = 0.0;
	curp = new plane(x,n);
	for (i = 0; i < 3; i ++) {
		curp->amb[i] = amb[i];
		curp->dif[i] = dif[i];
		curp->spec[i] = spec[i];
	}
	plast->next = curp;
	curp->specpow = 32;
	curp->refl = 0.0;
	curp->transp = 0.0;
	curp->index = 1.0;
	plast = curp;

	amb[0] = .2; amb[1] = .3; amb[2] = .2;
	dif[0] = .2; dif[1] = .4; dif[2] = .2; 
	spec[0] = .0; spec[1] = .0; spec[2] = .0;
	x[0] = 0; x[1] = 1.5; x[2] = 0.0;
	n[0] = 0.0; n[1] = -1.0; n[2] = 0.0;
	curp = new plane(x,n);
	for (i = 0; i < 3; i ++) {
		curp->amb[i] = amb[i];
		curp->dif[i] = dif[i];
		curp->spec[i] = spec[i];
	}
	plast->next = curp;
	curp->specpow = 32;
	curp->refl = 0.0;
	curp->transp = 0.0;
	curp->index = 1.0;
	plast = curp;

	The end of the colored planes for the room.
		*/ 

/* The lights for model1.
	amb[0] = 0.288675; amb[1] = 0.288675; amb[2] = 0.288675;
	dif[0] = 0.288675; dif[1] = 0.288675; dif[2] = 0.288675; 
	spec[0] = 0.288675; spec[1] = 0.288675; spec[2] = 0.288675;
	x[0] = 4.0; x[1] = 3.0; x[2] = 2.0;
	light_head = new light(x);
	for (i = 0; i < 3; i ++) {
		light_head->amb[i] = amb[i];
		light_head->dif[i] = dif[i];
		light_head->spec[i] = spec[i];
	}

	amb[0] = 0.288675; amb[1] = 0.288675; amb[2] = 0.288675;
	dif[0] = 0.288675; dif[1] = 0.288675; dif[2] = 0.288675; 
	spec[0] = 0.288675; spec[1] = 0.288675; spec[2] = 0.288675;
	x[0] = 1.0; x[1] = -4.0; x[2] = 4.0;
	llast = new light(x);
	for (i = 0; i < 3; i ++) {
		llast->amb[i] = amb[i];
		llast->dif[i] = dif[i];
		llast->spec[i] = spec[i];
	}
	light_head->next = llast;

	amb[0] = 0.288675; amb[1] = 0.288675; amb[2] = 0.288675;
	dif[0] = 0.288675; dif[1] = 0.288675; dif[2] = 0.288675; 
	spec[0] = 0.288675; spec[1] = 0.288675; spec[2] = 0.288675;
	x[0] = -3.0; x[1] = 1.0; x[2] = 5.0;
	curl = new light(x);
	for (i = 0; i < 3; i ++) {
		curl->amb[i] = amb[i];
		curl->dif[i] = dif[i];
		curl->spec[i] = spec[i];
	}
	llast->next = curl;
*/
/*
	amb[0] = 1.0; amb[1] = 1.0; amb[2] = 1.0;
	dif[0] = 1.0; dif[1] = 1.0; dif[2] = 1.0; 
	spec[0] = 1.0; spec[1] = 1.0; spec[2] = 1.0;
	x[0] = .45; x[1] = .7; x[2] = .3;
	light_head = new light(x);
	for (i = 0; i < 3; i ++) {
		light_head->amb[i] = amb[i];
		light_head->dif[i] = dif[i];
		light_head->spec[i] = spec[i];
	}

	amb[0] = .7; amb[1] = .7; amb[2] = .7;
	dif[0] = .7; dif[1] = .7; dif[2] = .7; 
	spec[0] = .7; spec[1] = .7; spec[2] = .7;
	x[0] = .4; x[1] = .1; x[2] = 1.2;
	llast = new light(x);
	for (i = 0; i < 3; i ++) {
		llast->amb[i] = amb[i];
		llast->dif[i] = dif[i];
		llast->spec[i] = spec[i];
	}
	light_head->next = llast;
*/

	amb[0] = .5; amb[1] = .5; amb[2] = .5;
	dif[0] = .5; dif[1] = .5; dif[2] = .5; 
	spec[0] = .5; spec[1] = .5; spec[2] = .5;
	x[0] = .05; x[1] = .3; x[2] = .2;
	light_head = new light(x);
	for (i = 0; i < 3; i ++) {
		light_head->amb[i] = amb[i];
		light_head->dif[i] = dif[i];
		light_head->spec[i] = spec[i];
	}

	amb[0] = .5; amb[1] = .5; amb[2] = .5;
	dif[0] = .5; dif[1] = .5; dif[2] = .5; 
	spec[0] = .5; spec[1] = .5; spec[2] = .5;
	x[0] = .05; x[1] = -.3; x[2] = .2;
	llast = new light(x);
	for (i = 0; i < 3; i ++) {
		llast->amb[i] = amb[i];
		llast->dif[i] = dif[i];
		llast->spec[i] = spec[i];
	}
	light_head->next = llast;

	amb[0] = .3; amb[1] = .3; amb[2] = .3;
	dif[0] = .3; dif[1] = .3; dif[2] = .3; 
	spec[0] = .3; spec[1] = .3; spec[2] = .3;
	x[0] = .98; x[1] = .98; x[2] = .98;
	curl = new light(x);
	for (i = 0; i < 3; i ++) {
		curl->amb[i] = amb[i];
		curl->dif[i] = dif[i];
		curl->spec[i] = spec[i];
	}
	llast->next = curl;
/*
	amb[0] = .4; amb[1] = .4; amb[2] = .4;
	dif[0] = .4; dif[1] = .4; dif[2] = .4; 
	spec[0] = .4; spec[1] = .4; spec[2] = .4;
	x[0] = .01; x[1] = .01; x[2] = .98;
	curl = new light(x);
	for (i = 0; i < 3; i ++) {
		curl->amb[i] = amb[i];
		curl->dif[i] = dif[i];
		curl->spec[i] = spec[i];
	}
	llast->next = curl;

*/


}

double Intersect_sphere(double e[], double v[], double c[], double r)
//Returns the extend along v to touch the sphere. Returns negative t if there's
//no real intersection, as no intersection is treated in the same way as an
//imaginary intersection.
{
	double temp1 = pow(v[0],2) + pow(v[1],2) + pow(v[2],2);
	double temp2 = v[0]*(e[0]-c[0]) + v[1]*(e[1]-c[1]) + v[2]*(e[2]-c[2]);
	temp2 *= 2.0;
	double temp3 = pow(e[0],2) + pow(e[1],2) + pow(e[2],2);
	temp3 += pow(c[0],2) + pow(c[1],2) + pow(c[2],2);
	temp3 = temp3 - 2.0*(e[0]*c[0] + e[1]*c[1] + e[2]*c[2]) - pow(r,2);


	double big = pow(temp2,2) - 4.0*temp1*temp3;
	if (big > 0) {
		big = sqrt(big)/(2.0*temp1);
		double opb = -temp2/(2.0*temp1);

		if (opb > big + 1.0e-10) {return (opb - big);}
		else {return (opb + big);}
	}else {
		return -1.0;
	}
}
	
double Intersect_plane(double e[], double v[], double p[], double n[])
//Returns the extent along v from e to touch the plane. Again, negative
//return value won't be viewed anyway, to if there isn't an intersection,
//negative return value.
{
	double D = n[0]*p[0] + n[1]*p[1] + n[2]*p[2];
	double numer = D - (n[0]*e[0] + n[1]*e[1] + n[2]*e[2]);
	double denom = n[0]*v[0] + n[1]*v[1] + n[2]*v[2];
	if (fabs(denom) < 1.0e-14) {
		return -1.0;
	}else {
		return (numer/denom);
	}
}

//debug 
//int q = 0, z = 0;

void Single_Raytrace(double v[], double p[], double color[], int level, double refract)
{
	//color[0] = 1.0; color[1] = 0.0;  color[2] = 1.0;
	double min_t = 1.0e90, t, sum_l[3], lightv[3];
	int s_p = 0, i, snp = 0; //-1 for sphere, 1 for plane.
	sphere *cl_s, *curs = sphere_head;
	plane *cl_p, *curp = plane_head;	//plane,sphere, and light pointers.
	light *curl = light_head;
	double am[3], dif[3], spe[3], ex;
	double rind = 0, tind = 0, newrefract = 0; //reflection and transparency info.
	double sum_spec[3];
	//I'm going to keep the specular information separate from other mat properties,
	//so that reflective and transparent surfaces can still have specular highlights.

	//First, intersect the ray with all the spheres.

	while (curs != NULL) {
		t = Intersect_sphere(p, v, curs->c, curs->rad);
		if ((t < min_t)&&(t > 1.0e-10)) {
			min_t = t;
			cl_s = curs;
			s_p = -1;
			//if (q == 0) printf("%lf\n",t);
			//q++;
		}
		curs = curs->next;
	}
	//Now, check for a closer plane.
	while (curp != NULL) {
		t = Intersect_plane(p, v, curp->p, curp->n);
		if ((t < min_t)&&(t > 1.0e-10)) {
			min_t = t;
			cl_p = curp;
			s_p = 1;
		}
		curp = curp->next;
	}

	//A separate case is if the ray hits a light before all of the others.

	light *cl_l;  //10 lights max will add some specular highlight to the surface.
	double vmag,lmag,ldot, mindot = .01;
	vmag = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
	while (curl != NULL) {
		//Use vmag as a ruler for determining how close a light is.  If the light
		//is closer and within (arbitrary) 30 degrees, then the light will be the
		//color that this ray returns.  This works best for lights outside of the 
		//scene volume.
		for (i=0;i<3;i++) {
			lightv[i] = curl->p[i] - p[i];
		}
		lmag = sqrt(lightv[0]*lightv[0]+lightv[1]*lightv[1]+lightv[2]*lightv[2]);
		ldot = v[0]*lightv[0] + v[1]*lightv[1] + v[2]*lightv[2];
		ldot = ldot/(lmag*vmag);
		lmag = lmag/vmag;
		if ((lmag < min_t)&&(lmag > 1.0e-10)&&(ldot >= .85)&&(ldot > mindot)) {
			//for angles between 0 and pi/6, cos varies from ~.85 to 1.
			mindot = ldot;
			cl_l = curl;
			snp = -121;
			//if (q == 0) printf("%lf\n",t);
			//q++;
		}
		curl = curl->next;
	}


	//Now do the light computation if there was a hit.
	if ((s_p ==  1)||(s_p == -1)) {
		sum_l[0] = 0.0; sum_l[1] = 0.0; sum_l[2] = 0.0;
		sum_spec[0] = 0.0; sum_spec[1] = 0.0; sum_spec[2] = 0.0; 
		//Pick up the material properties.
		if (s_p == -1) {
			am[0] = cl_s->amb[0]; am[1] = cl_s->amb[1]; am[2] = cl_s->amb[2];
			dif[0] = cl_s->dif[0]; dif[1] = cl_s->dif[1]; dif[2] = cl_s->dif[2];
			spe[0] = cl_s->spec[0]; spe[1] = cl_s->spec[1]; spe[2] = cl_s->spec[2];
			ex = cl_s->specpow;
		}else if (s_p == 1) {
			am[0] = cl_p->amb[0]; am[1] = cl_p->amb[1]; am[2] = cl_p->amb[2];
			dif[0] = cl_p->dif[0]; dif[1] = cl_p->dif[1]; dif[2] = cl_p->dif[2];
			spe[0] = cl_p->spec[0]; spe[1] = cl_p->spec[1]; spe[2] = cl_p->spec[2];
			ex = cl_p->specpow;
		}
		//Done picking up mat. props of the intersected object.

		double ip[3]; //Point of intersection of the ray and object.
		ip[0] = p[0]+min_t*v[0]; ip[1] = p[1]+min_t*v[1]; ip[2] = p[2]+min_t*v[2];
		
		double n[3];  //Normal at the point of intersection.  Additional modifications
					//will be needed for phong shaded triangles.
		if (s_p == -1) {
			n[0] = ip[0]-cl_s->c[0]; n[1] = ip[1]-cl_s->c[1]; n[2] = ip[2]-cl_s->c[2];
		}else if (s_p == 1) {
			n[0] = cl_p->n[0]; n[1] = cl_p->n[1]; n[2] = cl_p->n[2];
		}

		double norm = sqrt(pow(n[0],2) + pow(n[1],2) + pow(n[2],2));
		n[0] = n[0]/norm; n[1] = n[1]/norm; n[2] = n[2]/norm;

		double view[3];  //The view vector, from intersection to eye.
		view[0] = p[0] - ip[0]; view[1] = p[1] - ip[1]; view[2] = p[2] - ip[2];
		norm = sqrt(pow(view[0],2) + pow(view[1],2) + pow(view[2],2));
		view[0] = view[0]/norm; view[1] = view[1]/norm; view[2] = view[2]/norm; 

		//More variables, for the light calculations.
		double l[3], l_dif[3], l_spe[3], dot1, dot2, ref[3], add_l[3];
		int in_shadow = 0;
		
		curl = light_head;
		double shadind;  //If occluding surface is transparent.
		while (curl != NULL) {
			//The light values serve as a filter for the material properties of
			//the object intersected.
			in_shadow = 0;
			shadind = 1.0;
			for (i = 0; i < 3; i ++) {
				sum_l[i] += curl->amb[i];
				l[i] = curl->p[i] - ip[i];
				l_dif[i] = curl->dif[i];
				l_spe[i] = curl->spec[i];
			}

			norm = sqrt(pow(l[0],2) + pow(l[1],2) + pow(l[2],2));
			l[0] = l[0]/norm; l[1] = l[1]/norm; l[2] = l[2]/norm; 
			//Now, l is a unit light vector from the intersection point.


			//First we have to check if the point is in shadow for this light.
			curs = sphere_head;
			curp = plane_head;

			while ((curs != NULL)&&(!in_shadow)) {
				t = Intersect_sphere(ip, l, curs->c, curs->rad);
				if ((t < norm)&&(t > 1.0e-8)) {
					//if t is less than the magnitude of the vector from ip to the
					//light, then the sphere is closer than the light to ip.
					in_shadow = 1;
					shadind = shadind*curs->transp;
					//if (s_p == -1) {
					//	t = Intersect_sphere(ip, l, curs->c, curs->rad);
					//}
						
				}
				curs = curs->next;
			}
			//Now, check for a closer plane.
			while ((curp != NULL)&&(!in_shadow)) {
				t = Intersect_plane(ip, l, curp->p, curp->n);
				if ((t < norm)&&(t > 1.0e-8)) {
					in_shadow = 1;
					shadind = curp->transp;
				}
				curp = curp->next;
			}
	
			if (!in_shadow) {
				add_l[0] = 0.0; add_l[1] = 0.0; add_l[2] = 0.0;
				dot1 = l[0]*n[0] + l[1]*n[1] + l[2]*n[2];
				if (dot1 >= 0.0) {
					ref[0] = 2.0*dot1*n[0]-l[0]; ref[1] = 2.0*dot1*n[1]-l[1]; 
					ref[2] = 2.0*dot1*n[2]-l[2];
					dot2 = view[0]*ref[0] + view[1]*ref[1] + view[2]*ref[2];
					if (dot2 < 0.0) dot2 = 0.0;
	
					add_l[0] = dot1*dif[0]*l_dif[0] + pow(dot2,ex)*spe[0]*l_spe[0];
					add_l[1] = dot1*dif[1]*l_dif[1] + pow(dot2,ex)*spe[1]*l_spe[1];
					add_l[2] = dot1*dif[2]*l_dif[2] + pow(dot2,ex)*spe[2]*l_spe[2];

					sum_spec[0] += shadind*pow(dot2,ex)*spe[0]*l_spe[0];
					sum_spec[1] += shadind*pow(dot2,ex)*spe[1]*l_spe[1];
					sum_spec[2] += shadind*pow(dot2,ex)*spe[2]*l_spe[2];
				}
	
				//Now, update the color to be returned.
				color[0] += shadind*add_l[0]; color[1] += shadind*add_l[1]; 
				color[2] += shadind*add_l[2];
			}

			curl = curl->next; //Do the same for the next light.
		}

		//Add in the ambient terms and threshold the color value.
		for (i=0;i<3;i++) {
			if (sum_l[i] > 1.0) sum_l[i] = 1.0;
			color[i] += am[i]*sum_l[i];
			if (color[i] > 1.0) color[i] = 1.0;
		}

		//Now we can get the color from reflection and refraction rays.

		if (level < 20) {
			//color[0] = .7*color[0]; color[1] = .7*color[1]; color[2] = .7*color[2];

			double ref_col[3], trans_col[3];
			double refl[3];
			trans_col[0] = 0.0; trans_col[1] = 0.0; trans_col[2] = 0.0;
			ref_col[0] = 0.0; ref_col[1] = 0.0; ref_col[2] = 0.0;

			//Make the reflection ray.
			//dot1 =  view[0]*n[0] + view[1]*n[1] + view[2]*n[2];
			//refl[0] = -view[0]+2.0*dot1*n[0]; refl[1] = -view[1]+2.0*dot1*n[1];
			//refl[2] = -view[2]+2.0*dot1*n[2];


			double rind = 0, tind = 0, newrefract = 0; 
			if (s_p == -1) {
				rind = cl_s->refl;		//Weightings on transparencey and reflection.
				tind = cl_s->transp;
				newrefract = cl_s->index;
			}else if (s_p == 1) {
				rind = cl_p->refl;
				tind = cl_p->transp;
				newrefract = cl_p->index;
			}

			//Make the transparency ray.
			//double theta = acos(dot1), temp1[3], ratio = refract/newrefract;
			//dot1 = -1.0*view[0]*n[0] + -1.0*view[1]*n[1] + -1.0*view[2]*n[2];
			//temp1[0] = -view[0]-dot1*n[0]; temp1[1] = -view[1]-dot1*n[1]; 
			//temp1[2] = -view[2]-dot1*n[2]; 

			if (rind > 0.001) {
				//Make the reflection ray.
				dot1 =  view[0]*n[0] + view[1]*n[1] + view[2]*n[2];
				refl[0] = -view[0]+2.0*dot1*n[0]; refl[1] = -view[1]+2.0*dot1*n[1];
				refl[2] = -view[2]+2.0*dot1*n[2];

				Single_Raytrace(refl, ip, ref_col, level + 1, refract);
			}
			if (tind > 0.001) {
//				double tran[3];
				//Make the transparency ray.
				if (refract == 0 || newrefract == 0)
					printf ("Refraction component zero. uh-oh.\n");

				dot1 = -1.0*view[0]*n[0] + -1.0*view[1]*n[1] + -1.0*view[2]*n[2];
				double norm1 = sqrt(pow(view[0],2)+pow(view[1],2)+pow(view[2],2));
				norm1 *= sqrt(n[0]*n[0]+ n[1]*n[1]+ n[2]*n[2]);
				dot1 = dot1/norm1;
/*

				double dotnv, dotbv, norm2;
				dotnv = view[0]*n[0] + view[1]*n[1] + view[2]*n[2];
				double norm1 = sqrt(pow(view[0],2)+pow(view[1],2)+pow(view[2],2));
				norm1 *= sqrt(n[0]*n[0]+ n[1]*n[1]+ n[2]*n[2]);
				dotnv = dotnv/norm1;

				if ((fabs(fabs(dotnv)-1.0) < 1.0e-12)||(fabs(dotnv) < 1.0e-12)) {
					tran[0] = -view[0]; tran[1] = -view[1]; tran[2] = -view[2];
				}else {
					double ncview[3];
					ncview[0] = n[1]*view[2] - view[1]*n[2];
					ncview[1] = view[0]*n[2] - n[0]*view[2];
					ncview[2] = n[0]*view[1] - view[0]*n[1];
					double b[3];
					b[0] = n[1]*ncview[2] - ncview[1]*n[2];
					b[1] = ncview[0]*n[2] - n[0]*ncview[2];
					b[2] = n[0]*ncview[1] - ncview[0]*n[1];
					norm2 = sqrt(b[0]*b[0] + b[1]*b[1] + b[2]*b[2]);
					dotbv = b[0]*view[0] + b[1]*view[1] + b[2]*view[2];
					norm2 *= sqrt(view[0]*view[0] + view[1]*view[1] + view[2]*view[2]);
					dotbv = dotbv/norm2;
					//With knowledge of the view dotted with both the normal and the b
					//vector (defining which half of the refraction plane the result t
					//will be in), we can make the right choices and not be confused by
					//domain issues related to sine and cosine.

					double theta, sinthetaprim, costhetaprim, ratio;
					if (dotnv < 0) { 
						ratio = refract/newrefract;
					}
					else { 
						ratio = newrefract/1.0;
					}

					if (dotnv > 0) {
						if (dotbv > 0) {
							printf ("b dot view is positive?\n");
						}else {
							//theta should be < pi/2.
							z++;
							theta = acos(dotnv);
							//sinthetaprim = ratio*sin(theta);
							sinthetaprim = ratio*sqrt(1-dotnv*dotnv);
							//if (sinthetaprim > 1.0) sinthetaprim = 1.0;
							costhetaprim = sqrt(1 - sinthetaprim*sinthetaprim);
							tran[0] = sinthetaprim*b[0] - costhetaprim*n[0];
							tran[1] = sinthetaprim*b[1] - costhetaprim*n[1];
							tran[2] = sinthetaprim*b[2] - costhetaprim*n[2];
						}
					}else {
						if (dotbv > 0) {
							printf("2. b dot view is positive?\n");
						}else {
							q++;
							theta = acos(dotnv);
							sinthetaprim = ratio*sin(3.14159265 - theta);
							//if (sinthetaprim > 1.0) sinthetaprim = 1.0;
							costhetaprim = sqrt(1 - sinthetaprim*sinthetaprim);
							tran[0] = sinthetaprim*b[0] + costhetaprim*n[0];
							tran[1] = sinthetaprim*b[1] + costhetaprim*n[1];
							tran[2] = sinthetaprim*b[2] + costhetaprim*n[2];
						}
					}

					norm1 = sqrt(tran[0]*tran[0]+tran[1]*tran[1]+tran[2]*tran[2]);
					tran[0] = tran[0]/norm1; tran[1]=tran[1]/norm1; tran[2]=tran[2]/norm1; 

				}


*/




				double sintheta, ratio;
				if (dot1 < 1.0e-12) { 
					ratio = refract/newrefract;
				}
				else { 
					ratio = newrefract/1.0;
				}
				sintheta = sqrt(1.0 - dot1*dot1);
				double sinthetaprim = ratio*sintheta;
				double costhetaprim;
				if (dot1 < 0.0) 
					costhetaprim = -sqrt(1.0 - pow(sinthetaprim,2));
				else costhetaprim = sqrt(1.0 - pow(sinthetaprim,2));
				if (sinthetaprim > 1.0) {
					costhetaprim = .0; sinthetaprim = 1.0;
				}
				double b[3];
				b[0] = -view[0] + dot1*n[0]; b[1] = -view[1] + dot1*n[1]; 
				b[2] = -view[2] + dot1*n[2];
				if (fabs(sintheta) > 1.0e-10) {
					b[0] = b[0]/sintheta; b[1] = b[1]/sintheta; b[2] = b[2]/sintheta; 
				}
				double tran[3];
				if ((fabs(dot1) < 1.0e-10)||(fabs(sintheta) < 1.0e-10)) {
					tran[0] = -view[0]; tran[1] = -view[1]; tran[2] = -view[2];
				}else {
					tran[0] = sinthetaprim*b[0]-costhetaprim*n[0];
					tran[1] = sinthetaprim*b[1]-costhetaprim*n[1];
					tran[2] = sinthetaprim*b[2]-costhetaprim*n[2];
				}
/*
				double theta = acos(dot1), tran[3], ratio, ctemp;
				if (dot1 < 1.0e-12) ratio = refract/newrefract;
				else ratio = newrefract/1.0;
				tran[0] = -view[0]-dot1*n[0]; tran[1] = -view[1]-dot1*n[1]; 
				tran[2] = -view[2]-dot1*n[2]; 
				tran[0] = ratio*tran[0]; tran[1] = ratio*tran[1]; 
				tran[2] = ratio*tran[2]; 
				ctemp = 1 - theta*pow(dot1,2);
				ctemp *= pow(ratio, 2);
				ctemp = sqrt(1.0 - ctemp);
				tran[0] -= ctemp*n[0]; tran[1] -= ctemp*n[1]; tran[2] -= ctemp*n[2]; 
				//whew
*/				

				if (dot1 < 0)
					Single_Raytrace(tran, ip, trans_col, level + 1, newrefract);
				else Single_Raytrace(tran, ip, trans_col, level + 1, 1.0);
			}
			dot1 = 1.0 - tind - rind;
			color[0] = dot1*color[0] + rind*ref_col[0] + tind*trans_col[0]; 
			color[1] = dot1*color[1] + rind*ref_col[1] + tind*trans_col[1]; 
			color[2] = dot1*color[2] + rind*ref_col[2] + tind*trans_col[2];

			//Now add the specular highlights in.  This is simpler than the commented
			//out code below.

			color[0] += sum_spec[0];
			color[1] += sum_spec[1];
			color[2] += sum_spec[2];


			
/*  This was my light rendering hack.  I will just keep the specular highlights
untempered by the reflection and refraction components, as the diffuse and ambient
are not.

			if (snp == -121) {
				//the ray got close to a light. So return part of the color of the light.
				//printf("the ray hit a light\n");
				double temppwr;
				for (i=0;i<3;i++) {
					lightv[i] = cl_l->p[i] - p[i];
				}
				temppwr = sqrt(lightv[0]*lightv[0]+lightv[1]*lightv[1]+lightv[2]*lightv[2]);
				lightv[0]=lightv[0]/temppwr; lightv[1]=lightv[1]/temppwr; 
				lightv[2]=lightv[2]/temppwr; 
				if (s_p == -1) {
					t = Intersect_sphere(p, lightv, cl_s->c, cl_s->rad);
					if ((t < temppwr)&&(t > 1.0e-8)) temppwr = 0.0;
					else temppwr = pow(mindot,25);
				}else if (s_p == 1) {
					t = Intersect_plane(p,lightv, cl_p->p, cl_p->n);
					if ((t < temppwr)&&(t > 1.0e-8)) temppwr = 0.0;
					else temppwr = pow(mindot,25);
				}else {
					temppwr = pow(mindot,25);
				}
				color[0] += temppwr*(cl_l->amb[0] + cl_l->dif[0] + cl_l->spec[0]);
				color[1] += temppwr*(cl_l->amb[1] + cl_l->dif[1] + cl_l->spec[1]);
				color[2] += temppwr*(cl_l->amb[2] + cl_l->dif[2] + cl_l->spec[2]);
			}
*/
		}

		//printf(".");
		//color[0] = 1.0; color[1] = 1.0; color[2] = 1.0;
	}else {
		color[0] = background[0];
		color[1] = background[1];
		color[2] = background[2];

		if (snp == -121) {
			//the ray got close to a light. So return part of the color of the light.
			//printf("the ray hit a light\n");
			double temppwr;
				temppwr = pow(mindot,25);
				color[0] += temppwr*(cl_l->amb[0] + cl_l->dif[0] + cl_l->spec[0]);
				color[1] += temppwr*(cl_l->amb[1] + cl_l->dif[1] + cl_l->spec[1]);
				color[2] += temppwr*(cl_l->amb[2] + cl_l->dif[2] + cl_l->spec[2]);
		}

	}

	for(i=0;i<3;i++) 
	if (color[i] > 1.0) color[i] = 1.0;
}

/*
void Get_RandTri()
{
	//This function will put a random triangle in list.
	int q[3];

	q[0] = 255*(1.0*rand())/RAND_MAX;
	q[1] = 255*(1.0*rand())/RAND_MAX;
	q[2] = 255*(1.0*rand())/RAND_MAX;

	//cerr << RAND_MAX << endl;
	for (int i = 0; i < 3; i ++) {
		list[i].x = 1.0*rand()/RAND_MAX;
		list[i].y = 1.0*rand()/RAND_MAX;
		list[i].z = 1.0*rand()/RAND_MAX;
//cerr << list[i].x << "  " << list[i].y << "  " << list[i].z ;
		list[i].color[0] = (byte) q[i];
		list[i].color[1] = (byte) q[(i+1)%3];
		list[i].color[2] = (byte) q[(i+2)%3];
//cerr << " " << (int) list[i].color[0] << "  " << (int) list[i].color[1] << "  " << (int) list[i].color[2] << endl;
	}

}
*/
void FrameBuffer::Test()
{
	///srand( time(NULL) );
	if (1.0*rand()/RAND_MAX > .5)
		fb->Clear(red);
	else fb->Clear(blue);
}

/*
void Get_Triangles(ifstream &fin)
//Fills the triangle info array.
{
	double x, y, z;
	int num, i, temp[3];
	rgb color;

	
	fin >> num;
	num = 3*num;
	num_vert = num;
	//tl = new TriList(num);
	list = new vertex[num];	
	
	for (i = 1; i <= num; i++) {
		//fin.getline(line, sizeof(line));
		//get_xyzrgb(line, x, y, z, color);
		fin >> x >> y >> z;
//		cerr << x << " " << y << " " << z << " is x y z. ";
		fin >> temp[0] >> temp[1] >> temp[2];
		color[0] = (byte) temp[0]; color[1] = (byte) temp[1]; color[2] = (byte) temp[2];
//		cerr << (int) color[0] << " " << (int) color[1] << " " 
//			 << (int) color[2] << " are the r g b.\n";
		list[(i-1)].fill_vertex(x,y,z,color);
	}

}
*/

void myDisplay(void)
{
	// copy the framebuffer to the screen
	glDrawPixels(fb->width, fb->height, GL_RGB, GL_UNSIGNED_BYTE,  fb->pix);
	glFlush();
}

void myDraw()
{
	// I do my drawing into the framebuffer object here
	// this is just a hack example
	//fb->Clear(black);

	// then transfer to the screen
	// I can do this as often as I like, even after every pixel for debugging purposes.
	//for (int i = 0; i < num_vert; i += 3)
	//	fb->TriDraw(list, i);
	
	//fb->Test();
	//rgb p;
	//	p[0] = (byte) 255;
	//	p[1] = (byte) 0.04;
	//	p[2] = (byte) 255;
	//fb->Clear(p);
	
	/***Here's the overal raytracing algorithm.***/

	int i, j;
	double v[3], imp[3];
	//ray, center-of-proj, direction-of-view, focal-length, point on the image plane.
	double color[3];
	register byte *p;

	//File input parsing here will fill out the scene information.

	Parser_Test();


	//Clear the frame buffer to the background color.
	//rgb backcolor;
	//backcolor[0] = (byte) 255*background[0];
	//backcolor[1] = (byte) 255*background[1];
	//backcolor[2] = (byte) 255*background[2];
	//fb->Clear(backcolor);


	Pwfi = new Projection(dov, foc);

	
	for (i = 0; i < xres; i ++) {
		for (j = 0; j < yres; j ++) {
			v[0] = 1.0*i; v[1] = 1.0*j; v[2] = 1.0;
			Pwfi->Convert(v);
			imp[0]=cop[0]+v[0]; imp[1]=cop[1]+v[1]; imp[2]=cop[2]+v[2];
			color[0] = 0.0; color[1] = 0.0; color[2] = 0.0;
//			if ((i == 87)&&(j == 126)) {
//				int m = 0;
//			}
			if ((i > 256)&&(j > 256)) 
				int n = 0; 
			Single_Raytrace(v, imp, color, 1, 1.0);
			if ((byte) 255*color[0] > 10) {
				int m = 0;
			}
			
			//This does all the work with, using the global scene variables to
			//fill leave the correct color in color.
			
			p = fb->pix[j*xres + i];
			p[0] = (byte) 255*color[0]; p[1] = (byte) 255*color[1]; p[2] = (byte) 255*color[2];
		}
	}

	//printf("%i and %i are z and q\n",z,q);
	
	//for (i = 0; i < 3; i++) 
	//	printf("%lf\t%lf\t%lf\n", Pwfi->mat[i][0], Pwfi->mat[i][1], Pwfi->mat[i][2]);


	/***That's it though, with all the important stuff hidden.***/


	myDisplay();
}

void myInit()
{
	// initialize GL stuff for speed

//#if 1
	// this makes quite a difference on the PC with software GL
	glEnable(GL_DITHER);    
	glEnable(GL_BLEND);
//#endif

	// SGI claims some of these might produce speedups for some architectures
//#if 0
    glDisable(GL_ALPHA_TEST);
    glDisable(GL_DEPTH_TEST);
    glDisable(GL_FOG);
    glDisable(GL_LIGHTING);
    glDisable(GL_LOGIC_OP);
    glDisable(GL_STENCIL_TEST);
    glDisable(GL_TEXTURE_1D);
    glDisable(GL_TEXTURE_2D);
    glPixelTransferi(GL_MAP_COLOR, GL_FALSE);
    glPixelTransferi(GL_RED_SCALE, 1);
    glPixelTransferi(GL_RED_BIAS, 0);
    glPixelTransferi(GL_GREEN_SCALE, 1);
    glPixelTransferi(GL_GREEN_BIAS, 0);
    glPixelTransferi(GL_BLUE_SCALE, 1);
    glPixelTransferi(GL_BLUE_BIAS, 0);
    glPixelTransferi(GL_ALPHA_SCALE, 1);
    glPixelTransferi(GL_ALPHA_BIAS, 0);
	glShadeModel (GL_FLAT);
//#endif
	
}

void KeyInputHandler(unsigned char Key, int x, int y)
{
  switch(Key)
  {
    case 27: exit(0);
	default: 
			myDraw(); 
		break;
  };
}


void main(int argc, char* argv[])
{
	ifstream fin;
	char *filename = new char[80];

	fb = new FrameBuffer(xres, yres);

    glutInitDisplayMode(GLUT_RGB);
	glutInitWindowPosition(50, 50);
    glutInitWindowSize(xres, yres);


    glutCreateWindow("COMP 236 - ASSN2");

    glutKeyboardFunc(KeyInputHandler);

    glutDisplayFunc(myDisplay);

	myInit();

	printf("Type any character to compute an image.\n");

    glutMainLoop();
}