/*
 * CS5600 University of Utah
 * Charles McGarvey
 * mcgarvey@eng.utah.edu
 */

#ifndef _MAT_H_
#define _MAT_H_

#include "common.h"
#include "vec.h"


/*
 * A simple matrix class with column-major storage and notation.
 */
struct mat
{
    vec_t v[4];
};
typedef struct mat mat_t;

/*
 * Initialize a matrix with individual components, row by row.
 */
INLINE_MAYBE
void mat_init(mat_t* m, scal_t m11, scal_t m12, scal_t m13, scal_t m14,
                        scal_t m21, scal_t m22, scal_t m23, scal_t m24,
                        scal_t m31, scal_t m32, scal_t m33, scal_t m34,
                        scal_t m41, scal_t m42, scal_t m43, scal_t m44)
{
    m->v[0] = vec_new2(m11, m21, m31, m41);
    m->v[1] = vec_new2(m12, m22, m32, m42);
    m->v[2] = vec_new2(m13, m23, m33, m43);
    m->v[3] = vec_new2(m14, m24, m34, m44);
}


/*
 * Create a new matrix with individual components, row by row.
 */
INLINE_MAYBE
mat_t mat_new(scal_t m11, scal_t m12, scal_t m13, scal_t m14,
              scal_t m21, scal_t m22, scal_t m23, scal_t m24,
              scal_t m31, scal_t m32, scal_t m33, scal_t m34,
              scal_t m41, scal_t m42, scal_t m43, scal_t m44)
{
    mat_t m;
    mat_init(&m, m11, m12, m13, m14,
                 m21, m22, m23, m24,
                 m31, m32, m33, m34,
                 m41, m42, m43, m44);
    return m;
}

/*
 * Create a new matrix with four column vectors.
 */
INLINE_MAYBE
mat_t mat_new2(vec_t a, vec_t b, vec_t c, vec_t d)
{
    mat_t m;
    m.v[0] = a;
    m.v[1] = b;
    m.v[2] = c;
    m.v[3] = d;
    return m;
}

#define MAT_IDENTITY mat_new(S(1.0), S(0.0), S(0.0), S(0.0), \
                             S(0.0), S(1.0), S(0.0), S(0.0), \
                             S(0.0), S(0.0), S(1.0), S(0.0), \
                             S(0.0), S(0.0), S(0.0), S(1.0))


/*
 * Get a column vector (can also access the vector array directly).
 */
INLINE_MAYBE
vec_t mat_col(mat_t m, int i)
{
    return m.v[i];
}

/*
 * Get a row vector.
 */
INLINE_MAYBE
vec_t mat_row(mat_t m, int i)
{
    switch (i) {
    case 0:
        return vec_new2(m.v[0].x, m.v[1].x, m.v[2].x, m.v[3].x);
    case 1:
        return vec_new2(m.v[0].y, m.v[1].y, m.v[2].y, m.v[3].y);
    case 2:
        return vec_new2(m.v[0].z, m.v[1].z, m.v[2].z, m.v[3].z);
    case 3:
        return vec_new2(m.v[0].w, m.v[1].w, m.v[2].w, m.v[3].w);
    }
}


/*
 * Print the matrix to stdout.
 */
INLINE_MAYBE
void mat_print(mat_t m)
{
    printf("|");
    vec_print(mat_row(m, 0));
    printf("|\n|");
    vec_print(mat_row(m, 1));
    printf("|\n|");
    vec_print(mat_row(m, 2));
    printf("|\n|");
    vec_print(mat_row(m, 3));
    printf("|");
}


/*
 * Multiply two matrices together.
 */
INLINE_MAYBE
mat_t mat_mult(mat_t a, mat_t b)
{
#define _DOT(I,J) vec_dot2(mat_row(a,I), mat_col(b,J))
    return mat_new(_DOT(0,0), _DOT(0,1), _DOT(0,2), _DOT(0,3),
                   _DOT(1,0), _DOT(1,1), _DOT(1,2), _DOT(1,3),
                   _DOT(2,0), _DOT(2,1), _DOT(2,2), _DOT(2,3),
                   _DOT(3,0), _DOT(3,1), _DOT(3,2), _DOT(3,3));
#undef _DOT
}

