English: Given fixed angular momentum, for each value of kinetic energy, we have a family of curves for the possible value of angular velocity. The curves are plotted in black, on the orange ellipsoid of fixed angular momentum.

```python import numpy as np import matplotlib.pyplot as plt import scipy

1. Adjustable parameters

L = 1.0 L2 = L**2 I1, I2, I3 = 1.0, 2.0, 3.0 xmin, xmax, ymin, ymax, zmin, zmax = -L/I1, L/I1, -L/I2, L/I2, -L/I3, L/I3 E1, E2, E3 = 0.5*L2/I1, 0.5*L2/I2, 0.5*L2/I3

def parametric_plot(E, orbit_res=1000):

   thetas = np.linspace(0, 2*np.pi, orbit_res)
   if E3 < E and E < E2:
       rs = np.zeros(orbit_res)
       zs = np.zeros(orbit_res)
       for i, theta in enumerate(thetas):
           invmatrix = scipy.linalg.inv(np.array([[I1 * np.cos(theta)**2 + I2 * np.sin(theta)**2    , I3],
                                                 [(I1 * np.cos(theta))**2 + (I2 * np.sin(theta))**2, I3**2]]))
           r2z2 = invmatrix @ np.array([[2*E], [L2]])
           rs[i] = np.sqrt(r2z2[0,0])
           zs[i] = np.sqrt(r2z2[1,0])
       return rs * np.cos(thetas), rs * np.sin(thetas), zs
   if E2 < E and E < E1:
       xs = np.zeros(orbit_res)
       rs = np.zeros(orbit_res)
       for i, theta in enumerate(thetas):
           invmatrix = scipy.linalg.inv(np.array([[I1, I2 * np.cos(theta)**2 + I3 * np.sin(theta)**2],
                                                 [I1, (I2 * np.cos(theta))**2 + (I3 * np.sin(theta))**2]]))
           x2r2 = invmatrix @ np.array([[2*E], [L2]])
           xs[i] = np.sqrt(x2r2[0,0])
           rs[i] = np.sqrt(x2r2[1,0])
       return xs, rs * np.cos(thetas), rs * np.sin(thetas)

fig = plt.figure(figsize = (16,16)) ax = plt.axes(projection='3d')

u = np.linspace(0, 2 * np.pi, 100) v = np.linspace(0, np.pi, 100) x = np.outer(np.cos(u), np.sin(v)) y = np.outer(np.sin(u), np.sin(v)) z = np.outer(np.ones(np.size(u)), np.cos(v)) x *= L/I1 y *= L/I2 z *= L/I3

ax.plot_surface(x, y, z, color='orange', alpha=0.3)

epsilon=1e-3 for E in np.linspace(E3+epsilon, E1-epsilon, 20):

   xs, ys, zs = parametric_plot(E)
   xs *= 1+epsilon
   ys *= 1+epsilon
   zs *= 1+epsilon
   ax.plot3D(xs, ys, zs, linewidth=1, color='k')
   ax.plot3D(-xs, -ys, -zs, linewidth=1, color='k')

ax.axes.set_xlim3d(xmin, xmax) ax.axes.set_ylim3d(ymin, ymax) ax.axes.set_zlim3d(zmin, zmax)

ax.set_aspect('equal') ax.set_axis_off()

plt.show()

```