Python code, done in collaboration with Brett Hajdaj:

 
    # ________________________________________________________________________________________
# Chaotic Tumbling of Hyperion
# 
# Authors: Anna Totilca, Brett Hajdaj
# PHYS 225
#
# The chaotic tumbling of Hyperion, a moon of Saturn, refers to its unpredictable and 
# non-repeating rotation. Unlike most moons, which spin consistently, 
# Hyperion's irregular shape and highly elliptical orbit cause it to experience varying 
# gravitational forces from Saturn. These forces make its spin axis change erratically,
# resulting in a chaotic motion where small changes in initial conditions lead to 
# drastically different rotational behaviors, making its spin unpredictable and unstable.
#
# Tumbling visualized: https://youtu.be/YC1X_zgv3Po?si=TYnXG00RWXqbBgeJ
# ________________________________________________________________________________________

import numpy as np
import matplotlib.pyplot as plt
import random

# Variable Setup
s = 1000 # Number of steps
dt = 0.0001 # yrs
GM_SC = 4.0 * np.pi ** 2.0 # AU^3 / yr^2

# Array Initialization
t = np.arange(0, s + dt, dt)
N = len(t)

r_c = np.zeros(N) # Radial Position 
x = np.zeros(N) # X Position
y = np.zeros(N) # Y Position
v_x = np.zeros(N) # X Velocity
v_y = np.zeros(N) # Y Velocity
theta = np.zeros(N) # Angle 
omega = np.zeros(N) # Angular Frequency

# Initial Conditions
x[0] = 1.0  # + 0.01*random.randint(-5,5) #(HU)
y[0] = 0.0 # Initial y (HU)
v_x[0] = 0.0 # (HU/yr)
v_y[0] = 5.0  # + 0.01*random.randint(-5,5) #(HU/yr)
theta[0] = 0.0 # (rad)
omega[0] = 0.0 # (rad/yr)

# Calculate Euler's Method
for i in range(1,N):
    # Update Radial Position
    r_c[i-1] = np.sqrt(x[i-1]**2 + y[i-1]**2) # HU
    
    # Update Velocities
    v_x[i] = v_x[i-1] - (4 * np.pi**2 * x[i-1] * dt) / (r_c[i-1]**3)
    v_y[i] = v_y[i-1] - (4 * np.pi**2 * y[i-1] * dt) / (r_c[i-1]**3)

    # Use Velocities to Update Position
    x[i] = x[i-1] + dt * v_x[i]
    y[i] = y[i-1] + dt * v_y[i]
    r_c[i] = np.hypot(x[i], y[i])

    # Calculate Angular Velocity & Theta
    Term1 = GM_SC / (r_c[i-1]**5)
    Term2 = x[i-1] * np.sin(theta[i-1]) - y[i-1] * np.cos(theta[i-1])
    Term3 = x[i-1] * np.cos(theta[i-1]) + y[i-1] * np.sin(theta[i-1])

    omega[i] = omega[i-1] - Term1 * Term2 * Term3 * dt
    theta[i] = theta[i-1] + omega[i-1] * dt

    # Safeguard θ
    tmp_theta = theta[i-1] + omega[i] * dt
    theta[i] = (tmp_theta + np.pi) % (2 * np.pi) - np.pi

# Store Final Radius
r_c[-1] = np.hypot(x[-1], y[-1])
    
# Plotting
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 6), sharex=True)

# θ vs Time Plot
ax1.plot(t, theta, color='purple')
ax1.set_xlim(0, 10)          
ax1.set_ylim(-np.pi, np.pi) # −π to π
ax1.set_ylabel(r'$\theta\;\mathrm{(rad)}$')
ax1.set_title(r'$\theta$ versus time for Hyperion')

# Omega vs Time Plot
ax2.plot(t, omega, color='orange')
ax2.set_xlim(0, 10) 
ax2.set_ylim(-20, 60)
ax2.set_xlabel('time (yr)')
ax2.set_ylabel(r'$\omega\;\mathrm{(rad/yr)}$')
ax2.set_title(r'$\omega$ versus time for Hyperion')

fig.tight_layout()
plt.show()

# Conversions to Polar:
theta_orbit = np.arctan2(y,x) # radians; angle of orbit, found from the inverse tangent of position components
r_orbit = np.hypot(x,y) # radius of hyperion's orbit

# Polar Plot of Hyperion's orbit
fig_polar, ax_polar = plt.subplots(subplot_kw = {'projection': 'polar'}, figsize=(6,6))

ax_polar.plot(theta_orbit, r_orbit, "b-")
ax_polar.set_rticks([0.5, 1, 1.5, 2])  
ax_polar.grid(True)
ax_polar.set_title("Hyperion's 2D Orbit")
ax_polar.plot(0, 0, 'yo', markersize=10)  # Saturn for context

print("Mean radius of the orbit", np.mean(r_orbit), "HU")
print(r_orbit)
plt.show()

# Further work can focus on the 3-body problem, exploring the ways which Titan's gravitational impact on Hyperion,
# as they have 4:3 orbital resonance (for every 4 orbits Hyperion completes, Titan completes 4; see this video:
# https://youtube.com/shorts/l_1pWzOsVo8?si=Dd2ydIfUPXl7FyIM) 
# and Titan is almost 250,000x more massive than Hyperion. This will have an effect on Hyperion's tumbling rotation