mat-yjc-new.h

//////////////////////////////////////////////////////////////////////////////
//
//  --- mat-yjc-new.h (Modified by Yi-Jen Chiang) ---
//
// YJC Comments:
//  1. The matrices are all in *row order* (as opposed to *column order* in OpenGL).
//
//  2. Function for general rotation, mat4 Rotate(angle, x, y, z) is added, where
//     angle is in degrees, and vector (x, y, z) can have length != 1.0.
//
//  3. For mat2, mat3 and mat4, the transpose() functions are all *incorrect*.
//     The correct new functions "transpose1()" are added for each type.
//     ==> Use "transpose1()" instead.
//
//  3'. Originally the function LookAt() was *incorrect*. It has now been fixed.
//      The correct function is still called "LookAt()". (Search for it to see it.)
//
//  4. For matrix A, ust A[i][j] to get the item at i-th row and j-th column
//     (i, j both start from 0).
//   
//  5. mat4 can be specified as
//     i. mat4(a, b, c, d): a becomes the first row, b the 2nd row, etc., where
//                          a, b, c, d are each a vec4.
//    ii. mat4(m00, m10, m20, m30,..., m33): the 16 floats m00, ..., m33 are given
//        in *column order* but the matrix is stored in *row order*.
// (Note: Search for the part of mat4 and look at the comments marked with "YJC:".)
//
//  6. Function for obtaining the Normal Matrix,
//     mat3 NormalMatrix(mv, non_uniform_scale_flag) is added, where
//     mv is the Model-View matrix (mat4) and 
//     non_uniform_scale_flag == 1 if mv involves non-uniform scaling, and
//     non_uniform_scale_flag == 0 otherwise.
//                  
//     The following related functions are also added:
//
//     mat3 upperLeftMat3(m): return the upper-left 3x3 submatrix of mat4 m.
//     mat3 inverse(m): return the inverse of 3x3 matrix m.
//     mat4 mat4WithUpperLeftMat3(m): return the mat4 where the
//          upper-left 3x3 submatrix is m, the 4th column and the 4th row are
//          both (0, 0, 0, 1).
//                  
//////////////////////////////////////////////////////////////////////////////

#ifndef __ANGEL_MAT_H__
#define __ANGEL_MAT_H__

#include "vec.h"
#include <stdio.h>

// YJC: added the following for the general rotation function Rotate().
#define _USE_MATH_DEFINES  1 // Include constants defined in math.h
#include <math.h>

namespace Angel {

//----------------------------------------------------------------------------
//
//  mat2 - 2D square matrix
//

class mat2 {

    vec2  _m[2];

   public:
    //
    //  --- Constructors and Destructors ---
    //

    mat2( const GLfloat d = GLfloat(1.0) )  // Create a diagional matrix
	{ _m[0].x = d;  _m[1].y = d;   }

    mat2( const vec2& a, const vec2& b )
	{ _m[0] = a;  _m[1] = b;  }

    mat2( GLfloat m00, GLfloat m10, GLfloat m01, GLfloat m11 )   //YJC: These 4 items are given in *column order,
                                                                 //     but the matrix is stored in *row order*.
      { _m[0] = vec2( m00, m01 ); _m[1] = vec2( m10, m11 ); }    //YJC: This is in row order.

    mat2( const mat2& m ) {
	if ( *this != m ) {
	    _m[0] = m._m[0];
	    _m[1] = m._m[1];
	} 
    }

    //
    //  --- Indexing Operator ---
    //

    vec2& operator [] ( int i ) { return _m[i]; }
    const vec2& operator [] ( int i ) const { return _m[i]; }

    //
    //  --- (non-modifying) Arithmatic Operators ---
    //

    mat2 operator + ( const mat2& m ) const
	{ return mat2( _m[0]+m[0], _m[1]+m[1] ); }

    mat2 operator - ( const mat2& m ) const
	{ return mat2( _m[0]-m[0], _m[1]-m[1] ); }

    mat2 operator * ( const GLfloat s ) const 
	{ return mat2( s*_m[0], s*_m[1] ); }

    mat2 operator / ( const GLfloat s ) const {
#ifdef DEBUG
	if ( std::fabs(s) < DivideByZeroTolerance ) {
	    std::cerr << "[" << __FILE__ << ":" << __LINE__ << "] "
		      << "Division by zero" << std::endl;
	    return mat2();
	}
#endif // DEBUG
	
	GLfloat r = GLfloat(1.0) / s;
	return *this * r;
    }

