I am replacing my project's use of glRotatef because I need to be able to transform double matrices. glRotated is not an option because OpenGL does not guarantee the stored matrices or any operations performed to be double precision. However, my new implementation only rotates around the global axes, and does not give the same result as glRotatef.
I have looked at some implementations of glRotatef (like OpenGl rotate custom implementation) and don't see how they account for the initial transformation matrix's local axes when calculating the rotation matrix.
I have a generic rotate function, taken (with some changes) from https://community.khronos.org/t/implementing-rotation-function-like-glrotate/68603:
typedef double double_matrix_t[16];
void rotate_double_matrix(const double_matrix_t in, double angle, double x, double y, double z,
double_matrix_t out)
{
double sinAngle, cosAngle;
double mag = sqrt(x * x + y * y + z * z);
sinAngle = sin ( angle * M_PI / 180.0 );
cosAngle = cos ( angle * M_PI / 180.0 );
if ( mag > 0.0f )
{
double xx, yy, zz, xy, yz, zx, xs, ys, zs;
double oneMinusCos;
double_matrix_t rotMat;
x /= mag;
y /= mag;
z /= mag;
xx = x * x;
yy = y * y;
zz = z * z;
xy = x * y;
yz = y * z;
zx = z * x;
xs = x * sinAngle;
ys = y * sinAngle;
zs = z * sinAngle;
oneMinusCos = 1.0f - cosAngle;
rotMat[0] = (oneMinusCos * xx) + cosAngle;
rotMat[4] = (oneMinusCos * xy) - zs;
rotMat[8] = (oneMinusCos * zx) + ys;
rotMat[12] = 0.0F;
rotMat[1] = (oneMinusCos * xy) + zs;
rotMat[5] = (oneMinusCos * yy) + cosAngle;
rotMat[9] = (oneMinusCos * yz) - xs;
rotMat[13] = 0.0F;
rotMat[2] = (oneMinusCos * zx) - ys;
rotMat[6] = (oneMinusCos * yz) + xs;
rotMat[10] = (oneMinusCos * zz) + cosAngle;
rotMat[14] = 0.0F;
rotMat[3] = 0.0F;
rotMat[7] = 0.0F;
rotMat[11] = 0.0F;
rotMat[15] = 1.0F;
multiply_double_matrices(in, rotMat, out); // Generic matrix multiplication function.
}
}
I call this function with the same rotations I used to call glRotatef with and in the same order, but the result is different. All rotations are done around the global axes, while glRotatef would rotate around the local axis of in.
For example, I have a plane:

and I pitch up 90 degrees (this gives the expected result with both glRotatef and my rotation function) and persist the transformation:

If I bank 90 degrees with glRotatef (glRotatef(90, 0.0f, 0.0f, 1.0f)), the plane rotates around the transformation's local Z axis pointing out of the plane's nose, which is what I want:

But if I bank 90 degrees with my code (rotate_double_matrix(in, 90.0f, 0.0, 0.0, 1.0, out)), I get this:

The plane is still rotating around the global Z axis.
Similar issues happen if I change the order of rotations - the first rotation gives the expected result, but subsequent rotations still happen around the global axes.
How does glRotatef rotate around a matrix's local axes? What do I need to change in my code to get the same result? I assume rotate_double_matrix needs to modify the x, y, z values passed in based on the in matrix somehow, but I'm not sure.