I am like, 3 weeks into OpenGL programming.
My first little project was a first person spaceship
flying around a small solar system.
I used heading pitch roll and ran into exactly the
problems you’re having. I switched to quaternions
and all was better.
Here are the functions I use… when the user
presses on the joystick, I multiply the current
orientation (represented as a quaternion) with
a new quaternion representing the desired
oreintation change. When I go to render,
I turn the quat to a matrix and multiply it to the
modelview.
Caveat: as I said, I’m 3 weeks into OpenGL programming, so take everything “cum grano salis”
Hope this helps
/*
 Multiply two quaternions

 res is a*b

 Assumes a and b are normalised
 res is normalised on return
*/
void QUAT_Mult(const QUAT *a, const QUAT *b, QUAT res)
{
/ Result is [ ww’  v.v’, v X v’ + wv’ + v’w ] */
/* Temp */
float t[8];
float norm;
t[0] = (a>w + a>x) * (b>w + b>x);
t[1] = (a>z  a>y) * (b>y  b>z);
t[2] = (a>x  a>w) * (b>y + b>z);
t[3] = (a>y + a>z) * (b>x  b>w);
t[4] = (a>x + a>z) * (b>x + b>y);
t[5] = (a>x  a>z) * (b>x  b>y);
t[6] = (a>w + a>y) * (b>w  b>z);
t[7] = (a>w  a>y) * (b>w + b>z);
res>w = t[1] + ((t[4]  t[5] + t[6] + t[7]) * 0.5f);
res>x = t[0]  (( t[4] + t[5] + t[6] + t[7]) * 0.5f);
res>y = t[2] + (( t[4]  t[5] + t[6]  t[7]) * 0.5f);
res>z = t[3] + (( t[4]  t[5]  t[6] + t[7]) * 0.5f);
/* Normalise result */
norm = sqrt(res>w * res>w +
res>x * res>x +
res>y * res>y +
res>z * res>z);
res>w /= norm;
res>x /= norm;
res>y /= norm;
res>z /= norm;
}
/*
 Quaternion to matrix

 Converts quaternion q into a matrix returned in m.
 m is in GL format (i.e. rowmajor format)

 Assumes q is normalised.

*/
void QUAT_ToMatrix(const QUAT *q, float m[16])
{
float wx, wy, wz, xx, yy, yz, xy, xz, zz, x2, y2, z2;
x2 = q>x + q>x;
y2 = q>y + q>y;
z2 = q>z + q>z;
xx = q>x * x2; xy = q>x * y2; xz = q>x * z2;
yy = q>y * y2; yz = q>y * z2; zz = q>z * z2;
wx = q>w * x2; wy = q>w * y2; wz = q>w * z2;
/* The following commented out stuff returns the matrix in
columnmajor format /
/
m[0] = 1.0f  (yy + zz); m[1] = xy  wz;
m[2] = xz + wy; m[3] = 0.0f;
m[4] = xy + wz; m[5] = 1.0f  (xx + zz);
m[6] = yz  wx; m[7] = 0.0f;
m[8] = xz  wy; m[9] = yz + wx;
m[10] = 1.0f  (xx + yy); m[11] = 0.0f;
m[12] = 0.0f; m[13] = 0.0f;
m[14] = 0.0f; m[15] = 1.0f;
/
/ This code returns the matrix in rowmajor format */
m[0] = 1.0f  (yy + zz); m[4] = xy  wz;
m[8] = xz + wy; m[12] = 0.0f;
m[1] = xy + wz; m[5] = 1.0f  (xx + zz);
m[9] = yz  wx; m[13] = 0.0f;
m[2] = xz  wy; m[6] = yz + wx;
m[10] = 1.0f  (xx + yy); m[14] = 0.0f;
m[3] = 0.0f; m[7] = 0.0f;
m[11] = 0.0f; m[15] = 1.0f;
}