I am creating a very simple project on OpenGL and I'm stuck with rotations. I am trying to rotate an object indepentdently in all 3 axes: X, Y, and Z. I've had sleepless nights due to the "gimbal lock" problem after rotating about one axis. I've then learned that quaternions would solve my problem. I've researched about quaternions and implementd it, but I havent't been able to convert my rotations to quaternions. For example, if I want to rotate around Z axis 90 degrees, I just create the {0,0,1} vector for my quaternion and rotate it around that axis 90 degrees using the code here:

http://iphonedevelopment.blogspot.com/2009/06/opengl-es-from-ground-up-part-7_04.html (the most complicated matrix towards the bottom)

That's ok for one vector, but, say, I first want to rotate 90 degrees around Z, then 90 degrees around X (just as an example). What vector do I need to pass in? How do I calculate that vector. I am not good with matrices and trigonometry (I know the basics and the general rules, but I'm just not a whiz) but I need to get this done. There are LOTS of tutorials about quaternions, but I seem to understand none (or they don't answer my question). I just need to learn to construct the vector for rotations around more than one axis combined.

UPDATE: I've found this nice page about quaternions and decided to implement them this way: http://www.cprogramming.com/tutorial/3d/quaternions.html

Here is my code for quaternion multiplication:

void cube::quatmul(float* q1, float* q2, float* resultRef){
 float w = q1[0]*q2[0] - q1[1]*q2[1] - q1[2]*q2[2] - q1[3]*q2[3];
 float x = q1[0]*q2[1] + q1[1]*q2[0] + q1[2]*q2[3] - q1[3]*q2[2];
 float y = q1[0]*q2[2] - q1[1]*q2[3] + q1[2]*q2[0] + q1[3]*q2[1];
 float z = q1[0]*q2[3] + q1[1]*q2[2] - q1[2]*q2[1] + q1[3]*q2[0];

 resultRef[0] = w;
 resultRef[1] = x;
 resultRef[2] = y;
 resultRef[3] = z;
}

Here is my code for applying a quaternion to my modelview matrix (I have a tmodelview variable that is my target modelview matrix):

void cube::applyquat(){
 float& x = quaternion[1];
 float& y = quaternion[2];
 float& z = quaternion[3];
 float& w = quaternion[0];
 float magnitude = sqrtf(w * w + x * x + y * y + z * z);
 if(magnitude == 0){
     x = 1;
     w = y = z = 0;
 }else
 if(magnitude != 1){
     x /= magnitude;
     y /= magnitude;
     z /= magnitude;
     w /= magnitude;
  }

 tmodelview[0] = 1 - (2 * y * y) - (2 * z * z);
 tmodelview[1] = 2 * x * y + 2 * w * z;
 tmodelview[2] = 2 * x * z - 2 * w * y;
 tmodelview[3] = 0; 

 tmodelview[4] = 2 * x * y - 2 * w * z;
 tmodelview[5] = 1 - (2 * x * x) - (2 * z * z);
 tmodelview[6] = 2 * y * z - 2 * w * x;
 tmodelview[7] = 0;

 tmodelview[8] = 2 * x * z + 2 * w * y;
 tmodelview[9] = 2 * y * z + 2 * w * x;
 tmodelview[10] = 1 - (2 * x * x) - (2 * y * y);
 tmodelview[11] = 0;

 glMatrixMode(GL_MODELVIEW);
 glPushMatrix();
 glLoadMatrixf(tmodelview);
 glMultMatrixf(modelview);
 glGetFloatv(GL_MODELVIEW_MATRIX, tmodelview);
 glPopMatrix();
}

And my code for rotation (that I call externally), where quaternion is a class variable of the cube:

void cube::rotatex(int angle){
 float quat[4];
 float ang = angle * PI / 180.0;
 quat[0] = cosf(ang / 2);
 quat[1] = sinf(ang/2);
 quat[2] = 0;
 quat[3] = 0;
 quatmul(quat, quaternion, quaternion);
 applyquat();
}

void cube::rotatey(int angle){
 float quat[4];
 float ang = angle * PI / 180.0;
 quat[0] = cosf(ang / 2);
 quat[1] = 0;
 quat[2] = sinf(ang/2);
 quat[3] = 0;
 quatmul(quat, quaternion, quaternion);
 applyquat();
}

void cube::rotatez(int angle){
 float quat[4];
 float ang = angle * PI / 180.0;
 quat[0] = cosf(ang / 2);
 quat[1] = 0;
 quat[2] = 0;
 quat[3] = sinf(ang/2);
 quatmul(quat, quaternion, quaternion);
 applyquat();
}

I call, say rotatex, for 10-11 times for rotating only 1 degree, but my cube gets rotated almost 90 degrees after 10-11 times of 1 degree, which doesn't make sense. Also, after calling rotation functions in different axes, My cube gets skewed, gets 2 dimensional, and disappears (a column in modelview matrix becomes all zeros) irreversibly, which obviously shouldn't be happening with a correct implementation of the quaternions.