Quit[];

El lagrangiano, los impulsos generalizados y las velocidades escritas en términos de los impulsos:

L[θ_, θp_, ψp_, ϕp_]:= (I3/2)(((ψp + ϕp Cos[θ]))^2)+ (I1/2)(((θp)^2)+ ((ϕp)^2)((Sin[θ])^2)) - mgl Cos[θ];
sol = Collect[
Solve[{
pθ == D[L[θ, θp, ψp, ϕp], θp], 
pψ == D[L[θ, θp, ψp, ϕp], ψp],
pϕ == D[L[θ, θp, ψp, ϕp], ϕp]
},
{θp, ψp, ϕp}]
, {I1, I3}, Simplify]

Uso θp, ϕp y ψp para denotar  (OverDot[θ]), (OverDot[ϕ]), (OverDot[ψ]). Los impulsos generalizados con pθ, etc.

El hamiltoniano, construido como Σ (Subscript[(OverDot[q]), i])(Subscript[p, i])-L

Collect[θp pθ + ψp pψ + ϕp pϕ - L[θ, θp, ψp, ϕp] /. sol[[1]],{pθ,I1, I3}, TrigFactor]
Collect[%, {pθ,I1, I3}, Simplify]

H[θ_, pθ_, pψ_, pϕ_]:= (((pθ)^2)/(2 I1))+(((pψ)^2)/(2 I3))+(((((pϕ-pψ Cos[θ]))^2) ((Csc[θ])^2))/(2 I1))+mgl Cos[θ];

Podemos multiplicar el hamiltoniano por una constante sin afectar las ecuaciones de movimiento. Lo vamos a multiplicar por I1.

A la constante mgl I1 la llamamos α, y a la relación (Subscript[I, 1])/(Subscript[I, 3]) la llamamos β.

Hbis[θ_, pθ_, pψ_, pϕ_, α_, β_]:= (((pθ)^2)/2)+β(((pψ)^2)/2)+α Cos[θ]+(((((pϕ-pψ Cos[θ]))^2) ((Csc[θ])^2))/2);

Las coordenadas ψ y ϕ con cíclicas. Ignorando la constante aditiva ((β ((Subscript[p, ψ])^2))/2), podemos entonces definir un hamiltoniano efectivo para θ y (Subscript[p, θ]), que depende de un único parámetro, α, y de los valores de las constantes de movimiento (Subscript[p, ψ]) y (Subscript[p, ϕ]).

Hθ[θ_, pθ_, pψ_, pϕ_, α_]:=  (((pθ)^2)/2)+(((((pϕ-pψ Cos[θ]))^2) ((Csc[θ])^2))/2)+α Cos[θ];

Veamos el caso particular en que (Subscript[p, ψ])=(Subscript[p, ϕ])=p.

Collect[Hθ[θ, pθ, p, p, α], p, Simplify]

TrigReduce[((Tan[(θ/2)])^2)]

Hθ1[θ_, pθ_, p_, α_]:= (((pθ)^2)/2)+α Cos[θ]+(1/2) (p^2) ((1-Cos[θ])/(1+Cos[θ]))

Es un sistema conservativo, de modo que el punto en el espacio de fases θ-(Subscript[p, θ]) se mueve a lo largo de curvas de Hθ1 constante.
Dividamos por la constante p y redefinamos las coordenadas en el espacio de fases como θ y P = (Subscript[p, θ])/p. Hay que graficar las curvas de nivel de la función

h[θ_,P_, A_]:= ((P^2)/2)+A Cos[θ]+(1/2)  ((1-Cos[θ])/(1+Cos[θ]));

donde la constante A=α/p=mgl I1/(p^2) mide la importancia relativa de la gravedad respecto del efecto estabilizador del espín.

Será útil definir el potencial efectivo

Vef[θ_, A_]:= A Cos[θ]+(1/2)  ((1-Cos[θ])/(1+Cos[θ]));

Analicemos las curvas de nivel para algunos valores de (A^□):

