/* revision 5 */

/* @2009 Kamil Kaminski
 * 
 * this code is not yet endian aware
 * the style of the syntax is original k&r except there's \n
 * after the opening { and extra space after if statement,
 * and for/while loops
 */

#include "math3d.h"

void m3dFindNormal(M3DVector3f result, const M3DVector3f point1, 
                   const M3DVector3f point2, const M3DVector3f point3)
{
    M3DVector3f v1, v2;

    /* calculate two vectors from the three points, assumes 
     * counter clockwise winding
     */
    v1[0] = point1[0] - point2[0];
    v1[1] = point1[1] - point2[1];
    v1[2] = point1[2] - point2[2];

    v2[0] = point2[0] - point3[0];
    v2[1] = point2[1] - point3[1];
    v2[2] = point2[2] - point3[2];

    /* take the cross product of the two vectors to get the normal vector */
    m3dCrossProduct(result, v1, v2);
}

void m3dLoadIdentity44(M3DMatrix44f m) /* 4x4 float */
{
    /* don't be fooled, this is still column major */
    static M3DMatrix44f identity = { 1.0f, 0.0f, 0.0f, 0.0f,
                                     0.0f, 1.0f, 0.0f, 0.0f,
                                     0.0f, 0.0f, 1.0f, 0.0f,
                                     0.0f, 0.0f, 0.0f, 1.0f };

    memcpy(m, identity, sizeof(M3DMatrix44f));
}

/* creates a 4x4 rotation matrix, takes radians not degrees */
void m3dRotationMatrix44(M3DMatrix44f m, float angle, float x, float y, float z)
{
    float mag, s, c;
    float xx, yy, zz, xy, yz, zx, xs, ys, zs, one_c;

    s = (float) (sin(angle));
    c = (float) (cos(angle));

    mag = (float) (sqrt(x*x + y*y + z*z));

    /* identity matrix */
    if (mag == 0.0f)
    {
        m3dLoadIdentity44(m);
        return;
    }

    /* rotation matrix is normalized */
    x /= mag;
    y /= mag;
    z /= mag;

    #define M(row,col) m[col*4+row]

    xx = x * x;
    yy = y * y;
    zz = z * z;
    xy = x * y;
    yz = y * z;
    zx = z * x;

    xs = x * s;
    ys = y * s;
    zs = z * s;
    one_c = 1.0f - c;

    M(0,0) = (one_c * xx) + c;
    M(0,1) = (one_c * xy) - zs;
    M(0,2) = (one_c * zx) + ys;
    M(0,3) = 0.0f;

    M(1,0) = (one_c * xy) + zs;
    M(1,1) = (one_c * yy) + c;
    M(1,2) = (one_c * yz) - xs;
    M(1,3) = 0.0f;

    M(2,0) = (one_c * zx) - ys;
    M(2,1) = (one_c * yz) + xs;
    M(2,2) = (one_c * zz) + c;
    M(2,3) = 0.0f;

    M(3,0) = 0.0f;
    M(3,1) = 0.0f;
    M(3,2) = 0.0f;
    M(3,3) = 1.0f;

    #undef M
}

/* draw a torus (doughnut), using the current 1d texture for light shading */
/* this funct accepts 4x4 trans matrix to be applied to the vertices */
void gltDrawTorus(GLfloat majorRadius, GLfloat minorRadius, GLint numMajor, 
                  GLint numMinor)
{
    M3DVector3f vNormal;
    double majorStep = 2.0f * M3D_PI / numMajor;
    double minorStep = 2.0f * M3D_PI / numMinor;
    int i, j;

    for (i = 0; i < numMajor; ++i)
    {
        double a0 = i * majorStep;
        double a1 = a0 + majorStep;
        GLfloat x0 = (GLfloat) cos(a0);
        GLfloat y0 = (GLfloat) sin(a0);
        GLfloat x1 = (GLfloat) cos(a1);
        GLfloat y1 = (GLfloat) sin(a1);

        glBegin(GL_TRIANGLE_STRIP);
        for (j = 0; j <= numMinor; ++j)
        {
            double b = j * minorStep;
            GLfloat c = (GLfloat) cos(b);
            GLfloat r = minorRadius * c + majorRadius;
            GLfloat z = minorRadius * (GLfloat) sin(b);

            glTexCoord2f((float) (i) / (float) (numMajor), (float) (j) \
                         / (float) (numMinor));
            vNormal[0] = x0 * c;
            vNormal[1] = y0 * c;
            vNormal[2] = z / minorRadius;
            m3dNormalizeVector(vNormal);
            glNormal3fv(vNormal);
            glVertex3f(x0 * r, y0 * r, z);

            glTexCoord2f((float) (i + 1) / (float) (numMajor), (float) (j) \
                         / (float) (numMinor));
            vNormal[0] = x1 * c;
            vNormal[1] = y1 * c;
            vNormal[2] = z / minorRadius;
            m3dNormalizeVector(vNormal);
            glNormal3fv(vNormal);
            glVertex3f(x1 * r, y1 * r, z);
        }
        glEnd();
    }
}