    friend mat2 operator * ( const GLfloat s, const mat2& m )
	{ return m * s; }
	
    mat2 operator * ( const mat2& m ) const {
	mat2  a( 0.0 );

	for ( int i = 0; i < 2; ++i ) {
	    for ( int j = 0; j < 2; ++j ) {
		for ( int k = 0; k < 2; ++k ) {
		    a[i][j] += _m[i][k] * m[k][j];
		}
	    }
	}

	return a;
    }

    //
    //  --- (modifying) Arithmetic Operators ---
    //

    mat2& operator += ( const mat2& m ) {
	_m[0] += m[0];  _m[1] += m[1];  
	return *this;
    }

    mat2& operator -= ( const mat2& m ) {
	_m[0] -= m[0];  _m[1] -= m[1];  
	return *this;
    }

    mat2& operator *= ( const GLfloat s ) {
	_m[0] *= s;  _m[1] *= s;   
	return *this;
    }

    mat2& operator *= ( const mat2& m ) {
	mat2  a( 0.0 );

	for ( int i = 0; i < 2; ++i ) {
	    for ( int j = 0; j < 2; ++j ) {
		for ( int k = 0; k < 2; ++k ) {
		    a[i][j] += _m[i][k] * m[k][j];
		}
	    }
	}

	return *this = a;
    }
    
    mat2& operator /= ( const GLfloat s ) {
#ifdef DEBUG
	if ( std::fabs(s) < DivideByZeroTolerance ) {
	    std::cerr << "[" << __FILE__ << ":" << __LINE__ << "] "
		      << "Division by zero" << std::endl;
	    return mat2();
	}
#endif // DEBUG

	GLfloat r = GLfloat(1.0) / s;
	return *this *= r;
    }

    //
    //  --- Matrix / Vector operators ---
    //

    vec2 operator * ( const vec2& v ) const {  // m * v
	return vec2( _m[0][0]*v.x + _m[0][1]*v.y,
		     _m[1][0]*v.x + _m[1][1]*v.y );
    }
	
    //
    //  --- Insertion and Extraction Operators ---
    //
	
    friend std::ostream& operator << ( std::ostream& os, const mat2& m )
	{ return os << std::endl << m[0] << std::endl << m[1] << std::endl; }

    friend std::istream& operator >> ( std::istream& is, mat2& m )
	{ return is >> m._m[0] >> m._m[1] ; }

    //
    //  --- Conversion Operators ---
    //

    operator const GLfloat* () const
	{ return static_cast<const GLfloat*>( &_m[0].x ); }

    operator GLfloat* ()
	{ return static_cast<GLfloat*>( &_m[0].x ); }
};

//
//  --- Non-class mat2 Methods ---
//

inline
mat2 matrixCompMult( const mat2& A, const mat2& B ) {
    return mat2( A[0][0]*B[0][0], A[0][1]*B[0][1],
		 A[1][0]*B[1][0], A[1][1]*B[1][1] );
}

//YJC: The following is *Wrong* (the 4 items must be given in *column order*).
//     ==> Use transpose1() below instead.
inline
mat2 transpose( const mat2& A ) {
    return mat2( A[0][0], A[1][0],
		 A[0][1], A[1][1] );
}

//YJC: The above transpose() is *Wrong*.
//     ==> Use the following transpose1() instead.
inline
mat2 transpose1( const mat2& A ) {
    return mat2( A[0][0], A[0][1],
		 A[1][0], A[1][1] );  //YJC: Important!
                                      //     These 4 items must be given in *column order*,
                                      //     so in this way we get the transpose.
}

//----------------------------------------------------------------------------
//
//  mat3 - 3D square matrix 
//

class mat3 {

    vec3  _m[3];

   public:
    //
    //  --- Constructors and Destructors ---
    //

    mat3( const GLfloat d = GLfloat(1.0) )  // Create a diagional matrix
	{ _m[0].x = d;  _m[1].y = d;  _m[2].z = d;   }

    mat3( const vec3& a, const vec3& b, const vec3& c )
	{ _m[0] = a;  _m[1] = b;  _m[2] = c;  }

