/*
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 addVector(const Vector3 &);
    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;} 

void Vector3::addVector(const Vector3 & b){
    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:
    Quadmesh(){};
    virtual void draw();
    virtual ~Quadmesh(){};
protected:
    int segments_;
    int sectors_;
};

class Torus : public Quadmesh {
public:
    Torus(){};
    Torus(double, double, int, int);
    virtual void draw(){Quadmesh::draw();};
private:
    double majrad_;
    double minrad_; 
};

class Helitoroid : public Quadmesh {
public:
    Helitoroid(){};
    Helitoroid(double, double, double, int, int, int);
    virtual void draw(){Quadmesh::draw();};
private:
    double majrad_;
    double minrad_; 
    double tuberad_;
    int pitch_; 
};

/*
 Makes a quadmesh tessellation of a torus
*/
Torus::Torus(double majRad, double minRad, int segments, int sectors) :
majrad_(majRad), minrad_(minRad) {
    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;
        azvect.scalarMultiply( majrad_ );
        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;
            norm.scalarMultiply(minrad_);
            vertices_[k] = azpoint;
            vertices_[k].displace(norm);
        //Now we need to set up the face index lists.
            faces_[k] =  shared_array<int>(new int[4]);
            faces_[k][0] = k;
            faces_[k][1] = (j<sectors_-1    ?  k+1                           : i*sectors_ );
            faces_[k][2] = (i<segments_-1 ? faces_[k][1]+sectors_ : (j<sectors_-1 ? j+1 : 0));
            faces_[k][3] = (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

majrad_ is the radius of the directrix circle.
minrad_ is the constant distance from the directrix circle to the
  spiraling axis curve of the tube.
tubrad_ is the cross-sectional radius 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(
        double majorRadius, double minorRadius, double tubeRadius,
        int segments, int sectors, int pitch):
        majrad_(majorRadius),
        minrad_(minorRadius),
        tuberad_( tubeRadius),
        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
    ***/
      axis.x =  cosphi*(majrad_ + minrad_ * costheta);
      axis.y =  sinphi*(majrad_ + minrad_ * costheta);
      axis.z =  minrad_*sintheta;
    
    /*** 
    Now differentiate the above components with respect to i, using
    differentiation formulas of elementary calculus.
    ***/
      tanv.x = cosphi*(minrad_ * -dtheta*sintheta) +
              -dphi*sinphi*(majrad_ + minrad_ * costheta);
    
      tanv.y = sinphi*( minrad_ * -dtheta*sintheta) +
               dphi*cosphi*(majrad_ + minrad_ * costheta);
    
      tanv.z = minrad_ * dtheta * costheta;
    
    /*** 
    Differentiate again, same way. 
    ***/
    
      curv.x =  cosphi*(minrad_ * -dtheta*dtheta*costheta) +
               -dphi*sinphi*(minrad_ * -dtheta*sintheta) +
          -dphi*sinphi*( minrad_ * -dtheta*sintheta) +
          -dphi*dphi*cosphi*(majrad_ + minrad_ * costheta);
    
      curv.y =  sinphi*( minrad_ * -dtheta*dtheta*costheta) +
                dphi*cosphi*( minrad_ * -dtheta*sintheta) +
                dphi*cosphi*(minrad_ * -dtheta*sintheta)+
               -dphi*dphi*sinphi*(majrad_ + minrad_ * 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;
            scratch.scalarMultiply(tuberad_*cos(psi));
            point.displace(scratch);
            scratch = binorm ;
            scratch.scalarMultiply(tuberad_*sin(psi));
            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();
        scratch.addVector(scratch1);
   
        avect = bvect;
        bvect  = Vector3(vertices_[k], vertices_[k4]);// down vector 
        scratch1 = crossProduct(bvect, avect);
        scratch1.normalize();
        scratch.addVector(scratch1); 
        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[4]);
            faces_[k][0] = k;
            faces_[k][1] = (j<sectors_-1 ? k+1 : i*sectors_ );
            faces_[k][2] = (i<segments_-1 ? faces_[k][1]+sectors_ : (j<sectors_-1 ? j+1 : 0));
            faces_[k][3] = (i<segments_-1 ? k+sectors_ : j);
        }
    }
}

void Quadmesh::draw() {
    glBegin(GL_QUADS);
    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);
    glLoadIdentity();
    gluPerspective(60, 1, 1, 1000);        // setup a perspective projection

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();

    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);

    // glShadeModel(GL_FLAT);               // flat shading
    glShadeModel(GL_SMOOTH);                // smooth shading
}
//-----------------------------------------------------------------------
// 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; 
}