/* this function just specifically draws the jet */
/* FIXME needs to accepts parameters of location and lightning */
void DrawJet(int nShadow)
{
    M3DVector3f vNormal;

    /* nose cone, set material color, note we only have to set to black 
     * for the shadow once 
     */
    if (nShadow == 0)
        glColor3ub(128, 128, 128);
    else
        glColor3ub(0, 0, 0);

    /* nose cone, points straight down, set material color */
    /* follow few lines use manual approach */
    glBegin(GL_TRIANGLES);
    glNormal3f(0.0f, -1.0f, 0.0f);
    glNormal3f(0.0f, -1.0f, 0.0f);
    glVertex3f(0.0f, 0.0f, 60.0f);
    glVertex3f(-15.0f, 0.0f, 30.0f);
    glVertex3f(15.0f, 0.0f, 30.0f);

    /* verticies for this panel */
    {
    M3DVector3f vPoints[3] = { {15.0f, 0.0f, 30.0f}
    ,
    {0.0f, 15.0f, 30.0f}
    ,
    {0.0f, 0.0f, 60.0f}
    };

    /* calculate the normal for the plane */
    m3dFindNormal(vNormal, vPoints[0], vPoints[1], vPoints[2]);
    glNormal3fv(vNormal);
    glVertex3fv(vPoints[0]);
    glVertex3fv(vPoints[1]);
    glVertex3fv(vPoints[2]);
    }

    {
    M3DVector3f vPoints[3] = { {0.0f, 0.0f, 60.0f}
    ,
    {0.0f, 15.0f, 30.0f}
    ,
    {-15.0f, 0.0f, 30.0f}
    };

    m3dFindNormal(vNormal, vPoints[0], vPoints[1], vPoints[2]);
    glNormal3fv(vNormal);
    glVertex3fv(vPoints[0]);
    glVertex3fv(vPoints[1]);
    glVertex3fv(vPoints[2]);
    }

    /* body of the plane */
    {
    M3DVector3f vPoints[3] = { {-15.0f, 0.0f, 30.0f}
    ,
    {0.0f, 15.0f, 30.0f}
    ,
    {0.0f, 0.0f, -56.0f}
    };

    m3dFindNormal(vNormal, vPoints[0], vPoints[1], vPoints[2]);
    glNormal3fv(vNormal);
    glVertex3fv(vPoints[0]);
    glVertex3fv(vPoints[1]);
    glVertex3fv(vPoints[2]);
    }

    {
    M3DVector3f vPoints[3] = { {0.0f, 0.0f, -56.0f}
    ,
    {0.0f, 15.0f, 30.0f}
    ,
    {15.0f, 0.0f, 30.0f}
    };

    m3dFindNormal(vNormal, vPoints[0], vPoints[1], vPoints[2]);
    glNormal3fv(vNormal);
    glVertex3fv(vPoints[0]);
    glVertex3fv(vPoints[1]);
    glVertex3fv(vPoints[2]);
    }

    glNormal3f(0.0f, -1.0f, 0.0f);
    glVertex3f(15.0f, 0.0f, 30.0f);
    glVertex3f(-15.0f, 0.0f, 30.0f);
    glVertex3f(0.0f, 0.0f, -56.0f);

    /* left wing, large triangle for bottom of wing */
    {
    M3DVector3f vPoints[3] = { {0.0f, 2.0f, 27.0f}
    ,
    {-60.0f, 2.0f, -8.0f}
    ,
    {60.0f, 2.0f, -8.0f}
    };

    m3dFindNormal(vNormal, vPoints[0], vPoints[1], vPoints[2]);
    glNormal3fv(vNormal);
    glVertex3fv(vPoints[0]);
    glVertex3fv(vPoints[1]);
    glVertex3fv(vPoints[2]);
    }

    {
    M3DVector3f vPoints[3] = { {60.0f, 2.0f, -8.0f}
    ,
    {0.0f, 7.0f, -8.0f}
    ,
    {0.0f, 2.0f, 27.0f}
    };

    m3dFindNormal(vNormal, vPoints[0], vPoints[1], vPoints[2]);
    glNormal3fv(vNormal);
    glVertex3fv(vPoints[0]);
    glVertex3fv(vPoints[1]);
    glVertex3fv(vPoints[2]);
    }

    {
    M3DVector3f vPoints[3] = { {60.0f, 2.0f, -8.0f}
    ,
    {-60.0f, 2.0f, -8.0f}
    ,
    {0.0f, 7.0f, -8.0f}
    };

    m3dFindNormal(vNormal, vPoints[0], vPoints[1], vPoints[2]);
    glNormal3fv(vNormal);
    glVertex3fv(vPoints[0]);
    glVertex3fv(vPoints[1]);
    glVertex3fv(vPoints[2]);
    }

    {
    M3DVector3f vPoints[3] = { {0.0f, 2.0f, 27.0f}
    ,
    {0.0f, 7.0f, -8.0f}
    ,
    {-60.0f, 2.0f, -8.0f}
    };

    m3dFindNormal(vNormal, vPoints[0], vPoints[1], vPoints[2]);
    glNormal3fv(vNormal);
    glVertex3fv(vPoints[0]);
    glVertex3fv(vPoints[1]);
    glVertex3fv(vPoints[2]);
    }

    /* tail section */
    /* bottom of back fin */
    glNormal3f(0.0f, -1.0f, 0.0f);
    glVertex3f(-30.0f, -0.50f, -57.0f);
    glVertex3f(30.0f, -0.50f, -57.0f);
    glVertex3f(0.0f, -0.50f, -40.0f);

    {
    M3DVector3f vPoints[3] = { {0.0f, -0.5f, -40.0f}
    ,
    {30.0f, -0.5f, -57.0f}
    ,
    {0.0f, 4.0f, -57.0f}
    };

    m3dFindNormal(vNormal, vPoints[0], vPoints[1], vPoints[2]);
    glNormal3fv(vNormal);
    glVertex3fv(vPoints[0]);
    glVertex3fv(vPoints[1]);
    glVertex3fv(vPoints[2]);
    }

    {
    M3DVector3f vPoints[3] = { {0.0f, 4.0f, -57.0f}
    ,
    {-30.0f, -0.5f, -57.0f}
    ,
    {0.0f, -0.5f, -40.0f}
    };

    m3dFindNormal(vNormal, vPoints[0], vPoints[1], vPoints[2]);
    glNormal3fv(vNormal);
    glVertex3fv(vPoints[0]);
    glVertex3fv(vPoints[1]);
    glVertex3fv(vPoints[2]);
    }

    {
    M3DVector3f vPoints[3] = { {30.0f, -0.5f, -57.0f}
    ,
    {-30.0f, -0.5f, -57.0f}
    ,
    {0.0f, 4.0f, -57.0f}
    };

    m3dFindNormal(vNormal, vPoints[0], vPoints[1], vPoints[2]);
    glNormal3fv(vNormal);
    glVertex3fv(vPoints[0]);
    glVertex3fv(vPoints[1]);
    glVertex3fv(vPoints[2]);
    }

    {
    M3DVector3f vPoints[3] = { {0.0f, 0.5f, -40.0f}
    ,
    {3.0f, 0.5f, -57.0f}
    ,
    {0.0f, 25.0f, -65.0f}
    };

    m3dFindNormal(vNormal, vPoints[0], vPoints[1], vPoints[2]);
    glNormal3fv(vNormal);
    glVertex3fv(vPoints[0]);
    glVertex3fv(vPoints[1]);
    glVertex3fv(vPoints[2]);
    }

    {
    M3DVector3f vPoints[3] = { {0.0f, 25.0f, -65.0f}
    ,
    {-3.0f, 0.5f, -57.0f}
    ,
    {0.0f, 0.5f, -40.0f}
    };

    m3dFindNormal(vNormal, vPoints[0], vPoints[1], vPoints[2]);
    glNormal3fv(vNormal);
    glVertex3fv(vPoints[0]);
    glVertex3fv(vPoints[1]);
    glVertex3fv(vPoints[2]);
    }

    {
    M3DVector3f vPoints[3] = { {3.0f, 0.5f, -57.0f}
    ,
    {-3.0f, 0.5f, -57.0f}
    ,
    {0.0f, 25.0f, -65.0f}
    };

    m3dFindNormal(vNormal, vPoints[0], vPoints[1], vPoints[2]);
    glNormal3fv(vNormal);
    glVertex3fv(vPoints[0]);
    glVertex3fv(vPoints[1]);
    glVertex3fv(vPoints[2]);
    }

    glEnd();
}