    mat3( GLfloat m00, GLfloat m10, GLfloat m20,
	  GLfloat m01, GLfloat m11, GLfloat m21,
	  GLfloat m02, GLfloat m12, GLfloat m22 ) //YJC: These 9 items are given in *column order*,
                                                  //     but the matrix is stored in *row order*.
	{
	    _m[0] = vec3( m00, m01, m02 );        //YJC: This is in row order.
	    _m[1] = vec3( m10, m11, m12 );
	    _m[2] = vec3( m20, m21, m22 );
	}

    mat3( const mat3& m )
	{
	    if ( *this != m ) {
		_m[0] = m._m[0];
		_m[1] = m._m[1];
		_m[2] = m._m[2];
	    } 
	}

    //
    //  --- Indexing Operator ---
    //

    vec3& operator [] ( int i ) { return _m[i]; }
    const vec3& operator [] ( int i ) const { return _m[i]; }

    //
    //  --- (non-modifying) Arithmatic Operators ---
    //

    mat3 operator + ( const mat3& m ) const
	{ return mat3( _m[0]+m[0], _m[1]+m[1], _m[2]+m[2] ); }

    mat3 operator - ( const mat3& m ) const
	{ return mat3( _m[0]-m[0], _m[1]-m[1], _m[2]-m[2] ); }

    mat3 operator * ( const GLfloat s ) const 
	{ return mat3( s*_m[0], s*_m[1], s*_m[2] ); }

    mat3 operator / ( const GLfloat s ) const {
#ifdef DEBUG
	if ( std::fabs(s) < DivideByZeroTolerance ) {
	    std::cerr << "[" << __FILE__ << ":" << __LINE__ << "] "
		      << "Division by zero" << std::endl;
	    return mat3();
	}
#endif // DEBUG
	
	GLfloat r = GLfloat(1.0) / s;
	return *this * r;
    }

    friend mat3 operator * ( const GLfloat s, const mat3& m )
	{ return m * s; }
	
    mat3 operator * ( const mat3& m ) const {
	mat3  a( 0.0 );

	for ( int i = 0; i < 3; ++i ) {
	    for ( int j = 0; j < 3; ++j ) {
		for ( int k = 0; k < 3; ++k ) {
		    a[i][j] += _m[i][k] * m[k][j];
		}
	    }
	}

	return a;
    }

    //
    //  --- (modifying) Arithmetic Operators ---
    //

    mat3& operator += ( const mat3& m ) {
	_m[0] += m[0];  _m[1] += m[1];  _m[2] += m[2]; 
	return *this;
    }

    mat3& operator -= ( const mat3& m ) {
	_m[0] -= m[0];  _m[1] -= m[1];  _m[2] -= m[2]; 
	return *this;
    }

    mat3& operator *= ( const GLfloat s ) {
	_m[0] *= s;  _m[1] *= s;  _m[2] *= s; 
	return *this;
    }

    mat3& operator *= ( const mat3& m ) {
	mat3  a( 0.0 );

	for ( int i = 0; i < 3; ++i ) {
	    for ( int j = 0; j < 3; ++j ) {
		for ( int k = 0; k < 3; ++k ) {
		    a[i][j] += _m[i][k] * m[k][j];
		}
	    }
	}

	return *this = a;
    }

    mat3& operator /= ( const GLfloat s ) {
#ifdef DEBUG
	if ( std::fabs(s) < DivideByZeroTolerance ) {
	    std::cerr << "[" << __FILE__ << ":" << __LINE__ << "] "
		      << "Division by zero" << std::endl;
	    return mat3();
	}
#endif // DEBUG

	GLfloat r = GLfloat(1.0) / s;
	return *this *= r;
    }

    //
    //  --- Matrix / Vector operators ---
    //

    vec3 operator * ( const vec3& v ) const {  // m * v
	return vec3( _m[0][0]*v.x + _m[0][1]*v.y + _m[0][2]*v.z,
		     _m[1][0]*v.x + _m[1][1]*v.y + _m[1][2]*v.z,
		     _m[2][0]*v.x + _m[2][1]*v.y + _m[2][2]*v.z );
    }
	
    //
    //  --- Insertion and Extraction Operators ---
    //
	
    friend std::ostream& operator << ( std::ostream& os, const mat3& m ) {
	return os << std::endl 
		  << m[0] << std::endl
		  << m[1] << std::endl
		  << m[2] << std::endl;
    }

    friend std::istream& operator >> ( std::istream& is, mat3& m )
	{ return is >> m._m[0] >> m._m[1] >> m._m[2] ; }

    //
    //  --- Conversion Operators ---
    //

