Orbits in an Axisymmetric Potential

A flattened galaxy is well described by the Miyamoto-Nagai potential $\Phi(R,z) = -GM / \sqrt{R^2 + (a + \sqrt{z^2 + b^2})^2}$, with cylindrical radius $R$, height $z$, disk scale length $a$ and thickness $b$. Because this is not an inverse-square point-mass field, the radial and azimuthal oscillation frequencies are incommensurate, so a generic bound orbit does not close: it precesses, filling a rosette bounded by an inner and outer radius. Energy $E$ and the angular momentum component $L_z$ about the symmetry axis are conserved and pin the rosette, while a third (non-classical) integral confines the motion in the meridional plane. The playground integrates an orbit with a symplectic step and shows the rosette plus the conserved quantities holding flat.