void gltDrawUnitAxes(void)
{
    GLUquadricObj *pObj; /* temporary, used for quadrics */

    /* measurements */
    float fAxisRadius = 0.025f;
    float fAxisHeight = 1.0f;
    float fArrowRadius = 0.06f;
    float fArrowHeight = 0.1f;

    /* setup the quadric object */
    pObj = gluNewQuadric();
    gluQuadricDrawStyle(pObj, GLU_FILL);
    gluQuadricNormals(pObj, GLU_SMOOTH);
    gluQuadricOrientation(pObj, GLU_OUTSIDE);
    gluQuadricTexture(pObj, GLU_FALSE);

    /* draw the blue z axis first with arrowed head */
    glColor3f(0.0f, 0.0f, 1.0f);
    gluCylinder(pObj, fAxisRadius, fAxisRadius, fAxisHeight, 10, 1);
    glPushMatrix();
    glTranslatef(0.0f, 0.0f, 1.0f);
    gluCylinder(pObj, fArrowRadius, 0.0f, fArrowHeight, 10, 1);
    glRotatef(180.0f, 1.0f, 0.0f, 0.0f);
    gluDisk(pObj, fAxisRadius, fArrowRadius, 10, 1);
    glPopMatrix();

    /* draw the red x axis 2nd with arrowed head */
    glColor3f(1.0f, 0.0f, 0.0f);
    glPushMatrix();
    glRotatef(90.0f, 0.0f, 1.0f, 0.0f);
    gluCylinder(pObj, fAxisRadius, fAxisRadius, fAxisHeight, 10, 1);
    glPushMatrix();
    glTranslatef(0.0f, 0.0f, 1.0f);
    gluCylinder(pObj, fArrowRadius, 0.0f, fArrowHeight, 10, 1);
    glRotatef(180.0f, 0.0f, 1.0f, 0.0f);
    gluDisk(pObj, fAxisRadius, fArrowRadius, 10, 1);
    glPopMatrix();
    glPopMatrix();

    /* draw the green y axis 3rd with arrowed head */
    glColor3f(0.0f, 1.0f, 0.0f);
    glPushMatrix();
    glRotatef(-90.0f, 1.0f, 0.0f, 0.0f);
    gluCylinder(pObj, fAxisRadius, fAxisRadius, fAxisHeight, 10, 1);
    glPushMatrix();
    glTranslatef(0.0f, 0.0f, 1.0f);
    gluCylinder(pObj, fArrowRadius, 0.0f, fArrowHeight, 10, 1);
    glRotatef(180.0f, 1.0f, 0.0f, 0.0f);
    gluDisk(pObj, fAxisRadius, fArrowRadius, 10, 1);
    glPopMatrix();
    glPopMatrix();

    /* white sphere at origin */
    glColor3f(1.0f, 1.0f, 1.0f);
    gluSphere(pObj, 0.05f, 15, 15);

    /* delete the quadric */
    gluDeleteQuadric(pObj);
}