    operator const GLfloat* () const
	{ return static_cast<const GLfloat*>( &_m[0].x ); }

    operator GLfloat* ()
	{ return static_cast<GLfloat*>( &_m[0].x ); }
};

//
//  --- Non-class mat3 Methods ---
//

inline
mat3 matrixCompMult( const mat3& A, const mat3& B ) {
    return mat3( A[0][0]*B[0][0], A[0][1]*B[0][1], A[0][2]*B[0][2],
		 A[1][0]*B[1][0], A[1][1]*B[1][1], A[1][2]*B[1][2],
		 A[2][0]*B[2][0], A[2][1]*B[2][1], A[2][2]*B[2][2] );
}

//YJC: This transpose() is *incorrect*.
//     ==> Use transpose1() below instead.
inline
mat3 transpose( const mat3& A ) {
    return mat3( A[0][0], A[1][0], A[2][0],
		 A[0][1], A[1][1], A[2][1],
		 A[0][2], A[1][2], A[2][2] );
}

// YJC: The above transpose() function is *incorrect*.
//      ==> Use the following transpose1() instead.
inline
mat3 transpose1( const mat3& A ) {
    return mat3( A[0][0], A[0][1], A[0][2],
		 A[1][0], A[1][1], A[1][2],
		 A[2][0], A[2][1], A[2][2] );    //YJC: Important!
                                                 //     The 9 items must be given in *column order*,
                                                 //     so in this way we get the transpose.
}

///////////////////////////////////////////////////////////////
// YJC: Added the following:
//      inverse(): return the inverse of the given 3x3 matrix m
inline
mat3 inverse( const mat3& m ) {
  mat3 r;

  float det = m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1]) +
              m[0][1] * (m[1][2] * m[2][0] - m[1][0] * m[2][2]) +
              m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0]);


  // check if non-singular matrix
  if (std::abs(det) < (1e-8) * (1e-8))
    { printf("Error! Matrix Determinant is too close to 0!\n");
      exit(-1);
    }  

  // Using the following "assert()" command would need to include additional header:
  // #include <cassert>  // or #include <assert.h>
  // 
  // assert(std::abs(det) > (1e-8) * (1e-8));

  r[0][0] =  (m[1][1] * m[2][2] - m[1][2] * m[2][1]) / det;
  r[1][0] = -(m[1][0] * m[2][2] - m[1][2] * m[2][0]) / det;
  r[2][0] =  (m[1][0] * m[2][1] - m[1][1] * m[2][0]) / det;
  r[0][1] = -(m[0][1] * m[2][2] - m[0][2] * m[2][1]) / det;
  r[1][1] =  (m[0][0] * m[2][2] - m[0][2] * m[2][0]) / det;
  r[2][1] = -(m[0][0] * m[2][1] - m[0][1] * m[2][0]) / det;
  r[0][2] =  (m[0][1] * m[1][2] - m[0][2] * m[1][1]) / det;
  r[1][2] = -(m[0][0] * m[1][2] - m[0][2] * m[1][0]) / det;
  r[2][2] =  (m[0][0] * m[1][1] - m[0][1] * m[1][0]) / det;

  return r;
}

//----------------------------------------------------------------------------
//
//  mat4.h - 4D square matrix
//

class mat4 {

    vec4  _m[4];

   public:
    //
    //  --- Constructors and Destructors ---
    //

    mat4( const GLfloat d = GLfloat(1.0) )  // Create a diagional matrix
	{ _m[0].x = d;  _m[1].y = d;  _m[2].z = d;  _m[3].w = d; }

    mat4( const vec4& a, const vec4& b, const vec4& c, const vec4& d )
	{ _m[0] = a;  _m[1] = b;  _m[2] = c;  _m[3] = d; }
            //
           // YJC: a becomes the first row, b the 2nd row,
           //      c the 3rd row, d the 4th row.

    mat4( GLfloat m00, GLfloat m10, GLfloat m20, GLfloat m30,
	  GLfloat m01, GLfloat m11, GLfloat m21, GLfloat m31,
	  GLfloat m02, GLfloat m12, GLfloat m22, GLfloat m32,
	  GLfloat m03, GLfloat m13, GLfloat m23, GLfloat m33 )
            //
            //YJC: These 16 items are given in *column order*,
            //     but the matrix is stored in *row order*.
            //
	{
	    _m[0] = vec4( m00, m01, m02, m03 );  //YJC: This is in row order:
	    _m[1] = vec4( m10, m11, m12, m13 );  //     _m[0] is the first row,
	    _m[2] = vec4( m20, m21, m22, m23 );  //     _m[1] the 2nd row, etc.
	    _m[3] = vec4( m30, m31, m32, m33 );
	}

