Consider two bodies A and B. With respect to an inertial coordinate system with origin at point O, the coords of the particles in A are vectors $x_{a}\in V_{3}$ with $a=1,2,\ldots, N_{A}$ and similarly the coords of the particles of B are $x_{b}\in V_{3}$ with $b=1,2,\ldots,N_{B}$. The momenta wrt the inertial frame with origin at point O of the particles of A are $p_{a}\in V_{3}$ and the momenta of the particles of B are $p_{b}\in V_{3}$. The total angular momentum of the system wrt point O is,
$$
J=\sum_{a}x_{a}\times p_{a}+\sum_{b}x_{b}\times p_{b} \ .
$$
Let's introduce the centre of mass of body A as $X_{A}\in V_{3}$,
$$
X_{A}=\frac{\sum_{a}m_{a}x_{a}}{\sum_{a}m_{a}}=\frac{\sum_{a}m_{a}x_{a}}{M_{A}}
$$
and the centre of mass of body B as,
$$
X_{B}=\frac{\sum_{b}m_{b}x_{b}}{\sum_{b}m_{b}}=\frac{\sum_{b}m_{b}x_{b}}{M_{B}}
$$
Adding and subtracting the centre of mass coords,
$$
J=\sum_{a}(x_{a}-X_{A})\times p_{a}+\sum_{b}(x_{b}-X_{B})\times p_{b}+X_{A}\times \sum_{a}p_{a}+X_{B}\times \sum_{b}p_{b} \ .
$$
The first two terms on the RHS are the angular momenta of the bodies about their respective centres of mass $X_{A}$ and $X_{B}$. Let's write these contributions as $J_{A}\in V_{3}$ and $J_{B}\in V_{3}$. The total angular momentum is now,
$$
J=J_{A}+J_{B}+X_{A}\times \sum_{a}p_{a}+X_{B}\times \sum_{b}p_{b} \ .
$$
Let the linear momentum of the particles of body A be,
$$
P_{A}=\sum_{a}p_{a}
$$
with a similar formula for the sum of the momenta of the particles of body B. Put these formulae into the equation for the total angular momentum,
$$
J=J_{A}+J_{B}+X_{A}\times P_{A}+X_{B}\times P_{B} \ .
$$
This is precisely the splitting of the total angular momentum that Harold wrote in his question. $J_{A}$ and $J_{B}$ are the angular momenta of A and B about their own centres of mass. $X_{A}\times P_{A}$ is the orbital angular momentum of A about the origin O of the inertial coords and $X_{B}\times P_{B}$ is the orbital angular momentum of B about O. The point O could be taken as the centre of mass of the complete system, but it doesn't matter for the above result.
$\vec\omega = I^{-1} \vec L$, and $\vec L$ is constant in the absence of external forces. The bit that I think you're missing is that $I$ rotates with the rigid body, so it is not constant in general and neither is $\vec\omega$.
I played with your online example, and the angular velocity does seem to always remain constant when I'm not poking the block, which is consistent with an inertia tensor that's a scalar multiple of the identity. I've never used Unity, but it seems to support arbitrary tensors of inertia via Rigidbody.inertiaTensor and Rigidbody.inertiaTensorRotation, so maybe you just need to set those properly.
If you end up having to roll your own physics, here's some naive, untested, and potentially wrong pseudocode:
const double time_step = ...;
const Quaternion L = {0, Lx, Ly, Lz}; // angular momentum
const double K1 = 1/I1, K2 = 1/I2, K3 = 1/I3; // reciprocals of principal-axis moments of inertia; no idea if there's a standard letter for this
Quaternion orientation = {1, 0, 0, 0};
while (1) {
Quaternion transformed_L = conjugate(orientation) * L * orientation;
Quaternion transformed_omega = {0, K1 * transformed_L.i, K2 * transformed_L.j, K3 * transformed_L.k};
Quaternion omega = orientation * transformed_omega * conjugate(orientation);
orientation = quaternion_exp(time_step * omega) * orientation;
// ... or orientation = normalize((1 + time_step * omega) * orientation);
// ... or orientation = normalize((1 + time_step * omega + 0.5 * (time_step * omega)**2) * orientation);
output(orientation);
}
Best Answer
RKE is just a convenient way to describe KE in rotational systems. We have $1/2I\omega^2$ and $1/2mv^2$. These are equivalent, it just depends on what is convenient we have. Let's model the Moon orbiting the Earth as a point in orbit.
$I = mr^2$ for this orbit
We also note $v = \omega r$ therefore $\omega = \frac{v}{r}$
Taking $1/2I\omega^2$ We put in our expression for $I$
$1/2mr^2\omega^2$
and then put in our expression for omega
$1/2mv^2$
Which is the expression for KE. So it is apparent that KE and RKE are equivalent.
The same is true for the moon rotating about its own axis, however, it is much harder to show because the velocity is different for a piece of mass on the surface of the moon and a piece of mass closer to the center. You could integrate this and calculate the energy using $1/2mv^2$, however, it is much easier to look up $I$ and $\omega$. and calculate using $1/2I\omega^2$.
The moon therefore has two types of RKE or two types of KE or one KE and one RKE depending on how you want to look at it. All are equivalent, but some are easier to calculate.