/*
 * Transform a vector using a matrix.
 */
INLINE_MAYBE
vec_t mat_apply(mat_t m, vec_t v)
{
    return vec_new2(vec_dot2(v,mat_row(m,0)),
                    vec_dot2(v,mat_row(m,1)),
                    vec_dot2(v,mat_row(m,2)),
                    vec_dot2(v,mat_row(m,3)));
}


/*
 * Create a new translate matrix.
 */
INLINE_MAYBE
mat_t MAT_TRANSLATE(scal_t x, scal_t y, scal_t z)
{
    return mat_new(S(1.0), S(0.0), S(0.0),      x,
                   S(0.0), S(1.0), S(0.0),      y,
                   S(0.0), S(0.0), S(1.0),      z,
                   S(0.0), S(0.0), S(0.0), S(1.0));
}

/*
 * Create a new translate matrix from a vector.
 */
INLINE_MAYBE
mat_t MAT_TRANSLATE2(vec_t v)
{
    return MAT_TRANSLATE(v.x, v.y, v.z);
}

/*
 * Create a new scale matrix.
 */
INLINE_MAYBE
mat_t MAT_SCALE(scal_t x, scal_t y, scal_t z)
{
    return mat_new(x, S(0.0), S(0.0), S(0.0),
              S(0.0),      y, S(0.0), S(0.0),
              S(0.0), S(0.0),      z, S(0.0),
              S(0.0), S(0.0), S(0.0), S(1.0));
}

/*
 * Create a new scale matrix from a vector.
 */
INLINE_MAYBE
mat_t MAT_SCALE2(vec_t v)
{
    return MAT_SCALE(v.x, v.y, v.z);
}

/*
 * Create a rotation matrix (around the X axis).
 */
INLINE_MAYBE
mat_t MAT_ROTATE_X(scal_t theta)
{
    scal_t sin_a = scal_sin(theta);
    scal_t cos_a = scal_cos(theta);
    return mat_new(S(1.0), S(0.0), S(0.0), S(0.0),
                   S(0.0),  cos_a, -sin_a, S(0.0),
                   S(0.0),  sin_a,  cos_a, S(0.0),
                   S(0.0), S(0.0), S(0.0), S(1.0));
}

/*
 * Create a rotation matrix (around the Y axis).
 */
INLINE_MAYBE
mat_t MAT_ROTATE_Y(scal_t theta)
{
    scal_t sin_a = scal_sin(theta);
    scal_t cos_a = scal_cos(theta);
    return mat_new(cos_a, S(0.0),  sin_a, S(0.0),
                  S(0.0), S(1.0), S(0.0), S(0.0),
                  -sin_a, S(0.0),  cos_a, S(0.0),
                  S(0.0), S(0.0), S(0.0), S(1.0));
}

/*
 * Create a rotation matrix (around the Z axis).
 */
INLINE_MAYBE
mat_t MAT_ROTATE_Z(scal_t theta)
{
    scal_t sin_a = scal_sin(theta);
    scal_t cos_a = scal_cos(theta);
    return mat_new(cos_a, -sin_a, S(0.0), S(0.0),
                   sin_a,  cos_a, S(0.0), S(0.0),
                  S(0.0), S(0.0), S(1.0), S(0.0),
                  S(0.0), S(0.0), S(0.0), S(1.0));
}

/*
 * Create a rotation matrix (around an arbitrary axis).
 */
INLINE_MAYBE
mat_t MAT_ROTATE(scal_t theta, scal_t x, scal_t y, scal_t z)
{
    /*
     * This code is an implementation of an algorithm described by Glenn
     * Murray at http://inside.mines.edu/~gmurray/ArbitraryAxisRotation/
     */
    vec_t v = vec_normalize(vec_new(x, y, z));
    x = v.x;
    y = v.y;
    z = v.z;
    scal_t sin_a = scal_sin(theta);
    scal_t cos_a = scal_cos(theta);
    scal_t x2 = x * x;
    scal_t y2 = y * y;
    scal_t z2 = z * z;
    return mat_new(x2+(S(1.0)-x2)*cos_a, x*y*(S(1.0)-cos_a)-z*sin_a, x*z*(S(1.0)-cos_a)+y*sin_a, S(0.0),
                   x*y*(S(1.0)-cos_a)+z*sin_a, y2+(S(1.0)-y2)*cos_a, y*z*(S(1.0)-cos_a)-y*sin_a, S(0.0),
                   x*z*(S(1.0)-cos_a)-y*sin_a, y*z*(S(1.0)-cos_a)+x*sin_a, z2+(S(1.0)-z2)*cos_a, S(0.0),
                   S(0.0), S(0.0), S(0.0), S(1.0));
}

