Quaternion for rotation on X and Y

I’m trying to rotate a figure by global X and Y axes.

Rotate doesn’t work because it works over local axes.

Well I implemented quaternion, I do the nexts steps:

1.Create the quaternion representation of angles and axis, for each axe.
2.I want to rotate over X and Y axes, then mulpiple the vector quaternion of each axe.
3.The result converts to matrix form.
4.LoadMatrix / MultMatrix.

My Code:

[NOTE]gl.PushMatrix();
q = qxq(bx, 1, 0, 0, by, 0, 1, 0);
m = mat(q[0], q[1], q[2], q[3]);

        gl.LoadIdentity();
        gl.MultMatrix(m);

                                [INDENT]gl.Begin(OpenGL.GL_TRIANGLES);
                                gl.Color(1.0f, 0.0f, 0.0f);
                                gl.Vertex(0.0f, 1.0f, 0.0f);
                                gl.Vertex(-1.0f, -1.0f, 1.0f);
                                gl.Vertex(1.0f, -1.0f, 1.0f);

                                gl.Color(0.0f, 1.0f, 0.0f);
                                gl.Vertex(0.0f, 1.0f, 0.0f);
                                gl.Vertex(1.0f, -1.0f, 1.0f);
                                gl.Vertex(1.0f, -1.0f, -1.0f);

                                gl.Color(0.0f, 0.0f, 1.0f);
                                gl.Vertex(0.0f, 1.0f, 0.0f);
                                gl.Vertex(1.0f, -1.0f, -1.0f);
                                gl.Vertex(-1.0f, -1.0f, -1.0f);

                                gl.Color(1.0f, 1.0f, 1.0f);
                                gl.Vertex(0.0f, 1.0f, 0.0f);
                                gl.Vertex(-1.0f, -1.0f, -1.0f);
                                gl.Vertex(-1.0f, -1.0f, 1.0f);
                                gl.End();
                                //Axe X
                                gl.Begin(OpenGL.GL_LINES);
                                gl.Color(1.0, 0.698, 0.698);
                                gl.Vertex(-50.0, 0.0, 0.0);
                                gl.Vertex(50.0, 0.0, 0.0);
                                gl.End();

                                //Axe Y
                                gl.Begin(OpenGL.GL_LINES);
                                gl.Color(0.698, 1.0, 0.698);
                                gl.Vertex(0.0, -50.0, 0);
                                gl.Vertex(0.0, 50.0, 0);
                                gl.End();

                                //Axe Z
                                gl.Begin(OpenGL.GL_LINES);
                                gl.Color(0.698, 0.698, 1.0);
                                gl.Vertex(0.0, 0.0, -50.0);
                                gl.Vertex(0.0, 0.0, 50.0);
                                gl.End();   [/INDENT]                               

gl.PopMatrix();[/NOTE]

bx and by are the variables that I use to increment or decremente the angle for each axe.

The method qxq(), I use to read the angles of axes. Here I calculate the quanterion vector for each axe and I do the multiplicaction:

[NOTE]private double[] qxq(double w1, double x1, double y1, double z1, double w2, double x2, double y2, double z2)
{
double[] q = new double[4];
double qx1 = 0.0f;
double qy1 = 0.0f;
double qz1 = 0.0f;
double qw1 = 0.0f;

    double qx2 = 0.0f;
    double qy2 = 0.0f;
    double qz2 = 0.0f;
    double qw2 = 0.0f;

    qw1 = Math.Cos(w1 * (Math.PI / 180) / 2);
    qx1 = Math.Sin(w1 * (Math.PI / 180) / 2) * x1;
    qy1 = Math.Sin(w1 * (Math.PI / 180) / 2) * y1;
    qz1 = Math.Sin(w1 * (Math.PI / 180) / 2) * z1;

    qw2 = Math.Cos(w2 * (Math.PI / 180) / 2);
    qx2 = Math.Sin(w2 * (Math.PI / 180) / 2) * x2;
    qy2 = Math.Sin(w2 * (Math.PI / 180) / 2) * y2;
    qz2 = Math.Sin(w2 * (Math.PI / 180) / 2) * z2;

    q[0] = qw1*qw2 - qx1*qx2 - qy1*qy2 - qz1*qz2;
    q[1] = qw1*qx2 + qx1*qw2 + qy1*qz2 - qz1*qy2;
    q[2] = qw1*qy2 + qy1*qw2 + qz1*qx2 - qx1*qz2;
    q[3] = qw1*qz2 + qz1*qw2 + qx1*qy2 - qy1*qx2;

    return q;
}[/NOTE]

In the method mat() I transform my resulting vector of qxq() to a matrix form:

[NOTE]private double[] mat(double w, double x, double y, double z)
{
double[] m = new double[16];
double qx = 0.0f;
double qy = 0.0f;
double qz = 0.0f;
double qw = 0.0f;

    qw = Math.Cos(w / 2);
    qx = Math.Sin(w / 2) * x;
    qy = Math.Sin(w / 2) * y;
    qz = Math.Sin(w / 2) * z;

    m[0] = 1.0f - 2.0f * qy * qy - 2.0f * qz * qz;
    m[1] = 0.0f + 2.0f * qx * qy + 2.0f * qw * qz;
    m[2] = 0.0f + 2.0f * qx * qz - 2.0f * qw * qy;
    m[3] = 0.0f;

    m[4] = 0.0f + 2.0f * qx * qy - 2.0f * qw * qz;
    m[5] = 1.0f - 2.0f * qx * qx - 2.0f * qz * qz;
    m[6] = 0.0f + 2.0f * qy * qz + 2.0f * qw * qx;
    m[7] = 0.0f;

    m[8] = 0.0f + 2.0f * qx * qz + 2.0f * qw * qy;
    m[9] = 0.0f + 2.0f * qy * qz - 2.0f * qw * qx;
    m[10] = 1.0f - 2.0f * qx * qx - 2.0f * qy * qy;
    m[11] = 0.0f;

    m[12] = 0.0f;
    m[13] = 0.0f;
    m[14] = 0.0f;
    m[15] = 1.0f;

    return m;
}[/NOTE]

Well I think my compute calculations are ok, but I can’t do that this works in the 2 axes together.

If I don’t make the multiplicaction and I call mat() separate for each axe, and if I incremente for example angle by X, it works; but if after that, I increment the angle by Y, the figure reset to the original position (the rotation on x that I did disappear), and do the rotation by Y. And if after I increment angle by X, the figure returns to the position as before the rotation on Y.

This happens when:

[NOTE]m = mat(angle, x, y, z);
gl.LoadIdentity();
gl.LoadMatrix(m);[/NOTE]

I change the value of angle and (x,y,z) to indicate the angle for one axe. For example

m = mat(45, 0, 1, 0); -------> 45° by Y

m = mat(45, 1, 0, 0); -------> 45° by x

Well , I just want to rotate over the global axes. Thanks.