How can I simulate a pot lid rotating around an axis that is quickly rotating?

Question

How can I simulate the rotation of a pot lid about an axis that is quickly rotating as in the following image? I'm particular interested in the motion of the red dot.

Here is my attempt:

φ = π/6; r = 1;
line = Rotate[{Line[
     Table[{Cos[θ], Sin[θ], 0}, {θ, 0, 
       2 π, π/20}]], {Red, Pink, PointSize[Large], 
     Point[{1, 0, 0}]}}, φ, {0, 1, 0}];
Manipulate[
 Graphics3D[{Rotate[
    Rotate[line, θ, {0, 0, 
      1}], θ Cos[φ], {Sin[φ] Cos[θ], 
     Sin[φ] Sin[θ], Cos[φ]}, {0, 0, 0}], {Line[
     Table[{r Cos[φ] Cos[θ1], 
       r Cos[φ] Sin[θ1], -r Sin[φ]}, {θ1, 
       0, θ, π/20}]]}, {Blue, PointSize[Large], 
    Point[{0, 0, 0}]}}, Axes -> True, PlotRange -> 2, 
  AxesLabel -> {"x", "y", "z"}], {θ, 0, 2 π}]

Answer 1

In my opinion, this problem is the same as a rolling ellipsoid on a plane without slipping, which can solved by Euler equations for rigid body.

The Euler equations in body axes are

$$ I_1\dot\omega_1-\omega_2\omega_3(I_2-I_3)=0\\ I_2\dot\omega_2-\omega_3\omega_1(I_3-I_1)=0\\ I_3\dot\omega_3-\omega_1\omega_2(I_1-I_2)=0 $$

with initial conditions $$ \omega_1(0)=0,\omega_2(0)=0,\omega_3(0)=\omega_{30} $$

where $\omega_i$ and $I_i$ are the angular frequency and momentum of inertia in the three principle axes.

Then the three Euler angles $\theta, \psi, \phi$ are solved using the Euler angle equations:

\begin{eqnarray} \dot\phi\sin\theta\sin\psi+\dot\theta\cos\psi&=&\omega_1\\ \dot\phi\sin\theta\cos\psi-\dot\theta\sin\psi&=&\omega_2\\ \dot\phi\cos\theta+\dot\psi&=&\omega_3 \end{eqnarray}

After get the Euler angles function of time, we can easily make a movie. The line on the ellipsoid is called polhode (black curve) and the line on the plane is called herpolhode (red curve).

Clear["`*"]
a = 3; b = 4; c = 1;(*The Length of the cubic*)
I1 = M/12 (b^2 + c^2); I2 = M/12 (a^2 + c^2); I3 = M/12 (a^2 + b^2); M = 12;(*Three momentum of inertia in principle axes*)
w10= 1; w20 = 2; w30 = 1;(*The initial conditions*)
L = Sqrt[w10^2 I1^2 + w20^2 I2^2 + w30^2 I3^2];(*The magnitude of the angular momentum*)  
T = 1/2 w10^2 I1 + 1/2 w20^2 I2 + 1/2 w30^2 I3 ;(*The kenetic engergy*)
d = Sqrt[2 T]/L;(*The distance of the origin to the plane*)

solw = NDSolve[{I1 w1'[t] == (I2 - I3) w2[t] w3[t], 
     I2 w2'[t] == (I3 - I1) w1[t] w3[t], 
     I3 w3'[t] == (I1 - I2) w2[t] w1[t], w1[0] == w10, w2[0] == w20, 
     w3[0] == w30}, {w1, w2, w3}, {t, 0, 20}] // 
   Flatten;(*Solve Euler equations, in the body axes*)

solang = NDSolve[{\[Phi]'[
          t] Sin[\[Theta][t]] Sin[\[Psi][t]] + \[Theta]'[
          t] Cos[\[Psi][t]] == w1[t],
      \[Phi]'[t] Sin[\[Theta][t]] Cos[\[Psi][t]] - \[Theta]'[
          t] Sin[\[Psi][t]] == w2[t],
      \[Phi]'[t] Cos[\[Theta][t]] + \[Psi]'[t] == w3[t], \[Phi][0] == 
       0.01, \[Theta][0] == 0.01, \[Psi][0] == 0.01} /. 
     solw, {\[Phi], \[Theta], \[Psi]}, {t, 0, 20}] // 
   Flatten;(*Solve the equations of Euler angles*)

