/*
 
 Assignment 2: Modeling Seashells
 
 CS451: Advanced Computer Graphics
 Jason Lawrence
 Spring 2008
 
 Acknowledgement: starter code derived from the Trimesh2 library
 written by Szymon Rusinkiewicz.
 
 Note: this program makes use of the FLTK user-interface library
 available from http://www.fltk.org.
 
 */

#ifdef WIN32
#include <windows.h>
#endif

#include <FL/Fl.H> 
#include <FL/Fl_Window.H> 
#include <FL/Fl_Value_Slider.H> 
#include <FL/Fl_Button.H>
#include <FL/fl_file_chooser.h>
#include <FL/glut.h>
#include <OpenGL/gl.h> 
#include <OpenGL/glu.h> 

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

#include "Vec.h"
#include "XForm.h"
#include "GLCamera.h"

#include <vector>
#include <algorithm>

using namespace std;

struct MyMesh {
	vector<point> vertices;
	vector< Vec<3,int> > faces;
};

struct FrenetFrame
{
	Vec<3,double> u;
	Vec<3,double> v;
	Vec<3,double> w;
};

struct CurvePoint
{
	double _tval;
	double _arclength;
	int _Q;
};

// Globals
MyMesh mesh;
GLCamera camera;
xform xf;

vector<raster> spline_control_points;
int curr_sel=-1;

bool draw_edges = false;
bool draw_cs = true;
bool draw_ff = false;

#define MODE_SHELL 0
#define MODE_CURVE 1

int app_mode = MODE_SHELL;

vec center(0,0,5);
double size = 10.0;

// Shell parameters
double tmin=0.0;
double tmax=18.0;
double ff_tval=0.1;
int q_val = 0;
double r_0=2.0;
double eps_r=1.04;
double z_0=0.5;
double eps_z=1.2;
double eps_c=1.0;
int resolution = 15;

// FLTK objects
Fl_Window       *fltkMainWin; 
Fl_Value_Slider *tminSlider;
Fl_Value_Slider *tmaxSlider;
Fl_Value_Slider *rnaughtSlider;
Fl_Value_Slider *epsrSlider;
Fl_Value_Slider *znaughtSlider;
Fl_Value_Slider *epszSlider;
Fl_Value_Slider *epscSlider;
Fl_Value_Slider *resolutionSlider;
Fl_Button       *loadButton;
Fl_Button       *makeButton;
Fl_Button       *clearButton;
Fl_Button       *saveButton;
Fl_Button       *saveMeshButton;

// Forward declarations
raster compute_curve_point(double t, int j);
double scale(double t);

// Set up lights and materials
void setup_lighting()
{
	glEnable(GL_LIGHTING);
	GLfloat mat_front[4] = { 1.0, 0.95, 0.9, 1 };
	GLfloat mat_back[4] = { 0.7, 0.6, 0.6, 1 };
	GLfloat mat_specular[4] = { 0.18, 0.18, 0.18, 0.18 };
	GLfloat mat_shininess[] = { 64 };
	GLfloat global_ambient[] = { 0.02, 0.02, 0.05, 0.05 };
	GLfloat light0_ambient[] = { 0, 0, 0, 0 };
	GLfloat light0_diffuse[] = { 0.85, 0.85, 0.8, 0.85 };
	GLfloat light0_specular[] = { 0.85, 0.85, 0.85, 0.85 };
	glMaterialfv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, mat_front);
	glMaterialfv(GL_BACK, GL_AMBIENT_AND_DIFFUSE, mat_back);
	glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, mat_specular);
	glMaterialfv(GL_FRONT_AND_BACK, GL_SHININESS, mat_shininess);
	glLightfv(GL_LIGHT0, GL_AMBIENT, light0_ambient);
	glLightfv(GL_LIGHT0, GL_DIFFUSE, light0_diffuse);
	glLightfv(GL_LIGHT0, GL_SPECULAR, light0_specular);
	glLightModelfv(GL_LIGHT_MODEL_AMBIENT, global_ambient);
	glLightModeli(GL_LIGHT_MODEL_LOCAL_VIEWER, GL_FALSE);
	glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, GL_TRUE);
	glEnable(GL_LIGHT0);
}

