void

_math_matrix_rotate( GLmatrix *mat,*

GLfloat angle, GLfloat x, GLfloat y, GLfloat z )

{

/ This function contributed by Erich Boleyn (erich@uruk.org) */

GLfloat mag, s, c;

GLfloat xx, yy, zz, xy, yz, zx, xs, ys, zs, one_c;

GLfloat m[16];

s = sin( angle * DEG2RAD );

c = cos( angle * DEG2RAD );

mag = GL_SQRT( x*x + y*y + z*z );

if (mag <= 1.0e-4) {

/* generate an identity matrix and return */

MEMCPY(m, Identity, sizeof(GLfloat)*16);

return;

}

x /= mag;

y /= mag;

z /= mag;

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

/*

* Arbitrary axis rotation matrix.

*

* This is composed of 5 matrices, Rz, Ry, T, Ryâ€™, Rzâ€™, multiplied

* like so: Rz * Ry * T * Ryâ€™ * Rzâ€™. T is the final rotation

* (which is about the X-axis), and the two composite transforms

* Ryâ€™ * Rzâ€™ and Rz * Ry are (respectively) the rotations necessary

* from the arbitrary axis to the X-axis then back. They are

* all elementary rotations.

*

* Rzâ€™ is a rotation about the Z-axis, to bring the axis vector

* into the x-z plane. Then Ryâ€™ is applied, rotating about the

* Y-axis to bring the axis vector parallel with the X-axis. The

* rotation about the X-axis is then performed. Ry and Rz are

* simply the respective inverse transforms to bring the arbitrary

* axis back to itâ€™s original orientation. The first transforms

* Rzâ€™ and Ryâ€™ are considered inverses, since the data from the

* arbitrary axis gives you info on how to get to it, not how

* to get away from it, and an inverse must be applied.

*

* The basic calculation used is to recognize that the arbitrary

* axis vector (x, y, z), since it is of unit length, actually

* represents the sines and cosines of the angles to rotate the

* X-axis to the same orientation, with theta being the angle about

* Z and phi the angle about Y (in the order described above)

* as follows:

*

* cos ( theta ) = x / sqrt ( 1 - z^2 )

* sin ( theta ) = y / sqrt ( 1 - z^2 )

*

* cos ( phi ) = sqrt ( 1 - z^2 )

* sin ( phi ) = z

*

* Note that cos ( phi ) can further be inserted to the above

* formulas:

*

* cos ( theta ) = x / cos ( phi )

* sin ( theta ) = y / sin ( phi )

*

* â€¦etc. Because of those relations and the standard trigonometric

* relations, it is pssible to reduce the transforms down to what

* is used below. It may be that any primary axis chosen will give the

* same results (modulo a sign convention) using thie method.

*

* Particularly nice is to notice that all divisions that might

* have caused trouble when parallel to certain planes or

* axis go away with care paid to reducing the expressions.

* After checking, it does perform correctly under all cases, since

* in all the cases of division where the denominator would have

* been zero, the numerator would have been zero as well, giving

* the expected result.

*/

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

matrix_multf( mat, m, MAT_FLAG_ROTATION );

}