ContourPlot[h[θ, P, #], {θ, -π, π}, {P, -2, 2}, 
Contours->Table[#× (1.5^n), {n,-1,10}], 
PlotLabel->StringJoin["A=", ToString[#]], 
ImageSize->150] & /@ {0.05, 0.1, 0.25, 0.5, 0.75, 1}

Las curvas de nivel son en esencia el gráfico del potencial efectivo restringido a los tramos del eje θ en donde es mayor que cero.

En efecto: la curva de nivel correspondiente al valor h = ϵ está dada por


	P = ±2(Sqrt[ϵ-(A Cos[θ]+(1/2)  ((1-Cos[θ])/(1+Cos[θ])))]).
	
Los cambios cualitativos al aumentar A están ligados a un cambio cualitativo en la forma del potencial efectivo. Veámoslo graficando 

	ϵ-Vef[θ, A]

 para los mismos valores de A que usamos antes, tomando algún valor particular de ϵ, lo que no afecta a la forma de la curva.

Plot[2-Vef[θ, #], {θ, -π, π},
PlotRange-> {-0.5,2},  
 PlotLabel->StringJoin["A=", ToString[#]],
ImageSize->150, Frame-> True] & /@ {0.05, 0.1, 0.25, 0.5, 0.75, 1}

Las curvas de nivel en el espacio de fases no son otra cosa que +± la raíz cuadrada de estas funciones.

El punto crítico ocurre cuando  θ = 0 deja de ser un mínimo de (Subscript[V, ef]) para transformarse de un máximo. Para A mayores que este valor crítico, (Subscript[V, ef]) tiene un máximo en θ = 0 y dos mínimos, simétricos respecto de θ = 0. Los mínimos aparecen en θ = +0 y se mueven hacia ±π a medida que aumenta A, lo que tiene sentido, porque aumentar A equivale a aumentar g.

Veámoslo. Para calcular la derivada del Vef respecto de θ, definimos x = Cos[θ], derivamos respecto de x y, siguiendo la regla de la cadena,  luego multiplicamos por -(Sqrt[1-(x^2)]).

Simplify[Solve[D[A x+(1/2)  ((1-x)/(1+x)), x] × -(Sqrt[1-(x^2)])== 0, x]]

Aquí x es el Cos[θ], de modo que la solución Cos[θ] = -1-1/(Sqrt[A]) nunca corresponde a un ángulo real, puesto que es menor que -1. La solución Cos[θ] = -1+1/(Sqrt[A]) es realizable si 

	A> 1/4.

En tal caso, para Cos[θ] = -1+1/(Sqrt[A])existen dos ángulos posibles, simétricos respecto del cero, como es fácil convencerse a través de la siguiente figura:

Apr = 0.66;
Plot[{Cos[x], -1+(1/(Sqrt[Apr]))}, {x, -π, π}]

Resumiendo: (Subscript[V, ef]) tiene un sólo extremo para θ real si A < 1/2.  Es decir, para valores suficientemente pequeños de g el único punto de equilibrio es θ = 0.  Vemos por la forma del potencial que en ese caso todas las órbitas en θ son acotadas y giran alrededor del origen.

Si A> 1/4, entonces aparecen los otros dos extremos. La forma del potencial indica que hay dos clases de órbitas. Para energías bajas, θ se mueve alrededor de uno de los mínimos, mientras que para energías mayores que Vef[0, A], el ángulo θ se mueve en una órbita que rodea al mismo tiempo al origen y a los dos mínimos del potencial. 

Eso se ve claramente en el retrato de fases y en el gráfico del potencial efectivo:

Apr = .88;
{
ContourPlot[h[θ, P, Apr], {θ, -π, π}, {P, -2, 2},
 Contours->Table[Apr× (1.5^n), {n,-1,10}], 
PlotLabel->StringJoin["A=", ToString[Apr]], 
ImageSize->250],

Plot[Vef[θ, Apr], {θ, -π, π}, 
PlotRange-> {0, 2},
ImageSize->250,
 PlotLabel->  "(Subscript[V, ef])"]
}

La curva separatriz entre los dos tipos de órbitas es la que pasa por θ = 0, P = 0, y por lo tanto  corresponde a ϵ=A  . Su ecuación es

	P = ±2(Sqrt[(1- Cos[θ])(A-((1/2)/(1+Cos[θ]))  )]).
	
La intersección con el eje P= 0 se produce cuando

	Cos[θ] = (1/(2A))-1-> (Subscript[θ, +±])= ±ArcCos[(1/(2A))-1].
	
Esto será relevante cuando demos condiciones iniciales. Si inicialmente fijamos P= 0, entonces para θ comprendido entre los dos valores (Subscript[θ, +±]) tendremos órbitas alrededor de uno de los puntos de equilibrio estables. En cambio, si |θ| es mayor que (SubPlus[θ]) tendremos órbitas que encerrarán los tres puntos de equilibrio. (Ver la figura anterior.)

Apr = 1.5;
Plot[{2(Sqrt[(1- Cos[θ])(Apr-((1/2)/(1+Cos[θ]))  )]), -2(Sqrt[(1- Cos[θ])(Apr-((1/2)/(1+Cos[θ]))  )])}, {θ, -π, π}, 
PlotStyle-> Blue, 
Background->LightBlue, 
PlotLabel-> "Curva separatriz",ImageSize->300]

El movimiento de la peonza en el espacio: integración numérica de las ecuaciones

Partimos del hamiltoniano original

H[θ_, pθ_, pψ_, pϕ_]:= (((pθ)^2)/(2 I1))+(((pψ)^2)/(2 I3))+(((((pϕ-pψ Cos[θ]))^2) )/(2 I1 ((Sin[θ])^2)))+mgl Cos[θ];

Ya sabemos que (Subscript[p, ψ]) y (Subscript[p, ϕ]) son constantes. Las ecuaciones de movimiento relevantes son las que dan las velocidades en términos de los impulsos

	{θp->((pθ)/I1),ψp->((pψ)/I3)-(((pϕ-pψ Cos[θ]) Cos[θ])/(I1 ((Sin[θ])^2))),ϕp->(((pϕ-pψ Cos[θ]))/(I1 ((Sin[θ])^2)))}

y la que da (Subscript[(OverDot[p]), θ]), que se obtiene derivando H respecto de θ y multiplicando × -1.

Simplify[-D[H[θ, pθ, pψ, pϕ], θ]]

Es decir

	 {pθp->-((pψ (pϕ-pψ Cos[θ]) )/(I1 Sin[θ]))+(((((pϕ-pψ Cos[θ]))^2) Cos[θ])/(I1 ((Sin[θ])^2)))+mgl Sin[θ]}

Caso (Subscript[p, ϕ]) = (Subscript[p, ψ])= p.

Para fijar el problema usaremos (Subscript[I, 3]) = 2 (Subscript[I, 1]), como si la peonza fuera un disco atravesado por un eje sin masa. Dado p, podemos definir una velocidad angular característica del problema, 

	ω = p/(Subscript[I, 3]). 

Esto nos permite introducir una variable temporal adimensional 

	τ = ω t. 
	
Una variación de τ en una unidad, corresponde, a grandes rasgos, a una vuelta de la peonza sobre su eje.

Volviendo a las ecs. de movimiento, con (Subscript[p, ψ])= (Subscript[p, ϕ])= p,  (Subscript[I, 3]) = 2 (Subscript[I, 1]), y cambiando a la variable τ  (de modo que ahora θp es dθ/dτ, etc.), queda:

TrigReduce[Simplify[{θp->(I3/p)(((pθ)/I1)),ψp->((I3/p))(((pψ)/I3)-(((pϕ-pψ Cos[θ]) Cos[θ])/(I1 ((Sin[θ])^2)))),ϕp->((I3/p))((((pϕ-pψ Cos[θ]))/(I1 ((Sin[θ])^2))))} /. pϕ-> pψ /. I1-> (I3/2)/. pψ-> p  ]]

Collect[((I3/p))(((pψ (-pϕ+pψ Cos[θ]) Csc[θ])/I1)+(((((pϕ-pψ Cos[θ]))^2) Cot[θ] ((Csc[θ])^2))/I1)+mgl Sin[θ]) /. pϕ-> pψ /. I1-> (I3/2)/. pψ-> p  , mgl, FullSimplify]

Es decir

	{θp->2((pθ)/p),ψp->((1- Cos[θ])/(1+Cos[θ])),ϕp->(2/(1+Cos[θ])), pθp -> -8 p ((Csc[θ])^3) ((Sin[(θ/2)])^4)+((I3 mgl Sin[θ])/p)}

Ya casi estamos. Así como adimensionalizamos el tiempo, podemos adimensionalizar los impulsos, dividiéndolos por el impulso característico p.

Entonces, ahora escribiremos ((Subscript[p, θ])/p)-> P y ((pθp)/p)-> Pp. Al dividir  por p en el término asociado al potencial gravitatorio aparece la constante

	(Subscript[I, 3])mgl / (p^2)
	
que, según nuestras primeras definiciones, es igual a 2A, el parámetro que mide la importancia relativa de la gravedad respecto del espín.

Finalmente, las ecuaciones de movimiento adimensionalizadas son:

	{θp->2P,ψp->((1- Cos[θ])/(1+Cos[θ])),ϕp->(2/(1+Cos[θ])), Pp -> -8  ((Csc[θ])^3) ((Sin[(θ/2)])^4)+2A Sin[θ]}
	
Introduciendo algunas simplificaciones para las funciones trigonométricas resultan:

FullSimplify[{θp->2P,ψp->((1- Cos[θ])/(1+Cos[θ])),ϕp->(2/(1+Cos[θ])), Pp -> -8  ((Csc[θ])^3) ((Sin[(θ/2)])^4)+2A Sin[θ]}]

El verdedero problema aquí se reduce a integrar la ecuación para Pθ y θ. Vamos a hacerlo numéricamente.

Ya tenemos las ecs. diferenciales. Las condiciones iniciales para (Subscript[p, ϕ])= (Subscript[p, ψ])= p están contenidas en el parámetro A. 

Como ϕ y ψ son cíclicas, podemos fijar sus valores iniciales arbitrariamente. Uno estaría inclinado a usar ϕ = ψ = 0, pero ese punto corresponde a un punto singular de los ángulos de Euler, de modo que conviene elegir ϕ igual a cualquier otro valor.

Las condiciones iniciales para θ y pθ deben considerarse con mayor cuidado. Al analizar los retratos de fase vimos que θ puede hacer dos cosas bien diferentes: rodear el origen o rodear (si los hay) alguno de los puntos de equilibrio que aparecen cuando A es mayor que 1/4. En cualquier caso, algo es seguro, para todas las trayectorias hay puntos de retorno, en los cuales (Subscript[p, θ]) = 0. De manera que podemos tomar ésa como condición inicial para (Subscript[p, θ]).

Si estamos en el caso A< 1/4, para distintos valores iniciales de θ no obtendremos resultados cualitativamente diferentes. Siempre son órbitas alrededor del origen. 

Si estamos en el caso A> 1/4,  la forma de las órbitas dependerá de si el valor inicial de θ está dentro o fuera de la curva separatriz. 

	|θ|< (SubPlus[θ])= ArcCos[(1/(2A))-1]-> dentro de la separatriz

Probamos primero como marcha el asunto integrando un caso particular:

(Ustedes cambien el valor de Apr y el de θ0. Prueben con Apr mayor o menor que 1/4, y -en el primer caso- con θ0 mayor o menor que θmas)

Apr = Rationalize[0.33]; (*uso Rationalize para generar un número 'exacto'*)

(* Quiero tener a la vista el ángulo que separa las condiciones iniciales que para A > 1/4 caen de un lado o del otro de la separatriz*)
θmas = If[Apr > 1/4, ArcCos[(1/(2.Apr))-1], π] ;
Print[Style[StringJoin["(SubPlus[θ])= ", ToString[θmas]],Red, 16]]

ϕ0 =1;
ψ0 = 0;
θ0 =Rationalize[1];

Dsol = NDSolve[
{
θ'[t] == 2 P[t],
ψ'[t] == ((Tan[((θ[t])/2)])^2),
ϕ'[t] == ((Sec[((θ[t])/2)])^2),
 P'[t] == -8  ((Csc[θ[t]])^3) ((Sin[((θ[t])/2)])^4)+2Apr Sin[θ[t]],
P[0] == 0,
ϕ[0] == ϕ0,
ψ[0] == ψ0,
θ[0] == θ0},
{θ,P, ψ, ϕ}, {t, 0, Tmax = 50}]

theta  = Dsol[[1]][[1]][[2]];
ptheta  = Dsol[[1]][[2]][[2]];
psi= Dsol[[1]][[3]][[2]];
phi = Dsol[[1]][[4]][[2]];

ParametricPlot[{theta[t], ptheta[t]}, {t, 0, Tmax}, PlotRange-> All, AxesOrigin->{0,0}]

Print["θ+= ", θmas]

Para un análisis más sistemático, definimos una función que integra las ecuaciones de movimiento dados los valores de A, (Subscript[θ, 0]) y el tiempo máximo.

El intervalo de la solución va entre t = 0 y t = tmax, a menos que el Mathematica diga lo contrario.

solucion[A_, θ0_, tmax_]:= 
Module[{θ, ϕ, ψ, P,  ϕ0=1, ψ0=0, t, Dsol},
Dsol = NDSolve[
{
θ'[t] == 2 P[t],
ψ'[t] == ((Tan[((θ[t])/2)])^2),
ϕ'[t] == ((Sec[((θ[t])/2)])^2),
 P'[t] == -8  ((Csc[θ[t]])^3) ((Sin[((θ[t])/2)])^4)+2A Sin[θ[t]],
P[0] == 0,
ϕ[0] == ϕ0,
ψ[0] == ψ0,
θ[0] == θ0},
{θ,P, ψ, ϕ}, {t, 0, tmax}];
Dsol[[1]][[All, 2]]
];

Probemos nuestra función con un valor de A menor que 1/4 y valores iniciales de θ entre casi -π y casi π, esquivando el cero.

Apr = 0.1;
fraccion = Rationalize[0.85];
Monitor[tabla1 = Table[solucion[Rationalize[Apr], θ0, 20], {θ0, -fraccion π, fraccion π,  (π/10)}];, θ0]

g1 = ParametricPlot[Table[{tabla1[[i]][[1]][t], tabla1[[i]][[2]][t]}, {i, 1, Length[tabla1]}], {t, 0, 10}, PlotRange-> All, PlotStyle-> Red]

Hacemos la comparación con los retratos de fase hechos más arriba:

g1contour =ContourPlot[h[θ, P, Apr], {θ, -π, π}, {P, -2, 2},
 Contours->Table[h[θ0, 0, Apr], {θ0, -fraccion π, fraccion  π,  (π/20)}], 
PlotLabel->StringJoin["A=", ToString[Apr]], 
ImageSize->200]

Show[g1contour, g1, ImageSize-> 400]

Probemos ahora con un valor de A mayor que 1/4 y valores iniciales de θ entre casi -π y casi π, esquivando el cero.

Apr = 0.88;
fraccion = Rationalize[0.75];
Monitor[tabla2 = Table[solucion[Rationalize[Apr], θ0, 20], {θ0, -fraccion π, fraccion π,  (π/10)}];, θ0]

g2 = ParametricPlot[Table[{tabla2[[i]][[1]][t], tabla2[[i]][[2]][t]}, {i, 1, Length[tabla2]}], {t, 0, 10}, PlotRange-> All, PlotStyle-> Red]

La comparación con los retratos de fase hechos más arriba:

g2contour =ContourPlot[h[θ, P, Apr], {θ, -π, π}, {P, -2, 2},
 Contours->Table[h[θ0, 0, Apr], {θ0, -fraccion π, fraccion  π,  (π/20)}], 
PlotLabel->StringJoin["A=", ToString[Apr]], 
ImageSize->200]

Show[g2contour, g2, ImageSize-> 400]

Bárbaro!

El movimiento en el espacio 3D:

$HistoryLength = 0;
antialias[g_,n_: 1]:=ImageResize[Rasterize[g,"Image",ImageResolution->n 72, Background->None],Scaled[1/n]];
h[θ_,P_, A_]:= ((P^2)/2)+A Cos[θ]+(1/2)  ((1-Cos[θ])/(1+Cos[θ]));
solucion[A_, θ0_, tmax_]:= 
Module[{θ, ϕ, ψ, P,  ϕ0=1, ψ0=0, t, Dsol},
Dsol = NDSolve[
{
θ'[t] == 2 P[t],
ψ'[t] == ((Tan[((θ[t])/2)])^2),
ϕ'[t] == ((Sec[((θ[t])/2)])^2),
 P'[t] == -8  ((Csc[θ[t]])^3) ((Sin[((θ[t])/2)])^4)+2A Sin[θ[t]],
P[0] == 0,
ϕ[0] == ϕ0,
ψ[0] == ψ0,
θ[0] == θ0},
{θ,P, ψ, ϕ}, {t, 0, tmax}];
Dsol[[1]][[All, 2]]
];

Primero veamos cómo graficar la peonza para valores arbitrarios de los ángulos de Euler:

ez = {0,0,1};
eρ[ϕ_]:= { Cos[ϕ],  Sin[ϕ], 0};
eϕ[ϕ_]:= {-Sin[ϕ], Cos[ϕ], 0};

e1[θ_, ϕ_, ψ_]:= eρ[ϕ] Cos[ψ] + (Cos[θ] eϕ[ϕ]+ Sin[θ] ez) Sin[ψ];
e2[θ_, ϕ_, ψ_]:= -eρ[ϕ] Sin[ψ] +  (Cos[θ] eϕ[ϕ]+ Sin[θ] ez) Cos[ψ]; 
e3[θ_, ϕ_]:= Cos[θ] ez - Sin[θ] eϕ[ϕ];

peonza[θ_, ϕ_, ψ_, L_]:= Graphics3D[{
{Yellow, Tube[{{0,0,0},1.2L e3[θ, ϕ] }, 0.1]},
{Red, Tube[{{-e1[θ, ϕ, ψ]+L e3[θ, ϕ],+L e3[θ, ϕ] }}, 0.1]},
{Orange, Tube[{{L e3[θ, ϕ],+L e3[θ, ϕ]+e1[θ, ϕ, ψ] }}, 0.1]},
{Blue, Tube[{{-e2[θ, ϕ, ψ]+L e3[θ, ϕ],+L e3[θ, ϕ] }}, 0.1]},
{Green, Tube[{{L e3[θ, ϕ],+L e3[θ, ϕ]+e2[θ, ϕ, ψ] }}, 0.1]}
}];

Lpr = 3.;
Manipulate[
Show[
peonza[θ, ϕ, ψ, Lpr],
PlotRange-> 4{{-1, 1}, {-1, 1}, {-1, 1}},BoxRatios-> {1, 1, 1},Axes-> True, ViewAngle-> 20°],
{θ, 0, π, 0.01, Appearance-> "Open"}, {ϕ, 0, 8π, 0.01, Appearance-> "Open"},{ψ, 0, 8π, 0.01, Appearance-> "Open"}]

Ahora usamos las soluciones de las ecuaciones de movimiento: {θ,P, ψ, ϕ}.

Lpr = 3.;
Apr = 0.4;
 
(* Si A > 1/4, más-menos este valor son los ángulos de equilibrio *)
ArcCos[-1.+(1/(Sqrt[Apr]))]
(* +Si A > 1/4, más-menos esta cantidad marca el valor del ángulo inicial que separa las órbitas dentro y fuera de la separatriz *)
ArcCos[(1/(2Apr))-1.]

Calculamos la solución numérica de las ecs. de movimiento y graficamos algunas cosas antes de ir al gráfico 3D. Primero, la curva en el espacio de fases θ-(Subscript[p, θ]), para estar seguros de que estamos en el caso que queremos, y luego la trayectoria 3D que sigue el extremo de la peonza, para tener una idea general del movimiento.

θ0 =1.350099; (* Experimenten con este valor. Aquí está elegido a  ojo para tener una trayectoria casi cerrada *)
solpr = solucion[Rationalize[Apr],Rationalize[θ0],tmax =  100];  

orbitaθ = ParametricPlot[{solpr[[1]][t], solpr[[2]][t]}, {t, 0, 69.8611}, PlotStyle->Red];

Show[
Graphics3D[{Yellow, Opacity[0.5],Sphere[{0,0,0},1.2Lpr]}],
orbita3D = ParametricPlot3D[1.2Lpr e3[solpr[[1]][t], solpr[[4]][t]], {t, 0,100}, PlotStyle->{Blue,Thin}, PlotPoints->200], ImageSize-> 600,
PlotRange-> All
]

Cálculos auxiliares para ver dónde más o menos se cierra la trayectoria (lo que no siempre pasa)

Plot[{Log[((Norm[ e3[solpr[[1]][t], solpr[[4]][t]]-e3[solpr[[1]][0], solpr[[4]][0]]])^2)],Log[((Norm[ e3[solpr[[1]][t], solpr[[4]][t]]+e1[solpr[[1]][t], solpr[[4]][t], solpr[[3]][t]]-(e3[solpr[[1]][0], solpr[[4]][0]]+e1[solpr[[1]][0], solpr[[4]][0], solpr[[3]][0]])])^2)]}, {t, 68, 70}, PlotRange-> {-20, 3},PlotStyle->{Blue,Thin}, PlotPoints->200]

FindRoot[((Norm[ e3[solpr[[1]][t], solpr[[4]][t]]-e3[solpr[[1]][0], solpr[[4]][0]]])^2) == 0, {t, 68.87, 68.89}, MaxIterations->50000]

El retrato de fases θ-(Subscript[p, θ]) correspondiente al valor de A elegido, superpuesto con la órbita:

retrato = ContourPlot[h[θ, P, Apr], {θ, -3,3}, {P, -1,1},ContourShading->None,PlotPoints->100, ContourStyle->Thick,
 Contours->Join[{h0 = h[θ0, 0, Apr]},Table[(i/25) h0, {i, 1,40,2}]],FrameTicksStyle-> Directive[Black, 16],FrameLabel->{Style["θ", 20],Style["(Subscript[p, θ])", 20]} ,
ImageSize->625, PlotRange-> All, Background->White, ImageMargins->{{25,0},{0,0}}];
Show[retrato, orbitaθ, PlotRange-> All]

La cosa en sí

Auxiliar: las estelas de los extremos transversales de la peonza (ejecutar cada vez que cambian las condiciones iniciales)

θϕψtabla = Table[{solpr[[1]][τ], solpr[[4]][τ], solpr[[3]][τ]}, {τ, 0., tmax, 0.01}];
e1tabla = e1[#[[1]], #[[2]], #[[3]]]& /@ θϕψtabla ;
e2tabla = e2[#[[1]], #[[2]], #[[3]]]& /@ θϕψtabla ;
e3tabla = e3[#[[1]], #[[2]]]& /@ θϕψtabla ;
estela11 =Lpr e3tabla + e1tabla;
estela12 = Lpr e3tabla - e1tabla;
estela21 = Lpr e3tabla + e2tabla;
estela22 = Lpr e3tabla - e2tabla;

Finalmente, el movimiento 3D:

Lpr = 3.;
Δ = 3;
θΔ = solpr[[1]][Δ];
ψΔ = solpr[[3]][Δ];
ϕΔ = solpr[[4]][Δ];
peonzaΔ = peonza[θΔ, ϕΔ, ψΔ, Lpr];
Manipulate[
θ = solpr[[1]][t];
ψ = solpr[[3]][t];
ϕ = solpr[[4]][t];
pθ = solpr[[2]][t];
imin = Max[1, (piso =Floor[(t/0.01)])-coma];
imax = piso;
ipaso = 4;
psize = Table[((4.(i-imin) )/(imax-imin+0.001)), {i,imin, imax, ipaso}];
j = 0;

{Show[
Graphics3D[
Table[
++j;
{AbsolutePointSize[psize[[j]]],
{Orange, Point[estela11[[i]]]},
 {Red, Point[estela12[[i]]]},
{Green, Point[estela21[[i]]]},
{Blue, Point[estela22[[i]]]}}
, {i,imin, imax,ipaso}]

],
peonza[θ, ϕ, ψ, Lpr],
PlotRange-> 4{{-1, 1}, {-1, 1}, {-1.2, 1.2}},
Axes-> False,
ViewAngle->18°,
Background-> Black,
ImageSize-> 350
], 
Show[ retrato,Graphics[{Blue,AbsolutePointSize[8],Point[{θ,pθ}]}],ImageSize-> 360]
},

{{t,Δ}, Δ, tmax,0.05, Appearance-> "Open", AnimationRate-> 1},
{{coma, 60}, 0, 300, 10, Appearance-> "Open"},
TrackedSymbols:> {t, coma},
ControlPlacement->Above
]

Para exportar las imágenes:

t =3;
θ = solpr[[1]][t];
ψ = solpr[[3]][t];
ϕ = solpr[[4]][t];
pθ = solpr[[2]][t];
imin = Max[1, (piso =Floor[(t/0.01)])-200];
imax = piso;
ipaso = 1;
psize = Table[((6.(i-imin) )/(imax-imin+0.001)), {i,imin, imax, ipaso}];
j = 0;

im1 = antialias[Show[
Graphics3D[
Table[
++j;
{Opacity[(((((j/2)/((imax-imin)/ipaso))))^2)],AbsolutePointSize[psize[[j]]],
{Orange, Point[estela11[[i]]]},
 {Red, Point[estela12[[i]]]},
{Green, Point[estela21[[i]]]},
{Blue, Point[estela22[[i]]]}}
, {i,imin, imax,ipaso}]

],
peonza[θ, ϕ, ψ, Lpr],
PlotRange-> 4{{-1, 1}, {-1, 1}, {-1.2, 1.2}},
Axes-> False,
ViewAngle->18°,
Background-> Black,
ImageSize-> 500], 2];

dim1 = ImageDimensions[im1];

im2 =ImageResize[antialias[Show[ retrato,Graphics[{Blue,AbsolutePointSize[8],Point[{θ,pθ}]}],ImageSize-> 500],2], dim1];

imtemp =ImageAssemble[{{im1, im2}}]

dir = "C:\\users\\juan zanella\\desktop\\peonza\\";

nframe = 0;
Monitor[Table[
++nframe;
θ = solpr[[1]][t];
ψ = solpr[[3]][t];
ϕ = solpr[[4]][t];
pθ = solpr[[2]][t];
imin = Max[1, (piso =Floor[(t/0.01)])-200];
imax = piso;
ipaso = 1;
psize = Table[((6.(i-imin) )/(imax-imin+0.001)), {i,imin, imax, ipaso}];
j = 0;

im1 = antialias[Show[
Graphics3D[
Table[
++j;
{Opacity[(((((j/2)/((imax-imin)/ipaso))))^2)],AbsolutePointSize[psize[[j]]],
{Orange, Point[estela11[[i]]]},
 {Red, Point[estela12[[i]]]},
{Green, Point[estela21[[i]]]},
{Blue, Point[estela22[[i]]]}}
, {i,imin, imax,ipaso}]

],
peonza[θ, ϕ, ψ, Lpr],
PlotRange-> 4{{-1, 1}, {-1, 1}, {-1.2, 1.2}},
Axes-> False,
ViewAngle->18°,
Background-> Black,
ImageSize-> 500], 2];

dim1 = ImageDimensions[im1];

im2 =ImageResize[antialias[Show[ retrato,Graphics[{Blue,AbsolutePointSize[8],Point[{θ,pθ}]}],ImageSize-> 500],2], dim1];

imtemp =SetAlphaChannel[ImageAssemble[{{im1, im2}}]];

Export[StringJoin[dir, "frame_peonza_4_1_", If[nframe < 10, "000", If[nframe < 100, "00", "0"]], ToString[nframe], ".png"], imtemp, "PNG"],

{t, Δ, Δ+68.85,0.05}
];,{t, imtemp}]