// Draw the mesh
void draw_mesh(MyMesh &mesh)
{
	glPushMatrix();
	glMultMatrixd(xf);
	
	glDepthFunc(GL_LESS);
	glEnable(GL_DEPTH_TEST);
	glDisable(GL_CULL_FACE);
	glShadeModel(GL_FLAT);
	glEnable(GL_NORMALIZE);
	glShadeModel(GL_SMOOTH);
	
	// Draw mesh - first pass
	if (draw_edges) {
		glPolygonOffset(1, 1);
		glEnable(GL_POLYGON_OFFSET_FILL);
	}
	
	/*
	glPointSize(1.0f);
	glBegin(GL_POINTS);
	for (int i = 0; i < mesh.vertices.size(); i++)
	{
		glVertex3fv(mesh.vertices[i]);
	}*/
	glEnd();
	///*
	glBegin(GL_TRIANGLES);
	int nf = mesh.faces.size();
	for (int i = 0; i < nf; i++) {
		const point &v0 = mesh.vertices[mesh.faces[i][0]];
		const point &v1 = mesh.vertices[mesh.faces[i][1]];
		const point &v2 = mesh.vertices[mesh.faces[i][2]];
		
		glNormal3fv((v1 - v0) CROSS (v2 - v0));  // Needs GL_NORMALIZE
		glVertex3fv(v0);
		glVertex3fv(v1);
		glVertex3fv(v2);
	}
	glEnd();
	
	// Draw mesh - second pass to draw the edges on top
	if (draw_edges) {
		glDisable(GL_POLYGON_OFFSET_FILL);
		glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
		glLineWidth(1.0);
		GLfloat mat_diffuse[4] = { 0.0f, 0.0f, 1.0f, 1.0f };
		glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, mat_diffuse);
		glBegin(GL_TRIANGLES);
		for (int i = 0; i < nf; i++) {
			const point &v0 = mesh.vertices[mesh.faces[i][0]];
			const point &v1 = mesh.vertices[mesh.faces[i][1]];
			const point &v2 = mesh.vertices[mesh.faces[i][2]];
			
			glNormal3fv((v1 - v0) CROSS (v2 - v0));
			glVertex3fv(v0);
			glVertex3fv(v1);
			glVertex3fv(v2);
		}
		glEnd();
		glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
	}
	
	glPopMatrix();
	glDisable(GL_DEPTH_TEST);
	glDisable(GL_LIGHTING);
	// */
}

// Draw the coordinate system in which the helico spiral is defined
void draw_coordinate_system()
{
	glPushMatrix();
	glMultMatrixd(xf);
	
	glLineWidth(3);
	glColor3f(1,0,0);
	glBegin(GL_LINES);
	glVertex3f(0,0,0); glVertex3f(10,0,0);
	glEnd();
	glColor3f(0,1,0);
	glBegin(GL_LINES);
	glVertex3f(0,0,0); glVertex3f(0,10,0);
	glEnd();
	glColor3f(0,0,1);
	glBegin(GL_LINES);
	glVertex3f(0,0,0); glVertex3f(0,0,10);
	glEnd();
	
	glPopMatrix();
}

// Compute the Frenet Frame
FrenetFrame compute_ff(double t)
{
	FrenetFrame f;
	
	if (t == 0.0)
		t=1.0;
	
	// Compute differentials
	double dr = pow(eps_r, t) * r_0;
	
	double dx = dr * (cos(t) * log(eps_r) - sin(t));
	double dy = dr * (sin(t) * log(eps_r) + cos(t));;
	double dz = log(eps_z) * z_0 * pow(eps_z, t);
	
	double ler = log(eps_r);
	double ddx = dr * (cos(t) * ler * ler - 2 * ler * sin(t) - cos(t));
	double ddy = dr * (sin(t) * ler * ler + 2 * ler * cos(t) - sin(t));
	double ddz = log(eps_z) * log(eps_z) * z_0 * pow(eps_z, t);
	
	// Construct H' and H"
	Vec<3,double> dH(dx, dy, dz);
	Vec<3,double> ddH(ddx, ddy, ddz);
	
	// Compute w, v, and u
	f.w = dH; normalize(f.w);
	f.v = f.w CROSS ddH; normalize(f.v);
	f.u = f.v CROSS f.w;

	return f;
}