    mat4( const mat4& m )
	{
	    if ( *this != m ) {
		_m[0] = m._m[0];
		_m[1] = m._m[1];
		_m[2] = m._m[2];
		_m[3] = m._m[3];
	    } 
	}

    //
    //  --- Indexing Operator ---
    //

    vec4& operator [] ( int i ) { return _m[i]; }
    const vec4& operator [] ( int i ) const { return _m[i]; }

    //
    //  --- (non-modifying) Arithematic Operators ---
    //

    mat4 operator + ( const mat4& m ) const
	{ return mat4( _m[0]+m[0], _m[1]+m[1], _m[2]+m[2], _m[3]+m[3] ); }

    mat4 operator - ( const mat4& m ) const
	{ return mat4( _m[0]-m[0], _m[1]-m[1], _m[2]-m[2], _m[3]-m[3] ); }

    mat4 operator * ( const GLfloat s ) const 
	{ return mat4( s*_m[0], s*_m[1], s*_m[2], s*_m[3] ); }

    mat4 operator / ( const GLfloat s ) const {
#ifdef DEBUG
	if ( std::fabs(s) < DivideByZeroTolerance ) {
	    std::cerr << "[" << __FILE__ << ":" << __LINE__ << "] "
		      << "Division by zero" << std::endl;
	    return mat4();
	}
#endif // DEBUG
	
	GLfloat r = GLfloat(1.0) / s;
	return *this * r;
    }

    friend mat4 operator * ( const GLfloat s, const mat4& m )
	{ return m * s; }
	
    mat4 operator * ( const mat4& m ) const {
	mat4  a( 0.0 );

	for ( int i = 0; i < 4; ++i ) {
	    for ( int j = 0; j < 4; ++j ) {
		for ( int k = 0; k < 4; ++k ) {
		    a[i][j] += _m[i][k] * m[k][j];
		}
	    }
	}

	return a;
    }

    //
    //  --- (modifying) Arithematic Operators ---
    //

    mat4& operator += ( const mat4& m ) {
	_m[0] += m[0];  _m[1] += m[1];  _m[2] += m[2];  _m[3] += m[3];
	return *this;
    }

    mat4& operator -= ( const mat4& m ) {
	_m[0] -= m[0];  _m[1] -= m[1];  _m[2] -= m[2];  _m[3] -= m[3];
	return *this;
    }

    mat4& operator *= ( const GLfloat s ) {
	_m[0] *= s;  _m[1] *= s;  _m[2] *= s;  _m[3] *= s;
	return *this;
    }

    mat4& operator *= ( const mat4& m ) {
	mat4  a( 0.0 );

	for ( int i = 0; i < 4; ++i ) {
	    for ( int j = 0; j < 4; ++j ) {
		for ( int k = 0; k < 4; ++k ) {
		    a[i][j] += _m[i][k] * m[k][j];
		}
	    }
	}

	return *this = a;
    }

    mat4& operator /= ( const GLfloat s ) {
#ifdef DEBUG
	if ( std::fabs(s) < DivideByZeroTolerance ) {
	    std::cerr << "[" << __FILE__ << ":" << __LINE__ << "] "
		      << "Division by zero" << std::endl;
	    return mat4();
	}
#endif // DEBUG

	GLfloat r = GLfloat(1.0) / s;
	return *this *= r;
    }

    //
    //  --- Matrix / Vector operators ---
    //

    vec4 operator * ( const vec4& v ) const {  // m * v
	return vec4( _m[0][0]*v.x + _m[0][1]*v.y + _m[0][2]*v.z + _m[0][3]*v.w,
		     _m[1][0]*v.x + _m[1][1]*v.y + _m[1][2]*v.z + _m[1][3]*v.w,
		     _m[2][0]*v.x + _m[2][1]*v.y + _m[2][2]*v.z + _m[2][3]*v.w,
		     _m[3][0]*v.x + _m[3][1]*v.y + _m[3][2]*v.z + _m[3][3]*v.w
	    );
    }
	
    //
    //  --- Insertion and Extraction Operators ---
    //
	
