Airplane velocity matrix I'm working on a programming project that emulates an airplane flight. The airplane has a velocity vector that can change when the airplane yaws or pitches. I have recreated this matrix that rotates a vector around the axis, but an airplane rotates differently - for instance when an airplane yaws it rotates around the y-axis of the airplane and not the axis of the system. When the airplanes rotate, its axis will rotate with it.
So my question is, mathematically, how fo I preform an airplane rotation on its own axis when the roll, pitch, and yaw of the airplane is known and then transform it to the axis of the system?
If an airplane is starting when the velocity is on the airplanes Z-axis then:
Here is a visual explanation of the question.
 A: The simple answer is, you don't; you use versors instead.
A versor is a four-component value $\mathbf{q} = (q_r, q_i, q_j, q_k) = (w, x, y, z)$, that describes an orientation.  It is also known as an unit quaternion, $w^2 + x^2 + y^2 + z^2 = 1$.  When doing calculations with versors, you can normalize the versor by dividing each component by $\sqrt{w^2 + x^2 + y^2 + z^2}$.
For calculating actual transforms, you derive the 3D rotation matrix from versor $\mathbf{q}$ using this matrix, where $s = 1/\sqrt{q_r^2 + q_i^2 + q_j^2 + q_k^2}$.
The default orientation, or "no change in orientation", is $(1, 0, 0, 0)$.  (Technically, versors of form $(0, x, y, z)$ where $x^2+y^2+z^2=1$ are also "no change in orientation", except that there there is an axis, just no rotation around that axis; the $(1, 0, 0, 0)$ is roughly equivalent to "no axis, and no rotation".)
One way to define a versor is by an axis vector $\hat{a} = (x_a, y_a, z_a)$, and a rotation around that vector by angle $\varphi$. Then, using
$$L = \lVert \hat{a} \rVert = \sqrt{x_a^2 + y_a^2 + z_a^2}$$
the corresponding versor is
$$\mathbf{q} = \left( \cos(\varphi), ~ \sin(\varphi)\frac{x_a}{L}, ~ \sin(\varphi)\frac{y_a}{L}, ~ \sin(\varphi)\frac{z_a}{L} \right)$$
The "trick" is that you can use quaternion multiplication, the Hamilton product, to combine any number of rotations.  If $\mathbf{q}_1$ is your current orientation, and you modify it by (change of) orientation $\mathbf{q}_2$, the final orientation is $\mathbf{q}_3 = \mathbf{q}_2 \mathbf{q}_1$.  Note that the order is important: initial orientation is on the right, and the last rotation leftmost.
The most useful thing is that you can scale the rotation using some value $f$, between $0$ and $1$:
$$\mathbf{q}_{temp} = \mathbf{q}_{old} + f \mathbf{q}_{change} \mathbf{q}_{old}$$
followed by normalizing $\mathbf{q}_{temp}$ as explained above.
Note that when $f$ is just a normal scalar (number), $f \mathbf{q} = (f q_r, f q_i , f q_j , f q_k)$, i.e. you just multiply each component with that value.  Addition and subtraction is done component-wise.  It is only the multiplication that is really "odd" with versors.
A practical example works better.
Let's assume $\mathbf{q}_0 = (w_0, x_0, y_0, z_0)$ describes the current orientation of your plane.
Let $\mathbf{q}_R$ describe the maximum amount of yaw the rudder can produce in one time unit, and $-1 \le R \le +1$ describe the current position of the rudder.
Let $\mathbf{q}_E$ describe the maximum amount of pitch the elevators can produce in one time unit, and $-1 \le E \le +1$ describe the current position of the elevators.
Let $\mathbf{q}_A$ describe the maximum amount of roll the ailerons can produce in one time unit, and $-1 \le A \le +1$ describe their current position.
Finally, let $T$ be the amount of time elapsed with the controls $R$, $E$, and $A$ in their current positions, and we want to find $\mathbf{q}$ describing the ensuing orientation. For simplicity, we'll calculate $\mathbf{q}_{temp}$, which after normalization will give $\mathbf{q}$:
$$\mathbf{q}_{temp} = \mathbf{q}_0 + \frac{1}{T} \left( R \mathbf{q}_R \mathbf{q}_0 + E \mathbf{q}_E \mathbf{q}_0 + A \mathbf{q}_A \mathbf{q}_0 \right)$$
Note that unlike with Euler or Tait-Bryan angles, the order of the additions does not matter; the answer is always the same.  There is also no gimbal lock.
In practice, you'll also want to modify the scalar multipliers ($R$, $E$, and $A$) to depends on the airspeed, as the efficacy of the control surfaces depends on the airspeed (and possibly direction); I don't know, I'm not writing the simulator.
(As used above, they are more like thrusters (or thruster pairs in opposite directions) that apply rotation to the object, without preserving angular momentum.)
The reason why you cannot do the same with rotation matrices is that well, mathematically you can, but in practice, rounding errors creep in, and after a couple of dozen operations, your rotation matrix is no longer orthonormal, and space becomes taffy.  Versor normalization has no bias in direction, so it just keeps space "rigid".  Unfortunately, all matrix normalization techniques suitable for a simulator tend to have such bias directions, so a simulator using that will have controls that behave slightly oddly in certain directions.
It is best to use versors instead to describe the orientation and changes in orientation, and just turn it into a rotation matrix when needed.