void draw_frenet_frame()
{
	FrenetFrame F = compute_ff(ff_tval);
	
	glPushMatrix();
	glMultMatrixd(xf);
	
	// Find the coordinate on the helico spiral
	double t = ff_tval;
	double r = r_0*pow(eps_r,t);
	double z = z_0*pow(eps_z,t);
	double x = r*cos(t);
	double y = r*sin(t);	
	
	//printf("w: (%f,%f,%f)\n", F.w[0],F.w[1],F.w[2]);
	//printf("v: (%f,%f,%f)\n", F.v[0],F.v[1],F.v[2]);
	//printf("u: (%f,%f,%f)\n", F.u[0],F.u[1],F.u[2]);
	
	glLineWidth(3);
	// Red
	glColor3f(1,0,0);
	glBegin(GL_LINES);
		glVertex3f(x,y,z); glVertex3f(x+F.w[0],y+F.w[1],z+F.w[2]);
	glEnd();
	// Green
	glColor3f(0,1,0);
	glBegin(GL_LINES);
		glVertex3f(x,y,z); glVertex3f(x+F.v[0],y+F.v[1],z+F.v[2]);
	glEnd();
	// Blue
	glColor3f(0,0,1);
	glBegin(GL_LINES);
		glVertex3f(x,y,z); glVertex3f(x+F.u[0],y+F.u[1],z+F.u[2]);
	glEnd();
	
	vector<vec> spline_points;
	for (int i = 0; i < spline_control_points.size(); i++)
	{
		for (double t = 0.0f; t <= 1.0f; t += 0.1)
		{
			// Grab the current B-spline curve
			raster temp = compute_curve_point(t, i);
			// Scale by eps_c^t
			temp = temp * (float) scale(ff_tval);
			//printf("eps_c: %f\nff_tval: %f\neps_c^t: %f\n\n", eps_c, ff_tval, pow(eps_c, ff_tval));
			
			vec temp2;
			// Project onto Frenet Frame
			temp2[0] = x + temp[0] * F.u[0] + temp[1] * F.v[0];
			temp2[1] = y + temp[0] * F.u[1] + temp[1] * F.v[1];
			temp2[2] = z + temp[0] * F.u[2] + temp[1] * F.v[2];
			spline_points.push_back(temp2);
		}
	}
	
	// Draw the spline at the current helico spiral point
	glPointSize(1);
	glColor3f(0.0f, 1.0f, 0.0f);
	glBegin(GL_LINE_LOOP);
		for (int i=0; i<(int)spline_points.size(); i++)
		{
			glVertex3f(spline_points[i][0],spline_points[i][1],spline_points[i][2]);
		}
	glEnd();
	
	glColor3f(0.0f, 0.0f, 1.0f);
	
	glPopMatrix();
}

// Draw the helico spiral
void draw_helico_spiral()
{
	glPushMatrix();
	glMultMatrixd(xf);
	
	glLineWidth(1);
	glBegin(GL_LINE_STRIP);
		int steps = 100;
		for (int i=0; i<=steps; i++)
		{
			double t = tmin + (tmax-tmin)*(double(i)/double(steps));
			double r = r_0*pow(eps_r,t);
			double z = z_0*pow(eps_z,t);
			double x = r*cos(t);
			double y = r*sin(t);
			glVertex3d(x,y,z);
		}
	glEnd();
	
	glPopMatrix();
}

double compute_curve_length(double t_val = ff_tval, vector<CurvePoint>* f_s = NULL, FrenetFrame F = compute_ff(ff_tval))
{
	if (f_s != NULL)
		(*f_s).clear();
	
	double length = 0.0;
	
	double t = t_val;
	double delta_t = 0.1;
	double r = r_0*pow(eps_r,t);
	double z = z_0*pow(eps_z,t);
	double x = r*cos(t);
	double y = r*sin(t);	
	
	vec CO;
	vec Cm;
	for (int i = 0; i <= spline_control_points.size(); i++)
	{
		for (double t = 0.0f; t <= 1.0f; t += delta_t)
		{
			vec C;
			
			if (i == spline_control_points.size())
			{
				C = CO;
			} else
			{
				// Grab the current B-spline curve
				raster temp = compute_curve_point(t, i);
				// Scale by eps_c^t
				temp = temp * (float) scale(ff_tval);
				
				// Project onto Frenet Frame
				C[0] = x + temp[0] * F.u[0] + temp[1] * F.v[0];
				C[1] = y + temp[0] * F.u[1] + temp[1] * F.v[1];
				C[2] = z + temp[0] * F.u[2] + temp[1] * F.v[2];
			}
			
			if (i == 0)
			{
				Cm = C;
				CO = C;
				continue;
			}
			
			// Compute du/ds = (f[s+h][u] - f[s][u]) / h
			double du = (C[0] - Cm[0]) / delta_t;
			// Compute dv/ds = (f[s+h][v] - f[s][v]) / h
			double dv = (C[1] - Cm[1]) / delta_t;
			// Compute dw/ds = (f[s+h][w] - f[s][w]) / h
			double dw = (C[2] - Cm[2]) / delta_t;
			
			double term = sqrt(du*du + dv*dv + dw*dw);
			
			length += term * delta_t;
			Cm = C;
			
			// Save curve information if requested
			if (f_s != NULL)
			{
				CurvePoint P;
				P._tval = t;
				P._Q = i-1;
				P._arclength = length;
				(*f_s).push_back(P);
			}
		}
	}

	return length;
}