#define A(row,col)  a[(col<<2)+row]
#define B(row,col)  b[(col<<2)+row]
#define P(row,col)  product[(col<<2)+row]

void m3dMatrixMultiply44(M3DMatrix44f product, const M3DMatrix44f a, 
                         const M3DMatrix44f b)
{
    int i;
    for (i = 0; i < 4; i++)
    {
        float ai0 = A(i, 0), ai1 = A(i, 1), ai2 = A(i, 2), ai3 = A(i, 3);
        P(i, 0) = ai0 * B(0, 0) + ai1 * B(1, 0) + ai2 * B(2, 0) + ai3 * B(3, 0);
        P(i, 1) = ai0 * B(0, 1) + ai1 * B(1, 1) + ai2 * B(2, 1) + ai3 * B(3, 1);
        P(i, 2) = ai0 * B(0, 2) + ai1 * B(1, 2) + ai2 * B(2, 2) + ai3 * B(3, 2);
        P(i, 3) = ai0 * B(0, 3) + ai1 * B(1, 3) + ai2 * B(2, 3) + ai3 * B(3, 3);
    }
}

#undef A
#undef B
#undef P


/*************************************************************/
/* unrelated functions that do not have much to do with math */
/*                                                           */
/*************************************************************/

GLint gltWriteTGA(const char *szFileName)
{
    FILE *pFile;              /* file pointer */
    TGAHEADER tgaHeader;      /* tga file header */
    unsigned long lImageSize; /* size in bytes of image */
    GLbyte *pBits = NULL;     /* pointer to bits */
    GLint iViewport[4];       /* viewport in pixels */
    GLenum lastBuffer;        /* storage for the current read buffer setting */

    /* get the viewport dimensions */
    glGetIntegerv(GL_VIEWPORT, iViewport);

    /* how big is the image going to be (targas are tightly packed) */
    lImageSize = iViewport[2] * 3 * iViewport[3];

    /* allocate block, if this doesn't work, go home */
    pBits = (GLbyte *) malloc(lImageSize);
    if (pBits == NULL)
    {
        perror("malloc");
        return 0;
    }

    /* read bits from color buffer */
    glPixelStorei(GL_PACK_ALIGNMENT, 1);
    glPixelStorei(GL_PACK_ROW_LENGTH, 0);
    glPixelStorei(GL_PACK_SKIP_ROWS, 0);
    glPixelStorei(GL_PACK_SKIP_PIXELS, 0);

    /* get the current read buffer setting and save it, switch to
     * the front buffer and do the read operation, finally, restore
     * the read buffer state
     */
    glGetIntegerv(GL_READ_BUFFER, (GLint *) &lastBuffer);
    glReadBuffer(GL_FRONT);
    glReadPixels(0, 0, iViewport[2], iViewport[3], GL_BGR_EXT,
    GL_UNSIGNED_BYTE, pBits);
    glReadBuffer(lastBuffer);

    /* initialize the targa header */
    tgaHeader.identsize = 0;
    tgaHeader.colorMapType = 0;
    tgaHeader.imageType = 2;
    tgaHeader.colorMapStart = 0;
    tgaHeader.colorMapLength = 0;
    tgaHeader.colorMapBits = 0;
    tgaHeader.xstart = 0;
    tgaHeader.ystart = 0;
    tgaHeader.width = iViewport[2];
    tgaHeader.height = iViewport[3];
    tgaHeader.bits = 24;
    tgaHeader.descriptor = 0;

    /* attempt to open the file */
    pFile = fopen(szFileName, "wb");
    if (pFile == NULL)
    {
        perror("fopen");
        free(pBits); /* free buffer and return error */
        return 0;
    }
    
    /* write the header */
    fwrite(&tgaHeader, sizeof(TGAHEADER), 1, pFile);

    /* write the image data */
    fwrite(pBits, lImageSize, 1, pFile);

    /* free temporary buffer and close the file */
    free(pBits);
    fclose(pFile);

    return 1;
}

