/*
helitoroid.cpp   --

This is the C++ code that generates and draws the two surfaces in the
QuantComp Co logo, well-commented with a story about the geometry  --

The surfaces are the torus and the helitoroid wrapped around it.

More, this program allows interactive control of 3D position and orientation by keyboard,
(when compiled and run on a local workstation with support of the relevant OpenGL libraries).

"Helitoroid" is the name I have coined for this surface.
It's  a surface in the form of a closed tube whose central axis is
a closed curve that spirals around a directrix circle, at constant distance
with constant integer pitch, so closing on itself.

I like this surface because its construction uses the Frenet-Serret theory of space curves,
a very basic body of mathematical knowledge that's at the foundation of
classical differential geometry,
a very lovely subject.

The torus is the familiar donut, characterized by its major and minor radii.

The classes Torus and Helitoroid know how to compute (tessellations of) the respective surfaces in standard
position, with the directing circle in the xy-plane, centered at the origin.  The OpenGL modelview
transformation serves to place it anywhere, in any orientation.

The syntax highlighting in the html rendition of this file is produced by the Komodo editor.

Copyright 2007  Ronald Levine,  QuantComp Co.
*/

#include <cmath>
#include <vector>
#include <boost/shared_ptr.hpp>
#include <boost/shared_array.hpp>
#include <GL/glut.h>                    // GLUT
#include <GL/glu.h>                     // GLU
#include <GL/gl.h>                          // OpenGL

using namespace std;
using namespace boost;
static const double PI = 3.1415926535897931;
struct Vector3;  //forward decl for the benefit of Point3::displace()
struct Point3 {
double x, y, z;
Point3(){x = y = z = 0.;}
double *getArray(){return &x;}
void displace(const Vector3 & b);
};

struct Vector3 {
double x, y, z;
Vector3(){ x = y = z = 0.0; };
Vector3(double u, double v, double w):x(u), y(v), z(w){};
Vector3(const Point3 & u, const Point3 & v);
Vector3 operator+(Vector3 rhs){return Vector3(this->x + rhs.x, this->y + rhs.y, this->z + rhs.z);};
void normalize();
void scalarMultiply(double);
void subtractVector(const Vector3 &);
double *getArray(){return &x;};
};

//displace a point by a vector
void Point3::displace(const Vector3 & d){ x += d.x; y += d.y; z += d.z;}

//construct vector as difference between two points.
Vector3::Vector3(const Point3 & u, const Point3 & v){x = v.x -u.x; y = v.y-u.y; z=v.z-u.z;}

x += b.x;
y += b.y;
z += b.z;
}

void Vector3::subtractVector(const Vector3 & b){
x -= b.x;
y -= b.y;
z -= b.z;
}

void Vector3::scalarMultiply(double r){
x *= r;
y *= r;
z *= r;
}

void Vector3::normalize(){
this->scalarMultiply( 1.0/ sqrt(x*x+y*y+z*z));
}

Vector3  crossProduct(const Vector3 &a, const Vector3 &b){
return Vector3(
a.y*b.z - a.z*b.y,
a.z*b.x - a.x*b.z,
a.x*b.y - a.y*b.x
);
}

double dotProduct(const Vector3 & a, const Vector3 & b){
return a.x*b.x + a.y*b.y +a.z*b.z;
}

typedef   shared_array<int>   facetype;

class Polyhedron {  //An abstract class
public:
Polyhedron(){};
virtual void draw() = 0;
virtual ~Polyhedron(){};
protected:
shared_array<Point3>        vertices_;
shared_array<facetype>    faces_;
shared_array<Vector3>     faceNormals_;
shared_array<Vector3>     vertexNormals_;
};

class Quadmesh : public Polyhedron {
public:
virtual void draw();
protected:
int segments_;
int sectors_;
};

class Torus : public Quadmesh {
public:
Torus(){};
Torus(double, double, int, int);
private:
};

class Helitoroid : public Quadmesh {
public:
Helitoroid(){};
Helitoroid(double, double, double, int, int, int);
private:
int pitch_;
};