CurvePoint find_spline_parameter(vector<CurvePoint> f_s, double arclength)
{
	CurvePoint P;
	for (int i = 0; i < f_s.size(); i++)
	{
		if (f_s[i]._arclength >= arclength)
		{
			P = f_s[i];
			break;
		}
	}
	
	return P;
}

double scale(double t)
{
	return pow(eps_c, t)  * (t + 0.1) * eps_c;// * (1.28 + 0.15 * pow(sin(t) + cos(t), 2));
}

void construct_mesh()
{
	// Loop over the helico spiral
	for (double t = tmin; t <= tmax; t+= 0.1)
	{
		// Compute Frenet Frame
		FrenetFrame F = compute_ff(t);
		// Grab the arc length function
		vector<CurvePoint> f_s;
		double L = compute_curve_length(t, &f_s, F);
		
		// Number of integration steps
		int n = resolution;
		// Step length
		double ds = L / (double) n;
		double smax = (double) spline_control_points.size();
		double smin = 0.0;
		
		// Find the parameter values
		vector<CurvePoint> p;
		for (double s = smin; s <= smax; s+= ds)
		{
			CurvePoint P = find_spline_parameter(f_s, s);
			p.push_back(P);
		}
		
		// Save the mesh
		double r = r_0*pow(eps_r,t);
		double z = z_0*pow(eps_z,t);
		double x = r*cos(t);
		double y = r*sin(t);
		
		for (int i = 0; i <= p.size(); i++)
		{
			int back_vertex;
			int back_vertex_back_rim;
			int back_rim;
			int vertex;
			
			// Really we do the first one last
			if ( i == p.size() )
			{
				back_vertex =			mesh.vertices.size() - 1;
				back_vertex_back_rim =	mesh.vertices.size() - p.size() - 1;
				back_rim =				mesh.vertices.size() - 2*p.size();
				vertex =				mesh.vertices.size() - p.size();
			} else
			{
				// Grab the current B-spline curve
				raster temp = compute_curve_point(p[i]._tval, p[i]._Q);
				// Scale by eps_c^t
				temp = temp * (float) scale(t);
				
				vec temp2;
				// Project onto Frenet Frame
				temp2[0] = x + temp[0] * F.u[0] + temp[1] * F.v[0];
				temp2[1] = y + temp[0] * F.u[1] + temp[1] * F.v[1];
				temp2[2] = z + temp[0] * F.u[2] + temp[1] * F.v[2];
				mesh.vertices.push_back(temp2);
				
				// Construct faces
				back_vertex =			mesh.vertices.size() - 2;
				back_vertex_back_rim =	mesh.vertices.size() - p.size() - 2;
				back_rim =				mesh.vertices.size() - p.size() - 1;
				vertex =				mesh.vertices.size() - 1;
			}
			
			// Some boundary cases
			if (t == tmin || i == 0)
				continue;
			
			// Two faces are created
			
			Vec<3,int> face1;
			Vec<3,int> face2;
			
			face1[0] = vertex;
			face1[1] = back_vertex;
			face1[2] = back_vertex_back_rim;
			face2[0] = vertex;
			face2[1] = back_vertex_back_rim;
			face2[2] = back_rim;
			
			mesh.faces.push_back(face1);
			mesh.faces.push_back(face2);
		}
	}
}

