Solution (source code)

= Solution

Differentiation gives
$$
\boxed{p_\theta=I_1\dot\theta,\quad
p_\psi=I_3(\dot\psi+\dot\phi\cos\theta),\quad
p_\phi=I_1\dot\phi\sin^2\theta+p_\psi\cos\theta.}
$$
The angles $\phi$ and $\psi$ are cyclic, so their canonical momenta are conserved. The <Lagrangian> is autonomous, so the total <energy> is also conserved. The <momentum> $p_\theta$ is not generally constant because <gravity> and the angular kinetic terms depend on $\theta$.