/*
Makes a quadmesh tessellation of a torus
*/
sectors_ = sectors;
segments_ = segments;
int numPoints = segments_*sectors_;
vertices_  = shared_array<Point3>(new Point3[numPoints]);
vertexNormals_ = shared_array<Vector3>(new Vector3[numPoints]);
int numFaces = numPoints;  //for closed quad mesh
faces_ = shared_array<facetype> (new facetype[numFaces]);

for(int i=0; i< segments_; i++){
double phi = i*2.0*PI/segments_;
Vector3  nazvect(cos(phi), sin(phi) , 0.0);
Vector3 azvect = nazvect;
Point3  azpoint;
azpoint.displace(azvect);
for(int j=0; j<sectors_;j++){
int k = i*sectors_ + j;

double theta = j*2.0*PI/sectors_;
Vector3 norm = nazvect;
norm.scalarMultiply(cos(theta));
norm.z = sin(theta);
vertexNormals_[k] = norm;
vertices_[k] = azpoint;
vertices_[k].displace(norm);
//Now we need to set up the face index lists.
faces_[k] =  shared_array<int>(new int);
faces_[k] = k;
faces_[k] = (j<sectors_-1    ?  k+1                           : i*sectors_ );
faces_[k] = (i<segments_-1 ? faces_[k]+sectors_ : (j<sectors_-1 ? j+1 : 0));
faces_[k] = (i<segments_-1 ?  k+sectors_                  : j);
}
}
}

/***
Makes a quadmesh polyhedron in the form of a helical toroid.
That is, a surface in the form of a closed tube whose central axis is
a closed curve that spirals around a directrix circle, at constant distance
with constant pitch.

The directrix circle is in the xy-plane and centered on the origin

minrad_ is the constant distance from the directrix circle to the
spiraling axis curve of the tube.
pitch_ is the the number of times the tube winds around the directrix circle.

segments_ and sectors_ give the fineness of the tessellation of the tube,
segments_  along the axis of the tube, and
sectors_  around its circumference.
***/