raster compute_curve_point(double t, int j)
{
	raster Q;
	
	// Create point indices
	int p0 = j - 3;
	int p1 = j - 2;
	int p2 = j - 1;
	int p3 = j;
	
	// Determine whether we need to loop around spline points
	if (j < 3)
	{
		p0 = j - 3 < 0 ? spline_control_points.size() + j - 3 : j - 3;
		p1 = j - 2 < 0 ? spline_control_points.size() + j - 2 : j - 2;
		p2 = j - 1 < 0 ? spline_control_points.size() + j - 1 : j - 1;
	}
	
	// Extract points
	raster P0 = spline_control_points[p0];
	raster P1 = spline_control_points[p1];
	raster P2 = spline_control_points[p2];
	raster P3 = spline_control_points[p3];
	
	// Compute the spline point 
	Q =		P0 * (float) ( -1.0 / 6.0 * t * t * t + 0.5 * t * t - 0.5 * t + 1.0 / 6.0) +
			P1 * (float) ( 0.5 * t * t * t - t * t + 2.0 / 3.0) +
			P2 * (float) ( -0.5 * t * t * t + 0.5 * t * t + 0.5 * t + 1.0 / 6.0) +
			P3 * (float) ( 1.0 / 6.0 * t * t * t);
		
	return Q;
}

// Draw the spline that models the shell's "generating curve"
void draw_spline()
{
	// Draw the control points
	glPointSize(5);
	glBegin(GL_POINTS);
	for (int i=0; i<(int)spline_control_points.size(); i++)
    {
		glVertex2f(spline_control_points[i][0],spline_control_points[i][1]);
    }
	glEnd();
	
	vector<raster> spline_points;
	// Insert your code here for drawing the B-spline curve corresponding to
	// these control points.
	for (int i = 0; i < spline_control_points.size(); i++)
	{
		for (double t = 0.0f; t <= 1.0f; t += 0.1)
		{
			spline_points.push_back(compute_curve_point(t, i));
		}
	}
	
	// Draw the spline
	glPointSize(1);
	glColor3f(0.0f, 1.0f, 0.0f);
	glBegin(GL_LINE_LOOP);
	for (int i=0; i<(int)spline_points.size(); i++)
    {
		glVertex2f(spline_points[i][0],spline_points[i][1]);
    }
	glEnd();
	
	glColor3f(0.0f, 0.0f, 1.0f);
}

// Draw the complete scene
void redraw()
{
	switch (app_mode)
    {
		case MODE_SHELL:
		{
			camera.setupGL(xf * center, size);
			glClearColor(1, 1, 1, 0);
			glClearDepth(1);
			glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
			if (draw_cs)
				draw_coordinate_system();
			if (draw_ff)
				draw_frenet_frame();
			draw_helico_spiral();
			setup_lighting();
			draw_mesh(mesh);
			glutSwapBuffers();
		}
			break;
		case MODE_CURVE:
		{
			glMatrixMode(GL_MODELVIEW);
			glLoadIdentity();
			glMatrixMode(GL_PROJECTION);
			glLoadIdentity();
			glOrtho(-2,2,-2,2,-1,1);
			glClearColor(1, 1, 1, 0);
			glClearDepth(1);
			glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
			draw_spline();
			glutSwapBuffers();
		}
			break;
		default:
			break;
    }
}

// Handle mouse button and motion events
static unsigned buttonstate = 0;

// Reset the view...
void resetview()
{
	xf = xform::trans(0, 0, -5.0f * size) * xform::trans(-center);
	camera.stopspin();
	buttonstate = 0;
}

void mousemotionfunc(int x, int y)
{
	static const Mouse::button physical_to_logical_map[] = {
		Mouse::NONE, Mouse::ROTATE, Mouse::MOVEXY, Mouse::MOVEZ,
		Mouse::MOVEZ, Mouse::MOVEXY, Mouse::MOVEXY, Mouse::MOVEXY,
	};
	
	switch (app_mode)
    {
		case MODE_SHELL:
		{
			Mouse::button b = Mouse::NONE;
			if (buttonstate & (1 << 4))
				b = Mouse::WHEELUP;
			else if (buttonstate & (1 << 3))
				b = Mouse::WHEELDOWN;
			else
				b = physical_to_logical_map[buttonstate & 7];
			
			camera.mouse(x, y, b, xf * center, size, xf);
			if (b != Mouse::NONE)
				glutPostRedisplay();
		}
			break;
		case MODE_CURVE:
		{
			if (curr_sel >= 0)
			{
				spline_control_points[curr_sel] = Vec<2,float>(float(x-256)/128.0f,float((512-1-y)-256)/128.0f);
				glutPostRedisplay();
			}
		}
			break;
		default:
			break;
    }
}

