/* Simple C++ illustrating a 3D vector library with overloading 
   and 3D rotation about an arbitrary axis using Rodrigues's formula */

#include <math.h>
#include <GLUT/glut.h>	// GL/glut.h if not Mac

#define RAD 57.29577951
#define DegToRad(angle)	((angle)/RAD)

/* Minimal vector library, as needed here */

struct V3f	// declaration and initializatoin of the "V3f" structure
{
	float x, y, z;
	V3f(float x1, float y1, float z1)
		{x=x1; y=y1; z=z1;}
 	V3f()
 		{x=0; y=0; z=0;}
};

V3f operator+(V3f a, V3f b)			// vector addition
{
	return V3f(a.x+b.x, a.y+b.y, a.z+b.z);
}

V3f operator-(V3f a) 				// changing sign of a vector
{
	return V3f(-a.x, -a.y, -a.z);
}

V3f operator-(V3f a, V3f b)			// vectior subtraction
{
	return V3f(a.x-b.x, a.y-b.y, a.z-b.z);
}

V3f operator*(float c, V3f a)		// scalar multiplied by a vector
{
	return V3f(c*a.x, c*a.y, c*a.z);
}

V3f operator*(V3f a, float c)		// vector multiplied by a scalar
{
	return V3f(c*a.x, c*a.y, c*a.z);
}

float operator*(V3f a, V3f b)		// dot product
{
	return (a.x*b.x + a.y*b.y + a.z*b.z);
}

V3f operator%(V3f a, V3f b)			// cross product
{
	V3f c;
	c.x =   a.y * b.z - a.z * b.y;
	c.y = - a.x * b.z + a.z * b.x;
	c.z =   a.x * b.y - a.y * b.x;
	return c;
}

float VecLength(V3f v)				// vector length
{
	return sqrt (v*v);
}

V3f VecNormalize(V3f v)				// vector normalization
{
	return (1.0/VecLength(v)*v);
}

/* Define point to be rotated */

V3f P(0.8, 0.0, 0.0);
V3f O(0, 0, 0);

/* Define the axis of rotation n (arbitrary length) */

V3f axis(1,1,1);

/* Draw ccoordinate system and the axis of rotation */

void DrawAxes()
{
	glColor3f (0.0, 1.0, 0.0);
	glBegin(GL_LINES);	
	glVertex3f(0,0,0);
	glVertex3f(1,0,0);
	glVertex3f(0,0,0);
	glVertex3f(0,1,0);
	glVertex3f(0,0,0);
	glVertex3f(0,0,1);
	glEnd();

	V3f n = VecNormalize(axis);
	glColor3f (1.0, 1.0, 0.0);
	glBegin(GL_LINES);	
	glVertex3f(0,0,0);
	glVertex3f(n.x, n.y, n.z);
	glEnd();
}

/* Generate a sequence of point representing rotation of P around 
   the axis. This includes the essence of Rodrigues's algorithm */

void Generate()
{
	V3f n = VecNormalize(axis);
	V3f OP, OQ, QS, SP_prime, P_prime;

	glPointSize(5.0);
	glColor3f (1.0, 1.0, 1.0);
	glBegin(GL_POINTS);
	for (float phi=-180;phi<180;phi+=2) {
		OP = P-O;
		OQ = (OP*n)*n;
		QS = (OP-OQ)*cos(DegToRad(phi));
		SP_prime = (n%OP)*sin(DegToRad(phi));
		P_prime = O+(OQ+QS+SP_prime);
    	glVertex3f(P_prime.x, P_prime.y, P_prime.z);
	}
	glEnd();
	return;
}

void init (void) 						// initialize OpenGL
{
	glClearColor (0.0, 0.0, 0.0, 0.0);	// select clearing color

	glMatrixMode(GL_PROJECTION);		// set up viewing parameters
	glLoadIdentity();
 	gluPerspective(20, 1.0, 3.0, 7.0);
	gluLookAt(1, 1, 5, 0, 0.2, 0, 0, 1, 0);
 }

void display(void)						// display object
{
	glClear (GL_COLOR_BUFFER_BIT);		// cleaer all pixels
	DrawAxes();
	Generate();
	glFlush();
}

/* 
 * Declare initial window size, position, and display mode
 * (single buffer and RGBA).  Open window with program name
 * in its title bar.  Call initialization routines.
 * Register callback function to display graphics.
 * Enter main loop and process events.
 */
int main(int argc, char** argv)
{
   glutInit(&argc, argv);
   glutInitDisplayMode (GLUT_SINGLE | GLUT_RGB);
   glutInitWindowSize (750, 750); 
   glutInitWindowPosition (100, 100);
   glutCreateWindow (argv[1]);
   init ();
   glutDisplayFunc(display);
   glutMainLoop(); 
   return 0;   /* ANSI C requires main to return int. */
}