Helitoroid::Helitoroid(
int segments, int sectors, int pitch):
pitch_( pitch) {

segments_ = segments;
sectors_  = sectors;
int numPoints = segments*sectors;

vertices_  = shared_array<Point3>(new Point3[numPoints]);
vertexNormals_ = shared_array<Vector3>(new Vector3[numPoints]);
int numFaces = numPoints;  //for closed quad mesh
faces_ = shared_array<facetype> (new facetype[numFaces]);

Point3  axis;              // general point on spiral axis of helitoriod

//angle variables
double phi;                // angle of arc on directrix circle, from x-axis
double theta;             // angle between the vector from point on directrix
//   to point on spiral axis and xy-plane,
//    equal to pitch*phi
double psi;                // angle on the circular intersection of tube
//with a plane perpendicular to spiral axis,
//measured from curve normal vector

double       dphi = 2.0*PI/segments;/* constant increments for angle variables */
double        dtheta = pitch*dphi;
double         dpsi = -2.0*PI/sectors;

double  cosphi,sinphi,costheta,sintheta;
Vector3 tanv, curv, binorm;  /* vectors of Frenet frame of axis curve */
Vector3 scratch,scratch1,avect, bvect;  /* scratch vectors */

for(int i=0;i < segments;i++){

phi   = i*dphi;
theta = i*dtheta;

cosphi = cos(phi); sinphi = sin(phi);
costheta = cos(theta); sintheta = sin(theta);

/***
The following statements, computing axis, tanv and curv, are based on the
elementary differential geometry of space curves.  See any book on
elementary differential geometry, such as those by Dirk Struik
or Barrett O'Neill.

axis is a general point on the spiraling axis curve of the tube.

tanv is the tangent vector to the axis curve, i.e. the derivative of
the axis with respect to arclength.  Here I get it by differentiating
the axis function with respect to i and normalizing.

curv is the unit vector perpendicular to the tangent vector in
the osculating plane of the axis curve, pointing to the center of
of curvature.  It is the negative of the unit vector parallel to the
derivative of tanv with respect to arclength.  Here I get it by taking the
second derivative of the axis function with respect to i,
then taking the component perpendicular to tanv
(by the Gram-Schmidt formula), then normalizing.

binorm is the unit binormal; tanv, curv, and binorm form an orthonormal
Frenet frame for the axis curve.
***/

/***
First write down a parametric representation for axis, using
elementary analytic geometry.  Remember, theta = pitch * phi
***/

/***
Now differentiate the above components with respect to i, using
differentiation formulas of elementary calculus.
***/
tanv.x = cosphi*(minrad_ * -dtheta*sintheta) +

tanv.y = sinphi*( minrad_ * -dtheta*sintheta) +

tanv.z = minrad_ * dtheta * costheta;

/***
Differentiate again, same way.
***/

curv.x =  cosphi*(minrad_ * -dtheta*dtheta*costheta) +

curv.y =  sinphi*( minrad_ * -dtheta*dtheta*costheta) +

curv.z = -minrad_ * dtheta * dtheta * sintheta;

/***
Note: You could approximate tanv and curv, and avoid having to recall the
calculus formulas, by double-differencing axis.  But this would
substantially complicate the loop logic.
***/

/*** Gram-Schmidt formula to get orthonormal frame***/
tanv.normalize();
Vector3 scratch = tanv;
scratch.scalarMultiply(dotProduct(curv, tanv));
curv.subtractVector(scratch);
curv.normalize();    //curve normal vector

/*** complete the Frenet frame ***/
binorm = crossProduct(tanv, curv);         //curve binormal vector

/***
Now compute quadmesh vertices, uniformly spaced on circle of intersection
of tube with plane perpendicular to axis,  starting at the normal vector.
***/
for(int j=0; j< sectors_; j++){
Point3 point = axis;
psi = j*dpsi;
scratch = curv;
point.displace(scratch);
scratch = binorm ;
point.displace(scratch);
vertices_[i*sectors+j] = point;
};
};

/***
At this point, we have the vertices of the quad mesh.
For Gouraud shading, we also need vectors normal to the surface at
the vertices.  For the vertex normal at each vertex, I use the unit vector
parallel to the average of the geometric normals of the four facets
of the tessellation that share the vertex.  Each of these four facet
geometric normals is obtained by taking the unit vector parallel to
the cross product of the two vectors from the given vertex to the
two neighboring vertices of the respective facet.
***/

for(int i=0; i< segments_ ;i++)
for(int j=0; j< sectors_;j++)   {

int k = i* sectors + j;

/* in general */      /* on edges of quad mesh*/
int k1 = (j<sectors-1        ?        k+1                 :       k+1 - sectors                         );
int k2 = (i<segments-1        ?        k+sectors         :              j                                 );
int k3 = (j>0                        ?        k-1                 :           k+sectors-1                    );
int k4 = (i>0                        ?        k-sectors         :          j+(segments-1)*(sectors) );

avect = Vector3(vertices_[k], vertices_[k1]); // right vector
bvect = Vector3(vertices_[k], vertices_[k2]);        //up vector
scratch = crossProduct(bvect, avect);
scratch.normalize();

avect = bvect;
bvect = Vector3(vertices_[k], vertices_[k3]); //left vector
scratch1 = crossProduct(bvect, avect);
scratch1.normalize();

avect = bvect;
bvect  = Vector3(vertices_[k], vertices_[k4]);// down vector
scratch1 = crossProduct(bvect, avect);
scratch1.normalize();
scratch.normalize();
vertexNormals_[k] = scratch ;
};

//Now we need to set up the face index lists.
for(int i=0; i<segments_; i++){
for(int j=0; j<sectors_; j++){
int k = i*sectors_+j;
faces_[k] =  shared_array<int>(new int);
faces_[k] = k;
faces_[k] = (j<sectors_-1 ? k+1 : i*sectors_ );
faces_[k] = (i<segments_-1 ? faces_[k]+sectors_ : (j<sectors_-1 ? j+1 : 0));
faces_[k] = (i<segments_-1 ? k+sectors_ : j);
}
}
}