void mousebuttonfunc(int button, int state, int x, int y)
{
	switch (app_mode)
    {
		case MODE_SHELL:
		{
			if (state == GLUT_DOWN)
				buttonstate |= (1 << button);
			else
				buttonstate = 0;//buttonstate &= ~(1 << button);
			
			mousemotionfunc(x, y);
		}
			break;
		case MODE_CURVE:
		{
			
			if (state == GLUT_DOWN && button == GLUT_LEFT_BUTTON)
			{
				
				GLdouble modelview[16];
				glGetDoublev(GL_MODELVIEW_MATRIX, modelview);
				GLdouble projection[16];
				glGetDoublev(GL_PROJECTION_MATRIX, projection);
				GLint viewport[4];
				glGetIntegerv(GL_VIEWPORT, viewport);
				
				// Pick a point on the spline
				double minr = 100.0;
				for (int i=0; i<(int)spline_control_points.size(); i++)
				{
					GLdouble wx,wy,wz;
					gluProject(spline_control_points[i][0],spline_control_points[i][1],0.0,
							   modelview, projection, viewport,
							   &wx,&wy,&wz);
					double _r2 = (wx-x)*(wx-x)+(wy-(512-1-y))*(wy-(512-1-y));
					if (_r2 < 8.0 && _r2 < minr)
					{
						curr_sel = i;
						minr = _r2;
					}
				}
			}
			if (state == GLUT_UP && button == GLUT_LEFT_BUTTON)
			{
				curr_sel = -1;
			}
		}
			break;
		default:
			break;
    }
}

// File I/O routines

void load_from_file (char *filename)
{
	FILE *fid = fopen(filename,"rt");
	fscanf(fid,"%lg ", &tmin);
	fscanf(fid,"%lg ", &tmax);
	fscanf(fid,"%lg ", &r_0);
	fscanf(fid,"%lg ", &eps_r);
	fscanf(fid,"%lg ", &z_0);
	fscanf(fid,"%lg ", &eps_z);
	int cpts;
	fscanf(fid,"%d  ", &cpts);
	spline_control_points.clear();
	for (int i=0; i<cpts; i++)
    {
		float x,y;
		fscanf(fid,"%f %f",&x,&y);
		spline_control_points.push_back(Vec<2,float>(x,y));
    }
	fclose(fid);
}

void save_to_file (char *filename)
{
	FILE *fid = fopen(filename,"wt");
	fprintf(fid,"%f\n", tmin);
	fprintf(fid,"%f\n", tmax);
	fprintf(fid,"%f\n", r_0);
	fprintf(fid,"%f\n", eps_r);
	fprintf(fid,"%f\n", z_0);
	fprintf(fid,"%f\n", eps_z);
	fprintf(fid,"%d\n", spline_control_points.size());
	for (int i=0; i<(int)spline_control_points.size(); i++)
		fprintf(fid,"%f %f\n",spline_control_points[i][0],spline_control_points[i][1]);
	fclose(fid);
}

void save_mesh_to_file (char *filename)
{
	FILE *fid = fopen(filename,"wt");
	fprintf(fid,"OFF\n%d %d\n", mesh.vertices.size(), mesh.faces.size());
	for (int i=0; i<(int)mesh.vertices.size(); i++)
		fprintf(fid,"%f %f %f\n", mesh.vertices[i][0], mesh.vertices[i][1], mesh.vertices[i][2]);
	for (int i=0; i<(int)mesh.faces.size(); i++)
		fprintf(fid,"3 %d %d %d\n", mesh.faces[i][0], mesh.faces[i][1], mesh.faces[i][2]);
	fclose(fid);
}