GLbyte *gltLoadTGA(const char *szFileName, GLint *iWidth, GLint *iHeight, 
                   GLint *iComponents, GLenum *eFormat)
{
    FILE *pFile;              /* file pointer */
    TGAHEADER tgaHeader;      /* TGA file header */
    unsigned long lImageSize; /* size in bytes of image */
    short sDepth;             /* pixel depth; */
    GLbyte *pBits = NULL;     /* pointer to bits */

    /* default/failed values */
    *iWidth = 0;
    *iHeight = 0;
    *eFormat = GL_BGR_EXT;
    *iComponents = GL_RGB8;

    /* attempt to open the file */
    pFile = fopen(szFileName, "rb");
    if (pFile == NULL)
    {
        perror("fopen");
        return 0;
    }

    /* read in header (binary) */
    fread(&tgaHeader, 18 /* sizeof(TGAHEADER) */, 1, pFile);

    /* get width, height, and depth of texture */
    *iWidth = tgaHeader.width;
    *iHeight = tgaHeader.height;
    sDepth = tgaHeader.bits / 8;

    /* put some validity checks here, very simply, i only understand
     * or care about 8, 24, or 32 bit targas
     */
    if (tgaHeader.bits != 8 && tgaHeader.bits != 24 && tgaHeader.bits != 32)
        return NULL;

    /* calculate size of image buffer */
    lImageSize = tgaHeader.width * tgaHeader.height * sDepth;

    /* allocate memory and check for success */
    pBits = (GLbyte *) malloc(lImageSize * sizeof(GLbyte));
    if (pBits == NULL)
    {
        perror("malloc");
        return 0;
    }

    /* read in the bits */
    /* check for read error, this should catch rle or other  */
    /* weird formats that i don't want to recognize */
    if (fread(pBits, lImageSize, 1, pFile) != 1)
    {
        perror("fread");
        free(pBits);
        return NULL;
    }
    
    /* set opengl format expected */
    switch (sDepth)
    {
        case 3: /* most likely case */
            *eFormat = GL_BGR_EXT;
            *iComponents = GL_RGB8;
            break;
        case 4:
            *eFormat = GL_BGRA_EXT;
            *iComponents = GL_RGBA8;
            break;
        case 1:
            *eFormat = GL_LUMINANCE;
            *iComponents = GL_LUMINANCE8;
            break;
    }

    /* done with file */
    fclose(pFile);

    /* return pointer to image data */
    return pBits;
}