for(int i=0; i<segments_; i++){
for(int j=0; j<sectors_; j++){
int k = i*sectors_+j;
for (int l=0;l<4;l++){
glNormal3dv(vertexNormals_[faces_[k][l]].getArray());
glVertex3dv(vertices_[faces_[k][l]].getArray());
}
}
}
glEnd();
}

//Some global stuff
Helitoroid helitor(1., 0.2, 0.05, 500, 20, 10);
Torus         torus( 1.0, 0.1, 100, 10);
double rotAngle = 0;
double otrAngle = 0;
double dx = 0.;
double dy = 0.;
double dz = 0.;

void init()
{
glClearColor(0, 0, 0, 0);                // background color
glClearDepth(1.0);                        // background depth value

glMatrixMode(GL_PROJECTION);
gluPerspective(60, 1, 1, 1000);        // setup a perspective projection

glMatrixMode(GL_MODELVIEW);

gluLookAt(                                // set up the camera
0.0, 0.0, 5.0,                // eye position
0.0, 0.0, 0.0,                // lookat position
0.0, 1.0, 0.0);                // up direction

glEnable(GL_DEPTH_TEST);                // enable hidden surface removal

glEnable(GL_LIGHTING);                // enable lighting
glEnable(GL_LIGHT0);                // enable

float lpos[] = { 5, 5, 5, 0 };
glLightfv(GL_LIGHT0, GL_POSITION, lpos);

}
//-----------------------------------------------------------------------
// display callback function
//      This is called each time application needs to redraw itself.
//      Most of the opengl work is done through this function.
//-----------------------------------------------------------------------

void display()
{
glClear(
GL_COLOR_BUFFER_BIT |                // clear the frame buffer (color)
GL_DEPTH_BUFFER_BIT);                // clear the depth buffer (depths)

glPushMatrix();                        // save the current camera transform
glTranslated(dx, dy, dz);
glRotated(rotAngle, 0, 1, 0);        // rotate by rotAngle about y-axis
glRotated(otrAngle, 1, 0, 0);
glEnable(GL_COLOR_MATERIAL);        // specify object color
glColor3f(0.1, 0.1,1.0);         //bluish

helitor.draw();                        // draw the helitoroid
glColor3f(1.0, 0.1, 0.1);                // redish
torus.draw();
glPopMatrix();                        // restore the modelview matrix
glFlush();                                // force OpenGL to render now

glutSwapBuffers();                        // make the image visible
}

//-----------------------------------------------------------------------
// keyboard callback function
//      This is called whenever a keyboard key is hit.
//-----------------------------------------------------------------------

void keyboard(unsigned char k, int x, int y)
{
switch (k)
{
case 'a':
rotAngle += 5;                        // increase rotation by 5 degrees
break;
case 'A':
rotAngle -= 5;                        // decrease rotation by 5 degrees
break;
case 's':
otrAngle += 5;
break;
case 'S':
otrAngle -= 5;
break;
case 'x':
dx += 0.05;
break;
case 'X' :
dx -= 0.05;
break;
case 'y':
dy += 0.05;
break;
case 'Y' :
dy -= 0.05;
break;
case 'z':
dz += 0.05;
break;
case 'Z' :
dz -= 0.05;
break;
case 'q':
exit(0);                        // exit
}
glutPostRedisplay();                // redraw the image now
}

int main(int argc, char** argv){
//build geometry

glutInitDisplayMode(                // initialize GLUT
GLUT_DOUBLE |                // use double buffering
GLUT_DEPTH |                // request memory for z-buffer
GLUT_RGB );                // set RGB color mode

glutCreateWindow("Helitoroid/Torus Using OGL and GLUT");        // create the window

glutDisplayFunc(display);       // call display() to redraw window
glutKeyboardFunc(keyboard); // call keyboard() when key is hit

init();                                // our own initializations

glutMainLoop();                        // let GLUT take care of everything
return 0;
}