/*
 * Create a view matrix based on eye, spot, and an up vector.
 */
INLINE_MAYBE
mat_t MAT_LOOKAT(vec_t eye, vec_t spot, vec_t up)
{
    vec_t f = vec_normalize(vec_sub(spot, eye));
    vec_t s = vec_normalize(vec_cross(f, up));
    vec_t u = vec_cross(s, f);
    return mat_mult(mat_new(s.x,    s.y,    s.z, S(0.0),
                            u.x,    u.y,    u.z, S(0.0),
                           -f.x,   -f.y,   -f.z, S(0.0),
                         S(0.0), S(0.0), S(0.0), S(1.0)), MAT_TRANSLATE2(vec_neg(eye)));
}

/*
 * Create a 3D orthogonal projection matrix.
 */
INLINE_MAYBE
mat_t MAT_ORTHO(scal_t left, scal_t right,
                scal_t bottom, scal_t top,
                scal_t near, scal_t far)
{
    scal_t rml = right - left;
    scal_t rpl = right + left;
    scal_t tmb = top - bottom;
    scal_t tpb = top + bottom;
    scal_t fmn = far - near;
    scal_t fpn = far + near;
    return mat_new(S(2.0)/rml,     S(0.0),      S(0.0), -rpl / rml,
                       S(0.0), S(2.0)/tmb,      S(0.0), -tpb / tmb,
                       S(0.0),     S(0.0), S(-2.0)/fmn, -fpn / fmn,
                       S(0.0),     S(0.0),      S(0.0),     S(1.0));
}

/*
 * Create a frustum-based projection matrix.
 */
INLINE_MAYBE
mat_t MAT_FRUSTUM(scal_t left, scal_t right,
                  scal_t bottom, scal_t top,
                  scal_t near, scal_t far)
{
    scal_t rml = right - left;
    scal_t rpl = right + left;
    scal_t tmb = top - bottom;
    scal_t tpb = top + bottom;
    scal_t fmn = far - near;
    scal_t fpn = far + near;
    scal_t n2  = near * S(2.0);
    return mat_new(n2/rml, S(0.0),  rpl/rml,      S(0.0),
                   S(0.0), n2/tmb,  tpb/tmb,      S(0.0),
                   S(0.0), S(0.0), -fpn/fmn, -n2*far/fmn,
                   S(0.0), S(0.0),  S(-1.0),      S(0.0));
}

/*
 * Create a perspective projection matrix.
 */
INLINE_MAYBE
mat_t MAT_PERSPECTIVE(scal_t fovy, scal_t aspect, scal_t near, scal_t far)
{
    scal_t top    = near * scal_tan(fovy * S(0.5));
    scal_t bottom = -top;
    scal_t left   = bottom * aspect;
    scal_t right  = top * aspect;
    return MAT_FRUSTUM(left, right, bottom, top, near, far);
}

/*
 * Create a viewport matrix.
 */
INLINE_MAYBE
mat_t MAT_VIEWPORT(int x, int y, int w, int h)
{
    scal_t xs = (scal_t)x;
    scal_t ys = (scal_t)y;
    scal_t ws = (scal_t)w * S(0.5);
    scal_t hs = (scal_t)h * S(0.5);
    return mat_new(ws, S(0.0), S(0.0),  ws+xs,
               S(0.0),     hs, S(0.0),  hs+ys,
               S(0.0), S(0.0), S(1.0), S(0.0),
               S(0.0), S(0.0), S(0.0), S(1.0));
}


#endif // _MAT_H_