// Handle keyboard events
void keyboardfunc(unsigned char key, int x, int y)
{
	switch (key) {
		case ' ':
			resetview(); break;
		case 'c':
			draw_cs = !draw_cs; break;
		case 'e':
			draw_edges = !draw_edges; break;
		case 'f':
			draw_ff = !draw_ff; break;
		case 'm':
		case 'M':
			app_mode = (app_mode + 1) % 2; break;
		case '\033': // Esc
		case 'Q':
		case 'q':
			exit(0);
		case '=':
			ff_tval = ff_tval >= tmax ? ff_tval : ff_tval + 0.1; 
			break;
		case '-':
			ff_tval = ff_tval <= 0.1 ? ff_tval : ff_tval - 0.1; 
			break;
		case ']':
			q_val += q_val == spline_control_points.size()-1 ? 0 : 1; break;
		case '[':
			q_val -= q_val == 0 ? 0 : 1; break;
	}
	glutPostRedisplay();
}

// Callback functions for GUI support

void tminSliderCB(Fl_Value_Slider *pslider, long dum)
{
	tmin = pslider->value();
	glutPostRedisplay();
} 

void tmaxSliderCB(Fl_Value_Slider *pslider, long dum)
{
	tmax = pslider->value();
	glutPostRedisplay();
} 

void rnaughtSliderCB(Fl_Value_Slider *pslider, long dum)
{
	r_0 = pslider->value();
	glutPostRedisplay();
} 

void epsrSliderCB(Fl_Value_Slider *pslider, long dum)
{
	eps_r = pslider->value();
	glutPostRedisplay();
} 

void znaughtSliderCB(Fl_Value_Slider *pslider, long dum)
{
	z_0 = pslider->value();
	glutPostRedisplay();
} 

void epszSliderCB(Fl_Value_Slider *pslider, long dum)
{
	eps_z = pslider->value();
	glutPostRedisplay();
} 

void epscSliderCB(Fl_Value_Slider *pslider, long dum)
{
	eps_c = pslider->value();
	glutPostRedisplay();
} 

void resolutionSliderCB(Fl_Value_Slider *pslider, long dum)
{
	resolution = pslider->value();
	glutPostRedisplay();
}

void loadButtonCB(Fl_Button *button,void *)
{       
	char * filename = fl_file_chooser("Select filename","*.shl",0);
	if (filename)
		load_from_file(filename);
}

// Insert code here!
void makeButtonCB(Fl_Button *button,void *)
{   
    mesh.vertices.clear();
	mesh.faces.clear();
	construct_mesh();
	glutPostRedisplay();
}
void clearButtonCB(Fl_Button *button,void *)
{
	mesh.vertices.clear();
	mesh.faces.clear();
	glutPostRedisplay();
}

void saveButtonCB(Fl_Button *button,void *)
{       
	char * filename = fl_file_chooser("Select filename","*.shl",0);
	if (filename)
		save_to_file(filename);
}

void saveMeshButtonCB(Fl_Button *button,void *)
{       
	char * filename = fl_file_chooser("Select filename","*.off",0);
	if (filename)
		save_mesh_to_file(filename);
}

void usage(const char *myname)
{
	fprintf(stderr, "usage: %s\n", myname);
	exit(1);
}

