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
- 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()
```