w1[t_] = w1[t] /. solw; w2[t_] = w2[t] /. solw; w3[t_] = w3[t] /. solw;
\[Phi][t_] = \[Phi][t] /. solang; \[Theta][t_] = \[Theta][t] /. 
  solang; \[Psi][t_] = \[Psi][t] /. solang;
A[t_] = {{Cos[\[Psi][t]], Sin[\[Psi][t]], 0}, {-Sin[\[Psi][t]], 
     Cos[\[Psi][t]], 0}, {0, 0, 1}}.{{1, 0, 0}, {0, Cos[\[Theta][t]], 
     Sin[\[Theta][t]]}, {0, -Sin[\[Theta][t]], 
     Cos[\[Theta][t]]}}.{{Cos[\[Phi][t]], Sin[\[Phi][t]], 
     0}, {-Sin[\[Phi][t]], Cos[\[Phi][t]], 0}, {0, 0, 
     1}};(*Transformation matrix as a function of time*)

cpilherpolhode = 
  Compile[{{t, _Real}}, 
   Transpose[{{Cos[\[Psi][t]], Sin[\[Psi][t]], 0}, {-Sin[\[Psi][t]], 
        Cos[\[Psi][t]], 0}, {0, 0, 1}}.{{1, 0, 0}, {0, 
        Cos[\[Theta][t]], Sin[\[Theta][t]]}, {0, -Sin[\[Theta][t]], 
        Cos[\[Theta][t]]}}.{{Cos[\[Phi][t]], Sin[\[Phi][t]], 
        0}, {-Sin[\[Phi][t]], Cos[\[Phi][t]], 0}, {0, 0, 1}}].(1/Sqrt[
      2 T] {w1[t], w2[t], w3[t]})];(*Compile the herpolhode equations*)

polhode = 
  Table[Point[(1/Sqrt[2 T] {w1[t], w2[t], w3[t]})], {t, 0, 10, 
     0.1}] /. solw;(*polhode*) 
elps = Scale[{Sphere[], 
    Table[Rotate[
      Line@Table[{1.01 Cos[u], 1.01 Sin[u], 0}, {u, 0., 2 Pi, 
         2 Pi/30}], i, {0, 1, 0}], {i, 0, \[Pi], \[Pi]/3}], 
    Rotate[Line@
      Table[{1.01 Cos[u], 1.01 Sin[u], 0}, {u, 0., 2 Pi, 
        2 Pi/30}], \[Pi]/2, {1, 0, 0}]}, {1/Sqrt[I1], 1/Sqrt[I2] , 1/
    Sqrt[I3] }, {0, 0, 0}];(*ellipsoid of inertia with lines on it*)

elpspolhode = {elps, polhode};(*combine the ellipsoid and the polhode*)

plane = Plot3D[
   L/(w30 I3) (d - (w10 I1)/L x - (w20 I2)/L y), {x, -0.4, 
    0.4}, {y, -0.5, 0.5}, Mesh -> None, 
   PlotStyle -> Opacity[0.5]];(*plane*)

Manipulate[
 Show[{Graphics3D[
    Rotate[Rotate[
      Rotate[elpspolhode, \[Phi][t], {0, 0, 1}], \[Theta][
       t], {Cos[\[Phi][t]], Sin[\[Phi][t]], 0}], \[Psi][
      t], {Sin[\[Theta][t]] Sin[\[Phi][
         t]], -Sin[\[Theta][t]] Cos[\[Phi][t]], Cos[\[Theta][t]]}]], 
   plane}, ParametricPlot3D[cpilherpolhode[tt], {tt, 0, t}, 
   PlotStyle -> Red], AspectRatio -> 1, PlotRange -> 0.5, 
  Boxed -> False, ViewAngle -> 0.24252597899424388`, 
  ViewPoint -> {0.09410273512100081`, -1.389522619330309`, 
    3.083885813479511`}, 
  ViewVertical -> {0.4688415188450498`, 0.3712202655129461`, 
    0.8014880814332579`}], {t, 0.001, 20, 0.1}]