Heavy symmetric top Hamiltonian (source code)

= Heavy symmetric top Hamiltonian

With transverse and axial moments $I,J$, the Hamiltonian is
$$
H=\frac{p_\theta^2}{2I}
+\frac{(p_\phi-p_\psi\cos\theta)^2}{2I\sin^2\theta}
+\frac{p_\psi^2}{2J}+Mga\cos\theta.
$$