We use an angular momentum approach to study the Cassini states (CS) of large natural satellites such as the Galilean satellites and Titan. Unlike classical approaches where obliquity is the solution of a trigonometric equation, our approach allows us to identify not only the mean obliquity of satellites, but also their nutation in space as well as their polar motion (PM) with respect to the solid surface. Triaxiality of the satellite has a significant effect on the mean obliquities of CSI, CSII and CSIV. We assess the stability of the Cassini states over a wide range of free and forced precession frequency ratios and find that CSI and CSIII are always stable. Even if the different Galilean Moons are thought to occupy CSI, we therefore also analytically study CSIII. We here solve the dynamic equations governing CSI and CSIII up to order two in small quantities and without averaging the external torque over the mean anomaly to obtain the time-variable obliquity and polar motion (at long-period and short-period). By extending the system of equations governing CSI, including gravitational and pressure couplings between misaligned layers, we predict the orientation of the spin axes of the outer shell, internal ocean and solid interior for an ocean-bearing body.