int main(int argc, char *argv[])
{
	// Initialize spline points
	for (int i=0; i<10; i++)
    {
		double th = double(i)/double(10)*3.14159*2.0;
		spline_control_points.push_back(Vec<2,float>(cos(th),sin(th)));
    }
		
	// Create and setup main FLTK window
	fltkMainWin = new Fl_Window(100+512, 512, "Seashell Modeler");
	
	// Create a slider control for each of the set of parameters that
	// define the shell's surface
	
	tminSlider = new Fl_Value_Slider(0,0,100,20,"tmin");
	tminSlider->type(FL_HOR_NICE_SLIDER); 
	tminSlider->minimum(0.0); 
	tminSlider->maximum(6.0); 
	tminSlider->step(0.1);
	tminSlider->labelsize(12); 
	tminSlider->value(tmin);
	tminSlider->callback((Fl_Callback*)tminSliderCB,(void*)(0)); 
	
	tmaxSlider = new Fl_Value_Slider(0,40,100,20,"tmax");
	tmaxSlider->type(FL_HOR_NICE_SLIDER); 
	tmaxSlider->minimum(0.0); 
	tmaxSlider->maximum(30.0); 
	tmaxSlider->step(0.1);
	tmaxSlider->labelsize(12); 
	tmaxSlider->value(tmax);
	tmaxSlider->callback((Fl_Callback*)tmaxSliderCB,(void*)(0)); 
	
	rnaughtSlider = new Fl_Value_Slider(0,80,100,20,"r_0");
	rnaughtSlider->type(FL_HOR_NICE_SLIDER); 
	rnaughtSlider->minimum(0.0); 
	rnaughtSlider->maximum(10.0);
	rnaughtSlider->step(0.1);
	rnaughtSlider->labelsize(12); 
	rnaughtSlider->value(r_0);
	rnaughtSlider->callback((Fl_Callback*)rnaughtSliderCB,(void*)(0)); 
	
	epsrSlider = new Fl_Value_Slider(0,120,100,20,"eps_r");
	epsrSlider->type(FL_HOR_NICE_SLIDER); 
	epsrSlider->minimum(0.0); 
	epsrSlider->maximum(2.0);
	epsrSlider->step(0.001);
	epsrSlider->labelsize(12); 
	epsrSlider->value(eps_r);
	epsrSlider->callback((Fl_Callback*)epsrSliderCB,(void*)(0)); 
	
	znaughtSlider = new Fl_Value_Slider(0,160,100,20,"z_0");
	znaughtSlider->type(FL_HOR_NICE_SLIDER); 
	znaughtSlider->minimum(0.0); 
	znaughtSlider->maximum(10.0);
	znaughtSlider->step(0.1);
	znaughtSlider->labelsize(12); 
	znaughtSlider->value(z_0);
	znaughtSlider->callback((Fl_Callback*)znaughtSliderCB,(void*)(0)); 
	
	epszSlider = new Fl_Value_Slider(0,200,100,20,"eps_z");
	epszSlider->type(FL_HOR_NICE_SLIDER); 
	epszSlider->minimum(0.0); 
	epszSlider->maximum(2.0);
	epszSlider->step(0.001);
	epszSlider->labelsize(12); 
	epszSlider->value(eps_z);
	epszSlider->callback((Fl_Callback*)epszSliderCB,(void*)(0)); 
	
	epscSlider = new Fl_Value_Slider(0,240,100,20,"eps_c");
	epscSlider->type(FL_HOR_NICE_SLIDER); 
	epscSlider->minimum(0.0);
	epscSlider->maximum(3.0);
	epscSlider->step(0.01);
	epscSlider->labelsize(12); 
	epscSlider->value(eps_c);
	epscSlider->callback((Fl_Callback*)epscSliderCB,(void*)(0)); 
	
	resolutionSlider = new Fl_Value_Slider(0,280,100,20,"resolution");
	resolutionSlider->type(FL_HOR_NICE_SLIDER); 
	resolutionSlider->minimum(10);
	resolutionSlider->maximum(70);
	resolutionSlider->step(1);
	resolutionSlider->labelsize(12); 
	resolutionSlider->value(resolution);
	resolutionSlider->callback((Fl_Callback*)resolutionSliderCB,(void*)(0));
	
	// Create and setup buttons for file I/O
	
	makeButton = new Fl_Button(0,360,100,20,"Make Mesh");
	makeButton->labelsize(12);
	makeButton->callback((Fl_Callback*)makeButtonCB);
	
	clearButton = new Fl_Button(0,380,100,20,"Clear Mesh");
	clearButton->labelsize(12);
	clearButton->callback((Fl_Callback*)clearButtonCB);
	
	loadButton = new Fl_Button(0,420,100,20,"Open...");
	loadButton->labelsize(12);
	loadButton->callback((Fl_Callback*)loadButtonCB);
	
	saveButton = new Fl_Button(0,440,100,20,"Save...");
	saveButton->labelsize(12);
	saveButton->callback((Fl_Callback*)saveButtonCB);
	
	saveMeshButton = new Fl_Button(0,460,100,20,"Save Mesh...");
	saveMeshButton->labelsize(12);
	saveMeshButton->callback((Fl_Callback*)saveMeshButtonCB);
	
	fltkMainWin->show();
	
	// Create and setup GLUT pane within main window
	fltkMainWin->begin();
	glutInitWindowPosition(100, 0);
	glutInitWindowSize(512,512);
	glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH);
	glutCreateWindow("");
	fltkMainWin->end();
	
	// Register callback functions for drawing, mouse, and keyboard events
	glutDisplayFunc(redraw);
	glutMouseFunc(mousebuttonfunc);
	glutMotionFunc(mousemotionfunc);
	glutKeyboardFunc(keyboardfunc);
	
	// Reset camera parameters
	resetview();
	
	// Enter event loop
	glutMainLoop();
}