    friend std::ostream& operator << ( std::ostream& os, const mat4& m ) {
	return os << std::endl 
		  << m[0] << std::endl
		  << m[1] << std::endl
		  << m[2] << std::endl
		  << m[3] << std::endl;
    }

    friend std::istream& operator >> ( std::istream& is, mat4& m )
	{ return is >> m._m[0] >> m._m[1] >> m._m[2] >> m._m[3]; }

    //
    //  --- Conversion Operators ---
    //

    operator const GLfloat* () const
	{ return static_cast<const GLfloat*>( &_m[0].x ); }

    operator GLfloat* ()
	{ return static_cast<GLfloat*>( &_m[0].x ); }
};

//
//  --- Non-class mat4 Methods ---
//

inline
mat4 matrixCompMult( const mat4& A, const mat4& B ) {
    return mat4(
	A[0][0]*B[0][0], A[0][1]*B[0][1], A[0][2]*B[0][2], A[0][3]*B[0][3],
	A[1][0]*B[1][0], A[1][1]*B[1][1], A[1][2]*B[1][2], A[1][3]*B[1][3],
	A[2][0]*B[2][0], A[2][1]*B[2][1], A[2][2]*B[2][2], A[2][3]*B[2][3],
	A[3][0]*B[3][0], A[3][1]*B[3][1], A[3][2]*B[3][2], A[3][3]*B[3][3] );
}

//YJC: This transpose() function has No effect and is *incorrect*.
//     ==> Use transpose1() below instead.
inline
mat4 transpose( const mat4& A ) {
    return mat4( A[0][0], A[1][0], A[2][0], A[3][0],
		 A[0][1], A[1][1], A[2][1], A[3][1],
		 A[0][2], A[1][2], A[2][2], A[3][2],
		 A[0][3], A[1][3], A[2][3], A[3][3] );
}

// YJC: The above transpose() function has NO effect and is *incorrect*.
//      ==> Use the following transpose1() instead.
//          In particular this is to be used in the function Rotate().
inline
mat4 transpose1( const mat4& A ) {
    return mat4( A[0][0], A[0][1], A[0][2], A[0][3],
		 A[1][0], A[1][1], A[1][2], A[1][3],
		 A[2][0], A[2][1], A[2][2], A[2][3],
		 A[3][0], A[3][1], A[3][2], A[3][3] ); //YJC: Important!!
                                                       //     The 16 items must be given in *column order*,
                                                       //     so in this way we get the transpose.
}

//////////////////////////////////////////////////////////////////////////////
//
//  Helpful Matrix Methods
//
//////////////////////////////////////////////////////////////////////////////

#define Error( str ) do { std::cerr << "[" __FILE__ ":" << __LINE__ << "] " \
				    << str << std::endl; } while(0)

inline
vec4 mvmult( const mat4& a, const vec4& b )
{
    Error( "replace with vector matrix multiplcation operator" );

    vec4 c;
    int i, j;
    for(i=0; i<4; i++) {
	c[i] =0.0;
	for(j=0;j<4;j++) c[i]+=a[i][j]*b[j];
    }
    return c;
}

//----------------------------------------------------------------------------
//
//  Rotation matrix generators
//

inline
mat4 RotateX( const GLfloat theta )
{
    GLfloat angle = DegreesToRadians * theta;

    mat4 c;
    c[2][2] = c[1][1] = cos(angle);
    c[2][1] = sin(angle);
    c[1][2] = -c[2][1];
    return c;
}

inline
mat4 RotateY( const GLfloat theta )
{
    GLfloat angle = DegreesToRadians * theta;

    mat4 c;
    c[2][2] = c[0][0] = cos(angle);
    c[0][2] = sin(angle);
    c[2][0] = -c[0][2];
    return c;
}

inline
mat4 RotateZ( const GLfloat theta )
{
    GLfloat angle = DegreesToRadians * theta;

    mat4 c;
    c[0][0] = c[1][1] = cos(angle);
    c[1][0] = sin(angle);
    c[0][1] = -c[1][0];
    return c;
}

