This does not look right - or I really haven’t gotten what you are trying to do. I’d have to see the whole code to understand what’s going on.
What this looks like is that you are modifying the model matrix by a rotation about where the axis will end up. Which is the wrong way around. You should modify the model matrix by a rotation about where the axis is now.
If I do:
m = matrix() -- identity matrix
m = m:rotate(45,1,0,0) -- these are just to jumble up the matrix a little
m = m:rotate(45,0,1,0)
m = m:rotate(45,0,0,1)
m = m:scale(.5,1,2) -- and make sure it's not orthogonal any more
m = m:rotate(xa,m,m,m)
cube draws the unit cube and
xa is some parameter I can vary, then this rotates the cube about some seemingly random axis (not actually random, but may as well be).
However, if I put
m = m:rotate(xa,1,0,0) then the cube rotates about its natural x-axis.
Now where things get a bit confusing is if I do something like:
m = m:rotate(xa,1,0,0)
m = m:rotate(ya,0,1,0)
m = m:rotate(za,0,0,1)
because this rotates first about the z axis, then the y axis, and finally the x axis. If I vary xa, then ya, then za then it appears as though the cube rotates first about its natural x-axis, then about its natural y-axis, and finally about its natural z-axis. So this looks right. Until, that is, I go back and vary xa again. Now it looks as though it is no longer rotating about its natural x-axis but about the world’s x-axis. But that’s because I’m still applying the transformations in the same order. What I probably wanted to do was something more complicated which involves remembering the order in which I specify the operations at runtime not at start-up. So if I vary xa, then ya, and then vary xa again, I want it to apply the second variation of xa taking into account the ya rotation. The code snippet above can’t cope with that: it does all the xa’s then all the ya’s, and then all the za’s (actually, it does them in the opposite order … but that’s not the point).
So as you do, you need to keep track of the model matrix and accumulate the transformations. But each new transformation is with respect to the natural coordinate and so doesn’t need to know about the accumulated transformation up to that point.
Here’s some code that does that (there are better ways):
m = matrix()
ix,iy,iz = 0,0,0
if xa =~ ix then
m = m:rotate(xa - ix,1,0,0)
ix = xa
if ya =~ iy then
m = m:rotate(ya - iy,0,1,0)
iy = ya
if za =~ iz then
m = m:rotate(za - iz,0,0,1)
iz = za
This will always rotate the object about its natural axes.