Learning Scientific Programming with Python (2nd edition)

P8.2.7: The tumbling of Hyperion

Question P8.2.7

Hyperion is an irregularly-shaped moon of Saturn notable for its chaotic rotation. Its motion may be modelled as follows.

The orbit of Hyperion (H) about Saturn (S) is an ellipse with semi-major axis, $a$, and eccentricity, $e$. Let its point of closest approach (periapsis) be P. Its distance from the planet, SH, as a function of its true anomaly (orbital angle, $\phi$, measured from the line SP) is therefore $$ r = \frac{a(1-e^2)}{1+e\cos\phi}. $$

Define the angle $\theta$ to be that between the axis of the smallest principal moment of inertia (loosely, the longest axis of the moon) and SP, and the quantity $\Omega$ to be a scaled rate of change of $\theta$ with $\phi$ (i.e. the rate at which Hyperion spins as it orbits Saturn) as follows: $$ \Omega = \frac{a^2}{r^2}\frac{\mathrm{d}\theta}{\mathrm{d}\phi}. $$

The geometry of the orbit of Hyperion

The geometry of the orbit of Hyperion.

Now, it can be shown that $$ \frac{\mathrm{d}\Omega}{\mathrm{d}\phi} = - \frac{B-A}{C}\frac{3}{2(1-e^2)}\frac{a}{r}\sin [2(\theta-\phi)], $$ where $A$, $B$ and $C$ are the principal moments of inertia.

Use scipy.integrate.odeint to find and plot the spin rate, $\Omega$, as a function of $\phi$ for the initial conditions (a) $\theta = \Omega = 0$ at $\phi = 0$, and (b) $\theta = 0$, $\Omega = 2$ at $\phi=0$. Take $e=0.1$ and $(B-A)/C = 0.265$.