/*** YJC: Construct a matrix for general rotation
          This is taken and modified from "vmath.h" file of Red Book 8th Ed. in include/ dir
          Note: The rotation axis vector (x, y, z) can have length != 1.0.
*/
inline
mat4 Rotate(const GLfloat angle, const GLfloat x, const GLfloat y, const GLfloat z)
{
    mat4 result;

    //YJC: normalize (x, y, z) to a unit-length vector (x1, y1, z1) and use the latter.
    float len = sqrt(x * x + y * y + z * z);
    float x1, y1, z1;
    x1 = x; y1 = y; z1 = z;
    if (len < 0.00001)
      { printf("Error! Rotation axis vector is too close to (0,0,0)\n");
	exit(-1);
      }
    if (len != 1.0) { x1 = x1/len; y1 = y1/len; z1 = z1/len;}

    const GLfloat x2 = x1 * x1;
    const GLfloat y2 = y1 * y1;
    const GLfloat z2 = z1 * z1;

    float rads = float(angle) * 0.0174532925f;
    const float c = cosf(rads);
    const float s = sinf(rads);
    const float omc = 1.0f - c;

    result[0] = vec4(x2 * omc + c, y1 * x1 * omc + z1 * s, x1 * z1 * omc - y1 * s, 0.0);
    result[1] = vec4(x1 * y1 * omc - z1 * s, y2 * omc + c, y1 * z1 * omc + x1 * s, 0.0);
    result[2] = vec4(x1 * z1 * omc + y1 * s, y1 * z1 * omc - x1 * s, z2 * omc + c, 0.0);
    result[3] = vec4(0.0, 0.0, 0.0, 1.0);

    /*** YJC: compared with "glMatrixEA.js-YJC" mat4.rotate, the matrix "result" above is the same
              and is in *Column Order* ("vmath.h" file also indicates this, saying that the matrix
              is "Column primary data (essentially, array of vectors)").
         ===> Return transpose1(result) to make it *Row Order*, to be consistent with other functions here.
              (Note: we use transpose1() instead of transpose(); the latter is incorrect.)
    ***/
    return transpose1(result); 
    // return result;  /* Original */
}

//----------------------------------------------------------------------------
//
//  Translation matrix generators
//

inline
mat4 Translate( const GLfloat x, const GLfloat y, const GLfloat z )
{
    mat4 c;
    c[0][3] = x;
    c[1][3] = y;
    c[2][3] = z;
    return c;
}

inline
mat4 Translate( const vec3& v )
{
    return Translate( v.x, v.y, v.z );
}

inline
mat4 Translate( const vec4& v )
{
    return Translate( v.x, v.y, v.z );
}

//----------------------------------------------------------------------------
//
//  Scale matrix generators
//

inline
mat4 Scale( const GLfloat x, const GLfloat y, const GLfloat z )
{
    mat4 c;
    c[0][0] = x;
    c[1][1] = y;
    c[2][2] = z;
    return c;
}

inline
mat4 Scale( const vec3& v )
{
    return Scale( v.x, v.y, v.z );
}

//----------------------------------------------------------------------------
//
//  Projection transformation matrix geneartors
//
//    Note: Microsoft Windows (r) defines the keyword "far" in C/C++.  In
//          order to avoid any name conflicts, we use the variable names
//          "zNear" to reprsent "near", and "zFar" to reprsent "far".
//



inline
mat4 Ortho( const GLfloat left, const GLfloat right,
	    const GLfloat bottom, const GLfloat top,
	    const GLfloat zNear, const GLfloat zFar )
{
    mat4 c;
    c[0][0] = 2.0f / (right - left);
    c[1][1] = 2.0f / (top - bottom);
    c[2][2] = 2.0f / (zNear - zFar);
    c[3][3] = 1.0f;
    c[0][3] = -(right + left)/(right - left);
    c[1][3] = -(top + bottom)/(top - bottom);
    c[2][3] = -(zFar + zNear)/(zFar - zNear);
    return c;
}

inline
mat4 Ortho2D( const GLfloat left, const GLfloat right,
	      const GLfloat bottom, const GLfloat top )
{
    return Ortho( left, right, bottom, top, -1.0, 1.0 );
}

inline
mat4 Frustum( const GLfloat left, const GLfloat right,
	      const GLfloat bottom, const GLfloat top,
	      const GLfloat zNear, const GLfloat zFar )
{
    mat4 c;
    c[0][0] = 2.0f*zNear/(right - left);
    c[0][2] = (right + left)/(right - left);
    c[1][1] = 2.0f*zNear/(top - bottom);
    c[1][2] = (top + bottom)/(top - bottom);
    c[2][2] = -(zFar + zNear)/(zFar - zNear);
    c[2][3] = -2.0f*zFar*zNear/(zFar - zNear);
    c[3][2] = -1.0f;
    return c;
}

inline
mat4 Perspective( const GLfloat fovy, const GLfloat aspect,
		  const GLfloat zNear, const GLfloat zFar)
{
    GLfloat top   = tan(fovy*DegreesToRadians/2) * zNear;
    GLfloat right = top * aspect;

    mat4 c;
    c[0][0] = zNear/right;
    c[1][1] = zNear/top;
    c[2][2] = -(zFar + zNear)/(zFar - zNear);
    c[2][3] = -2.0f*zFar*zNear/(zFar - zNear);
    c[3][2] = -1.0f;
    return c;
}

//----------------------------------------------------------------------------
//
//  Viewing transformation matrix generation
//

inline
mat4 LookAt( const vec4& eye, const vec4& at, const vec4& up )
{
    vec4 n = normalize(eye - at);

    //YJC: Corrected below
    vec4 u = normalize( vec4(cross(up,n),0) );
    vec4 v = normalize( vec4(cross(n,u),0) );
    //
    //YJC: Original below --- Incorrect since vec4() construct will
    //     set the w component (the 4th component) to 1 (see "vec.h").
    //     But u, v are each a basis vector and w should be 0, as above.
    // vec4 u = normalize( cross(up,n) );
    // vec4 v = normalize( cross(n,u)  );
    
    vec4 t = vec4(0.0, 0.0, 0.0, 1.0);
    mat4 c = mat4(u, v, n, t);
    return c * Translate( -eye );
}

//---------------------------------------------------------------------------
// YJC: Added the following:
//      upperLeftMat3(m): Return the upper-left 3x3 submatrix of the given mat4 m
//
inline
mat3 upperLeftMat3( const mat4& m ) {
  
  return mat3(vec3(m[0][0], m[0][1], m[0][2]),
              vec3(m[1][0], m[1][1], m[1][2]),
              vec3(m[2][0], m[2][1], m[2][2]));
}        

// YJC: Added the following:
//      NormalMatrix(mv, non_uniform_scale_flag):
//
//      Return the Normal Matrix (mat3) given the Model-View matrix mv (mat4).
//      * non_uniform_scale_flag == 1 if mv involves non-uniform scaling
//        non_uniform_scale_flag == 0 otherwise      
inline
mat3 NormalMatrix( const mat4& mv, int non_uniform_scale_flag ) {
  
  // The upper-left 3x3 submatrix of mv
  mat3 m = upperLeftMat3( mv );

  if (non_uniform_scale_flag == 0) // No non-uniform scaling is involved
    return m;

  else // mv involves non-uniform scaling ==> return the transpose of inverse(m)
    return transpose1( inverse(m) );
}

// YJC: Added the following:
//      mat4WithUpperLeftMat3(m)
//
//      Return the 4x4 matrix where the upper-left 3x3 submatrix is m,
//      the 4th column is (0,0,0,1) --- i.e., No translation,
//      and the 4th row is (0,0,0,1).
//      (Basically this is a Model-View matrix with No translation and
//       the rotation and scaling parts are specified by m.)
inline
mat4 mat4WithUpperLeftMat3( const mat3& m){
     
  return mat4( vec4(m[0][0], m[0][1], m[0][2], 0.0),
               vec4(m[1][0], m[1][1], m[1][2], 0.0),
               vec4(m[2][0], m[2][1], m[2][2], 0.0),
               vec4(    0.0,     0.0,     0.0, 1.0) );
}

//----------------------------------------------------------------------------

inline
vec4 minus(const vec4& a, const vec4&  b )
{
    Error( "replace with vector subtraction" );
    return vec4(a[0]-b[0], a[1]-b[1], a[2]-b[2], 0.0);
}

inline
void printv(const vec4& a )
{
    Error( "replace with vector insertion operator" );
    printf("%f %f %f %f \n\n", a[0], a[1], a[2], a[3]);
}

inline
void printm(const mat4 a)
{
    Error( "replace with matrix insertion operator" );
    for(int i=0; i<4; i++) printf("%f %f %f %f \n", a[i][0], a[i][1], a[i][2], a[i][3]);
    printf("\n");
}

inline
mat4 identity()
{
    Error( "replace with either a matrix constructor or identity method" );
    mat4 c;
    for(int i=0; i<4; i++) for(int j=0; j<4; j++) c[i][j]=0.0;
    for(int i=0; i<4; i++) c[i][i] = 1.0;
    return c;
}


}  // namespace Angel

#endif // __ANGEL_MAT_H__